本文大部分内容轉載自部落格【congwulong】 (https://blog.csdn.net/congwulong/article/details/7576012)
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. 函數接口
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變換公式
![各種二分查找[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)