轉(zhuǎn)自CSDN博主：FFTW中文參考

一、 簡(jiǎn)介

先看一下使用FFTW編程的方法：

? #include <fftw3.h>

...

{

fftw_complex*in, *out;

fftw_planp;

...

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

// 輸入數(shù)據(jù)in賦值

p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(p);

// 執(zhí)行變換

...

fftw_destroy_plan(p);

fftw_free(in);

fftw_free(out);

}

大致是先用fftw_malloc分配輸入輸出內(nèi)存，然后輸入數(shù)據(jù)賦值，然后創(chuàng)建變換方案（fftw_plan），然后執(zhí)行變換（fftw_execute），最后釋放資源，還是比較簡(jiǎn)單的。

二、 一維復(fù)數(shù)據(jù)的DFT

1. 數(shù)據(jù)類(lèi)型

fftw_complex默認(rèn)由兩個(gè)double組成，在內(nèi)存中順序排列，實(shí)部在 前，虛部在后，即typedef double fftw_complex[2]。FFTW文檔指出如果有一個(gè)支持C99標(biāo)準(zhǔn)的C編譯器（如gcc），可以在#include 前加入#include ，這樣一來(lái)fftw_complex就被定義為本機(jī)復(fù)數(shù)類(lèi)型，而且與上述typedef二進(jìn)制兼容（指內(nèi)存排列），經(jīng) 測(cè)試不能用在Windows下。C++有一個(gè)復(fù)數(shù)模板類(lèi)complex，在頭文件下定義。C++標(biāo)準(zhǔn)委 員會(huì)最近同意該類(lèi)的存儲(chǔ)方式與C99二進(jìn)制兼容，即順序存儲(chǔ)，實(shí)部在前，虛部在后（見(jiàn)報(bào)告WG21/N1388），該解決方案在所有主流標(biāo)準(zhǔn)庫(kù)實(shí)現(xiàn)中都能正確工作。所以實(shí)際上可以用complex 來(lái)代替fftw_complex，比如有一個(gè)復(fù)數(shù)數(shù)組complex *x，則可以將其類(lèi)型轉(zhuǎn)換后作為參數(shù)傳遞給fftw：reinterpret_cast(x)。測(cè)試如下：開(kāi) 兩個(gè)數(shù)組fftw_complex x1[2]和complex x2[2]，然后賦相同值，在調(diào)試模式下可以看到它們的內(nèi)存排列是相同的。complex類(lèi)數(shù)據(jù)賦值的方式不是很直接，必須采用無(wú)名對(duì)象方式x2[i] = complex (1,2) 或成員函數(shù)方式x2[i].real(1);x2[i].imag(2);不能直接寫(xiě)x2[i].real=1;x2[i].imag=2。 fftw_complex賦值方式比較直接：x1[i][0]=1;x1[i][1]=2。最后，考慮到數(shù)據(jù)對(duì)齊（見(jiàn)后），最好使用 fftw_malloc分配內(nèi)存，所以可以將其返回的指針轉(zhuǎn)換為complex *類(lèi)型使用（比如賦值或讀取等），變換時(shí)再將其轉(zhuǎn)換為fftw_complex*。

2. 函數(shù)接口

fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);

n為數(shù)據(jù)個(gè)數(shù)，可以為任意正整數(shù)，但如果為一些小因子的乘積計(jì)算起來(lái)可以更有效，不過(guò)即使n為素?cái)?shù)算法仍然能夠達(dá)到O(nlogn)的復(fù)雜度。FFTW對(duì)N=2a?3b?5c?7d?11e?13f的變換處理得最好，其中e+f=0/1，其它冪指數(shù)可以為任意值。

如果in和out指針相同為原位運(yùn)算，否則為非原位運(yùn)算。

sign可以為正變換FFTW_FORWARD(-1)，也可以為逆變換FFTW_BACKWORD(+1)，實(shí)際上就是變換公式中指數(shù)項(xiàng)的符號(hào)。需注意FFTW的逆變換沒(méi)有除以N，即數(shù)據(jù)正變換再反變換后是原始數(shù)據(jù)的N倍。

flags參數(shù)一般情況下為FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW會(huì)先計(jì)算一些FFT并測(cè)量所用的時(shí)間，以便為大小為n的變換尋找最優(yōu)的計(jì)算方法。依據(jù) 機(jī)器配置和變換的大?。╪），這個(gè)過(guò)程耗費(fèi)約數(shù)秒（時(shí)鐘clock精度）。FFTW_ESTIMATE則相反，它直接構(gòu)造一個(gè)合理的但可能是次最優(yōu)的方 案?？傮w來(lái)說(shuō)，如果你的程序需要進(jìn)行大量相同大小的FFT，并且初始化時(shí)間不重要，可以使用FFTW_MEASURE，否則應(yīng)使用 FFTW_ESTIMATE。FFTW_MEASURE模式下in和out數(shù)組中的值會(huì)被覆蓋，所以該方式應(yīng)該在用戶初始化輸入數(shù)據(jù)in之前完成。

