请参考

2. FFT的DSP实现

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

[cpp] view plaincopyprint?

1. /* 
2.  * zx_fft.h 
3.  * 
4.  *  Created on: 2013-8-5 
5.  *      Author: monkeyzx 
6.  */  
7.   
8. #ifndef ZX_FFT_H_   
9. #define ZX_FFT_H_   
10.   
11. typedef float          FFT_TYPE;  
12.   
13. #ifndef PI   
14. #define PI             (3.14159265f)   
15. #endif   
16.   
17. typedef struct complex_st {  
18.     FFT_TYPE real;  
19.     FFT_TYPE img;  
20. } complex;  
21.   
22. int fft(complex *x, int N);  
23. int ifft(complex *x, int N);  
24. void zx_fft(void);  
25.   
26. #endif /* ZX_FFT_H_ */  

/*
 * 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 plaincopyprint?

1. /* 
2.  * zx_fft.c 
3.  * 
4.  * Implementation of Fast Fourier Transform(FFT) 
5.  * and reversal Fast Fourier Transform(IFFT) 
6.  * 
7.  *  Created on: 2013-8-5 
8.  *      Author: monkeyzx 
9.  */  
10.   
11. #include "zx_fft.h"   
12. #include <math.h>   
13. #include <stdlib.h>   
14.   
15. /* 
16.  * Bit Reverse 
17.  * === Input === 
18.  * x : complex numbers 
19.  * n : nodes of FFT. @N should be power of 2, that is 2^(*) 
20.  * l : count by bit of binary format, @l=CEIL{log2(n)} 
21.  * === Output === 
22.  * r : results after reversed. 
23.  * Note: I use a local variable @temp that result @r can be set 
24.  * to @x and won't overlap. 
25.  */  
26. static void BitReverse(complex *x, complex *r, int n, int l)  
27. {  
28.     int i = 0;  
29.     int j = 0;  
30.     short stk = 0;  
31.     static complex *temp = 0;  
32.   
33.     temp = (complex *)malloc(sizeof(complex) * n);  
34.     if (!temp) {  
35.         return;  
36.     }  
37.   
38.     for(i=0; i<n; i++) {  
39.         stk = 0;  
40.         j = 0;  
41.         do {  
42.             stk |= (i>>(j++)) & 0x01;  
43.             if(j<l)  
44.             {  
45.                 stk <<= 1;  
46.             }  
47.         }while(j<l);  
48.   
49.         if(stk < n) {             /* 满足倒位序输出 */  
50.             temp[stk] = x[i];  
51.         }  
52.     }  
53.     /* copy @temp to @r */  
54.     for (i=0; i<n; i++) {  
55.         r[i] = temp[i];  
56.     }  
57.     free(temp);  
58. }  
59.   
60. /* 
61.  * FFT Algorithm 
62.  * === Inputs === 
63.  * x : complex numbers 
64.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
65.  * === Output === 
66.  * the @x contains the result of FFT algorithm, so the original data 
67.  * in @x is destroyed, please store them before using FFT. 
68.  */  
69. int fft(complex *x, int N)  
70. {  
71.     int i,j,l,ip;  
72.     static int M = 0;  
73.     static int le,le2;  
74.     static FFT_TYPE sR,sI,tR,tI,uR,uI;  
75.   
76.     M = (int)(log(N) / log(2));  
77.   
78.     /* 
79.      * bit reversal sorting 
80.      */  
81.     BitReverse(x,x,N,M);  
82.   
83.     /* 
84.      * For Loops 
85.      */  
86.     for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */  
87.         le = (int)pow(2,l);  
88.         le2 = (int)(le / 2);  
89.         uR = 1;  
90.         uI = 0;  
91.         sR = cos(PI / le2);  
92.         sI = -sin(PI / le2);  
93.         for (j=1; j<=le2; j++) {   /* loop for each sub DFT */  
94.             //jm1 = j - 1;   
95.             for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */  
96.                 ip = i + le2;  
97.                 tR = x[ip].real * uR - x[ip].img * uI;  
98.                 tI = x[ip].real * uI + x[ip].img * uR;  
99.                 x[ip].real = x[i].real - tR;  
100.                 x[ip].img = x[i].img - tI;  
101.                 x[i].real += tR;  
102.                 x[i].img += tI;  
103.             }  /* Next i */  
104.             tR = uR;  
105.             uR = tR * sR - uI * sI;  
106.             uI = tR * sI + uI *sR;  
107.         } /* Next j */  
108.     } /* Next l */  
109.   
110.     return 0;  
111. }  
112.   
113. /* 
114.  * Inverse FFT Algorithm 
115.  * === Inputs === 
116.  * x : complex numbers 
117.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
118.  * === Output === 
119.  * the @x contains the result of FFT algorithm, so the original data 
120.  * in @x is destroyed, please store them before using FFT. 
121.  */  
122. int ifft(complex *x, int N)  
123. {  
124.     int k = 0;  
125.   
126.     for (k=0; k<=N-1; k++) {  
127.         x[k].img = -x[k].img;  
128.     }  
129.   
130.     fft(x, N);    /* using FFT */  
131.   
132.     for (k=0; k<=N-1; k++) {  
133.         x[k].real = x[k].real / N;  
134.         x[k].img = -x[k].img / N;  
135.     }  
136.   
137.     return 0;  
138. }  
139.   
140. /* 
141.  * Code below is an example of using FFT and IFFT. 
142.  */  
143. #define  SAMPLE_NODES              (128)   
144. complex x[SAMPLE_NODES];  
145. int INPUT[SAMPLE_NODES];  
146. int OUTPUT[SAMPLE_NODES];  
147.   
148. static void MakeInput()  
149. {  
150.     int i;  
151.   
152.     for ( i=0;i<SAMPLE_NODES;i++ )  
153.     {  
154.         x[i].real = sin(PI*2*i/SAMPLE_NODES);  
155.         x[i].img = 0.0f;  
156.         INPUT[i]=sin(PI*2*i/SAMPLE_NODES)*1024;  
157.     }  
158. }  
159.   
160. static void MakeOutput()  
161. {  
162.     int i;  
163.   
164.     for ( i=0;i<SAMPLE_NODES;i++ )  
165.     {  
166.         OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;  
167.     }  
168. }  
169.   
170. void zx_fft(void)  
171. {  
172.     MakeInput();  
173.   
174.     fft(x,128);  
175.     MakeOutput();  
176.   
177.     ifft(x,128);  
178.     MakeOutput();  
179. }