據說FFTW（Fastest Fourier Transform in the West）是世界上最快的FFT。為了詳細了解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的內容（大小不變）以用同一個方案執行新的變換。
三、 多維復數據的DFT
fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2, 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);
多維復數據的DFT同一維復數據的DFT用法類似，數組in/out為行優先方式 存儲。fftw_plan_dft是一個通用的復DFT函數，可以執行一維、二維或多維復DFT。比如對於圖像（2維數據），其變換為 fftw_plan_dft_2d（height,width,85）,因為原始圖像數據為height×width的矩陣，即第一維(n0)為行數 height。
四、 一維實數據的DFT
函數接口
fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex *in, double *out, unsigned flags);
r2c版本：實輸入數據，復Hermitian輸出，正變換。
c2r版本：復Hermitian輸入數據，實輸出數據，逆變換。
n：邏輯長度，不必為物理長度。由於實數據的DFT具有 Hermitian對稱性，所以只需要計算n/2+1（向下取整）個輸出就可以了。比如對於r2c，輸入in有n個數據，輸出out有floor(n /2)＋1個數據。對於原位運算，in和out為同一數組（out須強制類型轉換），所以其必須足夠大以容納所有數據，長度為2*(n/2+1)，in數 組的前n個數據為輸入數據，后面的數據用來保存輸出。
c2r逆變換在任何情況下（不管是否為原位運算）都破壞輸入數組in，如果有必要可以通過在flags中加入FFTW_PRESERVE_INPUT來阻止，但這會損失一些性能，而其這個標志位目前在多維實DFT中已不被支持。
比如對於n=4，in=[1 2 3 4]，out=[10 -2+2i -2 -2-2i]，out具有共軛對稱性，out的實際內存為10 0 -2 2 -2 0，共3個復數數據。對於n=5，in=[1 2 3 4 5]，out=[15 -2.5+3.44i -2.5+0.81i -2.5-0.81i -2.5-3.44i] ，out的實際內存為15 0 -2.5 3.44 -2.5 0.81，共3個復數數據。
實數據DFT中，首個變換數據為直流分量，為實數；如果n為偶數，由 Nyquist采樣定理，第n/2個變換數據也為實數；所以可以把這兩個數據組合在一起，即將第n/2個變換數據作為第0個變換數據的虛部，這樣一來輸入 數組就和輸出數組相等（n=n/2*2）。一些FFT的實現就是這么做的，但FFTW沒有這么做，因為它並不能很好地推廣到多維DFT中，而且存儲空間的 節省也是非常小以至於可以忽略不計。
一個一維c2r和r2c DFT的替代接口可以在r2r接口中找到，它利用了半復數輸出類型（即實部和虛部分開放在不通的區域里），使輸出數組具有和輸入數組同樣的長度和類型。該接口在多維變換中用處不大，但有時可能會有一些性能的提升。
五、 多維實數據的DFT
fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags);
這是r2c接口（正變換），c2r接口（逆變換）只是簡單的將輸入/輸出類型交換一下。用法大致同1d情況一樣，需要特別注意的是復數據的存放方式。對於n0×n1×n1×…×nd-1的實輸入數組（行優先），經過r2c正變換后，輸出為一個n0×n1×n1×…×(nd-1/2+1)的復數（fftw_complex）數組（行優先），其中除法向下取整。由於復數數據的總長度要大於實數據，所以如果需要進行原位運算則輸入實數組必須擴展以能夠容納所有復數據，即實數組的最后一維必須包含2(floor(nd-1/2)+1)個double元素。比如對於一個2維實正DFT，輸入數組為n0×n1大小，輸出復數組大小為n0×floor(n1/2+1)（由對稱性），其長度大於實輸入數組。所以對於原位運算，輸入數組要擴展到n0×2floor(n1/2+1)。如果n1為偶數，擴展為n0×（n1＋2）；如果n1為奇數，擴展為n0×（n1＋1）；這些擴展的內存不需要賦初值，它們只用來存放輸出數據。
比如對於3×3的實正DFT，in=[0 2 4;6 1 3;5 7 4]，out=[32 0.5+0.86i 0.5-0.86i;-7+5.2i -1-1.73i -8.5-6.06i;-7-5.2i -8.5+6.06i -1+1.73i]；out的實際內存為32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06；即為 3×2的復數組，換算成double類型為3×4，所以如果要進行原位運算，in數組大小必須為3×4，即最后一維（每行）擴展一個double元素。另 外，如果用這個out數組進行3×3的c2r逆變換，將得到實數據[0 18 36;54 9 27;45 63 36]，即原始數據的9（n0×n1）倍，這是因為FFTW的逆變換沒有歸一化。
六、 更多實數據的DFT
通過一個統一的r2r（real-to-real，實數－實數）接口，FFTW支持其它的一些變換類型，這些變換的輸入和輸出數組大小相同。這些r2r變換可以分為3個類型：DFT的實 數據輸入，complex-Hermitian（指復Hermitian對稱）以半復數格式的輸出；DCT/DST（離散正余弦變換）；DHT（離散 Hartley變換）。接口如下：
fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out, fftw_r2r_kind kind, unsigned flags); fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, unsigned flags); fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2, unsigned flags); fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags);
這里n為數組的物理尺寸。對於多維變換，數組按行優先方式存儲（與C++標准相同，與Fortran不同）。由於DFT是可分離變換，所以2維/3維/多維的變換是在每個維度上分別進行變換得到的，每個維度都可指定一個kind參數，指定該維的變換類型。
1. 半復數格式DFT（HalfComplex-format）
對於大小為n的1維DFT，輸出格式如下：
r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
上述(n+1)/2向下取整。rk是第k個輸出的實部，ik是 第k個輸出的虛部。對於一個半復數格式的數組hc[n]，第k個元素的實部為hc[k]，虛部為[n-k]；k==0或n/2（n為偶數）情況除外，這兩 種情況下虛部為0，不存儲。所以對於r2hc(real-half complex，正變換)變換，輸入輸出數組大小都為n，hc2r（half complex- real，逆變換）變換也是如此。除了數據格式的差異，r2hc和hc2r變換的結果與前述r2c和c2r的變換結果是相同的。
對於多維比如2維變換，由可分離性，可以先對行 變換，再對列變換，FFTW_R2HC方式行變換的結果是半復數行，如果采用FFTW_R2HC方式進行列變換，需要進行一些數據處理，r2r變換不會自 動處理，需要手動進行，所以對於多維實數據變換，推薦使用普通的r2c/c2r接口。
2. DCT/DST
DCT可以認為是實偶對稱數據DFT（REDFT,Real-Even DFT）, DST可以認為是實奇對稱數據DFT（RODFT,Real-Odd DFT）。REDFTab和RODFTab中的a,b為數據移位標志（1表示移位），這些構成了DCT和DST的I－IV類，比較常用的為DCT-II，FFTW支持所有這些類型的變換。
FFTW_REDFT00 (DCT-I): even around j=0 and even around j=n-1. FFTW_REDFT10 (DCT-II, the DCT): even around j=-0.5 and even around j=n-0.5. FFTW_REDFT01 (DCT-III, the IDCT): even around j=0 and odd around j=n. FFTW_REDFT11 (DCT-IV): even around j=-0.5 and odd around j=n-0.5. FFTW_RODFT00 (DST-I): odd around j=-1 and odd around j=n. FFTW_RODFT10 (DST-II): odd around j=-0.5 and odd around j=n-0.5. FFTW_RODFT01 (DST-III): odd around j=-1 and even around j=n-1. FFTW_RODFT11 (DST-IV): odd around j=-0.5 and even around j=n-0.5.
對稱性只是邏輯意義上的，對物理輸入數據沒有任何限制。比如對於n＝5的REDFT00 (DCT-I)，輸入數據為abcde，它對應n=8的abcdedcb的常規DFT；對於n＝4的REDFT10 (DCT-II)，輸入數據為abcd，它對應n=8的abcddcba的常規DFT。
所有這些變換都是可逆的。R*DFT00的逆變 換是R*DFT00，R*DFT10的逆變換是R*DFT01（即DCT和IDCT），R*DFT11的逆變換是R*DFT11。如同DFT一樣，這些變 換的結果都沒有歸一化，所以正變換再逆變換后數據變為原始數據的N倍，N為邏輯DFT大小。比如對於REDFT00變換，N=2(n-1)；對於 RODFT00，N=2n。
注意n=1的REDFT00對應與N=0的DFT，所以它是未定義的（返回值為NULL的fftw_plan）。
R*DFT01和R*DFT10要比R*DFT11稍微快一些，尤其對於奇數長度數據；而R*DFT00則要慢一些，尤其對於奇數長度數據。
比如對於in=[1 2 3 4]，n=4，N=2n=8。Matlab下dct變換的結果為[5 -2.2304 0 -0.15851]；FFTW的結果為（FFTW_REDFT10）out=[20 -6.3086 0 -0.4483415]，為matlab結果的√8（√N）倍；用out進行逆dct變換（FFTW_REDFT01）的結果為[8 16 24 32]，為原始數據的8（N）倍。
再比如對於in=[0 2 4;6 1 3;5 7 4]的二維DCT變換，n=3，N=2n=6。Matlab下dct2的變換結果為out_matlab=[10.667 0 0.4714;-4.0825 -2.5 1.4434;0.4714 -2.5981 -3.1667]；FFTW的結果為（fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT10, FFTW_REDFT10, FFTW_ESTIMATE）out_fftw=[128 0 4;-34.641 -15 8.66;4 -15.588 -19]，這與Matlab的結果同樣是有差別的。將Matlab的結果變換到FFTW結果的步驟為：
1. 將out_matlab乘以√6×√6（即√N×√N）；
2. 再將上一步得到的out_matlab的第一行和第一列都乘以√2，因此第一個元素（即左上角的元素）要乘以2。
第一個是歸一化的原因，第二個是FFTW為了將DCT變換與實偶對稱FFT相對應的結果。這些對於DCT變換的應用都影響不大。（見此）
最后對out_fftw進行IDCT變換 （fftw_plan_r2r_2d(3, 3, in, out, FFTW_REDFT01, FFTW_REDFT01, FFTW_ESTIMATE），將得到[0 72 144;216 36 108;180 252 144]；它是原始數據in的36（N×N，N=6）倍。
3. 其它
fftw_malloc考慮了數據對齊，以便使 用SIMD指令加速，所以最好不要用C函數malloc替代，而且不要將fftw_malloc、fftw_free和malloc、free、 delete等混用。盡量使用fftw_malloc分配數組，而不是下面的固定數組，因為固定數組是在棧上分配的，而棧空間較小；還因為這種方式沒有考 慮數據對齊，不便應用SIMD指令。
fftw_complex data[N0][N1][N2]; fftw_plan plan; ... plan = fftw_plan_dft_3d(N0, N1, N2, &data[0][0][0], &data[0][0][0], FFTW_FORWARD, FFTW_ESTIMATE); ...
對於多維數組也盡量使用fftw_malloc(n0*n1*n285*sizeof(double))，不要使用C的malloc。
fftw_complex *a_good_array; a_good_array = (fftw_complex*) fftw_malloc(5*12*27* sizeof(fftw_complex)); fftw_complex ***a_bad_array; /* another way to make a 5x12x27 array */ a_bad_array = (fftw_complex ***) malloc(5 * sizeof(fftw_complex **));
七、 函數參考
1. 復數DFT
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2, 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);
2. 實數DFT
fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags); fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags);
3. 實數－實數變換
fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out, fftw_r2r_kind kind, unsigned flags); fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, unsigned flags); fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2, unsigned flags); fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags);
4. 實數－實數變換類型
對於大小為n的下列變換，對應的邏輯DFT大小為N，N用來進行歸一化。FFTW的變換沒有歸一化，正變換后再逆變換為原數據的N倍（不是n倍），對於多維變換，為N的乘積(N0*N1*N285)。下表列出了變換類型及其對應的邏輯大小N及逆變換：
FFTW_R2HC computes a real-input DFT with output in halfcomplex format, i.e. real and imaginary parts for a transform of size n stored as:r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1 (Logical N=n, inverse is FFTW_HC2R.) FFTW_HC2R computes the reverse of FFTW_R2HC, above. (Logical N=n, inverse is FFTW_R2HC.) FFTW_DHT computes a discrete Hartley transform. (Logical N=n, inverse is FFTW_DHT.) FFTW_REDFT00 computes an REDFT00 transform, i.e. a DCT-I. (Logical N=2*(n-1), inverse is FFTW_REDFT00.) FFTW_REDFT10 computes an REDFT10 transform, i.e. a DCT-II (sometimes called the DCT). (Logical N=2*n, inverse is FFTW_REDFT01.) FFTW_REDFT01 computes an REDFT01 transform, i.e. a DCT-III (sometimes called the IDCT, being the inverse of DCT-II). (Logical N=2*n, inverse is FFTW_REDFT=10.) FFTW_REDFT11 computes an REDFT11 transform, i.e. a DCT-IV. (Logical N=2*n, inverse is FFTW_REDFT11.) FFTW_RODFT00 computes an RODFT00 transform, i.e. a DST-I. (Logical N=2*(n+1), inverse is FFTW_RODFT00.) FFTW_RODFT10 computes an RODFT10 transform, i.e. a DST-II. (Logical N=2*n, inverse is FFTW_RODFT01.) FFTW_RODFT01 computes an RODFT01 transform, i.e. a DST-III. (Logical N=2*n, inverse is FFTW_RODFT=10.) FFTW_RODFT11 computes an RODFT11 transform, i.e. a DST-IV. (Logical N=2*n, inverse is FFTW_RODFT11.)
八、 其它
1. 數據類型
FFTW有三個版本的數據類型：double、float和long double，使用方法如下：
- 鏈接對應的庫（比如libfftw3-3、libfftw3f-3、或ibfftw3l-3）
- 包含同樣的頭文件fftw3.h
- 將所有以小寫"fftw_"開頭的名字替換為"fftwf_"（float版本）或"fftwl_"（long double版本）。比如將fftw_complex替換為fftwf_complex，將fftw_execute替換為fftwf_execute等。
- 所有以大寫"FFTW_"開頭的名字不變
- 將函數參數中的double替換為float或long double
最后，雖然long double是C99的標准，但你的編譯器可能根本不支持該類型，或它並不能提供比double更高的精度。
2. 用同一個fftw_plan執行多個數據的變換
前面說過可以利用同一個fftw_plan通過對輸入數據賦不同值來實現不同的變換，實際上還可以利用同一個fftw_plan實現對不同輸入輸出數據的變換，也就是說可以有多個輸入輸出數據數組，各自進行變換，互不影響。當然這要滿足一定的條件：
- 輸入/輸出數據大小相等。
- 變換類型、是否原位運算不變。
- 對split數組（指實虛部分開），實部和虛部的分割方式與方案創建時相同。
- 數組的對齊方式相同。如果都是用fftw_malloc分配的則此項條件滿足，除非使用 FFTW_UNALIGNED標志。
如果想對新數組，比如大小相等的一批數組執行變換，可以使用以下接口：
void fftw_execute_dft( const fftw_plan p, fftw_complex *in, fftw_complex *out); void fftw_execute_split_dft( const fftw_plan p, double *ri, double *ii, double *ro, double *io); void fftw_execute_dft_r2c( const fftw_plan p, double *in, fftw_complex *out); void fftw_execute_split_dft_r2c( const fftw_plan p, double *in, double *ro, double *io); void fftw_execute_dft_c2r( const fftw_plan p, fftw_complex *in, double *out); void fftw_execute_split_dft_c2r( const fftw_plan p, double *ri, double *ii, double *out); void fftw_execute_r2r( const fftw_plan p, double *in, double *out);
這些函數的執行不會修改原始plan，並且可以和fftw_execute混合使用。
3. 多線程FFTW
FFTW可以多線程執行，但是多線程存在線程同步問題，這可能會降低性能。所以除非問題規模非常大，一般並不能從多線程中獲益。
4. FFTW變換公式