不知道上述說(shuō)法是不是這個(gè)意思：先用FFTW_MEASURE模式自動(dòng)選最優(yōu)方案，速度較慢；然后使用該模式變換數(shù)據(jù)就會(huì)較快。示例代碼為：

intlength = 50000;

fftw_complex* din? = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2);

fftw_complex* dout = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2);

fftw_planp? = fftw_plan_dft_1d(length, din, din, FFTW_FORWARD, FFTW_MEASURE);

fftw_execute(p);

// 輸入數(shù)據(jù)din賦值// ...

fftw_execute(p);

// 讀取變換結(jié)果// ...

fftw_destroy_plan(p);

fftw_free(din);

fftw_free(dout);

實(shí)驗(yàn)發(fā)現(xiàn)第一個(gè)fftw_execute耗費(fèi)了數(shù)秒，而第二個(gè)fftw_execute則瞬間完成，說(shuō)明上述猜想可能是對(duì)的。

創(chuàng)建完方案（fftw_plan）后，就可以用fftw_execute對(duì)指定的 數(shù)據(jù)in/out做任意次變換。如果想變換一個(gè)相同大小（N相等）但數(shù)據(jù)不同的另外一個(gè)數(shù)組in，可以創(chuàng)建一個(gè)新方案，F(xiàn)FTW會(huì)自動(dòng)重用上次方案的信 息。這一點(diǎn)其實(shí)是非常好的，比如你首先用FFTW_MEASURE模式創(chuàng)建了一個(gè)最優(yōu)的變換方案，只要變換數(shù)據(jù)的大小不變，你可以用 fftw_plan_dft_1d創(chuàng)建新的方案以對(duì)新數(shù)據(jù)執(zhí)行變換，同時(shí)新變換仍然是最優(yōu)的。一個(gè)fftw_plan只能對(duì)固定的in/out進(jìn)行變換， 但可以在變換后改變in的內(nèi)容（大小不變）以用同一個(gè)方案執(zhí)行新的變換。

三、 多維復(fù)數(shù)據(jù)的DFT

? ? fftw_plan fftw_plan_dft_2d(int n0, intn1,

fftw_complex *in, fftw_complex*out,

int sign, unsignedflags);

fftw_plan fftw_plan_dft_3d(int n0, int n1, intn2,

fftw_complex *in, fftw_complex*out,

int sign, unsignedflags);

fftw_plan fftw_plan_dft(int rank, const int*n,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

多維復(fù)數(shù)據(jù)的DFT同一維復(fù)數(shù)據(jù)的DFT用法類(lèi)似，數(shù)組in/out為行優(yōu)先方式 存儲(chǔ)。fftw_plan_dft是一個(gè)通用的復(fù)DFT函數(shù)，可以執(zhí)行一維、二維或多維復(fù)DFT。比如對(duì)于圖像（2維數(shù)據(jù)），其變換為 fftw_plan_dft_2d（height,width,85）,因?yàn)樵紙D像數(shù)據(jù)為height×width的矩陣，即第一維(n0)為行數(shù) height。

四、 一維實(shí)數(shù)據(jù)的DFT

函數(shù)接口

fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out, unsignedflags);

fftw_plan fftw_plan_dft_c2r_1d(int n, fftw_complex *in, double *out, unsigned flags);

r2c版本：實(shí)輸入數(shù)據(jù)，復(fù)Hermitian輸出，正變換。

c2r版本：復(fù)Hermitian輸入數(shù)據(jù)，實(shí)輸出數(shù)據(jù)，逆變換。

n：邏輯長(zhǎng)度，不必為物理長(zhǎng)度。由于實(shí)數(shù)據(jù)的DFT具有 Hermitian對(duì)稱(chēng)性，所以只需要計(jì)算n/2+1（向下取整）個(gè)輸出就可以了。比如對(duì)于r2c，輸入in有n個(gè)數(shù)據(jù)，輸出out有floor(n /2)＋1個(gè)數(shù)據(jù)。對(duì)于原位運(yùn)算，in和out為同一數(shù)組（out須強(qiáng)制類(lèi)型轉(zhuǎn)換），所以其必須足夠大以容納所有數(shù)據(jù)，長(zhǎng)度為2*(n/2+1)，in數(shù) 組的前n個(gè)數(shù)據(jù)為輸入數(shù)據(jù)，后面的數(shù)據(jù)用來(lái)保存輸出。

c2r逆變換在任何情況下（不管是否為原位運(yùn)算）都破壞輸入數(shù)組in，如果有必要可以通過(guò)在flags中加入FFTW_PRESERVE_INPUT來(lái)阻止，但這會(huì)損失一些性能，而其這個(gè)標(biāo)志位目前在多維實(shí)DFT中已不被支持。

