Linux下简单的FFT实现

Linux下简单的FFT实现

在Linux平台下，从音频接口采回PCM数据，进行FFT变换。最后通过ncurse库，在文本模式下粗糙得对频谱进行实时显示。

整个程序只是作为练手用，代码贴在下面，供自己查阅和大家参考。

Makefile

1. all: directshow.o fft.o   
2.     gcc -o spectrum directshow.o fft.o -lm -lcurses   
3. directshow.o: directshow.c fft.h   
4.     gcc -c directshow.c    
5. fft.o: fft.c   
6.     gcc -c fft.c    

fft.c

1. /*******************************************************************************  
2. * File:        fft.c  
3. * Author:      Jerome Wang   
4. *              wangquzhi@gmail.com  
5. * Last Edit:   2011/01/19  
6. * Description:  void fft(double *DataTimeDomain_r, double *DataTimeDomain_i);  
7. *     Basic funtion used for fft. Two parameters, for the real-part and j-part  
8. * of the input data. The results are also restored here.  
9. *     Chose the point of fft by define one of these macros:  
10. *     FFT_CONFIG_4  
11. *     FFT_CONFIG_8  
12. *     FFT_CONFIG_128  
13. *       
14. *     By define FFT_DEBUGE, we can get the debug message indicating the process  
15. * of the fft processing.  
16. *******************************************************************************/  
17. #include <math.h>   
18. #include <stdio.h>   
19. #define FFT_CONFIG_128   
20. /*#define FFT_DEBUGE*/  
21. #define PI 3.14   
22. inline void Butterfly(double *a_r, double *a_i, \   
23.                       double *b_r, double * b_i,\   
24.                       double wx, int N)   
25. {   
26.     double a_tmp_r = *a_r;   
27.     double a_tmp_i = *a_i;       
28.     double b_tmp_r = *b_r;   
29.     double b_tmp_i = *b_i;   
30.        
31.     double wx_tmp_r, wx_tmp_i;   
32.     wx_tmp_r =  cos(wx*2*PI/N);   
33.     wx_tmp_i = -sin(wx*2*PI/N);   
34.     *a_r = a_tmp_r + b_tmp_r*wx_tmp_r - b_tmp_i*wx_tmp_i;   
35.     *a_i = a_tmp_i + b_tmp_r*wx_tmp_i + b_tmp_i*wx_tmp_r;   
36.     *b_r = a_tmp_r - b_tmp_r*wx_tmp_r + b_tmp_i*wx_tmp_i;   
37.     *b_i = a_tmp_i - b_tmp_r*wx_tmp_i - b_tmp_i*wx_tmp_r;   
38. }   
39. inline int BitReverse(int n, int bits)   
40. {   
41.     int i = 0;   
42.     int ret = 0;   
43.     for(i = 0; i < bits; i++)   
44.     {   
45.         ret|=((n>>i)&1)<<(bits-i-1);   
46.     }   
47.        
48.     return ret;   
49. }   
50. void fft(double *DataTimeDomain_r, double *DataTimeDomain_i)   
51. {   
52.     #ifdef FFT_CONFIG_128   
53.     const int FFT_Pnt = 128;   
54.     /* For 2^7 = 128, 128 points' computation is divided to 7 stages */  
55.     const int FFT_Stg = 7;   
56.     /* In each stage, it contains BttrPerStg butterfly computations */  
57.     const int BttrPerStg = 64;   
58.     #endif   
59.     #ifdef FFT_CONFIG_8   
60.     const int FFT_Pnt = 8;   
61.     const int FFT_Stg = 3;   
62.     const int BttrPerStg = 4;   
63.     #endif   
64.     #ifdef FFT_CONFIG_4   
65.     const int FFT_Pnt = 4;   
66.     const int FFT_Stg = 2;   
67.     const int BttrPerStg = 2;   
68.     #endif   
69.     /* Current Computing Stage, ranges from 0 to (FFT_Stg - 1) */       
70.     int CrrtStg = 0;   
71.     /* Current Butterfly, ranges from 0 to (BttrPerStg - 1) */       
72.     int CrrtBttr = 0;   
73.        
74.     int TempPower = 0;   
75.     int BttrPara1 = 0;   
76.     int BttrPara2 = 0;   
77.     int BttrPara3 = 0;   
78.   
79.     for(CrrtStg = 0; CrrtStg < FFT_Stg; CrrtStg++)   
80.     {   
81.         #ifdef FFT_DEBUGE   
82.         printf("This is the %dnd Stage:\n", CrrtStg+1);   
83.         #endif   
84.         for(CrrtBttr = 0; CrrtBttr < BttrPerStg; CrrtBttr++)   
85.         {   
86.             TempPower = 1 << CrrtStg;   
87.             BttrPara1 = (CrrtBttr/TempPower)*(TempPower<<1)+CrrtBttr%TempPower;   
88.             BttrPara2 = BttrPara1 + TempPower;           
89.             BttrPara3 = (CrrtBttr % TempPower) * \   
90.                             (1 << (BttrPerStg - CrrtStg - 1)) >> 1;       
91.             #ifdef FFT_DEBUGE   
92.             printf("\tNow we perform the butterfly: %d,%d,%d\n", \   
93.                                         BitReverse(BttrPara1, FFT_Stg),\   
94.                                         BitReverse(BttrPara2, FFT_Stg),\   
95.                                         BttrPara3);   
96.             #endif   
97.             Butterfly(DataTimeDomain_r + BitReverse(BttrPara1, FFT_Stg),\   
98.                         DataTimeDomain_i + BitReverse(BttrPara1, FFT_Stg),\   
99.                         DataTimeDomain_r + BitReverse(BttrPara2, FFT_Stg),\   
100.                         DataTimeDomain_i + BitReverse(BttrPara2, FFT_Stg),\   
101.                         BttrPara3,\   
102.                         FFT_Pnt);   
103.         }   
104.     }   
105. }  

fft.h

1. #ifndef FFT_H   
2. #define FFT_H   
3. void fft(double *DataTimeDomain_r, double *DataTimeDomain_i);   
4. #endif