FFT的DSP实现

fish001

FFT的DSP实现 [复制链接]

使用C语言实现的FFT及IFFT算法实例，能计算任意以2为对数底的采样点数的FFT，算法参考上面给的流程图。

[cpp] view plain copy
print?

• /*
• * zx_fft.h
• *
• *  Created on: 2013-8-5
• *      Author: monkeyzx
• */  
•   
• #ifndef ZX_FFT_H_  
• #define ZX_FFT_H_  
•   
• typedef float          FFT_TYPE;  
•   
• #ifndef PI  
• #define PI             (3.14159265f)  
• #endif  
•   
• typedef struct complex_st {  
•     FFT_TYPE real;  
•     FFT_TYPE img;  
• } complex;  
•   
• int fft(complex *x, int N);  
• int ifft(complex *x, int N);  
• void zx_fft(void);  
•   
• #endif /* ZX_FFT_H_ */  
  

[cpp] view plain copy
print?

• /*
• * zx_fft.c
• *
• * Implementation of Fast Fourier Transform(FFT)
• * and reversal Fast Fourier Transform(IFFT)
• *
• *  Created on: 2013-8-5
• *      Author: monkeyzx
• */  
•   
• #include "zx_fft.h"  
• #include   
• #include   
•   
• /*
• * Bit Reverse
• * === Input ===
• * x : complex numbers
• * n : nodes of FFT. @N should be power of 2, that is 2^(*)
• * l : count by bit of binary format, @l=CEIL{log2(n)}
• * === Output ===
• * r : results after reversed.
• * Note: I use a local variable @temp that result @r can be set
• * to @x and won't overlap.
• */  
• static void BitReverse(complex *x, complex *r, int n, int l)  
• {  
•     int i = 0;  
•     int j = 0;  
•     short stk = 0;  
•     static complex *temp = 0;  
•   
•     temp = (complex *)malloc(sizeof(complex) * n);  
•     if (!temp) {  
•         return;  
•     }  
•   
•     for(i=0; i
•         stk = 0;  
•         j = 0;  
•         do {  
•             stk |= (i>>(j++)) & 0x01;  
•             if(j
•             {  
•                 stk <<= 1;  
•             }  
•         }while(j
•   
•         if(stk < n) {             /* 满足倒位序输出 */  
•             temp[stk] = x;  
•         }  
•     }  
•     /* copy @temp to @r */  
•     for (i=0; i
•         r = temp;  
•     }  
•     free(temp);  
• }  
•   
• /*
• * FFT Algorithm
• * === Inputs ===
• * x : complex numbers
• * N : nodes of FFT. @N should be power of 2, that is 2^(*)
• * === Output ===
• * the @x contains the result of FFT algorithm, so the original data
• * in @x is destroyed, please store them before using FFT.
• */  
• int fft(complex *x, int N)  
• {  
•     int i,j,l,ip;  
•     static int M = 0;  
•     static int le,le2;  
•     static FFT_TYPE sR,sI,tR,tI,uR,uI;  
•   
•     M = (int)(log(N) / log(2));  
•   
•     /*
•      * bit reversal sorting
•      */  
•     BitReverse(x,x,N,M);  
•   
•     /*
•      * For Loops
•      */  
•     for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */  
•         le = (int)pow(2,l);  
•         le2 = (int)(le / 2);  
•         uR = 1;  
•         uI = 0;  
•         sR = cos(PI / le2);  
•         sI = -sin(PI / le2);  
•         for (j=1; j<=le2; j++) {   /* loop for each sub DFT */  
•             //jm1 = j - 1;  
•             for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */  
•                 ip = i + le2;  
•                 tR = x[ip].real * uR - x[ip].img * uI;  
•                 tI = x[ip].real * uI + x[ip].img * uR;  
•                 x[ip].real = x.real - tR;  
•                 x[ip].img = x.img - tI;  
•                 x.real += tR;  
•                 x.img += tI;  
•             }  /* Next i */  
•             tR = uR;  
•             uR = tR * sR - uI * sI;  
•             uI = tR * sI + uI *sR;  
•         } /* Next j */  
•     } /* Next l */  
•   
•     return 0;  
• }  
•   
• /*
• * Inverse FFT Algorithm
• * === Inputs ===
• * x : complex numbers
• * N : nodes of FFT. @N should be power of 2, that is 2^(*)
• * === Output ===
• * the @x contains the result of FFT algorithm, so the original data
• * in @x is destroyed, please store them before using FFT.
• */  
• int ifft(complex *x, int N)  
• {  
•     int k = 0;  
•   
•     for (k=0; k<=N-1; k++) {  
•         x[k].img = -x[k].img;  
•     }  
•   
•     fft(x, N);    /* using FFT */  
•   
•     for (k=0; k<=N-1; k++) {  
•         x[k].real = x[k].real / N;  
•         x[k].img = -x[k].img / N;  
•     }  
•   
•     return 0;  
• }  
•   
• /*
• * Code below is an example of using FFT and IFFT.
• */  
• #define  SAMPLE_NODES              (128)  
• complex x[SAMPLE_NODES];  
• int INPUT[SAMPLE_NODES];  
• int OUTPUT[SAMPLE_NODES];  
•   
• static void MakeInput()  
• {  
•     int i;  
•   
•     for ( i=0;i
•     {  
•         x.real = sin(PI*2*i/SAMPLE_NODES);  
•         x.img = 0.0f;  
•         INPUT=sin(PI*2*i/SAMPLE_NODES)*1024;  
•     }  
• }  
•   
• static void MakeOutput()  
• {  
•     int i;  
•   
•     for ( i=0;i
•     {  
•         OUTPUT = sqrt(x.real*x.real + x.img*x.img)*1024;  
•     }  
• }  
•   
• void zx_fft(void)  
• {  
•     MakeInput();  
•   
•     fft(x,128);  
•     MakeOutput();  
•   
•     ifft(x,128);  
•     MakeOutput();  
• }