按時間抽取的基2FFT算法分析

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

2、傅里葉變換(DFT)的快速算法。DFT的定狡式為N-IX伙)=工心)呼心在所有復指數值Wf的值全部已算好的情況下，要計算一個X伙)需要N次 復數乘法和N-1次復數加法。算出全部N點X伙)共需N?次復數乘法和 N(N-I)次復數加法。即計算量是與N?成正比的。FFT的基本思想：將大點數的DFT分解為若干個小點數DFT的組合，從 而城少運算量。WN因子具有以下兩個特性，可使DFT運算量盡量分解為小點數的DFT 運算：(1) 周期性：Wn = WN)k(2) 對稱性：WN,2 =利用這兩個性質，可以使DFT運算中有些項合并，以減少乘法次數。例子:求當N=4時，X(2)的值3X(2)=工 Xs)WJ

3、= x(O)VV4o + (1)W42 + (2)VV44 + x(3)46H-O=-r(0) + (2)VV40 +Wl) + x(3)VV42 (周期性)=(0) + x(2)-x(l) + X 網 (對稱性)通過合并，使乘法次數由4次減少到1次，運算量減少。FFT的算法形式有很多種，但基本上可以分為兩大類：按時間抽取(DlT)和按頻率抽取(DlF)O按時間抽取(DlT)的FTT為了將大點數的DFT分解為小點數的DFT運算，要求序列的長度N為復 合數，最常用的是N = 2“的情況(M為正整數)。該情況下的變換稱為基 2FFTo下面討論基2情況的算法。先將序列X (n)按奇偶項分解為兩組x(

4、2r) = XI (r)N 、x(2r + ) = x2(r)廠_，''勺"_將DFT運算也相應分為兩組JV-IX(k) = DFTx(n)=x(n)W0NTNT=2>S)W(+;?=0H=OH為偶數n為奇數JV/2-1Af/2-1=Yjx(2r)Vrk ÷ £心 + 1)W嚴r-()r-()JV2-1W/2-1(r)W +; x2(r)r rOr-0/2-12-lX"叱2+叫2(M農(因為叱嚴=叱2)r-()r-0= X") + WjX2 伙)其中XI伙）、X2（Ar）分別是x1（n）s勺5）的"2點的DFTN

5、/2-172-lNX1W = X1 (r)VZ2 =工v(2r)VV2,0k-r-()r-0LN/2-1NJ 27NX2(k)=jx2(r)WZ2 = (2r + l)2,0Arl-l r-0r-0L至此，一個N點DFT被分解為兩個N/2點的DFT。上面是否將全部N點的X伙）求解出來了N分析：X伙）和X?伙）只有N/2個點（&=0丄,一一1）,則由 2X伙）=Xl伙）+ WX2伙）只能求出X伙）的前N/2個點的DFT,要求出全部N點的X伙），需要找出Xl W > X2（k）和X伙+ N2）的關系，其中 =0,l, -,-Io 由式子 X伙)=Xl(k) + WX2(k)可得2X伙

