FFT快速傅里葉變換蝶形算法詳解

本章目錄直接計(jì)算DFT的問題及改進(jìn)的途徑按時(shí)間抽取的基2-FFT算法

按頻率抽取的基2-FFT算法

快速傅里葉逆變換(IFFT)算法

Matlab實(shí)現(xiàn)1現(xiàn)在是1頁(yè)\一共有52頁(yè)\編輯于星期日5.1引言

DFT在實(shí)際應(yīng)用中很重要:可以計(jì)算信號(hào)的頻譜、功率譜和線性卷積等。直接按DFT變換進(jìn)行計(jì)算，當(dāng)序列長(zhǎng)度N很大時(shí)，計(jì)算量非常大，所需時(shí)間會(huì)很長(zhǎng)。FFT并不是一種與DFT不同的變換，而是DFT的一種快速計(jì)算的算法。

2現(xiàn)在是2頁(yè)\一共有52頁(yè)\編輯于星期日5.2直接計(jì)算DFT的問題及改進(jìn)的途徑

DFT的運(yùn)算量

設(shè)復(fù)序列x(n)長(zhǎng)度為N點(diǎn)，其DFT為k=0，，…，N-1（1）計(jì)算一個(gè)X(k)值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：N復(fù)數(shù)加法次數(shù)：N－13現(xiàn)在是3頁(yè)\一共有52頁(yè)\編輯于星期日5.2.1DFT的運(yùn)算量（2）計(jì)算全部N個(gè)X(k)值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：N2復(fù)數(shù)加法次數(shù)：N(N－1)（3）對(duì)應(yīng)的實(shí)數(shù)運(yùn)算量4現(xiàn)在是4頁(yè)\一共有52頁(yè)\編輯于星期日一次復(fù)數(shù)乘法：4次實(shí)數(shù)乘法2次實(shí)數(shù)加法＋一個(gè)X(k)：4N次實(shí)數(shù)乘法＋2N+2(N-1)=2(2N-1)次實(shí)數(shù)加法所以整個(gè)N點(diǎn)DFT運(yùn)算共需要：N×2(2N-1)=2N(2N-1)實(shí)數(shù)乘法次數(shù)：4N2實(shí)數(shù)加法次數(shù)：5現(xiàn)在是5頁(yè)\一共有52頁(yè)\編輯于星期日DFT運(yùn)算量的結(jié)論N點(diǎn)DFT的復(fù)數(shù)乘法次數(shù)舉例NN2NN22464404941612816384864256655361625651226214432102810241048576結(jié)論：當(dāng)N很大時(shí)，其運(yùn)算量很大，對(duì)實(shí)時(shí)性很強(qiáng)的信號(hào)處理來(lái)說(shuō)，要求計(jì)算速度快，因此需要改進(jìn)DFT的計(jì)算方法，以大大減少運(yùn)算次數(shù)。6現(xiàn)在是6頁(yè)\一共有52頁(yè)\編輯于星期日

5.2.2減少運(yùn)算工作量的途徑

主要原理是利用系數(shù)的以下特性對(duì)DFT進(jìn)行分解：（1）對(duì)稱性（2）周期性（3）可約性另外，7現(xiàn)在是7頁(yè)\一共有52頁(yè)\編輯于星期日5.3按時(shí)間抽取的基2-FFT算法

算法原理按時(shí)間抽取基-2FFT算法與直接計(jì)算DFT運(yùn)算量的比較按時(shí)間抽取的FFT算法的特點(diǎn)按時(shí)間抽取FFT算法的其它形式流程圖8現(xiàn)在是8頁(yè)\一共有52頁(yè)\編輯于星期日5.3.1算法原理

設(shè)N＝2L，將x(n)按n的奇偶分為兩組：

