FFTW(Fastest Fourier Transform in the West)是世界上最快的FFT, 實測計算長度為10000的double數組, 單次運行時間在2ms左右。為了詳細了解FFTW以及為編程方便,特將用戶手冊看了一下,并結合手冊制作了以下FFTW中文參考。其中大部分是原文重點內容的翻譯,并加入了一些注解。
??先看一下使用FFTW編程的方法:
????? #include <fftw3.h> ???? ... ???? { ???????? fftw_complex *in, *out; ???????? fftw_plan p; ???????? ... ???????? in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); ???????? out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); ????? ???// 輸入數據in賦值 ???????? p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); ???????? fftw_execute(p); // 執行變換 ???????? ... ???????? fftw_destroy_plan(p); ???????? fftw_free(in); ???????? fftw_free(out); ???? } |
大致是先用fftw_malloc分配輸入輸出內存,然后輸入數據賦值,然后創建變換方案(fftw_plan),然后執行變換(fftw_execute),最后釋放資源,還是比較簡單的。
一維復數據的DFT
? 1. 數據類型
? ?fftw_complex默認由兩個double組成,在內存中順序排列,實部在 前,虛部在后,即typedef double fftw_complex[2]。FFTW文檔指出如果有一個支持C99標準的C編譯器(如gcc),可以在#include <fftw3.h>前加入#include <complex.h>,這樣一來fftw_complex就被定義為本機復數類型,而且與上述typedef二進制兼容(指內存排列),經 測試不能用在Windows下。C++有一個復數模板類complex<T>,在頭文件<complex>下定義。C++標準委 員會最近同意該類的存儲方式與C99二進制兼容,即順序存儲,實部在前,虛部在后(見報告WG21/N1388),該解決方案在所有主流標準庫實現中都能正確工作。所以實際上可以用complex <double> 來代替fftw_complex,比如有一個復數數組complex<double> *x,則可以將其類型轉換后作為參數傳遞給fftw:reinterpret_cast<fftw_complex*>(x)。測試如下:開 兩個數組fftw_complex x1[2]和complex<double> x2[2],然后賦相同值,在調試模式下可以看到它們的內存排列是相同的。complex<T>類數據賦值的方式不是很直接,必須采用無名對象方式x2[i] = complex <double>(1,2) 或成員函數方式x2[i].real(1);x2[i].imag(2);不能直接寫x2[i].real=1;x2[i].imag=2。 fftw_complex賦值方式比較直接:x1[i][0]=1;x1[i][1]=2。最后,考慮到數據對齊(見后),最好使用 fftw_malloc分配內存,所以可以將其返回的指針轉換為complex <double> *類型使用(比如賦值或讀取等),變換時再將其轉換為fftw_complex*。
? 2. 函數接口
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); |
?n為數據個數,可以為任意正整數,但如果為一些小因子的乘積計算起來可以更有效,不過即使n為素數算法仍然能夠達到O(nlogn)的復雜度。FFTW對N=2a 3b 5c 7d 11e 13f的變換處理得最好,其中e+f=0/1,其它冪指數可以為任意值。
? ?如果in和out指針相同為原位運算,否則為非原位運算。
? ?sign可以為正變換FFTW_FORWARD(-1),也可以為逆變換FFTW_BACKWORD(+1),實際上就是變換公式中指數項的符號。需注意FFTW的逆變換沒有除以N,即數據正變換再反變換后是原始數據的N倍。
? ?flags參數一般情況下為FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW會先計算一些FFT并測量所用的時間,以便為大小為n的變換尋找最優的計算方法。依據 機器配置和變換的大小(n),這個過程耗費約數秒(時鐘clock精度)。FFTW_ESTIMATE則相反,它直接構造一個合理的但可能是次最優的方案。總體來說,如果你的程序需要進行大量相同大小的FFT,并且初始化時間不重要,可以使用FFTW_MEASURE,否則應使用 FFTW_ESTIMATE。FFTW_MEASURE模式下in和out數組中的值會被覆蓋,所以該方式應該在用戶初始化輸入數據in之前完成。
? ?不知道上述說法是不是這個意思:先用FFTW_MEASURE模式自動選最優方案,速度較慢;然后使用該模式變換數據就會較快。示例代碼為:
? int length = 50000; ? fftw_complex* din? = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); ? fftw_complex* dout = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); ? fftw_plan p?? = fftw_plan_dft_1d(length, din, din, FFTW_FORWARD, FFTW_MEASURE); ? fftw_execute(p); ? // 輸入數據din賦值 ? // ... ? fftw_execute(p); ? // 讀取變換結果 ? // ... ? fftw_destroy_plan(p); ? fftw_free(din); ? fftw_free(dout); |
實驗發現第一個fftw_execute耗費了數秒,而第二個fftw_execute則瞬間完成,說明上述猜想可能是對的。
? ? 創建完方案(fftw_plan)后,就可以用fftw_execute對指定的 數據in/out做任意次變換。如果想變換一個相同大小(N相等)但數據不同的另外一個數組in,可以創建一個新方案,FFTW會自動重用上次方案的信息。這一點其實是非常好的,比如你首先用FFTW_MEASURE模式創建了一個最優的變換方案,只要變換數據的大小不變,你可以用 fftw_plan_dft_1d創建新的方案以對新數據執行變換,同時新變換仍然是最優的。
一個fftw_plan只能對固定的in/out進行變換, 但可以在變換后改變in的內容(大小不變)以用同一個方案執行新的變換。
?軟件整體界面如下:
1、可設置正弦信號的頻率、采樣頻率、采樣點數、直流信號、是否加噪聲。
2、點擊“計算FFT”按鈕可生成頻譜。
3、還可以進行自動計算。信號頻譜按10Hz依次增加,再依次計算這個頻率對應的頻譜。運行效果如下:
將原始信號修改為 sin(f0) + sin(f1) + sin(f2) + sin(f3)? 模擬4個不同頻率的信號合成信號。
再看它的FFT變換結果。
源碼出售,完整代碼,包教會?