FFTW - 計算 FFT 的函式庫 (C 語言)

FFTW 是用 C 語言寫成的函式庫，用以計算 DFT 轉換。輸入的訊號可以是任意長度，可以是實數型式或是複數型式的訊號。當然，從它的名稱知道它是以 FFT 來實作。本來我已經有修改自 NUMERICAL RECIPES 的 fft 轉換程式，不過 NUMERICAL RECIPES 的授權限制滿大的，相較之下 FFTW 則是採用 GNU General Public License 授權，比較可以放心使用。

要在 VC 下使用 FFTW ，可以從這裡下載 zip 檔。

安裝︰

解壓縮之後，開啟 命令提示字元 ，轉換目錄到解壓縮的目錄，再執行以下三個指令︰

lib /machine:i386 /def:libfftw3-3.def
lib /machine:i386 /def:libfftw3f-3.def
lib /machine:i386 /def:libfftw3l-3.def

這會在目錄下建立三個 lib 檔。

將 libfftw3l-3.dll, libfftw3f-3.dll, libfftw3-3.dll 這3個 dll 檔複製到 system32 目錄。
在 VC 專案中指定 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個lib檔及 fftw3.h 檔所在的目錄。

在 VC 專案中加入 libfftw3l-3.lib, libfftw3f-3.lib, libfftw3-3.lib 這3個 lib。

使用︰
首先要引入 fftw3.h 這個標頭檔，如果要做的是一維 DFT 轉換，可參考以下的程式碼(取自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);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
...
fftw_execute(p); /* repeat as needed */
...
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
}

fftw_complex 是 FFTW 自訂的複數型態，是一個2個元素的浮點陣列 double[2]，若 in 為 fftw_complex 所宣告的一個陣列

fftw_complex * in = new fftw_complex[10];

，則第一個元素的實部為 in[0][0], 虛部為 in[0][1]。
fftw_malloc 是類似於 malloc 的函數，之間的差別我還不太清楚，不過文件上有提到，使用其他開記憶體的方式也是可以的(例如: new)，使用 fftw_malloc 的好處是，開啟的陣列元素(包括 fftw_complex 的實虛部)會初始化為0(並非如此)。

fftw_plan_dft_1d - 從函式名稱大概可以知道，這個函式是用來產生一個準備的工作，它並不會真的執行 DFT 轉換，真正的 DFT 轉換要再用 fftw_exectue 來執行。

多維度號的轉換只需要將 fftw_plan_dft_1d 轉為以下任一個函式。

fftw_plan fftw_plan_dft_2d(int nx, int ny,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

fftw_plan fftw_plan_dft_3d(int nx, int ny, int nz,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

fftw_plan fftw_plan_dft(int rank, const int *n,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);
輸出訊號格式︰
fftw_plan_dft_1d :
若輸入訊號的長度為 32 ，則輸出的FFT訊號從第17個元素(從0起算)開始重複。
out[17][0] == out[15][0], out[17][1] == out[15][1]
out[18][0] == out[14][0], out[18][1] == out[14][1]
                       .
                       .
out[31][0] == out[2][0], out[31][1] == out[2][1]


fftw_plan_dft_r2c_1d :
這個函式的傳入陣列為實數陣列，若輸入訊號長度為32，則輸出的FFT訊號會在
第17個元素之後均為0