詳解FFT快速傅里葉變換FFT

1、第四章快速傅里葉變換有限長(zhǎng)序列可以通過(guò)離散傅里葉變換(DFT) 將其頻域也離散化成有限長(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 的多種算法，基本算法是基2 DIT和基2 DIF。FFT在離散傅里葉反變換、線性卷積和線性相關(guān)等方面也有重要應(yīng)用?？焖俑道锶~變換(FFT

2、)是計(jì)算離散傅里葉變換(DFT)的快速算法。DFT 的定義式為N- 1knX (k )=匯 x(n)W rn (k)n=0kn在所有復(fù)指數(shù)值 WN的值全部已算好的情況下，要計(jì)算一個(gè)X (k )需要N2次復(fù)數(shù)乘法和N 1 次復(fù)數(shù)加法。算出全部N 點(diǎn) X (k ) 共需 N 次復(fù)數(shù)乘法2和 N ( N - 1) 次復(fù)數(shù)加法。即計(jì)算量是與N 2 成正比的。FFT 的基本思想：將大點(diǎn)數(shù)的DFT 分解為若干個(gè)小點(diǎn)數(shù)DFT 的組合，從而減少運(yùn)算量。WN 因子具有以下兩個(gè)特性，可使 DFT 運(yùn)算量盡量分解為小點(diǎn)數(shù)的DFT運(yùn)算：=Wkn = W (n+N )k( 1) 周期性：W (k+N )n N NN(

3、2) 對(duì)稱(chēng)性：W (k+N / 2) = - W kNN利用這兩個(gè)性質(zhì)，可以使 DFT 運(yùn)算中有些項(xiàng)合并，以減少乘法次數(shù)。例子：求當(dāng)n = 4時(shí)，X(2)的值3匯 44444X (2) =x(n)W2n = x(0)W 0 + x(1W 2 + x(2)W 4 + x(3)W 6n=0= x(0) +x(2)W 0 +x(1) + x(3)W2(周期性)44=x(0) + x(2)-x(1) + x(3)W(4(對(duì)稱(chēng)性)通過(guò)合并，使乘法次數(shù)由4次減少到1次，運(yùn)算量減少。FFT的算法形式有很多種，但基本上可以分為兩大類(lèi)：按時(shí)間抽取(DIT)和按頻率抽取(DIF)。4.1按時(shí)間抽取(DIT )的F

4、TT為了將大點(diǎn)數(shù)的 DFT分解為小點(diǎn)數(shù)的復(fù)合數(shù)，最常用的是 N = 2 M的情況(M基2FFT。下面討論基2情況的算法。先將序列x(n)按奇偶項(xiàng)分解為兩組?x(2r ) = x(r)DFT運(yùn)算，要求序列的長(zhǎng)度N為為正整數(shù))。該情況下的變換稱(chēng)為?N?x(2r +1) = x2(r )r = 0，1L，2- 1將DFT運(yùn)算也相應(yīng)分為兩組N - 1Xz knX (k ) = DFT x(n) =x(n)WNn=0N- 1N- 1匯 x(n)Wkn 匯 x(n)Wkn+n=0n為偶數(shù)n=0n為奇數(shù)N / 2- 1E (2 )2rk二 xr=0r WnN / 2- 1E(2 x r=0N / 2- 12

5、rkN / 2- 1+ 1)(2 r+1 )kWn2rkx2(r)WN=Xi (r)WNWnr=0N / 2- 1k N / 2- 1Nrk/ 2 +WNrk (因?yàn)?W 2 rk W rk )=2_j X12_J X2N / 2WN - WN / 2)(r)W(r)Wr =0r=0= X 1(k ) + WNkX 2(k )其中 Xi(k)、X 2(k)分別是 Xi(n)、X2 (n)的 N/2 點(diǎn)的 dftN / 2- 1N / 2- 1 rkrkXi(k )= E Xi(r)WN / 2 = E x(2r)WN / 2 ,0 <k <- 1r=0r=0N / 2- 1N /

6、2- 1rkrkX 2(k)= E X2 (r)WN / 2 = E x(2r +1)Wn / 2 ,0 <k<1r=0r=0至此，一個(gè)N點(diǎn)DFT被分解為兩個(gè)N/2 點(diǎn)的 DFT。上面是否將全部 N點(diǎn)的X (k )求解出來(lái)了 ?分析：X1(k )和JN- 1)，則由X 2 (k )只有 N/2 個(gè)點(diǎn)(k = 0,1,L ,2X (k) = X 1 (k ) +W ：X 2(k )只能求出X (k )的前N/2個(gè)點(diǎn)的DFT,要求出全部N點(diǎn)的X (k )，需要找出X(k )、X 2 (k )和X (k + N / 2)的關(guān)系，其kN(k ) +WNX2(k )可得中 k = 0,1L

7、, - 1。由式子 X (k ) = X 1X (k + N / 2) = X (k + N / 2) +W k+N / 2 X 2(k + N / 2)化簡(jiǎn)得1NkX (k + N / 2) = = X 1(k ) - Wn X 2 (k ) ， k = 0,1,L這樣N點(diǎn)DFT可全部由下式確定出來(lái):kk = 0,1L N2*)?X (k) = X 1 (k)+WN X 2(k)?k?X (k + N / 2) = X 1(k) - W X2 (k)上式可用一個(gè)專(zhuān)用的碟形符號(hào)來(lái)表示，這個(gè)符號(hào)對(duì)應(yīng)一次復(fù)乘和兩次復(fù)加運(yùn)N圖蝶形運(yùn)算符號(hào)通過(guò)這樣的分解以后,每一個(gè)N/2點(diǎn)的DFT只需要(N)22N2