r=0，1，…，則9現(xiàn)在是9頁(yè)\一共有52頁(yè)\編輯于星期日式中，X1(k)和X2(k)分別是x1(n)和x2(n)的N/2的DFT。另外，式中k的取值范圍是：0，1，…，N/2－1。10現(xiàn)在是10頁(yè)\一共有52頁(yè)\編輯于星期日因此，只能計(jì)算出X(k)的前一半值。后一半X(k)值，N/2，N/2＋1，…，N？利用可得到同理可得11現(xiàn)在是11頁(yè)\一共有52頁(yè)\編輯于星期日考慮到因此可得后半部分X(k)及前半部分X(k)k=0，1，…，N/2－1k=0，1，…，N/2－112現(xiàn)在是12頁(yè)\一共有52頁(yè)\編輯于星期日蝶形運(yùn)算蝶形運(yùn)算式蝶形運(yùn)算信號(hào)流圖符號(hào)

因此，只要求出2個(gè)N/2點(diǎn)的DFT，即X1(k)和X2(k)，再經(jīng)過(guò)蝶形運(yùn)算就可求出全部X(k)的值，運(yùn)算量大大減少。13現(xiàn)在是13頁(yè)\一共有52頁(yè)\編輯于星期日以8點(diǎn)為例第一次按奇偶分解以N=8為例，分解為2個(gè)4點(diǎn)的DFT，然后做8/2=4次蝶形運(yùn)算即可求出所有8點(diǎn)X(k)的值。14現(xiàn)在是14頁(yè)\一共有52頁(yè)\編輯于星期日蝶形運(yùn)算量比較復(fù)數(shù)乘法次數(shù)：

N2復(fù)數(shù)加法次數(shù)：

N(N－1)復(fù)數(shù)乘法次數(shù)：

2*(N/2)2+N/2=N2/2+N/2復(fù)數(shù)加法次數(shù)：

2*(N/2)(N/2－1)+2*N/2=N2/2N點(diǎn)DFT的運(yùn)算量

分解一次后所需的運(yùn)算量＝2個(gè)N/2的DFT＋N/2蝶形：因此通過(guò)一次分解后，運(yùn)算工作量減少了差不多一半。

15現(xiàn)在是15頁(yè)\一共有52頁(yè)\編輯于星期日進(jìn)一步按奇偶分解

由于N＝2L，因而N/2仍是偶數(shù)，可以進(jìn)一步把每個(gè)N/2點(diǎn)子序列再按其奇偶部分分解為兩個(gè)N/4點(diǎn)的子序列。以N/2點(diǎn)序列x1(r)為例

則有

k=0,1,…,

16現(xiàn)在是16頁(yè)\一共有52頁(yè)\編輯于星期日且k=0,1,…,由此可見，一個(gè)N/2點(diǎn)DFT可分解成兩個(gè)N/4點(diǎn)DFT。同理，也可對(duì)x2(n)進(jìn)行同樣的分解，求出X2(k)。17現(xiàn)在是17頁(yè)\一共有52頁(yè)\編輯于星期日以8點(diǎn)為例第二次按奇偶分解18現(xiàn)在是18頁(yè)\一共有52頁(yè)\編輯于星期日算法原理

對(duì)此例N=8，最后剩下的是4個(gè)N/4=2點(diǎn)的DFT，2點(diǎn)DFT也可以由蝶形運(yùn)算來(lái)完成。以X3(k)為例。k=0,1即這說(shuō)明，N=2M的DFT可全部由蝶形運(yùn)算來(lái)完成。19現(xiàn)在是19頁(yè)\一共有52頁(yè)\編輯于星期日以8點(diǎn)為例第三次按奇偶分解N=8按時(shí)間抽取法FFT信號(hào)流圖20現(xiàn)在是20頁(yè)\一共有52頁(yè)\編輯于星期日5.3.2按時(shí)間抽取基2-FFT算法與直接計(jì)算DFT運(yùn)算量的比較

由按時(shí)間抽取法FFT的信號(hào)流圖可知，當(dāng)N=2L時(shí)，共有

級(jí)蝶形運(yùn)算；每級(jí)都由

個(gè)蝶形運(yùn)算組成，而每個(gè)蝶形有

次復(fù)乘、

次復(fù)加，因此每級(jí)運(yùn)算都需

次復(fù)乘和

次復(fù)加。

LN/2

N/2

12N21現(xiàn)在是21頁(yè)\一共有52頁(yè)\編輯于星期日這樣

