第五章 快速傅立葉變換(FFT)

1、第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT) 5.1 引言引言5.2 直接計算直接計算DFT的特點及減少運算量的基本途徑的特點及減少運算量的基本途徑5.3 基基2FFT算法算法習(xí)題與上機題習(xí)題與上機題第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.1 引言引言影響數(shù)字信號處理發(fā)展的最主要因素之一是處理速度。DFT使計算機在頻域處理信號成為可能，但是當(dāng)N很大時，直接計算N點DFT的計算量非常大?？焖俑盗⑷~變換(FFT： Fast Fourier Transform)可使實現(xiàn)DFT的運算量下降幾個數(shù)量級，從

2、而使數(shù)字信號處理的速度大大提高。自從1965年第一篇DFT快速算法的論文發(fā)表以來，人們已經(jīng)研究出多種FFT算法，它們的復(fù)雜度和運算效率各不相同。本章主要介紹最基本的基2FFT算法及其編程方法。第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.2 直接計算直接計算DFT的特點及減少運算量的基的特點及減少運算量的基本途徑本途徑長度為N的序列x(n)的N點DFT為WmN的周期性：10)()()(NnknNnWnxnxDFTkX點k=0,1,N-1mNmNjlNmNjlNmNWeeW2)(2(5.2.1)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) WmN的對稱性：

3、 (WN-mN) *= WmN (5.2.2)(5.2.3)mNNmNWW2第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3 基基2FFT算法算法5.3.1 DIT-FFT算法序列x(n)的N(N=2-M)點DFT為knNNnNWnxnxDFTkX10)()()(點k=0, 1, , N-1第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 將上面的和式按n的奇偶性分解為令x1(l)=x(2l), x2(l)=x(2l+1)。因為W2klN=WklN/2, 所以上式可寫成 )12(12/212/) 12()2()()()(lkNNnlkNNnknNnknNnW

4、lxWlxWnxWnxkX奇數(shù)偶數(shù)奇數(shù)偶數(shù))12(2/122/12/012)()()(lkNMnkNklNNlWlxWWlxkX奇數(shù)(5.3.1) 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) (5.3.1)式說明，按n的奇偶性將x(n)分解為兩個N/2長的序列x1(l)和x2(l)，則N點DFT可分解為兩個N/2點DFT來計算。用X1(k)和X2(k)分別表示將(5.3.2)式和(5.3.3)式代入(5.3.1)式，并利用和X1(k)、 X2(k)的隱含周期性可得到：12,.,1 , 0 )()()(12,.,1 , 0 )()()(12/02/22/2212/02/12/

5、11NkWlxlxDFTkXNkWlxlxDFTkXNlklNNNlklNN點點kNNkNWW2第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 這樣，就將N點DFT的計算分解為計算兩個N/2點離散傅立葉變換X1(k)和X2(k)，再計算(5.3.4)式。為了將如上分解過程用運算流圖表示，以便估計其運算量，觀察運算規(guī)律，總結(jié)編程方法，先介紹一種表示(5.3.4)式的蝶形圖。蝶形圖及其運算功能如圖5.3.1所示。12,.1 , 0)()()2()()()(2121NkkXWkXNkXkXWkXkXkNkN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.1

6、蝶形運算圖ABCA CBA CB第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 根據(jù)圖5.3.2可以求得第一次分解后的運算量。圖5.3.2包括兩個N/2點DFT和N/2個蝶形，每個 點DFT需要 次復(fù)數(shù)乘法和 次復(fù)數(shù)加法運算，每個蝶形只有一次復(fù)數(shù)乘法運算和兩次復(fù)數(shù)加法運算。所以，總的復(fù)數(shù)乘法次數(shù)為總的復(fù)數(shù)加法次數(shù)為2N2)2(N2) 12(NN2| ) 1(222)2(212NNNNNN22222) 12(2NNNN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.2 8點DFT一次時域抽取分解運算流圖N/2點DFTWN0N/2點DFTWN1WN2WN

7、3x(0)X1(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) N=8點DIT-FFT的運算流圖如圖5.3.3(a)所示。根據(jù)WkN/m=WkmN，將圖5.3.3(a)轉(zhuǎn)換成如圖5.3.3(b)所示的標(biāo)準(zhǔn)形式的運算流圖。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.3 8點DIT-FFT運算流圖WN0WN1WN2WN3WN04x(0)x(4)x(2)x(6)

8、x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1

9、)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)WN04WN04WN0412NW02NW12NW02NW(a)(b)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3.2 DIT-FFT的運算效率DIT-FFT的運算效率指直接計算DFT的運算量與DIT-FFT的運算量之比。由圖5.3.3可見，N=2M時，其DIT-FFT運算流圖由M級蝶形構(gòu)成，每級有N/2個蝶形。因此，每級需要N/2次復(fù)數(shù)乘法運算和N次復(fù)數(shù)加法運算，M級形共需復(fù)數(shù)乘法次數(shù)CM(2)和復(fù)數(shù)加法次數(shù)CA(2)分別為 (5.3.5) CA(2) =

