基2FFT算法分析的介紹

1、快速傅里葉變換有限長序列可以通過離散傅里葉變換(DFT)將其頻域也離散化成有限長序列.但其計算量太大，很難實時地處理問題，因此引出了快速傅里葉變換(FFT).1965年，Cooley和Tukey提出了計算離散傅里葉變換(DFT)的快速算法，將DFT的運(yùn)算量減少了幾個數(shù)量級。從此，對快速傅里葉變換(FFT)算法的研究便不斷深入，數(shù)字信號處理這門新興學(xué)科也隨FFT的出現(xiàn)和發(fā)展而迅速發(fā)展。根據(jù)對序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，基本算法是基2DIT和基2DIF。FFT在離散傅里葉反變換、線性卷積和線性相關(guān)等方面也有重要應(yīng)用?？焖俑道锶~變換(FFT)是計算離散傅里葉變換(DFT)的快速

2、算法。DFT的定義式為NknX(k)=xx(n)WNRN(k)n=0在所有復(fù)指數(shù)值W：n的值全部已算好的情況下，要計算一個X(k)需要N次復(fù)數(shù)乘法和N1次復(fù)數(shù)加法。算出全部N點(diǎn)X(k)共需N2次復(fù)數(shù)乘法和N(N-1)次復(fù)數(shù)加法。即計算量是與N2成正比的。FFT的基本思想：將大點(diǎn)數(shù)的DFT分解為若干個小點(diǎn)數(shù)DFT的組合，從而減少運(yùn)算量。Wn因子具有以下兩個特性，可使DFT運(yùn)算量盡量分解為小點(diǎn)數(shù)的DFT運(yùn)算：(1)周期性：W:咻1)n=W；n=W：nH)k對稱性：wFn/2)_W：利用這兩個性質(zhì)，可以使DFT運(yùn)算中有些項合并，以減少乘法次數(shù)。例子:求當(dāng)N=4時，X(2)的值3X(2)=x(n)W：

3、n=x(0)W40xW42x(2)w44x(3)W：nz0=x(0)x(2)W40x(1)x(3)W42(周期性)=x(0)+x(2)x(1)+x(3)W：(對稱性)通過合并，使乘法次數(shù)由4次減少到1次，運(yùn)算量減少。FFT的算法形式有很多種，但基本上可以分為兩大類：按時間抽取(DIT)和按頻率抽取(DIF)。4.1按時間抽取(DIT)的FTT為了將大點(diǎn)數(shù)的DFT分解為小點(diǎn)數(shù)的DFT運(yùn)算，要求序列的長度N為復(fù)合數(shù)，最常用的是N=2M的情況(M為正整數(shù))。該情況下的變換稱為基2FFT。下面討論基2情況的算法。先將序列x(n)按奇偶項分解為兩組,x(2r)=x1(r)x(2r+1)=x2(r)N川,

4、丁將DFT運(yùn)算也相應(yīng)分為兩組N-1X(k)=DFTx(n)=、'x(n)W；nn=0N4N4=x(n)W；n+、x(n)W；nn=0n=0n為偶數(shù)n為奇數(shù)N/24N/24=''x(2r)W：k'x(2r1)WN2r伊r=0r=0N/2N/22rkk2rk=,、Xi(r)WnWn-X2()Wnr0r0N/2N/2=ZXi(r)wNk/2十w；工X2(r)wNk/2(因為W；rk=W：k/2)rOr4=Xi(k)W；X2(k)其中X1(k)、X2(k)分別是x1(n)、x2(n)的N/2點(diǎn)的dftN/2N/2rk一rkNX1(k)=ZXi(r)WN/2=£

5、x(2)Wn/2,0<k<-1r0r02N/2JN/2rk一rkNX2(k)=X2(r)WN/2=X(2r1)WN/2,0<k_-1r=0r=02至此，一個N點(diǎn)DFT被分解為兩個N/2點(diǎn)的DFT。上面是否將全部N點(diǎn)的X(k)求解出來了？N分析：X(k)和X2(k)只有N/2個點(diǎn)(k=0,1，一一1),則由2kX(k)=X1(k)+WnX2«)只能求出X(k)的前N/2個點(diǎn)的DFT,要求出全部N點(diǎn)的X(k),需要找出X1(k)、X2(k)和X(k+N/2)的關(guān)系，其中k=0,1，£1。由式子X(k)=X1(k)+W；Xz(k)可得X(k+N/2)=X1(k+

6、N/2)+W；氏/2X2(k+N/2)化簡得kNX(k+N/2)=X1(k)WNX2(k),k=0,1:，一12這樣N點(diǎn)DFT可全部由下式確定出來：X(k)=X1(k)+W；Xz(k)N.kk=0,1,-1(*)X(k+N/2)=X1(k)WNX2(k)2上式可用一個專用的碟形符號來表示，這個符號對應(yīng)一次復(fù)乘和兩次復(fù)加運(yùn)a-.j-bWNkaW；ba-W；b圖蝶形運(yùn)算符號2N2N2通過這樣的分解以后，每一個N/2點(diǎn)的DFT只需要(一)2=次復(fù)數(shù)乘24NoN法，兩個N/2點(diǎn)的DFT需要2()2=次復(fù)乘，再加上將兩個N/2點(diǎn)22N2DFT合并成為N點(diǎn)DFT時有N/2次與W因子相乘，一共需要NN2一,

