DSP 第4章_0905

1、第4章 快速傅里葉變換(FFT) 第第4章章 快速傅里葉變換快速傅里葉變換(FFT) 4.1 引言引言 4.2 基基2FFT算法算法 4.3 進一步減少運算量的措施進一步減少運算量的措施 4.4 其他快速算法簡介其他快速算法簡介習題與上機題習題與上機題第4章 快速傅里葉變換(FFT) 4.1 引引 言言 在快速傅里葉變換FFT(Fast Fourier Transform)出現(xiàn)以前，直接用DFT算法進行譜分析和信號的實時處理是不切實際的。算法改進的途徑：算法改進的途徑：1、把長序列分解成短序列，然后組合出結果。、把長序列分解成短序列，然后組合出結果。2、利用計算中的對稱性，減少計算量。、利用計

2、算中的對稱性，減少計算量。旋轉因子nkNW第4章 快速傅里葉變換(FFT) 根據(jù)第3章所講內容，可知有以下幾個特點:(1) 共軛對稱性： (2) 周期性：。 (3) 可約性： 。 nkNWnkNWnkNWnkNNnkNnkNnkNWWWW2)(nNkNNnkNnkNWWW)()(nkNnkmNmWWnkNW以N=4為例，利用周期性，利用對稱性，54144404)()(,WWWWWWWnNkNNnkNnkN022,)(NNnkNNnkNnkNnkNWWWWWW第4章 快速傅里葉變換(FFT) WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWNNNN

3、NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN101000001010000012302020321000009630642032100000第4章 快速傅里葉變換(FFT) 求出四點DFT實際上只需要一次復數(shù)乘法。111111)3() 1 ()2()0() 3()3() 1 ()2()0()2()3() 1 ()2()0() 1 ()3() 1 ()2()0()0() 3()2() 1 ()0(111111111111) 3()2() 1 ()0(WxxxxXxxxxXWxxxxXxxxxXxxxxWWWWXXXX第4章 快速傅里葉變換(FFT) 1

4、、將時域序列按奇偶劃分，稱之為時間抽取(Decimation in Time，DIT)算法或基2時分算法;2、將頻域DFT按奇偶劃分，稱之為頻率抽取(Decimation in Frequency，DIF)算法或基2頻分算法。第4章 快速傅里葉變換(FFT) 4.2.2 時域抽取法時域抽取法(基基2FFT)基本原理基本原理設序列x(n)的長度為N，且滿足N=2M，M為自然數(shù)。按n的奇偶把x(n)分解為兩個N/2點的子序列12( )(2 )0112( )(21)0112Nx rxrrNx rxrr， ， ，第4章 快速傅里葉變換(FFT) 則x(n)的DFT為/2 1/2 12(21)00/2

5、1/2 1221200( )( )( )(2 )(21)( )( )knknNNnnNNkrkrNNrrNNkrkkrNNNrrX kx n Wx n Wxr WxrWx r WWx r W偶數(shù)奇數(shù)第4章 快速傅里葉變換(FFT) (4.2.7) (4.2.8) 運算可用圖4.2.1所示的流圖符號表示，稱為蝶形運算符號。 1210)()()(21NkkXWkXkXkN，1210)()()2(21NkkXWkXNkXkN，第4章 快速傅里葉變換(FFT) 圖4.2.1 DIT - FFT蝶式運算信號流圖(a) 蝶形流圖; (b)蝶形流圖簡化形式第4章 快速傅里葉變換(FFT) 經過一次奇偶抽取分

6、解后，N點DFT運算圖可以用圖4.2.2表示。 經過一次分解，可使運算量減少一半。 可以對N/2點DFT再作進一步分解，直到分解成N個1點DFT。共有M級蝶形運算1點DFT就是時域序列本身。第4章 快速傅里葉變換(FFT) 【例1】 將8點序列的DFT用基2時分FFT算法進行分解。解解 (1) 第一次分解。將8點序列x(n)分解為兩個4點序列x1(n)和x2(n)，x1(n)為偶數(shù)序列，x2(n)為奇數(shù)序列，即x1(n)=x(0)，x(2)，x(4)，x(6)x2(n)=x(1)，x(3)，x(5)，x(7)第4章 快速傅里葉變換(FFT) 4點序列x1(n)和x2(n)分別做4點DFT得到X

7、1(k)和X2(k)，由X1(k)和X2(k)通過蝶形構造獲得8點DFTX(k)。第一次分解的旋轉因子為 (k=0，1，2，3)，運算流圖如圖4.3.2所示。kW8第4章 快速傅里葉變換(FFT) 圖4.2.2 N點DFT基2時分第一次分解運算流圖(N=8)第4章 快速傅里葉變換(FFT) (2) 第二次分解。將兩個4點序列x1(n)和x2(n)分解為四個2點序列x3(n)、x4(n)、 和。)(3nx)(4nx)4(),0()2(),0()(113xxxxnx)6(),2()3(),1 ()(114xxxxnx)5(),1 ()2(),0()(223xxxxnx)7(),3()3(),1 (

8、)(224xxxxnx第4章 快速傅里葉變換(FFT) 圖4.2.3 N點DFT基2時分第二次分解運算流圖(N=8)第4章 快速傅里葉變換(FFT) (3) 第三次分解。將四個2點序列分解為8個1點序列，其中)4()1 ()(),0()0()(3635xxnxxxnx)6()1 ()(),2()0()(4847xxnxxxnx)5()1 ()(),1 ()0()(3635xxnxxxnx)7()1 ()(),3()0()(4847xxnxxxnx第4章 快速傅里葉變換(FFT) 圖4.2.4 N點DFT基2時分FFT運算流圖(N=8)?)4()0(04/xWxN) 7 () 3 (04/xWx

9、aN)5()1 (04/xWxaNaWabN02/第4章 快速傅里葉變換(FFT) 4.2.3 DIT-FFT算法與直接計算算法與直接計算DFT運算量的比較運算量的比較當N=2M 時，其運算流圖應有M級蝶形，M級運算總共需要的復數(shù)乘次數(shù)為lb22MNNCMNlbACN MNN復數(shù)加次數(shù)為第4章 快速傅里葉變換(FFT) 當N1時，N2(N/2) lbN，所以，DIT-FFT算法比直接計算DFT的運算次數(shù)大大減少。例如，N=210=1024時，8 .2045120576 048 1lb22NNN第4章 快速傅里葉變換(FFT) 圖4.2.5 DIT-FFT算法與直接計算DFT所需復數(shù)乘法次數(shù)的比

10、較曲線第4章 快速傅里葉變換(FFT) 4.2.4 DIT-FFT的運算規(guī)律及編程思想的運算規(guī)律及編程思想1 原位計算原位計算DIT-FFT的運算過程的規(guī)律。1、N=2M點的FFT共進行M級運算;2、同一級中，每個蝶形的兩個輸入數(shù)據(jù)只對計算本蝶形有用；3、在計算完一個蝶形后，所得輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲單元(數(shù)組元素)。 利用同一存儲單元存儲蝶形計算輸入、輸出數(shù)據(jù)的方法稱為原位(址)計算。第4章 快速傅里葉變換(FFT) 2 旋轉因子的變化規(guī)律旋轉因子的變化規(guī)律在N點DIT-FFT運算流圖中，每個蝶形都要乘以因子，稱其為旋轉因子，p為旋轉因子的指數(shù)。找出旋轉因子與運算級數(shù)的關系

11、。pNW個不同的旋轉因子。級共有數(shù)。第表示從左到右的運算級用1 -L2LL第4章 快速傅里葉變換(FFT) 對N=2M的一般情況，第L級的旋轉因子為3, 2, 1, 0 31, 0 20 1222/24/JWWWLJWWWLJWWWLJJNpNJJNpNJJNpNLLL時時時第4章 快速傅里葉變換(FFT) 因為所以(4.2.12) (4.2.13)用(4.2.12)和(4.2.13)式確定第L級運算的旋轉因子(實際編程序時，L為最外層循環(huán)變量)。L-12,0,1,2,21LpJNWWJMLMLMLN222212, 2, 1, 0122LJNJNpNJWWWLMMLLMJp2第4章 快速傅里葉

12、變換(FFT) 3 蝶形運算規(guī)律蝶形運算規(guī)律設序列x(n)經時域抽選(倒序)后，按圖4.2.4所示的次序(倒序)存入數(shù)組A中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應用原位計算，則蝶形運算可表示成如下形式：式中下標L表示第L級運算，AL(J)則表示第L級運算后的數(shù)組元素A(J)的值（即第L級蝶形的輸出數(shù)據(jù)）。而AL1(J)表示第L級運算前A(J)的值（即第L級蝶形的輸入數(shù)據(jù)）。11( )( )()pLLLNA JAJAJB W11()( )()pLLLNA JBAJAJB W12 0,1,21; 1,2,MLLpJJLM第4章 快速傅里葉變換(FFT) 4. 編程思想及程序框圖編程思想及程序框

13、圖 觀察圖4.2.4，歸納出對編程有用的運算規(guī)律：第L級中，每個蝶形的兩個輸入數(shù)據(jù)相距B=2L1個點；每級有B個不同的旋轉因子；同一旋轉因子對應著間隔為2L點的2ML個蝶形。三重循環(huán)程序實現(xiàn)DIT-FFT運算，程序框圖如圖4.2.6所示。 第4章 快速傅里葉變換(FFT) 圖4.2.6 DIT-FFT運算和程序框圖第4章 快速傅里葉變換(FFT) DIT-FFT算法運算流圖的輸出X(k)為自然順序。輸入序列不是按x(n)的自然順序排列，排序稱為序列x(n)的倒序(倒位)。在運算M級蝶形之前應先對序列x(n)進行倒序。5 序列的倒序序列的倒序規(guī)律如圖4.2.7所示。第4章 快速傅里葉變換(FFT

14、) 圖4.2.7 形成例序的樹狀圖（N=23）第4章 快速傅里葉變換(FFT) 表4.2.1 順序和倒序二進制數(shù)對照表第4章 快速傅里葉變換(FFT) 圖4.2.9所示的倒序的程序框圖中的虛線框內就是完成計算倒序值的運算流程圖。 第4章 快速傅里葉變換(FFT) 圖4.2.9 倒序程序框圖第4章 快速傅里葉變換(FFT) 第3章介紹的MATLAB函數(shù)fft是一個計算DFT的智能程序，如果計算點數(shù)N=2M，則自動按DIT-FFT快速算法計算。第4章 快速傅里葉變換(FFT) 4.2.5 頻域抽取法頻域抽取法FFT(DIF-FFT)在基2FFT算法中，頻域抽取法FFT也是一種常用的快速算法，簡稱D

15、TF- FFT。設序列x(n)長度為N=2M，首先將x(n)前后對半分開，得到兩個子序列，其DFT可表示為如下形式：10/2 110/2/2 1/2 1(/2)00/2 1/20( )DFT ( )( )( )( )( )2( )2NknNnNNknknNNnn NNNknk n NNNnnNkNknNNnX kx nx n Wx n Wx n WNx n Wx nWNx nWx nW第4章 快速傅里葉變換(FFT) 式中 將X(k)分解成偶數(shù)組與奇數(shù)組，當k取偶數(shù)(k=2m, m=0, 1, , N/21)時/21( 1)1kNkNkWk 偶數(shù)奇數(shù)，/220/2 1/20(2 )( )2(

16、)(4.2.14)2NmnNnNmnNnNXmx nx nWNx nx nW第4章 快速傅里葉變換(FFT) 當k取奇數(shù)(k=2m+1, m=0, 1, , N/21)時，/2 1(21)0/2 1/20(21)( )2( )2 (4.2.15)NnmNnNnnmNNnNXmx nx nWNx nx nWW令122102)()(2)()(21NnWNnxnxnxNnxnxnxnN，第4章 快速傅里葉變換(FFT) 將x1(n)和x2(n)分別代入(4.2.14)和(4.2.15)式，可得 (4.2.16)(4.2.16)式表明，X(k)按奇偶k值分為兩組，其偶數(shù)組是x1(n)的N/2點DFT，

17、奇數(shù)組則是x2(n)的N/2點DFT。x1(n)、x2(n)和x(n)之間的關系也可用圖4.2.10所示的蝶形運算流圖符號表示。圖4.2.11表示N=8時第一次分解的運算流圖。/2 11/20/2 12/20(2 )( )(21)( )NmnNnNmnNnXmx n WXmx n W第4章 快速傅里葉變換(FFT) 圖4.2.10 DTFFFT蝶形運算流圖符號第4章 快速傅里葉變換(FFT) 圖4.2.11 DIF-FFT第一次分解運算流圖（N=8）第4章 快速傅里葉變換(FFT) 由于N=2M，N/2仍然是偶數(shù)，繼續(xù)將N/2點DFT分成偶數(shù)組合奇數(shù)組。圖4.2.12表示N=8時第二次分解運算

18、流圖。 經過M1次分解，最后分解為2M1個兩點DFT，兩點DFT就是一個基本蝶形運算流圖。 當N=8，經兩次分解，便分解為四個兩點DFT。 N = 8的完整DIF-FFT運算流圖如圖4.2.13所示。第4章 快速傅里葉變換(FFT) 圖4.2.12 DIF-FFT第二次分解運算流圖（N = 8）第4章 快速傅里葉變換(FFT) 圖4.2.13 DIF-FFT運算流圖（N =8）第4章 快速傅里葉變換(FFT) 這種算法是對X(k)進行奇偶抽取分解的結果，所以稱之為頻域抽取法FFT。觀察圖4.2.13可知，DIF-FFT算法與DIT-TTF算法類似。不同的是DIF-FFT算法輸入序列為自然順序，

19、而輸出為倒序排列。 要對輸出數(shù)據(jù)進行倒序才能得到自然順序的X(k)。 另外，蝶形運算略有不同，DIT-FFT蝶形先乘后加(減)，而DIF-FFT蝶形先加(減)后相乘。第4章 快速傅里葉變換(FFT) 4.2.6 IDFT的高效算法的高效算法上述FFT算法流圖也可以用于計算IDFT。比較DFT和IDFT的運算公式：1010( )DFT ( )( )1( )IDFT ( )( )NknNnNknNkX kx nx n Wx nx nX k WN第4章 快速傅里葉變換(FFT) 只要將DFT運算式中的系數(shù)改變?yōu)椋詈蟪艘?/N，就是IDFT運算公式。 在DIT-FFT與DIF-FFT算法中的旋轉因子