10、NM=N lb N (5.3.6)lbNNMNCM22)2(第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 式中，lb N=log2 N。直接計算N點DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(N-1)N。當(dāng)N1時, N2/CM(2)1，所以N越大，DIT-FFT運算效率越高。DIT-FFT算法與DFT所需乘法次數(shù)與N的關(guān)系曲線如圖5.3.4所示。例如，N=210=1024時，DIT-FFT的運算效率為而當(dāng)N=211=2048時， 8 .20410210241024)2(22MCNFFTDITDFT的乘法次數(shù)的乘法次數(shù)37.372112048222)2(22MNMNNCNM第五

11、章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.4 DIT-FFT與DFT所需乘法次數(shù)比較曲線N(取樣點數(shù))1024512256128640直接計算DFTDIT-FFT算法64128 2565121024乘法次數(shù)(1000)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3.3 DIT-FFT的運算規(guī)律及編程思想1. 原位計算 由圖5.3.3可以看出，DIT-FFT的運算過程很有規(guī)律。2. 旋轉(zhuǎn)因子的變化規(guī)律觀察圖5.3.3(a)不難發(fā)現(xiàn)，第L級共有2L-1個不同的旋轉(zhuǎn)因子。N=23=8時的各級旋轉(zhuǎn)因子表示如下：L=1時， WpN=WJN/4=WJ2L

12、 J=0 L=2時， WpN =WJN/2=WJ2L J=0, 1L=3時， WpN = WJN=WJL 2J=0, 1, 2, 3第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 對N=2M的一般情況，第L級的旋轉(zhuǎn)因子為WpN=WJL 2J=0, 1, 2, , 2L-1-1由于 2L=2M2L-M=N2L-M所以 WpN =WJN2L-M=WJ2N M-L NJ=0, 1, 2, , 2L-1-1 (5.3.7) p=J2M-L (5.3.8) 這樣，就可按(5.3.7)式和(5.3.8)式確定第L級運算的旋轉(zhuǎn)因子(實際編程序時，L為最外層循環(huán)變量)。 第五章第五章 快速傅立

13、葉變換快速傅立葉變換(FFT)(FFT) 3. 蝶形運算規(guī)律設(shè)序列x(n)經(jīng)時域抽選(倒序)后，存入數(shù)組X中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應(yīng)用原位計算，則蝶形運算可表示成如下形式： XL(J) XL-1 (J)+XL-1 (J+B)WpN XL(J+B) XL-1 (J)-XL-1 (J+B)WpN式中 p=J2M-L; J=0, 1, , 2L-1-1; L=1, 2, , M第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 設(shè)T=XL-1 (J+B)WpN=TR+jTI式中，下標(biāo)R表示取實部，I表示取虛部，)()()(11JXjJXJXRLIIIRRRIIIRRRI

14、RLRIIIRRTJXBJXTJXBJXTJXJXTJXJXJjXJXJXpNBJXpNBJXTpNBJXpNBJXT)()()()()()()()()()()(2sin)(2cos)(2sin)(2cos)(第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 4. 編程思想及程序框圖仔細(xì)觀察圖 5.3.3，還可歸納出一些編程序有用的運算規(guī)律： 第L級中，每個蝶形的兩個輸入數(shù)據(jù)相距B=2L-1個點； 同一旋轉(zhuǎn)因子對應(yīng)著間隔為2L點的2M-L個蝶形?？偨Y(jié)上述運算規(guī)律，便可采用下述運算方法。先從輸入端(第1級)開始，逐級進(jìn)行，共進(jìn)行M級運算。在進(jìn)行第L級運算時，依次求出2L-1個不同的

15、旋轉(zhuǎn)因子，每求出一個旋轉(zhuǎn)因子，就計算完它對應(yīng)的所有2M-L個蝶形。這樣，我們可用三重循環(huán)程序?qū)崿F(xiàn)DIT-FFT運算，程序框圖如圖 5.3.5 所示。第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.5 DIT-FFT運算和程序框圖開 始送入 x(n)， MN2 M倒 序L1 , MJ0 , B 1P2 M LJk J , N1 , 2LTkXWBkXkXBkXWBkXkXTpNpN)( )()()()()( 輸 出結(jié) 束12LB第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5. 序列的倒序DIT-FFT算法的輸入序列的排序看起來似乎很亂，但仔細(xì)分析就會

16、發(fā)現(xiàn)這種倒序是很有規(guī)律的。由于N=2M，因此順序數(shù)可用M位二進(jìn)制數(shù)(nM-1 nM-2n1n0)表示。M次偶奇時域抽選過程如圖5.3.6所示。第一次按最低位n0的0和1將x(n)分解為偶奇兩組，第二次又按次低位n1的0、 1值分別對偶奇組分解； 依次類推，第M次按nM-1位分解，最后所得二進(jìn)制倒序數(shù)如圖5.3.6所示。表5.3.1列出了N=8時以二進(jìn)制數(shù)表示的順序數(shù)和倒序數(shù)。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.6 形成例序的樹狀圖(N=23)01010101010101(n2n1n0)200004261537100010110001101011111第五

17、章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 表 5.3.1 順序和倒序二進(jìn)制數(shù)對照表 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 形成倒序J后，將原存儲器中存放的輸入序列重新按倒序排列。設(shè)原輸入序列x(n)先按自然順序存入數(shù)組A中。例如，對N=8, A(0)，A(1)，A(2)，A(7)中依次存放著x(0),x(1), , x(7)。對x(n)的重新排序(倒序)規(guī)律如圖5.3.7所示。倒序的程序框圖如圖5.3.8所示，圖5.3.8中的虛線框內(nèi)是完成計算倒序值的運算流程圖。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.7 倒序規(guī)

18、律x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.8 倒序程序框圖221NNLHJNLHI1 , N1I J？T)J(A)J(X)I(A)I(XTJ K？LHK KJJ2KKKJJYNNY第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3.4 DIF-FFT在基2快速算法中，頻域抽取法FFT也是

19、一種常用的快速算法，簡稱DIF-FFT。 設(shè)序列x(n)長度為N=2M, 首先將x(n)前后對半分開，得到兩個子序列，其DFT可表示為如下形式：12/012/02/12/0)2/(12/012/12/010)2()()2()()()()()()(NnknNNnkNNNnNnkNNnknNNNnknNNnknNNnknNWNnxWnxWNnxWnxWnxWnxWnxnxDFTkX第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 式中 WkN/2N=(-1) k=將X(k)分解成偶數(shù)組與奇數(shù)組，當(dāng)k取偶數(shù)(k=2r, r=0, 1, , N/2-1)時，1 k=偶數(shù)-1 k=奇數(shù)rn

20、NNnrnNNnWNnxnxWNnxnxrX22/12/0212/0)2()()2()()2(5.3.9) 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 當(dāng)k取奇數(shù)(k=2r+1, r=0, 1, , N/2-1)時令nrNnNNnrnNNnWWNnxnxWNnxnxrX2/212/0)12(12/0)2()()2()() 12(5.3.10)12/,.2 , 1 , 0)2()()()2()()(221NnWNnxnxnxNnxnxnxnN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 將x1(n)和x2(n)分別代入(5.3.9)式和(5.3.10)式，可

21、得 x1(n) 、 x2(n)和x (n)的關(guān)系也可用圖5.3.9所示的蝶形運算流圖符號表示。圖5.3.10表示N=8時一次分解的運算流圖。rnNNnrnNNnWnxrXWnxrX2/12/022/12/01)() 12()()2(5.3.11)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.9 DIF-FFT蝶形運算流圖符號x(n)2(Nnx)2()()(1NnxnxnxnNWNnxnxnx)2()()(2nNW第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.10 DIF-FFT一次分解運算流圖(N=8)N/2點DFTWN0N/2點DFT

22、WN1WN2WN3X(0)x1(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.11示出了N=8時二次分解運算流圖。這樣繼續(xù)分解下去，經(jīng)過M-1次分解，最后分解為2M-1個2點DFT，2點DFT就是一個基本蝶形運算流圖。當(dāng)N=8時，經(jīng)兩次分解，便分解為四個2點DFT, 如圖5.3.11所示。N=8的完整DIF-FFT運算流圖如圖5.3.12所示。第五章第五章 快速傅立葉變換快速

23、傅立葉變換(FFT)(FFT) 圖5.3.11 DIF-FFT二次分解運算流圖(N=8)N/4點DFTWN0WN1WN2WN3x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)WN0WN2WN0WN2N/4點DFTN/4點DFTN/4點DFT第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.12 DIF-FFT運算流圖(N=8)WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)x(0)x(1)x(2)x(3)

