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

2、而使數字信號處理的速度大大提高。自從1965年第一篇DFT快速算法的論文發表以來，人們已經研究出多種FFT算法，它們的復雜度和運算效率各不相同。本章主要介紹最基本的基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奇數偶數奇數偶數)12(2/122/12/012)()()(lkNMnkNklNNlWlxWWlxkX奇數(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)式。為了將如上分解過程用運算流圖表示，以便估計其運算量，觀察運算規律，總結編程方法，先介紹一種表示(5.3.4)式的蝶形圖。蝶形圖及其運算功能如圖5.3.1所示。12,.1 , 0)()()2()()()(2121NkkXWkXNkXkXWkXkXkNkN第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.1

6、蝶形運算圖ABCA CBA CB第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 根據圖5.3.2可以求得第一次分解后的運算量。圖5.3.2包括兩個N/2點DFT和N/2個蝶形，每個 點DFT需要 次復數乘法和 次復數加法運算，每個蝶形只有一次復數乘法運算和兩次復數加法運算。所以，總的復數乘法次數為總的復數加法次數為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)所示。根據WkN/m=WkmN，將圖5.3.3(a)轉換成如圖5.3.3(b)所示的標準形式的運算流圖。 第五章第五章 快速傅立葉變換快速傅立葉變換(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級蝶形構成，每級有N/2個蝶形。因此，每級需要N/2次復數乘法運算和N次復數加法運算，M級形共需復數乘法次數CM(2)和復數加法次數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的復數乘法次數為N2,復數加法次數為(N-1)N。當N1時, N2/CM(2)1，所以N越大，DIT-FFT運算效率越高。DIT-FFT算法與DFT所需乘法次數與N的關系曲線如圖5.3.4所示。例如，N=210=1024時，DIT-FFT的運算效率為而當N=211=2048時， 8 .20410210241024)2(22MCNFFTDITDFT的乘法次數的乘法次數37.372112048222)2(22MNMNNCNM第五

11、章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 圖5.3.4 DIT-FFT與DFT所需乘法次數比較曲線N(取樣點數)1024512256128640直接計算DFTDIT-FFT算法64128 2565121024乘法次數(1000)第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5.3.3 DIT-FFT的運算規律及編程思想1. 原位計算 由圖5.3.3可以看出，DIT-FFT的運算過程很有規律。2. 旋轉因子的變化規律觀察圖5.3.3(a)不難發現，第L級共有2L-1個不同的旋轉因子。N=23=8時的各級旋轉因子表示如下：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級的旋轉因子為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級運算的旋轉因子(實際編程序時，L為最外層循環變量)。 第五章第五章 快速傅立

13、葉變換快速傅立葉變換(FFT)(FFT) 3. 蝶形運算規律設序列x(n)經時域抽選(倒序)后，存入數組X中。如果蝶形運算的兩個輸入數據相距B個點，應用原位計算，則蝶形運算可表示成如下形式： 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) 設T=XL-1 (J+B)WpN=TR+jTI式中，下標R表示取實部，I表示取虛部，)()()(11JXjJXJXRLIIIRRRIIIRRRI

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

15、旋轉因子，每求出一個旋轉因子，就計算完它對應的所有2M-L個蝶形。這樣，我們可用三重循環程序實現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)( )()()()()( 輸 出結 束12LB第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 5. 序列的倒序DIT-FFT算法的輸入序列的排序看起來似乎很亂，但仔細分析就會

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

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

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。 設序列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)分解成偶數組與奇數組，當k取偶數(k=2r, r=0, 1, , N/2-1)時，1 k=偶數-1 k=奇數rn

20、NNnrnNNnWNnxnxWNnxnxrX22/12/0212/0)2()()2()()2(5.3.9) 第五章第五章 快速傅立葉變換快速傅立葉變換(FFT)(FFT) 當k取奇數(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)的關系也可用圖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時二次分解運算流圖。這樣繼續分解下去，經過M-1次分解，最后分解為2M-1個2點DFT，2點DFT就是一個基本蝶形運算流圖。當N=8時，經兩次分解，便分解為四個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) 在實際中，有時為了防止運算過程中發生溢出，將1/N分配到每一級蝶形運算中。由于1/N=(1/2) M, 因此每級的每個蝶形輸出支路均有一相乘因子1/2，這種運算的蝶形流圖如圖5.3.14所示。由圖可知，乘法次數比圖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) 如果希望直接調用FFT子程序計算IFFT，則可用下面的方法：對上式兩邊同時取共軛，得knNNkknNNkWkXNnxWkXNnx1010)(1)()(1)(由于 因此 *)(1)(1)(*10kXDFTNWkXNnxknNNk第五章第五章