快速傅里葉變換

第四章快速傅里葉變換有限長(zhǎng)序列可以通過(guò)離散傅里葉變換①FT)將其頻域也離散化成有限長(zhǎng)序列.但其計(jì)算量太大,很難實(shí)時(shí)地處理問(wèn)題,因此引出了快速傅里葉變換(FFT).1965年，Cooley和Tukey提出了計(jì)算離散傅里葉變換(DFT)的快速算法，將DFT的運(yùn)算量減少了幾個(gè)數(shù)量級(jí)。從此，對(duì)快速傅里葉變換(FFT)算法的研究便不斷深入，數(shù)字信號(hào)處理這門(mén)新興學(xué)科也隨FFT的出現(xiàn)和發(fā)展而迅速發(fā)展。根據(jù)對(duì)序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，基本算法是基2DIT和基2DIF。FFT在離散傅里葉反變換、線(xiàn)性卷積和線(xiàn)性相關(guān)等方面也有重要應(yīng)用。快速傅里葉變換(FFT)是計(jì)算離散傅里葉變換(DFT)的快速算法。DFT的定義式為X(k)=婚x(n)WknR(k)NNnn在所有復(fù)指數(shù)值Wkn的值全部已算好的情況下，要計(jì)算一個(gè)X(k)需要NN次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法。算出全部N點(diǎn)X(k)共需N2次復(fù)數(shù)乘法和N(N■1)次復(fù)數(shù)加法。即計(jì)算量是與N2成正比的。FFT的基本思想：將大點(diǎn)數(shù)的DFT分解為若干個(gè)小點(diǎn)數(shù)DFT的組合，從而減少運(yùn)算量。WN因子具有以下兩個(gè)特性，可使DFT運(yùn)算量盡量分解為小點(diǎn)數(shù)的DFT運(yùn)算：(1)周期性：W(kBN)n■Wkn■W(nWkN NN(2)對(duì)稱(chēng)性：W(k?N/2)BIWkNN利用這兩個(gè)性質(zhì)，可以使DFT運(yùn)算中有些項(xiàng)合并，以減少乘法次數(shù)。例子:求當(dāng)N=4時(shí)，X(2)的值X(2)■?X(n)W2n■x(0)W0■X⑴W2■x(2)W4■x(3)W64 44 44nID■[X(0)■X(2)]W0■[X(1)■X(3)]W2(周期性)44=[x(0)■x⑵]一[x⑴■x⑶]W0(對(duì)稱(chēng)性)4通過(guò)合并，使乘法次數(shù)由4次減少到1次，運(yùn)算量減少。FFT的算法形式有很多種，但基本上可以分為兩大類(lèi)：按時(shí)間抽取(DIT)和按頻率抽取(DIF)。按時(shí)間抽取(DIT)的FTT為了將大點(diǎn)數(shù)的DFT分解為小點(diǎn)數(shù)的DFT運(yùn)算，要求序列的長(zhǎng)度N為復(fù)合數(shù)，最常用的是N■2M的情況(M為正整數(shù))。該情況下的變換稱(chēng)為基2FFT。下面討論基2情況的算法。先將序列x(n)按奇偶項(xiàng)分解為兩組x(2r)■x(r) N1 r■0,1,■，—■1x(2r■1)■x(r) 22將DFT運(yùn)算也相應(yīng)分為兩組X(k)■DFT[x(n)]■ x(n)WknNnIDx(n)Wkn＋x(n)Wkn

NNnnn為偶數(shù)nnnn為偶數(shù)n為奇數(shù)=.x(2r)W2rk■.x(2r■1)W(2r■)kNNrn rn="x(="x(r)W2r■WkNI/^^^x1NrID(r)W2rk2N=/^x(=/^x(r)Wrk1N/2rn■WkNIH^'x(r)WkN2N/2rID(因?yàn)閃2r■Wrk)

N N/2X(k)■WkX(k)

1 N2其中X(k)、X(k)分別是x(n)、x(n)的其中X(k)、X(k)分別是x(n)、x(n)的N/2點(diǎn)的DFTX(k)=N^HIjc(r)Wk1 1 N/2rnN^H1c(2r)Wk,0■k■N■1N/2 2rnX(k)=*x(r)Wrk2 2N/2r?N2x(2r■1)Wrk,0■k■NnN/2 2r?至此，一個(gè)N點(diǎn)DFT被分解為兩個(gè)N/2點(diǎn)的DFT。上面是否將全部N點(diǎn)的X(k)求解出來(lái)了？N分析：X1(k)和X2(k)只有N/2個(gè)點(diǎn)(k■0,1，…3?1),貝|由X(k)■X(k)^WkX(k)只能求出X(k)的前N/2個(gè)點(diǎn)的DFT，要求出N2全部N點(diǎn)的X(k)，需要找出X(k)、X(k)和X(k■N/2)的關(guān)系，其N(xiāo)中k■0,1，…3?1。由式子X(jué)(k)?X1(k)?“X2(k)可得X(k■N/2)■X(k■N/2)BWkz/2X(k■N/2)化簡(jiǎn)得X(k■N/2)=■X(k)BWkX(k)N2Nk■0,1,???，一■12這樣N點(diǎn)DFT可全部由下式確定出來(lái):IX(k)■X(k)■WkX(k)N2IX(k■N/2)■X(k)■WkX(k)1 N2Nk■0,1,…，一■1 (%)2上式可用一個(gè)專(zhuān)用的碟形符號(hào)來(lái)表示，這個(gè)符號(hào)對(duì)應(yīng)一次復(fù)乘和兩次復(fù)加運(yùn)算。a^WkbNWk -1a^WkbNWk -1N圖蝶形運(yùn)算符號(hào)aIWkbNNN2通過(guò)這樣的分解以后，每一個(gè)N/2點(diǎn)的DFT只需要(萬(wàn))2■4次復(fù)數(shù)乘法，兩個(gè)N/2點(diǎn)的DFT需要2(NN)2■N次復(fù)乘，再加上將兩個(gè)N/2點(diǎn)DFT合并成為N點(diǎn)DFT時(shí)有N／2次與W因子相乘，一共需要N2NN2—"y■一廠次復(fù)乘?？梢?jiàn)，通過(guò)這樣的分解，運(yùn)算量節(jié)省了近一半。乙 乙 乙因?yàn)镹■2m,N/2仍然是偶數(shù)，因此可以對(duì)兩個(gè)N/2點(diǎn)的DFT再分別作進(jìn)一步的分解，將兩個(gè)N/2點(diǎn)的DFT分解成兩個(gè)N/4點(diǎn)的DFT。例如對(duì)“一),可以在按其偶數(shù)部分及奇數(shù)部分進(jìn)行分解:N1■0,1，…，一■14x(2N1■0,1，…，一■141 3(21■1)■x(l)14則的運(yùn)算可相應(yīng)分為兩組：X(k)=.x(21)W2ik■.x(21■1)W(21■)kTOC\o"1-5"\h\z1 1 N/2 1 N/21n 1n=(1)W1k3N/41ID=(1)W1k3N/41IDN/2 4N/41n■X(k)■X(k)■WkX(k)N/2 4k■0,1,…，一■14將系數(shù)統(tǒng)一為以N為周期，即Wk ■W2k，可得N/2NNk■0,1,…，一■14IX(k)■X(k)■W2kNk■0,1,…，一■14■1 3 N4IX(k■N/4)■X(k)■W2kX(k)1 3 N4同樣,對(duì)X2(k)也可進(jìn)行類(lèi)似的分解。一直分解下去,最后是2點(diǎn)的DFT,2點(diǎn)DFT的運(yùn)算也可用碟形符號(hào)來(lái)表示。這樣，對(duì)于一個(gè)N■23■8的DFT運(yùn)算，其按時(shí)間抽取的分解過(guò)程及完整流圖如下圖所示。

x[0x[2x[6x[1x[3x[5x[7x[0x[2x[6x[1x[3x[5x[7這種方法，由于每一步分解都是按輸入序列在時(shí)域上的次序是屬于偶數(shù)還是奇數(shù)來(lái)抽取的，故稱(chēng)為“時(shí)間抽取法”。

分析上面的流圖，N■2M，一共要進(jìn)行M次分解，構(gòu)成了從x(n)到X(k)的M級(jí)運(yùn)算過(guò)程。每一級(jí)運(yùn)算都是由N/2個(gè)蝶形運(yùn)算構(gòu)成，因此每一級(jí)運(yùn)算都需要N/2次復(fù)乘和N次復(fù)加，則按時(shí)間抽取的M級(jí)運(yùn)算后總共需要NN復(fù)數(shù)乘法次數(shù)：mF-2■M■210g2N復(fù)數(shù)加法次數(shù)：〃■N■M■NlogNF2根據(jù)上面的流圖，分析FFT算法的兩個(gè)特點(diǎn)，它們對(duì)FFT的軟硬件構(gòu)成產(chǎn)生很大的影響。(1)原位運(yùn)算也稱(chēng)為同址運(yùn)算，當(dāng)數(shù)據(jù)輸入到存儲(chǔ)器中以后，每一級(jí)運(yùn)算的結(jié)果仍然存儲(chǔ)在原來(lái)的存儲(chǔ)器中，直到最后輸出，中間無(wú)需其它的存儲(chǔ)器。根據(jù)運(yùn)算流圖分析原位運(yùn)算是如何進(jìn)行的。原位運(yùn)算的結(jié)構(gòu)可以節(jié)省存儲(chǔ)單元，降低設(shè)備成本。(2)變址分析運(yùn)算流圖中的輸入輸出序列的順序，輸出按順序，輸入是“碼位倒置”的順序。見(jiàn)圖。碼位倒置的變址處理碼位倒置的變址處理自然順序二進(jìn)制表示碼位倒置碼位倒置順序0000000010011004201001023011110641000011510110156110011371111117碼位倒置順序在實(shí)際運(yùn)算中，直接將輸入數(shù)據(jù)x(n)按碼位倒置的順序排好輸入很不方便，一般總是先按自然順序輸入存儲(chǔ)單元，然后通過(guò)變址運(yùn)算將自然順序的存儲(chǔ)換成碼位倒置順序的存儲(chǔ)，這樣就可以進(jìn)行FFT的原位運(yùn)算。變質(zhì)的功能如圖所示。用軟件實(shí)現(xiàn)是通用采用雷德(Rader)算法，算出I的倒序J以后立即將輸入數(shù)據(jù)X(I)和XJ)對(duì)換。盡管變址運(yùn)算所占運(yùn)算量的比例很小，但對(duì)某些高要求的應(yīng)用(尤其在實(shí)時(shí)信號(hào)處理中)，也可設(shè)法用適當(dāng)?shù)碾娐方Y(jié)構(gòu)直接實(shí)現(xiàn)變址。例如單片數(shù)字信號(hào)處理器TMS320C25就有專(zhuān)用于FFT的二進(jìn)制碼變址模式。4.2按頻率抽取(DIF)的FTT除時(shí)間抽取法外，另外一種普遍使用的FFT結(jié)構(gòu)是頻率抽取法。頻率抽取法將輸入序列不是按奇、偶分組，而是將N點(diǎn)DFT寫(xiě)成前后兩部分：X(k)■DFT[x(n)]■■x(n)WknNnm■(N1/2l，;x(n)Wkn^■x(n)WknNNn，0 n，N/2=.x(n)Wnk■.x(n■N/2)Wa⑵kNNn，0 n，0=NRx(n)■W(N/2)kx(n，N/2)]WnkNNn，0因?yàn)閃N/2，，1,W(N/2)k■(■)k,k為偶數(shù)時(shí)(■)k■1,k為奇數(shù)時(shí)NN(■)k■■,由此可將X(k)分解為偶數(shù)組和奇數(shù)組：X(k尸NRx(n)■(m)kx(n■N/2)]WnNn，0X(2r尸N?[x(n)■x(n■N/2)]W2nrNn，0，N，/2，[1x(n)，x(n，N/2)]WnrN/2n，0

X(2r■1)="Wx(n)■x(n■N/2)]W(2r■)nNnID■NRx(n)■x(n■N/2)WnWnrNN/2nn■(n)■x(n)■x(n■N/2)令.1 n■0,1,■,N/2■1.(n)■[x(n)■x(n■N/2)]Wn2N這兩個(gè)序列都是N/2點(diǎn)的序列，對(duì)應(yīng)的是兩個(gè)N/2點(diǎn)的DFT運(yùn)算:X(2r尸.x(n)W〃1N/2nnX(2r■1)=.x(n)Wrn2N/2nn這樣，同樣是將一個(gè)N點(diǎn)的DFT分解為兩個(gè)N/2點(diǎn)的DFT了。頻率抽選法對(duì)應(yīng)的碟形運(yùn)算關(guān)系圖如下：對(duì)于N=8時(shí)頻率抽取法的FFT流圖如下:WnN(aBb)W對(duì)于N=8時(shí)頻率抽取法的FFT流圖如下:WnN3^-567[([[(XXXXX23^-567[([[(XXXXX2A2A2A2A這種分組的辦法由于每次都是按輸出X(k)在頻域的順序上是屬于偶數(shù)還是奇數(shù)來(lái)分組的，稱(chēng)為頻率抽取法。與前面按時(shí)間抽取的方法相比，相同點(diǎn)問(wèn)題：如何利用快速算法計(jì)算IDFT？分析IDFT的公式：x(n)■IDFT[X(k)]■—■X(k)W~k,n■0,1，…,NHlNNkn比較DFT的公式：X(k)■DFT[x(n)]■■x(n)Wnk,k■0,1，…，N■1Nnn得知可用兩種方法來(lái)實(shí)現(xiàn)IDFT的快速算法：(1)只要把DFT運(yùn)算中的每1一個(gè)系數(shù)Wnk該為W皿，并且最后再乘以常數(shù)人，，就可以用時(shí)間抽取法NN N或頻率抽取的FFT算法來(lái)直接計(jì)算IDFT。這種方法需要對(duì)FFT的程序和參數(shù)稍加改動(dòng)才能實(shí)現(xiàn)。(2 )因?yàn)閤(n)■1-N1X*(k)Wnk]*■—{DFT[X*(k)]},n■0,1,■,N■1，也就N NNk是說(shuō)，可先將X(k)取共軛變換，即將X(k)的虛部乘以一i,就可直接調(diào)用FFT的程序，最后再對(duì)運(yùn)算結(jié)果取一次共軛變換并乘以常數(shù)1/N即可得到x(n)的值。這種方法中，F(xiàn)FT運(yùn)算和IFFT運(yùn)算都可以共用一個(gè)子程序塊，在使用通用計(jì)算機(jī)或用硬件實(shí)現(xiàn)時(shí)比較方便。4.1.3混合基FFT算法以上討論的是基2的FFT算法，即N■2m的情況，這種情況實(shí)際上使用得最多，這種FFT運(yùn)算，程序簡(jiǎn)單，效率很高，用起來(lái)很方便。另外，在實(shí)際應(yīng)用時(shí)，有限長(zhǎng)序列的長(zhǎng)度N到底是多少在很大程度時(shí)是由人為因素確定的，因此，大多數(shù)場(chǎng)合人們可以將N選定為N■2m，從而可以直接調(diào)用以2為基數(shù)的FFT運(yùn)算程序。如果長(zhǎng)度N不能認(rèn)為確定，而N的數(shù)值又不是以2為基數(shù)的整數(shù)次方，一般可有以下兩種處理方法：將x(n)用補(bǔ)零的方法延長(zhǎng)，使N增長(zhǎng)到最鄰近的一個(gè)N■2m數(shù)值。例如，N=30,可以在序列x(n)中補(bǔ)進(jìn)x(30)=x(31)=O兩個(gè)零值點(diǎn)，使N=32。如果計(jì)算FFT的目的是為了了解整個(gè)頻譜，而不是特定頻率點(diǎn)，則此法可行。因?yàn)橛邢揲L(zhǎng)序列補(bǔ)零以后并不影響其頻譜X(ejw),只是頻譜的采樣點(diǎn)數(shù)增加而已。如果要求特定頻率點(diǎn)的頻譜，則N不能改變。如果N為復(fù)合數(shù)，則可以用以任意數(shù)為基數(shù)的FFT算法來(lái)計(jì)算?？焖俑道锶~變換的基本思想就是要將DFT的運(yùn)算盡量分小。例如，N=6時(shí)，可以按照N=3X2分解，將6點(diǎn)DFT分解為3組2點(diǎn)DFT。舉例：N=9時(shí)的快速算法。4.2快速傅里葉變換的應(yīng)用凡是可以利用傅里葉變換來(lái)進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法及運(yùn)用數(shù)字計(jì)算技術(shù)加以實(shí)現(xiàn)。FFT在數(shù)字通信、語(yǔ)音信號(hào)處理、圖像處理、匹配濾波以及功率譜估計(jì)、仿真、系統(tǒng)分析等各個(gè)領(lǐng)域都得到了廣泛的應(yīng)用。但不管FFT在哪里應(yīng)用，一般都以卷積積分或相關(guān)積分的具體處理為依據(jù)，或者以用FFT作為連續(xù)傅里葉變換的近似為基礎(chǔ)。利用FFT求線(xiàn)性卷積一快速卷積在實(shí)際中常常遇到要求兩個(gè)序列的線(xiàn)性卷積。如一個(gè)信號(hào)序列x(n)通過(guò)FIR濾波器時(shí)，其輸出y(n)應(yīng)是x(n)與h(n)的卷積：y(n)■x(n)■h(n)■■■x(m)h(n■m)m■■■有限長(zhǎng)序列x(n)與h(n)的卷積的結(jié)果y(n)也是一個(gè)有限長(zhǎng)序列。假設(shè)x(n)與h(n)的長(zhǎng)度分別為N1和N2,則y(n)的長(zhǎng)度為N=N1+N2-1。若通過(guò)補(bǔ)零使x(n)與h(n)都加長(zhǎng)到N點(diǎn)，就可以用圓周卷積來(lái)計(jì)算線(xiàn)性卷積。這樣得到用FFT運(yùn)算來(lái)求y(n)值(快速卷積)的步驟如下：(1)對(duì)序列x(n)與h(n)補(bǔ)零至長(zhǎng)為N，使N三N1+N2-1，并且N■2M(M為整數(shù))，即/(n),n■0,1，…,N1■1x(n)■■取n■N1,N1■1,….,N■1*(n),n■0,1，…，N2■1h(n)■■ …取n■N2,N2?1,…,N■1(2)用FFT計(jì)算x(n)與h(n)的離散傅里葉變換x(n)■(FFT)■X(k) (N點(diǎn))h(n)■(FFT)■H(k) (N點(diǎn))(3)計(jì)算Y(k)=X(k)H(k)(4)用IFFT計(jì)算Y(k)的離散傅里葉反變換得:y(n)=IFFT[Y(k)](N點(diǎn))利用FFT求相關(guān)一快速相關(guān)互相關(guān)及自相關(guān)的運(yùn)算已廣泛的應(yīng)用于信號(hào)分析與統(tǒng)計(jì)分析，應(yīng)用于連續(xù)時(shí)間系統(tǒng)也用于離散事件系統(tǒng)。用FFT計(jì)算相關(guān)函數(shù)稱(chēng)為快速相關(guān)，它與快速卷積完全類(lèi)似，不同的是一個(gè)應(yīng)用離散相關(guān)定理，另一個(gè)應(yīng)用離散卷積定理。同樣都要注意到離散傅里葉變換固有的周期性，也同樣用補(bǔ)零的方法來(lái)繞過(guò)這個(gè)障礙。設(shè)兩個(gè)離散時(shí)間信號(hào)x(n)與丫6)為已知，離散互相關(guān)函數(shù)記作R(n)，xy定義為R(n)■■■x(m)y(n■m)xym■■■如果x(n)與y(n)的序列長(zhǎng)度分別為N1和N2,則用FFT求相關(guān)的計(jì)算步驟如下：(1)對(duì)序列x(n)與y(n)補(bǔ)零至長(zhǎng)為N,使N三N1+N2-1,并且N■2M(M為整數(shù))，即/(n),n■0,1,…,N1■1x(n)■■取n■N1,N1?1,…,N■1/(n),n■0,1，…,N2■1y(n)■■取n■N2,N2?1,…,N■1(2)用FFT計(jì)算x(n)與y(n)的離散傅里葉變換x(n)■(FFT)■X(k) (N點(diǎn))y(n)■(FFT)■y(k) (N點(diǎn))(3)將X(k)的虛部Im[X(k)]改變符號(hào)，求得其共軛X