20、改為，最后的輸出再乘以1/N就可以用來計算IDFT。流圖的輸入是X(k)，輸出就是x(n)。 原來的DIT-FFT改為IFFT后，稱為DIF-IFFT更合適；DIF-FFT改為IFFT后, 應稱為DIT-IFFT。knNWpNWpNWknNW第4章 快速傅里葉變換(FFT) 教材第教材第4章習題與上機題解答章習題與上機題解答快速傅里葉變換（FFT）是DFT的快速算法， 沒有新的物理概念。 1 如果某通用單片計算機的速度為平均每次復數(shù)乘需要4 s， 每次復數(shù)加需要1 s， 用來計算N=1024點DFT， 問直接計算需要多少時間。 用FFT計算呢？照這樣計算， 用FFT進行快速卷積對信號進行處理時

21、， 估計可實現(xiàn)實時處理的信號最高頻率。 第4章 快速傅里葉變換(FFT) 解解: 當N=1024=210時， 直接計算DFT的復數(shù)乘法運算次數(shù)為N2=10241024=1 048 576次復數(shù)加法運算次數(shù)為N（N1）=10241023=1 047 552次直接計算所用計算時間TD為TD=410610242+1 047 552106=5.241 856 s用FFT計算1024點DFT所需計算時間TF為第4章 快速傅里葉變換(FFT) 66F665 10lblb10210245 1010 1024 10 10230.72 msNTNN N 快速卷積時， 需要計算一次N點FFT（考慮到H(k)=DF