7、+也次復(fù)乘?？梢姡ㄟ^這樣的分解，運(yùn)算量節(jié)省了近一半。22因為N=2M,N/2仍然是偶數(shù)，因此可以對兩個N/2點(diǎn)的DFT再分別作進(jìn)一步的分解，將兩個N/2點(diǎn)的DFT分解成兩個N/4點(diǎn)的DFT。例如對x1(r),可以在按其偶數(shù)部分及奇數(shù)部分進(jìn)行分解:xi(2l)=X3(l)_NI。1,一內(nèi)(2l+1)=X4(l)4則的運(yùn)算可相應(yīng)分為兩組：N/44Xi(k)=Zxi(2l)WN2l/k2l=0N/44xi(2l1)wN2l21)kl-0N/44N/44=x3(l)Wn/4WN/2x4(l)WN/4l=0l=0kN=X3(k)W：/2X4(k)k=0,1,14將系數(shù)統(tǒng)一為以N為周期，即W；/2=W；

8、k,可得,Xi(k)=X3(k)+WnX4(k)N.2kk=0,1,一-1Xi(kN/4)=X3(k)-Wr2kX4(k)4同樣，對X2（k）也可進(jìn)行類似的分解。一直分解下去，最后是2點(diǎn)的3一DFT,2點(diǎn)DFT的運(yùn)算也可用碟形符號來表不。這樣，對于一個N=2=8產(chǎn)*Mh.*1-出網(wǎng)化enp(3Hx(x(x(xlx(x(xX31)才司可可劃刀xfxl(xxlx(xlx1x(的DFT運(yùn)算，其按時間抽取的分解過程及完整流圖如下圖所示。x0bxHL8點(diǎn)DFTX(011Jux(?lQVroiQT口x30ftVfQlijy1xMl門OVfjfx510門f£lxMl»口x網(wǎng)沏°

9、X(G)°X(7)0)可©s'1)5)國71flflttfl-_l_tflXXXXXXXX*T*T，*TXXXXXXXX這種方法，由于每一步分解都是按輸入序列在時域上的次序是屬于偶數(shù)還是奇數(shù)來抽取的，故稱為“時間抽取法”。分析上面的流圖，N=2M,一共要進(jìn)行M次分解，構(gòu)成了從x(n)到X(k)的M級運(yùn)算過程。每一級運(yùn)算都是由N/2個蝶形運(yùn)算構(gòu)成，因此每一級運(yùn)算都需要N/2次復(fù)乘和N次復(fù)加，則按時間抽取的M級運(yùn)算后總共需要復(fù)數(shù)乘法次數(shù)：mF=M=；log2N復(fù)數(shù)加法次數(shù)：aF=NM=Nlog2N根據(jù)上面的流圖，分析FFT算法的兩個特點(diǎn)，它們對FFT的軟硬件構(gòu)成產(chǎn)生很大

10、的影響。(1)原位運(yùn)算也稱為同址運(yùn)算，當(dāng)數(shù)據(jù)輸入到存儲器中以后，每一級運(yùn)算的結(jié)果仍然存儲在原來的存儲器中，直到最后輸出，中間無需其它的存儲器。根據(jù)運(yùn)算流圖分析原位運(yùn)算是如何進(jìn)行的。原位運(yùn)算的結(jié)構(gòu)可以節(jié)省存儲單元，降低設(shè)備成本。(2)變址分析運(yùn)算流圖中的輸入輸出序列的順序，輸出按順序，輸入是“碼位倒置”的順序。見圖。自然順序一進(jìn)制表小碼位倒置碼位倒置順序0000000010011004201001023011110641000011510110156110011371111117X(0)123456X(4)X(2)X(6)X(1)X(5)X(3)碼位倒置的變址處理X(7)碼位倒置順序在實際運(yùn)算中

11、，直接將輸入數(shù)據(jù)x(n)按碼位倒置的順序排好輸入很不方便，一般總是先按自然順序輸入存儲單元，然后通過變址運(yùn)算將自然順序的存儲換成碼位倒置順序的存儲，這樣就可以進(jìn)行FFT的原位運(yùn)算。變質(zhì)的功能如圖所示。用軟件實現(xiàn)是通用采用雷德(Rader)算法，算出I的倒序J以后立即將輸入數(shù)據(jù)X(I)和X(J)對換。盡管變址運(yùn)算所占運(yùn)算量的比例很小，但對某些高要求的應(yīng)用(尤其在實時信號處理中)，也可設(shè)法用適當(dāng)?shù)碾娐方Y(jié)構(gòu)直接實現(xiàn)變址。例如單片數(shù)字信號處理器TMS320C25就有專用于FFT的二進(jìn)制碼變址模式。4.2按頻率抽?。―IF）的FTT除時間抽取法外，另外一種普遍使用的FFT結(jié)構(gòu)是頻率抽取法。頻率抽取法將輸

12、入序列不是按奇、偶分組，而是將N點(diǎn)DFT寫成前后兩部分：N1knX(k)=DFTx(n)='、x(n)WNn=0(N/2)1N1=x(n)W；n+x(n)W；nnHn-N/2N/2N/2='、x(n)WN1k、x(nN/2)WNnN/2)kn=0n=SN/2="x(n)WNN/2)kx(nN/2)WN1kn=e因為W；/2=1,WNN/2)k=(1)k,k為偶數(shù)時(1)k=1,k為奇數(shù)時(-1)k=-1,由此可將X(k)分解為偶數(shù)組和奇數(shù)組：N/24X(k)=£x(n)+(-1)kx(n+N/2)W:n=0N/2X(2r)=vx(n)x(nN/2)W；nrn

13、=0N/24='x(n)x(nN/2)WNn；2n=0N/24X(2r1)="x(n)x(nN/2)WN2r1)nn=0N/24=、x(n)-x(nN/2)W；W；/2nz00,1,N/2-1x1(n)=x(n)+x(n+N/2)令<nn=x2(n)=x(n)-x(nN/2)Wn這兩個序列都是N/2點(diǎn)的序列，對應(yīng)的是兩個N/2點(diǎn)的DFT運(yùn)算:N/2JX(2r尸'、Xi(n)W%n衛(wèi)N/2dX(2r1)='、X2(n)W；/2n=0這樣，同樣是將一個N點(diǎn)的DFT分解為兩個N/2點(diǎn)的DFT。頻率抽選法對應(yīng)的碟形運(yùn)算關(guān)系圖如下:<X(0XlUJx1)xZ

14、)X(1)XR)x3oxMloMWxX團(tuán)x(3)oxffl*x5)ox(6)oM小8點(diǎn)DFT網(wǎng)oX(4X(5)X(6>4點(diǎn)DFT4占DFTX0)X2)AX*oX6)+。Xl*。X3)zX5>ox(7)X0)XIX(2)oX7)X6)x(05xlh劃x30x(5)ox6)0刈J。這種分組的辦法由于每次都是按輸出X(k)在頻域的順序上是屬于偶數(shù)還是奇數(shù)來分組的，稱為頻率抽取法。與前面按時間抽取的方法相比，相同點(diǎn)問題：如何利用快速算法計算IDFT?分析IDFT的公式:1N1x(n)=IDFTX(k)二一'、X(k)WNTn=0,1,N-1N百比較DFT的公式：N-4X(k)=DF