8、次復(fù)數(shù)乘法，兩個(gè)N/2點(diǎn)的DFT需要2(N ) 2 = N一次復(fù)乘，再加上將兩個(gè) N/2點(diǎn)22DFT合并成 為N點(diǎn)DFT時(shí)有 N /2次與 W 因子相乘，一共需要N 2一次復(fù)乘?？梢?jiàn)，通過(guò)這樣的分解，運(yùn)算量節(jié)省了近一半。2因?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ì)xi(r),可以在按其偶數(shù)部分及奇數(shù)部分進(jìn)行分解:?Xi(2l ) = X3 (l )?Xi (2l +1) = X4(l )則的運(yùn)算可相應(yīng)分為兩組：1=0,1L，4N / 4- 121kX1(k )=三 X1 (21)WN /

9、2l=0N / 4- 1+ E xQ1W1=0(2l+1)kN / 2N / 4- 1E 3()lk/ 4x 1 Wn X 1=0N / 4 1E 4()x l Wl=0lkN / 4一k將系數(shù)統(tǒng)一為以N為周期，即 W N / 2k=X 3(k ) +Wn /2 X 4(k )k = 0,1L ,4 - 1=W 2 k,可得N2kk = 0,1L-N - 1?X1 (k) = X 3 (k) +Wn X 4 (k)22kXi (k + N / 4) = X3 (k) - Wn X 4 (k)同樣，對(duì)X2 (k)也可進(jìn)行類(lèi)似的分解。一直分解下去，最后是2點(diǎn)的DFT , 2點(diǎn)DFT的運(yùn)算也可用碟形

