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

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

2、而使數(shù)字信號處理的速度大大提高。自從1965年第一篇DFT快速算法的論文發(fā)表以來，人們已經(jīng)研究出多種FFT算法，它們的復(fù)雜度和運(yùn)算效率各不相同。本章主要介紹最基本的基2FFT算法及其編程方法。第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.2 直接計(jì)算直接計(jì)算DFT的特點(diǎn)及減少運(yùn)算量的基的特點(diǎn)及減少運(yùn)算量的基本途徑本途徑長度為N的序列x(n)的N點(diǎn)DFT為WmN的周期性：10)()()(NnknNnWnxnxDFTkX點(diǎn)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)點(diǎn)DFT為knNNnNWnxnxDFTkX10)()()(點(diǎn)k=0, 1, , N-1第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 將上面的和式按n的奇偶性分解為令x1(l)=x(2l), x2(l)=x(2l+1)。因?yàn)閃2klN=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)分解為兩個(gè)N/2長的序列x1(l)和x2(l)，則N點(diǎn)DFT可分解為兩個(gè)N/2點(diǎn)DFT來計(jì)算。用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點(diǎn)點(diǎn)kNNkNWW2第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 這樣，就將N點(diǎn)DFT的計(jì)算分解為計(jì)算兩個(gè)N/2點(diǎn)離散傅立葉變換X1(k)和X2(k)，再計(jì)算(5.3.4)式。為了將如上分解過程用運(yùn)算流圖表示，以便估計(jì)其運(yùn)算量，觀察運(yùn)算規(guī)律，總結(jié)編程方法，先介紹一種表示(5.3.4)式的蝶形圖。蝶形圖及其運(yùn)算功能如圖5.3.1所示。12,.1 , 0)()()2()()()(2121NkkXWkXNkXkXWkXkXkNkN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.1

6、蝶形運(yùn)算圖ABCA CBA CB第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 根據(jù)圖5.3.2可以求得第一次分解后的運(yùn)算量。圖5.3.2包括兩個(gè)N/2點(diǎn)DFT和N/2個(gè)蝶形，每個(gè) 點(diǎn)DFT需要 次復(fù)數(shù)乘法和 次復(fù)數(shù)加法運(yùn)算，每個(gè)蝶形只有一次復(fù)數(shù)乘法運(yùn)算和兩次復(fù)數(shù)加法運(yùn)算。所以，總的復(fù)數(shù)乘法次數(shù)為總的復(fù)數(shù)加法次數(shù)為2N2)2(N2) 12(NN2| ) 1(222)2(212NNNNNN22222) 12(2NNNN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.2 8點(diǎn)DFT一次時(shí)域抽取分解運(yùn)算流圖N/2點(diǎn)DFTWN0N/2點(diǎn)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點(diǎn)DIT-FFT的運(yùn)算流圖如圖5.3.3(a)所示。根據(jù)WkN/m=WkmN，將圖5.3.3(a)轉(zhuǎn)換成如圖5.3.3(b)所示的標(biāo)準(zhǔn)形式的運(yùn)算流圖。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.3 8點(diǎn)DIT-FFT運(yùn)算流圖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的運(yùn)算效率DIT-FFT的運(yùn)算效率指直接計(jì)算DFT的運(yùn)算量與DIT-FFT的運(yùn)算量之比。由圖5.3.3可見，N=2M時(shí)，其DIT-FFT運(yùn)算流圖由M級蝶形構(gòu)成，每級有N/2個(gè)蝶形。因此，每級需要N/2次復(fù)數(shù)乘法運(yùn)算和N次復(fù)數(shù)加法運(yùn)算，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。直接計(jì)算N點(diǎn)DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(N-1)N。當(dāng)N1時(shí), N2/CM(2)1，所以N越大，DIT-FFT運(yùn)算效率越高。DIT-FFT算法與DFT所需乘法次數(shù)與N的關(guān)系曲線如圖5.3.4所示。例如，N=210=1024時(shí)，DIT-FFT的運(yùn)算效率為而當(dāng)N=211=2048時(shí)， 8 .20410210241024)2(22MCNFFTDITDFT的乘法次數(shù)的乘法次數(shù)37.372112048222)2(22MNMNNCNM第五

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

12、 J=0 L=2時(shí)， WpN =WJN/2=WJ2L J=0, 1L=3時(shí)， 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級運(yùn)算的旋轉(zhuǎn)因子(實(shí)際編程序時(shí)，L為最外層循環(huán)變量)。 第五章第五章 快速傅立

13、葉變換快速傅立葉變換(FFT)(FFT) 3. 蝶形運(yùn)算規(guī)律設(shè)序列x(n)經(jīng)時(shí)域抽選(倒序)后，存入數(shù)組X中。如果蝶形運(yùn)算的兩個(gè)輸入數(shù)據(jù)相距B個(gè)點(diǎn)，應(yīng)用原位計(jì)算，則蝶形運(yùn)算可表示成如下形式： 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表示取實(shí)部，I表示取虛部，)()()(11JXjJXJXRLIIIRRRIIIRRRI

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

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

16、發(fā)現(xiàn)這種倒序是很有規(guī)律的。由于N=2M，因此順序數(shù)可用M位二進(jìn)制數(shù)(nM-1 nM-2n1n0)表示。M次偶奇時(shí)域抽選過程如圖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時(shí)以二進(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后，將原存儲(chǔ)器中存放的輸入序列重新按倒序排列。設(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)是完成計(jì)算倒序值的運(yùn)算流程圖。 第五章第五章 快速傅立葉變換快速傅立葉變換(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)前后對半分開，得到兩個(gè)子序列，其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)時(shí)，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)時(shí)令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所示的蝶形運(yùn)算流圖符號表示。圖5.3.10表示N=8時(shí)一次分解的運(yùn)算流圖。rnNNnrnNNnWnxrXWnxrX2/12/022/12/01)() 12()()2(5.3.11)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.9 DIF-FFT蝶形運(yùn)算流圖符號x(n)2(Nnx)2()()(1NnxnxnxnNWNnxnxnx)2()()(2nNW第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.10 DIF-FFT一次分解運(yùn)算流圖(N=8)N/2點(diǎn)DFTWN0N/2點(diǎn)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時(shí)二次分解運(yùn)算流圖。這樣繼續(xù)分解下去，經(jīng)過M-1次分解，最后分解為2M-1個(gè)2點(diǎn)DFT，2點(diǎn)DFT就是一個(gè)基本蝶形運(yùn)算流圖。當(dāng)N=8時(shí)，經(jīng)兩次分解，便分解為四個(gè)2點(diǎn)DFT, 如圖5.3.11所示。N=8的完整DIF-FFT運(yùn)算流圖如圖5.3.12所示。第五章第五章 快速傅立葉變換快速

23、傅立葉變換(FFT)(FFT) 圖5.3.11 DIF-FFT二次分解運(yùn)算流圖(N=8)N/4點(diǎn)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點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFT第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.12 DIF-FFT運(yùn)算流圖(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的運(yùn)算公式：由DIF-FFT運(yùn)算流圖改成的DIT-IFFT運(yùn)算流圖如圖5.3.13所示。rnNNkrnNNnWkXNnxIDEFnxWnxnxDEFkX1012/0)(1)()()()()(第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖 5.3.13 DIT-IFFT運(yùn)算流圖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) 在實(shí)際中，有時(shí)為了防止運(yùn)算過程中發(fā)生溢出，將1/N分配到每一級蝶形運(yùn)算中。由于1/N=(1/2) M, 因此每級的每個(gè)蝶形輸出支路均有一相乘因子1/2，這種運(yùn)算的蝶形流圖如圖5.3.14所示。由圖可知，乘法次數(shù)比圖5.3.13增加了(N/2)(M-1)次。 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.14 DIT-IFFT運(yùn)

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子程序計(jì)算IFFT，則可用下面的方法：對上式兩邊同時(shí)取共軛，得knNNkknNNkWkXNnxWkXNnx1010)(1)()(1)(由于 因此 *)(1)(1)(*10kXDFTNWkXNnxknNNk第五章第五章