15、Tx(n)='x(n)wNnk,k=0,1,N-1n=0得知可用兩種方法來實現(xiàn)IDFT的快速算法：(1)只要把DFT運(yùn)算中的每一個系數(shù)W:該為Wn,并且最后再乘以常數(shù),就可以用時間抽取法N或頻率抽取的FFT算法來直接計算IDFT。這種方法需要對FFT的程序和參數(shù)稍加改動才能實現(xiàn)。(2)因為1*nk*1._*、x(n)=£X(k)WN=DFTX(k),n=0,1,，N-1,也就Nk衛(wèi)N是說，可先將X(k)取共軻變換，即將X(k)的虛部乘以1,就可直接調(diào)用FFT的程序，最后再對運(yùn)算結(jié)果取一次共軻變換并乘以常數(shù)1/N即可得到x(n)的值。這種方法中，F(xiàn)FT運(yùn)算和IFFT運(yùn)算都可以共

16、用一個子程序塊，在使用通用計算機(jī)或用硬件實現(xiàn)時比較方便。4.1.3混合基FFT算法以上討論的是基2的FFT算法，即N=2M的情況，這種情況實際上使用得最多，這種FFT運(yùn)算，程序簡單，效率很高，用起來很方便。另外，在實際應(yīng)用時，有限長序列的長度N到底是多少在很大程度時是由人為因素確定的，因此，大多數(shù)場合人們可以將N選定為N=2M,從而可以直接調(diào)用以2為基數(shù)的FFT運(yùn)算程序。如果長度N不能認(rèn)為確定，而N的數(shù)值又不是以2為基數(shù)的整數(shù)次方，一般可有以下兩種處理方法：(1)將x(n)用補(bǔ)零的方法延長，使N增長到最鄰近的一個N=2M數(shù)值。例如，N=30,可以在序列x(n)中補(bǔ)進(jìn)x(30)=x(31)=0兩