級(jí)運(yùn)算總共需要：L復(fù)數(shù)乘法：

復(fù)數(shù)加法：

直接DFT算法運(yùn)算量復(fù)數(shù)乘法：

復(fù)數(shù)加法：

N2N(N－1)直接計(jì)算DFT與FFT算法的計(jì)算量之比為M22現(xiàn)在是22頁(yè)\一共有52頁(yè)\編輯于星期日FFT算法與直接DFT算法運(yùn)算量的比較NN2計(jì)算量之比MNN2計(jì)算量之比M2414.01281638444836.641644.025665536102464.0864125.45122621442304113.816256328.0102410485765120204.83210288012.82048419430411264372.464404919221.423現(xiàn)在是23頁(yè)\一共有52頁(yè)\編輯于星期日5.3.3按時(shí)間抽取的FFT算法的特點(diǎn)序列的逆序排列同址運(yùn)算（原位運(yùn)算）蝶形運(yùn)算兩節(jié)點(diǎn)間的距離的確定24現(xiàn)在是24頁(yè)\一共有52頁(yè)\編輯于星期日序列的逆序排列

由于x(n)被反復(fù)地按奇、偶分組，所以流圖輸入端的排列不再是順序的，但仍有規(guī)律可循：因?yàn)镹=2M，

對(duì)于任意n（0≤n≤N-1)，可以用M個(gè)二進(jìn)制碼表示為：

n反復(fù)按奇、偶分解時(shí)，即按二進(jìn)制碼的“0”“1”分解。序列的逆序排列25現(xiàn)在是25頁(yè)\一共有52頁(yè)\編輯于星期日倒位序的樹狀圖（N=8）

26現(xiàn)在是26頁(yè)\一共有52頁(yè)\編輯于星期日碼位的倒位序(N=8)

自然順序n二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序順序數(shù)000000001001100420100102301111064100001151011015611001137111111727現(xiàn)在是27頁(yè)\一共有52頁(yè)\編輯于星期日倒位序的變址處理（N=8）28現(xiàn)在是28頁(yè)\一共有52頁(yè)\編輯于星期日同址運(yùn)算（原位運(yùn)算）

某一列任何兩個(gè)節(jié)點(diǎn)k和j的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算后，得到結(jié)果為下一列k、j兩節(jié)點(diǎn)的節(jié)點(diǎn)變量，而和其他節(jié)點(diǎn)變量無(wú)關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲(chǔ)單元，降低設(shè)備成本。運(yùn)算前運(yùn)算后例同址運(yùn)算（原位運(yùn)算）29現(xiàn)在是29頁(yè)\一共有52頁(yè)\編輯于星期日觀察原位運(yùn)算規(guī)律30現(xiàn)在是30頁(yè)\一共有52頁(yè)\編輯于星期日蝶形運(yùn)算兩節(jié)點(diǎn)間的距離

以N=8為例：第一級(jí)蝶形，距離為：第二級(jí)蝶形，距離為：第三級(jí)蝶形，距離為：規(guī)律：對(duì)于共L級(jí)的蝶形而言，其m級(jí)蝶形運(yùn)算的節(jié)點(diǎn)間的距離為124蝶形運(yùn)算兩節(jié)點(diǎn)間的距離

31現(xiàn)在是31頁(yè)\一共有52頁(yè)\編輯于星期日

的確定

以N=8為例：

的確定

32現(xiàn)在是32頁(yè)\一共有52頁(yè)\編輯于星期日5.4按頻率抽取的基2-FFT算法

算法原理再把輸出X(k)按k的奇偶分組先把輸入按n的順序分成前后兩半設(shè)序列長(zhǎng)度為N=2L，L為整數(shù)前半子序列x(n)后半子序列