22、Th(n)已計算好存入內存）、 N次頻域復數(shù)乘法和一次N點IFFT。 所以， 計算1024點快速卷積的計算時間Tc約為第4章 快速傅里葉變換(FFT) cF2102471680 s4 1024 s65536 sTT 次復數(shù)乘計算時間所以， 每秒鐘處理的采樣點數(shù)（即采樣速率）s6102415 625 /65536 10F次 秒由采樣定理知， 可實時處理的信號最高頻率為smax156257.8125 kHz22Ff第4章 快速傅里葉變換(FFT) 應當說明， 實際實現(xiàn)時， fmax還要小一些。 這是由于實際中要求采樣頻率高于奈奎斯特速率， 而且在采用重疊相加法時， 重疊部分要計算兩次。 重疊部分長

23、度與h(n)長度有關， 而且還有存取數(shù)據(jù)和指令周期等消耗的時間。 2 如果將通用單片機換成數(shù)字信號處理專用單片機TMS320系列， 計算復數(shù)乘和復數(shù)加各需要10 ns。 請重復做上題。 解解: 與第1題同理。 直接計算1024點DFT所需計算時間TD為TD=1010910242+101091 047 552=20.961 28 ms第4章 快速傅里葉變換(FFT) 用FFT計算1024點DFT所需計算時間TF為99F8810 10lb10 10lb210241010 101024 1020.1536 msNTNNN快速卷積計算時間Tc約為cF3921024 2 0.1536 1010 1010