17、個零值點(diǎn)，使N=32。如果計算FFT的目的是為了了解整個頻譜，而不是特定頻率點(diǎn)，則此法可行。因為有限長序列補(bǔ)零以后并不影響其頻譜X(ejw),只是頻譜的采樣點(diǎn)數(shù)增加而已。(2)如果要求特定頻率點(diǎn)的頻譜，則N不能改變。如果N為復(fù)合數(shù)，則可以用以任意數(shù)為基數(shù)的FFT算法來計算?？焖俑道锶~變換的基本思想就是要將DFT的運(yùn)算盡量分小。例如，N=6時，可以按照N=3X2分解，將6點(diǎn)DFT分解為3組2點(diǎn)DFT。舉例：N=9時的快速算法。4.2快速傅里葉變換的應(yīng)用凡是可以利用傅里葉變換來進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法及運(yùn)用數(shù)字計算技術(shù)加以實現(xiàn)。FFT在數(shù)字通信、語音信號處理、圖像處理、匹配

18、濾波以及功率譜估計、仿真、系統(tǒng)分析等各個領(lǐng)域都得到了廣泛的應(yīng)用。但不管FFT在哪里應(yīng)用，一般都以卷積積分或相關(guān)積分的具體處理為依據(jù)，或者以用FFT作為連續(xù)傅里葉變換的近似為基礎(chǔ)。4.2.1利用FFT求線性卷積一快速卷積在實際中常常遇到要求兩個序列的線性卷積。如一個信號序列x(n)通過FIR濾波器時，其輸出y(n)應(yīng)是x(n)與h(n)的卷積：cdy(n)=x(n)h(n)=x(m)h(n-m)m二:二有限長序列x(n)與h(n)的卷積的結(jié)果y(n)也是一個有限長序列。假設(shè)x(n)與h(n)的長度分別為N1和N2,則y(n)的長度為N=N1+N2-1。若通過補(bǔ)零使x(n)與h(n)都加長到N點(diǎn)，

19、就可以用圓周卷積來計算線性卷積。這樣得到用FFT運(yùn)算來求y(n)值(快速卷積)的步驟如下：(1)對序列x(n)與h(n)補(bǔ)零至長為N,使N>N1+N2-1,并且N=2M(M為整數(shù))，即x(n)=«x(n),n=0,1，N1-1。n=N1,N1+1,，N-1、|h(n),n=0,1，N2-1h(n)=«0,n=N2,N2+1，N1(2)用FFT計算x(n)與h(n)的離散傅里葉變換x(n)u(FFT)uX(k)(N點(diǎn))h(n)u(FFT)uH(k)(N點(diǎn))(3)計算丫(k)=X(k)H(k)(4)用IFFT計算Y(k)的離散傅里葉反變換得：y(n)=IFFTY(k)(N

20、點(diǎn))4.2.2利用FFT求相關(guān)一快速相關(guān)互相關(guān)及自相關(guān)的運(yùn)算已廣泛的應(yīng)用于信號分析與統(tǒng)計分析，應(yīng)用于連續(xù)時間系統(tǒng)也用于離散事件系統(tǒng)。用FFT計算相關(guān)函數(shù)稱為快速相關(guān)，它與快速卷積完全類似，不同的是一個應(yīng)用離散相關(guān)定理，另一個應(yīng)用離散卷積定理。同樣都要注意到離散傅里葉變換固有的周期性，也同樣用補(bǔ)零的方法來繞過這個障礙。設(shè)兩個離散時間信號x(n)與y(n)為已知，離散互相關(guān)函數(shù)記作Rxy(n),定義為QORxy(n)='x(m)y(nm)m二：二如果x(n)與y(n)的序列長度分別為N1和N2,則用FFT求相關(guān)的計算步驟如下：(1)對序列x(n)與y(n)補(bǔ)零至長為N,使N>N1+N2-1,并且N=2m(M為整數(shù))，即x(n),n=0,1，N1-1x(n)=«0,n=N1,N1+1,，N1,、y(n),n=0,1，N2-1y(n)=30,n=N2,N2+1，N-1(2)用FFT計算x(n)與y(n)的離散傅里葉變換x(n)u(FFT)uX(k)(N點(diǎn))y(n)u(FFT)uY(k)(N點(diǎn))(3)將X(k)的虛部ImX(k)改變符號，求得其共軻X*(k)(4)計算Rxy(k)