0≤n≤0≤n≤33現(xiàn)在是33頁(yè)\一共有52頁(yè)\編輯于星期日5.4.1算法原理由DFT定義得k=0，1，…，N34現(xiàn)在是34頁(yè)\一共有52頁(yè)\編輯于星期日由于所以則k=0，1，…，N35現(xiàn)在是35頁(yè)\一共有52頁(yè)\編輯于星期日然后按k的奇偶可將X(k)分為兩部分r=0，1，…，則式可轉(zhuǎn)化為36現(xiàn)在是36頁(yè)\一共有52頁(yè)\編輯于星期日令n=0，1，…，代入r=0，1，…，可得為2個(gè)N/2點(diǎn)的DFT，合起來(lái)正好是N點(diǎn)X(k)的值。37現(xiàn)在是37頁(yè)\一共有52頁(yè)\編輯于星期日蝶形運(yùn)算將稱為蝶形運(yùn)算與時(shí)間抽選基2FFT算法中的蝶形運(yùn)算符號(hào)略有不同。38現(xiàn)在是38頁(yè)\一共有52頁(yè)\編輯于星期日例按頻率抽取(N=8)

例按頻率抽取，將N點(diǎn)DFT分解為兩個(gè)N/2點(diǎn)DFT的組合(N=8)39現(xiàn)在是39頁(yè)\一共有52頁(yè)\編輯于星期日

與時(shí)間抽取法的推導(dǎo)過(guò)程一樣，由于N=2L，N/2仍然是一個(gè)偶數(shù)，因而可以將每個(gè)N/2點(diǎn)DFT的輸出再分解為偶數(shù)組與奇數(shù)組，這就將N/2點(diǎn)DFT進(jìn)一步分解為兩個(gè)N/4點(diǎn)DFT。N=840現(xiàn)在是40頁(yè)\一共有52頁(yè)\編輯于星期日5.4.2頻率抽取法與時(shí)間抽取法的異同

頻率抽取法輸入是自然順序，輸出是倒位序的；時(shí)間抽取法正好相反。頻率抽取法的基本蝶形與時(shí)間抽取法的基本蝶形有所不同。頻率抽取法運(yùn)算量與時(shí)間抽取法相同。頻率抽取法與時(shí)間抽取法的基本蝶形是互為轉(zhuǎn)置的。

41現(xiàn)在是41頁(yè)\一共有52頁(yè)\編輯于星期日5.5快速傅里葉逆變換(IFFT)算法IDFT公式DFT公式比較可以看出，IDFT多出M個(gè)1/2可分解到M級(jí)蝶形運(yùn)算中。42現(xiàn)在是42頁(yè)\一共有52頁(yè)\編輯于星期日例頻率抽取IFFT流圖(N=8)

43現(xiàn)在是43頁(yè)\一共有52頁(yè)\編輯于星期日快速傅里葉逆變換另一種算法44現(xiàn)在是44頁(yè)\一共有52頁(yè)\編輯于星期日5.8Matlab實(shí)現(xiàn)用FFT進(jìn)行譜分析的Matlab實(shí)現(xiàn)用CZT進(jìn)行譜分析的Matlab實(shí)現(xiàn)在Matlab中使用的線性調(diào)頻z變換函數(shù)為czt，其調(diào)用格式為>>X=czt(x,M,W,A)其中，x是待變換的時(shí)域信號(hào)x(n)，其長(zhǎng)度為N，M是變換的長(zhǎng)度，W確定變換的步長(zhǎng)，A確定變換的起點(diǎn)。若M=N，A=1，則CZT變成DFT。45現(xiàn)在是45頁(yè)\一共有52頁(yè)\編輯于星期日5.8.1用FFT進(jìn)行譜分析的Matlab實(shí)現(xiàn)例5.1設(shè)模擬信號(hào)，以t=0.01n(n=0:N-1)進(jìn)行取樣，試用fft函數(shù)對(duì)其做頻譜分析。N分別為：(1)N=45；(2)N=50；(3)N=55；(2)N=60。程序清單如下%計(jì)算N=45的FFT并繪出其幅頻曲線N=45;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y))title('FFTN=45')46現(xiàn)在是46頁(yè)\一共有52頁(yè)\編輯于星期日例5.1程序清單%計(jì)算N=50的FFT并繪出其幅頻曲線N=50;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y))title('FFTN=50')47現(xiàn)在是47頁(yè)\一共有52頁(yè)\編輯于星期日%計(jì)算