24、24 0.317 44 msTT次復數(shù)乘計算時間第4章 快速傅里葉變換(FFT) 可實時處理的信號最高頻率fmax為maxsc1110241 = 3.1158 MHz=1.6129 MHz222fFT由此可見， 用DSP專用單片機可大大提高信號處理速度。 所以， DSP在數(shù)字信號處理領域得到廣泛應用。 機器周期小于1 ns的DSP產品已上市， 其處理速度更高。 第4章 快速傅里葉變換(FFT) 3 已知X(k)和Y(k)是兩個N點實序列x(n)和y(n)的DFT， 希望從X(k)和Y(k)求x(n)和y(n)， 為提高運算效率， 試設計用一次N點IFFT來完成的算法。 解解: 因為x(n)和y

25、(n)均為實序列， 所以， X(k)和Y(k)為共軛對稱序列， jY(k)為共軛反對稱序列。 可令X(k)和jY(k)分別作為復序列F(k)的共軛對稱分量和共軛反對稱分量， 即F(k)=X(k)+jY(k)=Fep(k)+Fop(k)計算一次N點IFFT得到f(n)=IFFTF(k)=Ref(n)+j Imf(n)第4章 快速傅里葉變換(FFT) 由DFT的共軛對稱性可知Ref(n)=IDFTFep(k)=IDFTX(k)=x(n)j Imf(n)=IDFTFop(k)=IDFTjY(k)=jy(n)故1( ) ( )( )2x nf nfn1( ) ( )( )2jy nf nfn4 設x(

26、n)是長度為2N的有限長實序列， X(k)為x(n)的2N點DFT。 (1) 試設計用一次N點FFT完成計算X(k)的高效算法。 (2) 若已知X(k) ，試設計用一次N點IFFT實現(xiàn)求X(k)的2N點IDFT運算。第4章 快速傅里葉變換(FFT) 解解: 本題的解題思路就是DIT-FFT思想。（1） 在時域分別抽取偶數(shù)和奇數(shù)點x(n)， 得到兩個N點實序列x1(n)和x2(n)： x1(n)=x(2n) n=0, 1, , N1x2(n)=x(2n+1) n=0, 1, , N1根據(jù)DIT-FFT的思想， 只要求得x1(n)和x2(n)的N點DFT， 再經過簡單的一級蝶形運算就可得到x(n)