10、符號(hào)來(lái)表示。這樣，對(duì)于一個(gè)N = 23 = 8DFT 運(yùn)算，其按時(shí)間抽取的分解過(guò)程及完整流圖如下圖所示。X0)KX1*oX2聲。X3&ft4點(diǎn)T/小DFT-*ft.F0F 號(hào) X(7)4點(diǎn)OFT0)制|)印 引力 xlHMSTHxfRxf01234567-Il fl -I ( ( ( -Il -Il xxxxxxxx1J 1J. u 1J 0 4 2 6X X X XJ u 口出p R X X X Xx0)oxnb1日點(diǎn)DFTX(0 ill 111 JuX 2n°xnjn V (9-1I1-J x310°xkJo 乂 fOlx 4。r 41x|51o0A 1*1rt

11、W 不 1| "JI UxfGlo°XPIft VJF £!xQo“x網(wǎng)(71x0<> 乂網(wǎng)。*團(tuán) 乂呵 xl) 乂冏 x3) x(7)”XD Xl) XRl X網(wǎng)X同 ?X5 小X同 M兇7 -1這種方法，由于每一步分解都是按輸入序列在時(shí)域上的次序是屬于偶數(shù) 還是奇數(shù)來(lái)抽取的，故稱(chēng)為“時(shí)間抽取法”分析上面的流圖， N = 2 M , 一共要進(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)算后總共需 要復(fù)數(shù)乘法次數(shù)：mF = N ?

12、= N log n M 222復(fù)數(shù)加法次數(shù)：aF = N ? M = N 10g2 N根據(jù)上面的流圖，分析F FT算法的兩個(gè)特點(diǎn)，它們對(duì)F FT的軟硬件構(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)制表小碼位倒置碼位倒置順序00000000100110042010010230111106

13、41000011510110156110011371111117碼位倒置順序X(0) X(4) X(2)X(6) X(1) X(5)X(3)X(7)碼位倒置的變址處理在實(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)和X(J)對(duì)換。盡管變址運(yùn)算所占運(yùn)算量的比例很小，但對(duì)某些高要求的應(yīng)用(尤其在實(shí)時(shí)信號(hào)處理中)，也可設(shè)法用適當(dāng)?shù)碾娐方Y(jié)構(gòu)直接實(shí)現(xiàn)變址。例

14、如單片數(shù)字信號(hào)處理器TMS320C25 就有專(zhuān)用于 FFT 的二進(jìn)制碼變址模式。4. 2 按頻率抽取(DIF )的FTT除時(shí)間抽取法外，另外一種普遍使用的FFT 結(jié)構(gòu)是頻率抽取法。頻率抽取法將輸入序列不是按奇、偶分組，而是將N點(diǎn)DFT寫(xiě)成前后兩部分：N -1knx(n)WNn =0N -1( N / 2)-1匯 x(n)Wkn+ 匯 x(n)Wknn=0Nn=N / 2N / 2- 1( )N / 2- 1nk +(n+N / 2)kE ( +/ 2)x n WN x n NWNn=0n=0nkN / 2- 1(N / 2)k x(n + N / 2)WN=E x(n) +Wnn =0W N

15、/ 2 = - 1,W( N / 2) k = (- 1)k ， k 為 偶 數(shù)時(shí) (- 1)k =1 ， k 為奇數(shù)時(shí) NN(- 1)k = - 1 ，由此可將X(k) 分解為偶數(shù)組和奇數(shù)組：N / 2- 1X (k)= n=0 x(n) + (- 1)k x(n + N / 2)WnkN / 2- 1n=0X (2r)=N / 2- 1Ex(n) + x(n + N / 2)W 2nrx(n) + x(n + N / 2)W n=0N / 2 nrN / 2-1X (2r +1>n=)x(n) - x(n + N / 2)W(2r+1)nN / 2- 1x(n) - x(n + N

16、/ 2)W nWnr n=0n = 0,1L , N / 2 - 1N?x1(n) = x(n) + x(n + N / 2) 令?X2 (n) = x(n) - x(n + N / 2)W n這兩個(gè)序列都是 N/2點(diǎn)的序列，對(duì)應(yīng)的是兩個(gè) N/2點(diǎn)的DFT運(yùn)算:N / 2- 1二1N / 2n=0nrX (2r) x (n)WN / 2- 1k 2 N / 2rnX (2r+1)= n=0 x (n)W這樣，同樣是將一個(gè)N點(diǎn)的DFT分解為兩個(gè)N/2點(diǎn)的DFT 了。頻率抽選法 對(duì)應(yīng)的碟形運(yùn)算關(guān)系圖如下：-1Wn對(duì)于N=8時(shí)頻率抽取法的 FFT流圖如下:x(Dj0xnio8 點(diǎn) DFTr V fQ

17、lI?！眡3L"X ZrtVniX門(mén)q 八 1JJnV Mlx51oX°x網(wǎng)乂C(7)X X(Q)*。X|2) xM) x|6| XU)a。X|3)"X(5)x(0). x(1)o 刈4 x(30x(5)0 x(6)0 孫2點(diǎn)DFT2點(diǎn)DFT2占DFTX(0) or口>o2占 匚八% DFTX X(7)這種分組的辦法由于每次都是按輸出X(k)在頻域的順序上是屬于偶數(shù)還是奇數(shù)來(lái)分組的，稱(chēng)為頻率抽取法。與前面按時(shí)間抽取的方法相比，相同點(diǎn)問(wèn)題：如何利用快速算法計(jì)算IDFT ?分析IDFT的公式：N - 1X=IDFTX (k) 二三 X (k)WN-nk,n =

18、0,1L ,N -N k=o比較DFT的公式:N - 1X (k) = DFTx(n) = £； x(n)WNnk ,k = 0,1L , N - 1n=0得知可用兩種方法來(lái)實(shí)現(xiàn) IDFT的快速算法：(1)只要把 DFT運(yùn)算中的每個(gè)系數(shù)Wnk該為W-nk,并且最后再乘以常數(shù)NN1、，一，、一.一，就可以用時(shí)間抽取法N或頻率抽取的FFT算法來(lái)直接計(jì)算IDFT。這種方法需要對(duì) FFT的程序和參數(shù)稍加 改動(dòng)才能實(shí)現(xiàn)。 (2) 因 為/、_ 1 N-11*x(n) # X * (k)W/nk * >DFTX (k), n = 0,1L ,N - 1,也就N k=0N是說(shuō)，可先將 X(k

19、)取共軻變換，即將 X(k)的虛部乘以1 ,就可直接調(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 = 2 M的情況，這種情況實(shí)際上使用得最多，這種 FFT運(yùn)算，程序簡(jiǎn)單，效率很高，用起來(lái)很方便。另 外，在實(shí)際應(yīng)用時(shí)，有限長(zhǎng)序列的長(zhǎng)度N到底是多少在很大程度時(shí)是由人為因素確定的，因此，大多數(shù)場(chǎng)合人們可以將N選定為N = 2 M ,從而可以直接調(diào)用以2為基數(shù)的 FFT運(yùn)算程序。如果長(zhǎng)度N不能認(rèn)為確

20、定，而 N的數(shù)值又不是以2為基數(shù)的整數(shù)次 方，一般可有以下兩種處理方法：(1) 將x(n)用補(bǔ)零的方法延長(zhǎng)，使 N增長(zhǎng)到最鄰近的一個(gè)N = 2 M數(shù)值。仞0口， N=30,可以在序列 x(n)中補(bǔ)進(jìn)x(30) = x(31) = 0兩個(gè)零值點(diǎn)， 使N=32。如果af算 FFT 的目的是為了了解整個(gè)頻譜， 而不是特定頻率點(diǎn)， 則此 法可行。因?yàn)橛邢揲L(zhǎng)序列補(bǔ)零以后并不影響其頻譜X (ejw ),只是頻譜的采樣點(diǎn)數(shù)增加而已。(2) 如果要求特定頻率點(diǎn)的頻譜，則 N不能改變。如果 N為復(fù)合數(shù)，則可以用以任意數(shù)為基數(shù)的FFT算法來(lái)計(jì)DFT的運(yùn)算N=3X 2分解，算?？焖俑道锶~變換的基本思想就是要將 盡量

21、分小。例如，N=6時(shí)，可以按照DFT。將6點(diǎn)DFT分解為3組2點(diǎn)舉例： 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ǔ)。4.2.1 利用 FFT 求線性卷積快速卷積在實(shí)際中常常遇到要求兩個(gè)序列的線性卷積。如一個(gè)信號(hào)序列x(n)通過(guò)FIR濾波器時(shí)，其輸出y(n)

22、應(yīng)是x(n)與h(n)的卷積：y(n) = x(n) ? h(n) =x(m)h(n - 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ì)算線性卷積。這樣得到用FFT 運(yùn)算來(lái)求y(n) 值(快速卷積)的步驟如下：(1) 對(duì)序列 x(n)與h(n)補(bǔ)零至長(zhǎng)為 N，使 N > N1+N2-1 ,并且 N =皆(M為整數(shù))，即?x(n), n = 01,L , N1 - 1 x(n) = ?0, n = N1, N

23、1 +1,L , N - 1?h(n), n = 01,L , N 2 - 1 h(n) = ?0, n = N2, N 2 +1,L , 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)4.2 快速傅里葉變換的應(yīng)用(3 ) 計(jì)算 Y(k)=X(k)H(k)(4) 用IFFT計(jì)算Y(k)的離散傅里葉反變換得：y(n)=IFFTY(k) ( N點(diǎn))4.2.2 利用 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)。

24、用 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)與y(n)為已知，離散互相關(guān)函數(shù)記作Rxy (n),定義為Rxy (n)=匯 x(m)y(n + m)m=- °°如果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 =2" (M為整數(shù))，即?x(n), n = 01,L , N1 - 1 x(n) = ?0, n = N1, N1 +1,L , N - 1?y(n), n =01,L , N 2 - 1 y(n) = ?0, n = N2, N 2 +1,L , N - 1(2 )用FFT計(jì)算x(n)與y(n)的離散傅里葉變換x(n) ? (FFT ) ?X (k )( N 點(diǎn))y(n) ? ( FF