24、x(4)x(5)x(6)x(7)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3.5 IDFT的高效算法上述FFT算法流圖也可以用于離散傅立葉逆變換(IDFT: Inverse Discrete Fourier Transform)。比較DFT和IDFT的運算公式：由DIF-FFT運算流圖改成的DIT-IFFT運算流圖如圖5.3.13所示。rnNNkrnNNnWkXNnxIDEFnxWnxnxDEFkX1012/0)(1)()()()()(第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.13 DIT-IFFT運算流圖WN0WN1WN2WN3WN

25、0WN0N1x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN2WN2N1N1N1N1N1N1N1第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 在實際中，有時為了防止運算過程中發(fā)生溢出，將1/N分配到每一級蝶形運算中。由于1/N=(1/2) M, 因此每級的每個蝶形輸出支路均有一相乘因子1/2，這種運算的蝶形流圖如圖5.3.14所示。由圖可知，乘法次數(shù)比圖5.3.13增加了(N/2)(M-1)次。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.14 DIT-IFFT運

26、算流圖(防止溢出)WN02121x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)212121WN121WN221WN3212121WN021WN2212121WN021WN22121WN02121WN0212121WN021WN021第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方法：對上式兩邊同時取共軛，得knNNkknNNkWkXNnxWkXNnx1010)(1)()(1)(由于 因此 *)(1)(1)(*10kXDFTNWkXNnxknNNk第五章第五章