27、的2N點DFT。 因為x1(n)和x2(n)均為實序列， 所以根據(jù)DFT的共軛對稱性， 可用一次N點FFT求得X1(k)和X2(k)。 具體方法如下：第4章 快速傅里葉變換(FFT) 令 y(n)=x1(n)+jx2(n)Y(k)=DFTy(n) k=0, 1, , N1則*11ep*22ep1( )DFT ( )( ) ( )()21j( )DFTj( )( ) ( )()2X kx nYkY kYNkXkx nYkY kYNk2N點DFTx（n）=X（k）可由X1(k)和X2(k)得到122122( )( )( ) 0,1,1()( )( )kNkNX kX kWXkkNX kNX kWX

28、k第4章 快速傅里葉變換(FFT) 這樣， 通過一次N點IFFT計算就完成了計算2N點DFT。 當然還要進行由Y(k)求X1(k)、 X2(k)和X(k)的運算(運算量相對很少)。 （2） 與（1）相同， 設x1(n)=x(2n) n=0, 1, , N1x2(n)=x(2n+1) n=0, 1, , N1X1(k)=DFTx1(n) k=0, 1, , N1X2(k)=DFTx2(n) k=0, 1, , N1則應滿足關系式122122( )( )( ) 0,1,1()( )( )kNkNX kX kWXkkNX kNX kWXk第4章 快速傅里葉變換(FFT) 由上式可解出1221( )(

29、 )()2 0,1,2,11( )( )()2kNX kX kX kNkNXkX kX kN W由以上分析可得出運算過程如下： 由X(k)計算出X1(k)和X2(k): 1221( )( )()21( )( )()2kNX kX kX kNXkX kX kN W第4章 快速傅里葉變換(FFT) 由X1(k)和X2(k)構成N點頻域序列Y(k): Y(k)=X1(k)+jX2(k)=Yep(k)+Yop(k)其中， Yep(k)=X1(k), Yop(k)=jX2(k)， 進行N點IFFT， 得到y(tǒng)(n)=IFFTY(k)=Rey(n)+j Imy(n) n=0, 1, , N1由DFT的共軛對

