快速傅裡葉變換（FFT）源代碼

為了看明白那堆積分變換，不得不把複變函數掃了一遍，可看完了，才發現原來這堆變換說白了隻是一些數字遊戲，也沒用到啥複變函數的知識。最後，用C++程式實作了下FFT，也算告一段落，代碼如下：

1. #include <iostream>
2. #include <fstream>
3. #include <math.h>
4. using namespace std;
5. const double PI = 3.14159265358979323846;
6. int n;  // 資料個數 = 2的logn次方
7. int logn;
8. /// 複數結構體
9. struct stCompNum 
10. {
11.     double re;
12.     double im;
13. };
14. stCompNum* pData1 = NULL;
15. stCompNum* pData2 = NULL;
16. /// 正整數位逆序後輸出
17. int reverseBits(int value, int bitCnt)
18. {
19.     int i;
20.     int ret = 0;
21.     for(i=0; i<bitCnt; i++)
22.     {
23.         ret |= (value & 0x1) << (bitCnt - 1 - i);
24.         value >>= 1;
25.     }
26.     return ret;
27. }
28. void main()
29. {
30.     ifstream fin("data.txt");
31.     int i,j,k;
32.     // input logn
33.     fin>>logn;
34.     // calculate n
35.     for(i=0, n=1; i<logn; i++) n *= 2;
36.     // malloc memory space
37.     pData1 = new stCompNum[n];
38.     pData2 = new stCompNum[n];
39.     // input raw data
40.     for(i=0; i<n; i++) fin>>pData1[i].re;
41.     for(i=0; i<n; i++) fin>>pData1[i].im;
42.     // FFT transform
43.     int cnt = 1;
44.     for(k=0; k<logn; k++)
45.     {
46.         for(j=0; j<cnt; j++)
47.         {
48.             int len = n / cnt;
49.             double c = - 2 * PI / len;
50.             for(i=0; i<len/2; i++)
51.             {
52.                 int idx = len * j + i;
53.                 pData2[idx].re = pData1[idx].re + pData1[idx + len/2].re;
54.                 pData2[idx].im = pData1[idx].im + pData1[idx + len/2].im;
55.             }
56.             for(i=len/2; i<len; i++)
57.             {
58.                 double wcos = cos(c * (i - len/2));
59.                 double wsin = sin(c * (i - len/2));
60.                 int idx = len * j + i;
61.                 stCompNum tmp;
62.                 tmp.re = pData1[idx - len/2].re - pData1[idx].re;
63.                 tmp.im = pData1[idx - len/2].im - pData1[idx].im;
64.                 pData2[idx].re = tmp.re * wcos - tmp.im * wsin;
65.                 pData2[idx].im = tmp.re * wsin + tmp.im * wcos;
66.             }
67.         }
68.         cnt <<= 1;
69.         stCompNum* pTmp = NULL;
70.         pTmp   = pData1;
71.         pData1 = pData2;
72.         pData2 = pTmp;
73.     }
74.     // resort
75.     for(i=0; i<n; i++)
76.     {
77.         int rev = reverseBits(i, logn);
78.         stCompNum tmp;
79.         if(rev > i)
80.         {
81.             tmp         = pData1[i];
82.             pData1[i]   = pData1[rev];
83.             pData1[rev] = tmp;
84.         }
85.     }
86.     // output result data
87.     for(i=0; i<n; i++) cout<<pData1[i].re<<"/t";
88.     cout<<endl;
89.     for(i=0; i<n; i++) cout<<pData1[i].im<<"/t";
90.     cout<<endl;
91.     // free memory space
92.     delete []pData1;
93.     delete []pData2;
94.     fin.close();
95.     system("pause");
96. }

輸入檔案data.txt的内容如下：

4

2.2 4.5 6.7 8.5 10.2 12.3 14.5 16.2 19.3 21.2 25.2 29.4 36.4 39.2 45.2 50.1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0