比如對(duì)于n=4，in=[1 2 3 4]，out=[10 -2+2i -2 -2-2i]，out具有共軛對(duì)稱(chēng)性，out的實(shí)際內(nèi)存為10 0 -2 2 -2 0，共3個(gè)復(fù)數(shù)數(shù)據(jù)。對(duì)于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的實(shí)際內(nèi)存為15 0 -2.5 3.44 -2.5 0.81，共3個(gè)復(fù)數(shù)數(shù)據(jù)。

實(shí)數(shù)據(jù)DFT中，首個(gè)變換數(shù)據(jù)為直流分量，為實(shí)數(shù)；如果n為偶數(shù)，由 Nyquist采樣定理，第n/2個(gè)變換數(shù)據(jù)也為實(shí)數(shù)；所以可以把這兩個(gè)數(shù)據(jù)組合在一起，即將第n/2個(gè)變換數(shù)據(jù)作為第0個(gè)變換數(shù)據(jù)的虛部，這樣一來(lái)輸入 數(shù)組就和輸出數(shù)組相等（n=n/2*2）。一些FFT的實(shí)現(xiàn)就是這么做的，但FFTW沒(méi)有這么做，因?yàn)樗⒉荒芎芎玫赝茝V到多維DFT中，而且存儲(chǔ)空間的 節(jié)省也是非常小以至于可以忽略不計(jì)。

一個(gè)一維c2r和r2c DFT的替代接口可以在r2r接口中找到，它利用了半復(fù)數(shù)輸出類(lèi)型（即實(shí)部和虛部分開(kāi)放在不通的區(qū)域里），使輸出數(shù)組具有和輸入數(shù)組同樣的長(zhǎng)度和類(lèi)型。該接口在多維變換中用處不大，但有時(shí)可能會(huì)有一些性能的提升。

五、 多維實(shí)數(shù)據(jù)的DFT

? ? fftw_plan fftw_plan_dft_r2c_2d(int n0, intn1,

double *in, fftw_complex*out,

unsignedflags);

fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, intn2,

double *in, fftw_complex*out,

unsignedflags);

fftw_plan fftw_plan_dft_r2c(int rank, const int*n,

double *in, fftw_complex*out,

unsigned flags);

這是r2c接口（正變換），c2r接口（逆變換）只是簡(jiǎn)單的將輸入/輸出類(lèi)型交換一下。用法大致同1d情況一樣，需要特別注意的是復(fù)數(shù)據(jù)的存放方式。對(duì)于n0×n1×n1×…×nd-1的實(shí)輸入數(shù)組（行優(yōu)先），經(jīng)過(guò)r2c正變換后，輸出為一個(gè)n0×n1×n1×…×(nd-1/2+1)的復(fù)數(shù)（fftw_complex）數(shù)組（行優(yōu)先），其中除法向下取整。由于復(fù)數(shù)數(shù)據(jù)的總長(zhǎng)度要大于實(shí)數(shù)據(jù)，所以如果需要進(jìn)行原位運(yùn)算則輸入實(shí)數(shù)組必須擴(kuò)展以能夠容納所有復(fù)數(shù)據(jù)，即實(shí)數(shù)組的最后一維必須包含2(floor(nd-1/2)+1)個(gè)double元素。比如對(duì)于一個(gè)2維實(shí)正DFT，輸入數(shù)組為n0×n1大小，輸出復(fù)數(shù)組大小為n0×floor(n1/2+1)（由對(duì)稱(chēng)性），其長(zhǎng)度大于實(shí)輸入數(shù)組。所以對(duì)于原位運(yùn)算，輸入數(shù)組要擴(kuò)展到n0×2floor(n1/2+1)。如果n1為偶數(shù)，擴(kuò)展為n0×（n1＋2）；如果n1為奇數(shù)，擴(kuò)展為n0×（n1＋1）；這些擴(kuò)展的內(nèi)存不需要賦初值，它們只用來(lái)存放輸出數(shù)據(jù)。

比如對(duì)于3×3的實(shí)正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的實(shí)際內(nèi)存為32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06；即為 3×2的復(fù)數(shù)組，換算成double類(lèi)型為3×4，所以如果要進(jìn)行原位運(yùn)算，in數(shù)組大小必須為3×4，即最后一維（每行）擴(kuò)展一個(gè)double元素。另 外，如果用這個(gè)out數(shù)組進(jìn)行3×3的c2r逆變換，將得到實(shí)數(shù)據(jù)[0 18 36;54 9 27;45 63 36]，即原始數(shù)據(jù)的9（n0×n1）倍，這是因?yàn)镕FTW的逆變換沒(méi)有歸一化。