30、稱性知*ep1*op21Re ( ) ( )( )DFT( )( )21jIm ( ) ( )( )DFT( )j( )2y ny ny nYkx ny ny ny nYkx n第4章 快速傅里葉變換(FFT) 由x1(n)和x2(n)合成x(n):12 2( )1 2nxnx nnxn偶數(shù)奇數(shù)，0n2N1在編程序實現(xiàn)時， 只要將存放x1(n)和x2(n)的兩個數(shù)組的元素分別依次放入存放x(n)的數(shù)組的偶數(shù)和奇數(shù)數(shù)組元素中即可。第4章 快速傅里葉變換(FFT) 5 分別畫出16點基2DIT-FFT和DIF-FFT運算流圖， 并計算其復數(shù)乘次數(shù)， 如果考慮三類碟形的乘法計算， 試計算復乘次數(shù)。

31、解解： 本題比較簡單， 仿照教材中的8點基2DIT-FFT和DIF-FFT運算流圖很容易畫出16點基2DIT-FFT和DIF-FFT運算流圖。 但畫圖占篇幅較大， 這里省略本題解答， 請讀者自己完成。第4章 快速傅里葉變換(FFT) 6* 按照下面的IDFT算法編寫MATLAB語言 IFFT程序， 其中的FFT部分不用寫出清單， 可調用fft函數(shù)。 并分別對單位脈沖序列、 矩形序列、 三角序列和正弦序列進行FFT和IFFT變換， 驗證所編程序。 *1( )IDFT( )DFT( )x nX kXkN第4章 快速傅里葉變換(FFT) 解解： 為了使用靈活方便， 將本題所給算法公式作為函數(shù)編寫ifft46.m如下： %函數(shù)ifft46.m%按照所給算法公式計算IFETfunction xn=ifft46(Xk，