6、 + N/2) = X1 (k + N2) + Wf2X2伙+ N2)化簡得NX(k + N2) = = Xl(R)-W(X2(Q, R = (M,一_12這樣N點DFT可全部由下式確定出來:X(k + N2) = Xk)-W.X2(k)2上式可用一個專用的碟形符號來表示，這個符號對應一次復乘和兩次復加運算。a + Wba-W通過這樣的分解以后，每一個N/2點、的DFT只需要()2 =次復數乘2NN法，兩個N/2點的DFT需要2()2 =次復乘，再加上將兩個N/2點 22NN NDFT合并成為N點DFT時有N/2次與W因子相乘，一共需要 + 222次復乘?？梢姡ㄟ^這樣的分解，運算量節(jié)省了近一

7、半。因為N = 2'". N/2仍然是偶數，因此可以對兩個N/2點的DFT再分別作進一步的分解，將兩個N/2點的DFT分解成兩個N/4點的DFTo例如對x1(r),可以在按其偶數部分及奇數部分進行分解：Fg心0,1,上1g + l) = x4()4則的運算可相應分為兩組：yv4-lf4-1X")=工 m(2)Vc2 + > + 1)叱於“/-0/-()jV/4-1N/4-1=>3(W為+叱爲工兀叱為/-0Z-OX3(k) + W,2X4 伙)jt=0,l,-l一48 點 DFT×(0 (I) ×(2) 必(3) <×(

8、4) 啦5 嘆6ZXXXXXXX(XXXXXXXX (0)M)(21問 £ x(×(×(×(x(×(×(X2dftXXXXXXXXOOOOoOoOO-OOOQO OoX1W = X3伙) + w,i伙),Al N I<k = O _ 1X、(k + N4) = X3(k)-WtkX4 伙) 4同樣，對X2(Zr)也可進行類似的分解。一直分解下去，最后是2點的DFT, 2點DFT的運算也可用碟形符號來表示。這樣，對于一個N = 23 =S 的DFT運算，其按時間抽取的分解過程及完整流圖如下圖所示。x(0- XrIb ×(

9、2kX (3Io- × M)O- x(5J0.X (6卩 ×Rk(O)同(2)問(1)(5)Pl ×(x(x(×(×(x(x(X這種方法，由于每一步分解都是按輸入序列在時域上的次序是屬于偶數 還是奇數來抽取的，故稱為"時間抽取法S分析上面的流圖，N = 2“，一共要進行M次分解，構成了從x(n)到 X(k)的M級運算過程。毎一級運算都是由N/2個蝶形運算構成，因此每一級 運算都需要N/2次復乘和N次復加，則按時間抽取的M級運算后總共需要NN復數乘法次數：HIF =M=-IosNf 22 復數加法次數：UIt =N-M= /Vlog2

10、N根據上面的流圖，分析FFT算法的兩個特點，它們對FFT的軟硬件 構成產生很大的影響。(1 )原位運算也稱為同址運算，當數據輸入到存儲器中以后，毎一級運算的結果仍然存儲 在原來的存儲器中，直到我后輸出，中間無需其它的存儲器。根據運算流圖 分析原位運算是如何進行的。原位運算的結構可以節(jié)省存儲單元，降低設備 成本。(2) 變址分析運算流圖中的輸入輸出序列的順序，輸出按順序，輸入是“碼位倒置” 的順序。見圖。自然順序二進制表示碼位倒置碼位倒置順序OOOO000010011004201001023011110641000011510110156110011371111117在實際運算中，直接舟輸入數據

11、X (n)按碼位倒置的順序排好輸入很 不方便，一般總是先按自然順序輸入存儲單元，然后通過變址運算將自然順 序的存儲換成碼位倒置順序的存儲，這樣就可以進行FFT的原位運算。變質 的功能如圖所示。用軼件實現是通用釆用雷徳(Rader)算法，算出I的倒 序J以后立即將輸入數據X(I)和X(J)對換。盡管變址運算所占運算董的比 例很小，但對某些高要求的應用(尤其在實時信號處理中)，也可設法用適 當的電路結構直接實現變址。例如單片數字信號處理器TMS320C25就有專用 于FFT的二進制碼變址模式。4. 2按頻率抽取(DlF)的FTT除吋間抽取法外，另外一種普遍使用的FFT結構是頻率抽取法。頻率抽取法將

12、輸入序列不是按奇、偶分組，而是將N點DFT寫成前后兩部分:N-IX(k) = DFnx(n)= x(n)WnH-O(T2)-1N-I= DS)W+ 2>(")wj0/r-4V2JV/2-1N/2-1=+ 2>(n + N2)Wyz2M-0r0N/2-1=yjx(n) +JV/2)V(H + / 2)W-O因為Wf" =_1,Wf 2)* =(_i)x , k為偶數時(_1)女=1, k為奇數時(一1)*=一1,由此可將X(k)分解為偶數紐和奇數紐：72-lXg= 2>() + (-1)飯+ N/2)卩曾n-0jV/2-1X(2r)=工x()+心+ N2)W

13、"H-Of/2-1=x(n) + x(n + N 2)V2rU)JV/2-1X (2r +1)=工x() - x(n + N / 2) Wi"-0N/2-1 =2>U)7G + N2)W.W篇 2片-0Xl (；0 = X(II) + x(zz + yv2)令 <H = O 丄，N/2 12() = x() 一 x(n + N2) ;這兩個序列都是N/2點的序列，對應的是兩個N/2點的DFT運算：2-lX(2J= x,(nM'2H-O2-lX（2 廠+ 1）=J>2陀;2-4）這樣，同樣是將一個N點的DFT分解為兩個N/2點的DFT 了。頻率抽選法

14、對應的碟形運算關系圖如下：對于N=8吋頻率抽取法的FFT流圖如下:* IJ -IJ IJ J -IJ -IJ IJ 卩 1 234 567 fl¾ fl¾ fl¾ ( ( H ( ××××××××OOOAVO : : 廠 ； 二O歹，矽OAVOo-Or*cb77-.2占C-八“DFT送2占DFT-1廠F>o2占C-八"DFTo-12占厶八、DFTr U 一 一/ “ 卩g尼他 P R XXXX XXXX 06OOO600 kkk-IJ -IJ U -IJ -IJ IJ

15、-IJ 0426 1537«11 ll «11 «11 H »11 <l Xxxx Xxxx這種分組的辦法由于每次都是按輸出X(k)在頻域的順序上是屬于偶數還是 奇數來分組的，稱為頻率抽取法。與前面按時間抽取的方法相比，相同點 問題：如何利用快速算法計算IDFT分析IDFT的公式：1 AMX(H) = IDFT X 伙)=若 X 伙)叱卩","=0,1,N 1比較DFT的公式：,-lXa) = DFTx(n) = YXS)W化 & =0,屮一1n0得知可用兩種方法來實現IDFT的快速算法：(1 )只要把DFT運算中的每一

16、 個系數Wy該為Wy 并且最后再乘以常數丄，就可以用時間抽取法或N頻率抽取的FFT算法來直接計算IDFTo這種方法需要對FFT的程序和參數稍 加 改 動 才 能 實 現。 (2) 因 為1 N-1W0 = 77工(QW嚴=77MT(燈M = O,1,N-1,也就 N A：-oN是說，可先將X(k)取共純變換，即將X(k)的虛部乘以一 1 ,就可直接調用 FFT的程序，罠后再對運算結果取一次共純變換并乘以常數1/N即可得到X (n) 的值。這種方法中，FFT運算和IFFT運算都可以共用一個子程序塊，在使 用通用計算機或用硬件實現吋比較方便。4.1.3混合基FFT算法以上討論的是基2的FFT算法，

17、即N = 2M的情況，這種情況實際上 使用得最多，這種FFT運算，程序簡單，效率很鬲，用起來很方便。另外， 在實際應用時，有限長序列的長度N到底是多少在很大程度時是由人為因素 確定的，因此，大多數場合人們可以將N選定為N = 2M 9從而可以直接調 用以2為基數的FFT運算程序。如果長度N不能認為確定，而N的數值又不是以2為基數的整數次方， 一般可有以下兩種處理方法：(1 ) 將x(n)用補零的方法延長，使N增長到最鄰近的一個 N = 2M數值。例如,N=30,可以在序列x(n)中補進x(30)=x(31) =O兩個零值點，使N=32。如果計算FFT的目的是為了 了解整 個頻譜，而不是特定頻率

18、點，則此法可行。因為有限長序列補 零以后并不影響其頻譜X(ejw),只是頻譜的采樣點數增加而 已。(2)如果要求特定頻率點的頻譜，則N不能改變。如果N為復合數， 則可以用以任意數為基數的FFT #法來計算?？焖俑道锶~變換的基本思想就是要將DFT的運算盡量分小。例如，N=6吋，可以按照N=3X2分解，將6點DFT分解為3組2點DFTo舉例：N=9時的快速算法。4.2快速傅里葉變換的應用凡是可以利用傅里葉變換來進行分析、綜合、變換的地方，都可以利 用FFT算法及運用數字計算技術加以實現。FFT在數字通信、語音信號處理、 圖像處理、匹配濾波以及功率譜仕計.仿真.系統(tǒng)分析等各個領域都得到了 廣泛的應用

19、。但不管FFT在哪里應用，一般都以卷積積分或相關積分的具體 處理為依據，或者以用FFT作為連續(xù)傅里葉變換的近似為基礎。利用FFT求線性卷積一快速卷積在實際中常常遇到要求兩個序列的線性卷積。如一個信號序列x(n)通 過FlR濾波器時，其輸出y(n)應是x(n)與h(n)的卷積：Xy(n) = x(n) * h(n)= 工 x(ni)h(n 一 In)“10C有限長序列x(n)與h(n)的卷積的結果y(n)也是一個有限長序列。假設x (n) 與h(n)0長度分別為Nl和N2,則y(n)的長度為N=N1+N2T°若通過補零使 x(n)與h(n)都加長到N點，就可以用圓周卷積來計算線性卷枳。

20、這樣得到 用FFT運算來求y(n)值(快速卷積)的步驟如下：(1 )對序列x(n)與h(n)補零至長為N,使NN1+N2-1 ,并且N = 2M (M為整數),即Lr(n), H = 0,1,9 /Vl -1X02)=O, H = NhNz ，N-If/?(/?), n = O,1,7V2-1 Il(H)= <O, n = V2,N2 + l,N-1(2) 用FFT計算x(n)與h(n)的離散傅里葉變換X(U) <=> (FFT) o X伙)(N 點)h(n) <=> (FFT) o H伙)(N 點)(3) 計算 Y(k)=X(k)H(k)(4) 用IFFT計算Y

21、(k)的離散傅里葉反變換得：y(n) = IFFTY(k) (N 點)4. 2. 2利用FFT求相關一快速相關互相關及自相關的運算已廣泛的應用于信號分析與統(tǒng)計分析，應用于連 續(xù)時間系統(tǒng)也用于離散事件系統(tǒng)。用FFT計算相關函數稱為快速相關，它與快速卷積完全類似，不同的是 一個應用離散相關定理，另一個應用離散卷積定理。同樣都要注意到離散傅 里葉變換固有的周期性，也同樣用補零的方法來繞過這個障礙。設兩個離散時間信號x(n)與y(n)為已知，離散互相關函數記作&、.(”)，定義為OCRan= YJX(If j)y(n + m)J?!-X如果x(n)與y(n)的序列長慶分別為NI和N2,則用FFT求相關的計算步驟 如下：(1 )對序列x(n)與y(n)補零至長為N,使NMNI+N2T,并且N = 2"(M為整數),即x(n),2 = 0丄，Nl-IXW =O, n = Nl,M + l,，N-1b(n), =0,，N2-1y(n)=, n = N2N2 + X 、N-(2) 用FFT計算x(n)與y(n)的離散傅里葉變換X(H) O (FFT) o X伙)(N 點)y(n) 0 (FFT) 0 Yg(N 點)(3) 將X(k)的虛部lmX