第四章 快速傅里葉變換n

1、2n直接計算直接計算DFT的問題及改進(jìn)的途徑的問題及改進(jìn)的途徑 3 4 設(shè)復(fù)序列設(shè)復(fù)序列x(n) 長度為長度為N點(diǎn)，其點(diǎn)，其DFT為為10( )( )NnkNnX kx n Wk=0，N-1 （1）計算一個）計算一個X(k) 值的運(yùn)算量值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：復(fù)數(shù)乘法次數(shù)： N復(fù)數(shù)加法次數(shù)：復(fù)數(shù)加法次數(shù)： N15（2）計算全部）計算全部N個個X(k) 值的運(yùn)算量值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：復(fù)數(shù)乘法次數(shù)： N2復(fù)數(shù)加法次數(shù)：復(fù)數(shù)加法次數(shù)： N(N1)（3）對應(yīng)的實數(shù)運(yùn)算量）對應(yīng)的實數(shù)運(yùn)算量1100( )( )Re ( )Im ( )ReImNNnknknkNNNnnX kx n Wx njx nWj

2、W10Re ( ) ReIm ( ) ImNnknkNNnx nWx nWRe ( ) ImIm( ) RenknkNNjx nWx nW6一次復(fù)數(shù)乘法一次復(fù)數(shù)乘法： 4次實數(shù)乘法次實數(shù)乘法 2次實數(shù)加法次實數(shù)加法 一個一個X(k)數(shù)值數(shù)值 ： 4N次實數(shù)乘法次實數(shù)乘法2N+2(N-1)= 2(2N-1)次實數(shù)加法次實數(shù)加法 所以所以 整個整個N點(diǎn)點(diǎn)DFT運(yùn)算共需要：運(yùn)算共需要：N2(2N-1)= 2N(2N-1)實數(shù)乘法次數(shù)：實數(shù)乘法次數(shù)：4 N2實數(shù)加法次數(shù)：實數(shù)加法次數(shù)：7N點(diǎn)點(diǎn)DFT的復(fù)數(shù)乘法次數(shù)舉例的復(fù)數(shù)乘法次數(shù)舉例NN2NN22464404941612816384864256 65

3、 536 16256512 262 144 3210281024 1 048 576 結(jié)論結(jié)論：當(dāng)：當(dāng)N很大時，其運(yùn)算量很大，對實時性很強(qiáng)的信號處理很大時，其運(yùn)算量很大，對實時性很強(qiáng)的信號處理來說，要求計算速度快，因此需要改進(jìn)來說，要求計算速度快，因此需要改進(jìn)DFT的計算方法，以的計算方法，以大大減少運(yùn)算次數(shù)。大大減少運(yùn)算次數(shù)。 8 nkNW主要原理是利用系數(shù)主要原理是利用系數(shù) 的以下特性對的以下特性對DFT進(jìn)行分解：進(jìn)行分解： nkNW（1）對稱性）對稱性 ()nkNW()k N nNW（2）周期性）周期性 ()()nN kn kNnkNNNWWW（3）可約性）可約性 mnknkmNNWW/

4、nknk mNNmWW另外，另外，12/NNWkNNkNWW)2/(9101(2 )( )xrx r設(shè)設(shè)N2L，將，將x(n)按按 n 的奇偶分為兩組：的奇偶分為兩組： 2(21)( )xrx rr =0，1， 12N10( ) ( )( )NnkNnX kDFT x nx n W則則1010)()(NnnnkNNnnnkNWnxWnx為奇數(shù)為偶數(shù)11120)12(1202) 12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1010)()(NnnnkNNnnnkNWnxWnx為奇數(shù)為偶數(shù)式中，式中，X1(k

5、)和和X2(k)分別是分別是x1(n)和和x2(n)的的N/2的的DFT。另外，式中另外，式中k的取值范圍是：的取值范圍是：0，1， ，N/21 。12因此，因此， 只能計算出只能計算出X(k)的前一半值。的前一半值。12( )( )( )kNX kXkW Xk后一半后一半X(k) 值，值， N/2 ， N/2 1， ，N ？rkNW2(2)2r NkNW利用利用可得到可得到 1()2NXk2 1(2)120( )Nr NkNrx r W2 1120( )NrkNrx r W)(1kX同理可得同理可得22()( )2NXkXk13考慮到考慮到 kNkNNNkNNWWWW2)2(因此可得后半部分

6、因此可得后半部分X(k) )2()2()2(221NkXWNkXNkXNkN12( )( )( )kNX kX kW Xk由前半部分由前半部分X(k) )()(21kXWkXkNk=0，1， ，N/21k=0，1， ，N/2114222 N2N蝶形運(yùn)算式蝶形運(yùn)算式蝶形運(yùn)算信號蝶形運(yùn)算信號流圖符號流圖符號 因此，只要求出因此，只要求出2個個N/2點(diǎn)的點(diǎn)的DFT，即，即X1(k)和和X2(k)，再經(jīng)過，再經(jīng)過蝶形運(yùn)算就可求出全部蝶形運(yùn)算就可求出全部X(k)的值。運(yùn)算量由次復(fù)數(shù)乘減為的值。運(yùn)算量由次復(fù)數(shù)乘減為運(yùn)算量減少一半。運(yùn)算量減少一半。)()()2()()()(2121kXWkXNkXkXWkX

7、kXkNkN150NW1NW2NW3NW以以N=8為例，為例，分解為分解為2個個4點(diǎn)的點(diǎn)的DFT，然后做然后做8/2=4次蝶形運(yùn)次蝶形運(yùn)算即可求出所有算即可求出所有8點(diǎn)點(diǎn)X(k)的值。的值。16復(fù)數(shù)乘法次數(shù)： N2復(fù)數(shù)加法次數(shù)： N(N1)復(fù)數(shù)乘法次數(shù)： 2*(N/2)2+N/2=N2/2+N/2復(fù)數(shù)加法次數(shù)： 2*(N/2)(N/21)+2*N/2=N2/2nN點(diǎn) 17 由于N2L，因而N/2仍是偶數(shù) ，可以進(jìn)一步把每個N/2點(diǎn)子序列再按其奇偶部分分解為兩個N/4點(diǎn)的子序列。 以N/2點(diǎn)序列x1(r)為例 1314(2 )( )0,1,1(21)( )4xlx lNlxlxl則有 rkNNr

8、WrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401) 12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(42/3kXWkXkNk=0,1, 14N18且且13/24( )( )4kNNXkXkWXkk=0,1, 14N由此可見，一個由此可見，一個N/2點(diǎn)點(diǎn)DFT可分解成兩個可分解成兩個N/4點(diǎn)點(diǎn)DFT。 同理，也可對同理，也可對x2(n)進(jìn)行同樣的分解，求出進(jìn)行同樣的分解，求出X2(k)。 192013/40( )lkNlx l W02(0)(4)xW x0(0)(4)NxW x 對此例對此例N=8，最后剩下的是，

9、最后剩下的是4個個N/4= 2點(diǎn)的點(diǎn)的DFT，2點(diǎn)點(diǎn)DFT也可以由蝶形運(yùn)算來完成。以也可以由蝶形運(yùn)算來完成。以X3(k)為例。為例。 /4 133/40( )( )NlkNlXkx l Wk=0, 1即即03323(0)(0)(1)XxW x13323(1)(0)(1)XxW x12(0)(4)xW x0(0)(4)NxW x這說明，這說明，N=2M的的DFT可全部由蝶形運(yùn)算來完成?？扇坑傻芜\(yùn)算來完成。21N=8按時間抽取法按時間抽取法FFT信號流圖信號流圖 22由按時間抽取法FFT的信號流圖可知，當(dāng)N=2L時，共有 級蝶形運(yùn)算；每級都由 個蝶形運(yùn)算組成，而每個蝶形有 次復(fù)乘、 次復(fù)加，因

10、此每級運(yùn)算都需 次復(fù)乘和 次復(fù)加。 LN/2 N/2 12N23這樣這樣 級運(yùn)算總共需要：級運(yùn)算總共需要： L復(fù)數(shù)乘法： NNLN2log22復(fù)數(shù)加法： NNLN2log直接直接DFT算法運(yùn)算量算法運(yùn)算量 復(fù)數(shù)乘法： 復(fù)數(shù)加法： N2N(N1)直接計算直接計算DFT與與FFT算法的計算量之比為算法的計算量之比為MNNNNNM222log2log224NN2計算量之比M NN2計算量之比M 2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 048 5765 1202

11、04.83210288012.820484 194 30411 264372.464404919221.4NN2log2NN2log225n序列的逆序排列n同址運(yùn)算（原位運(yùn)算）n蝶形運(yùn)算兩節(jié)點(diǎn)間的距離n 的確定rNW26返回返回27)(01221)()(BINMMDECnnnnnn 由于由于 x(n) 被反復(fù)地按奇、偶分組，所以流圖輸被反復(fù)地按奇、偶分組，所以流圖輸入端的入端的排列不再是順序的，但仍有規(guī)律可循：排列不再是順序的，但仍有規(guī)律可循： 因為因為 N=2M ， 對于任意對于任意 n（0n N-1)，可以用可以用M個二進(jìn)制個二進(jìn)制碼表示為：碼表示為：10,01221nnnnnMM n 反

12、復(fù)按奇、偶分解時，即按二進(jìn)制碼的反復(fù)按奇、偶分解時，即按二進(jìn)制碼的“0” “1” 分解。分解。n序列的逆序排列2829自然順序自然順序 n二進(jìn)制數(shù)二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序順序數(shù)倒位序順序數(shù)0000000010011004201001023011110641000011510110156110011371111117 n30為目前計算機(jī)已具有的功能為目前計算機(jī)已具有的功能31返回返回32 某一列任何兩個節(jié)點(diǎn)某一列任何兩個節(jié)點(diǎn)k 和和j 的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算后，的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算后，得到結(jié)果為下一列得到結(jié)果為下一列k、j兩節(jié)點(diǎn)的節(jié)點(diǎn)變量，而和其他節(jié)點(diǎn)變量兩節(jié)點(diǎn)的節(jié)點(diǎn)變量，而

13、和其他節(jié)點(diǎn)變量無關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲單元，降低設(shè)備成本。無關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲單元，降低設(shè)備成本。)(kX)2(NkX)(1kX)(2kXkNW運(yùn)算前運(yùn)算前運(yùn)算后運(yùn)算后)2(NkA)(kA)2(NkA)(kA例例n同址運(yùn)算（原位運(yùn)算）33返回返回34以以N=8為例：為例：第一級蝶形，距離為：第一級蝶形，距離為：第二級蝶形，距離為：第二級蝶形，距離為：第三級蝶形，距離為：第三級蝶形，距離為：規(guī)律規(guī)律：對于共：對于共L級的蝶形而言，其級的蝶形而言，其m級蝶形運(yùn)算的節(jié)級蝶形運(yùn)算的節(jié) 點(diǎn)間的距離為點(diǎn)間的距離為12412mn蝶形運(yùn)算兩節(jié)點(diǎn)間的距離 35rNWrNWrNW返回返回r

14、NWrNW36 rNW以N=8為例：0,10224/lWWWWmllNrNm時，1 , 0,2422/lWWWWmlllNrNm時，3 , 2 , 1 , 0,382lWWWWmlllNrNm時，級：第LNM,212 , 2 , 1 , 0,12LlrNlWWLMLMLMLN2222LMLMMLMLlNlNjlNjlNrNWeeWW222222rNWn 的確定 37n算法原理算法原理 再把輸出再把輸出X(k)按按k的奇偶分組的奇偶分組先把輸入按先把輸入按n的順序分成前后兩半的順序分成前后兩半設(shè)序列長度為設(shè)序列長度為N=2L，L為整數(shù)為整數(shù) 前半子序列前半子序列x(n) 后半子序列后半子序列)2

15、(Nnx 0n12N0n12N3810( )( )NnkNnX kx n W由由DFT定義得定義得/2 110/2( )( )NNnknkNNnn Nx n Wx n W12/0)2(12/0)2()(NnkNnNNnnkNWNnxWnx12/02)2()(NnnkNkNNWWNnxnxk=0，1， ，N39由于由于 1222jNNjNNeeW/2 120( )( )()2NNknkNNnNX kx nx nWW所以所以kkNNW) 1(2則則 12/0)2() 1()()(NnnkNkWNnxnxkXk=0，1， ，N40然后按然后按k的奇偶可將的奇偶可將X(k)分為兩部分分為兩部分 221

16、krkrr=0，1， ，12N則式則式 12/0)2() 1()()(NnnkNkWNnxnxkX可轉(zhuǎn)化為可轉(zhuǎn)化為 nrNNnWNnxnxrX212/0)2()()2(12/02/)2()(NnnrNWNnxnx)12(12/0)2()() 12(rnNNnWNnxnxrXnrNnNNnWWNnxnx212/0)2()(41/2 1/20(2 )( )()2NnrNnNXrx nx nW令令 nNWNnxnxnxNnxnxnx)2()()()2()()(21n=0，1， ，12N代入代入 /2 120(21) ( )()2NnnrNNnNXrx nx nWWnrNNnnrNNnWnxrWnxr

17、2120221201)() 12()()2(r=0，1， ，12N可得可得為為2個個N/2點(diǎn)的點(diǎn)的DFT，合起來正好是，合起來正好是N點(diǎn)點(diǎn)X(k)的值。的值。42nNWNnxnxnxNnxnxnx)2()()()2()()(21將將稱為蝶形運(yùn)算稱為蝶形運(yùn)算與時間抽選基與時間抽選基2FFT算法中的蝶形運(yùn)算符號略有不同。算法中的蝶形運(yùn)算符號略有不同。43例例 按頻率抽取，將按頻率抽取，將N點(diǎn)點(diǎn)DFT分解為兩個分解為兩個N/2點(diǎn)點(diǎn)DFT的組合的組合(N=8) 44 與時間抽取法的推導(dǎo)過程一樣，由于與時間抽取法的推導(dǎo)過程一樣，由于 N=2L，N/2仍然是仍然是一個偶數(shù)，因而可以將每個一個偶數(shù)，因而可以

18、將每個N/2點(diǎn)點(diǎn)DFT的輸出再分解為偶數(shù)組的輸出再分解為偶數(shù)組與奇數(shù)組，這就將與奇數(shù)組，這就將N/2點(diǎn)點(diǎn)DFT進(jìn)一步分解為兩個進(jìn)一步分解為兩個N/4點(diǎn)點(diǎn)DFT。 N=845按頻率抽取的按頻率抽取的FFT（N=8）信號流圖）信號流圖 返回返回 1 10NW2NWx(0)x(1)x(2)x(3) 1 1x(4)x(5)x(6)x(7)0NW1NW2NW3NWX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW 1 1 1 1 1 1 1 10NW0NW0NW0NW46 由由P31圖圖 與與P42圖圖 相比較，初看起來，相比較，初看起來，DIF法與法與DIT法的區(qū)別是：法的區(qū)

19、別是： DIF輸入是自然順序，輸入是自然順序， 輸出是倒位序的，輸出是倒位序的， 這與這與DIT法正好相反。法正好相反。 但但這不是實質(zhì)性的區(qū)別，因為這不是實質(zhì)性的區(qū)別，因為DIF法與法與DIT法一樣，都可將輸入或輸出法一樣，都可將輸入或輸出進(jìn)行重排，使二者的輸入或輸出順序變成自然順序或倒位序順序。進(jìn)行重排，使二者的輸入或輸出順序變成自然順序或倒位序順序。 DIF的基本蝶形（的基本蝶形（ P42圖圖 ）與）與DIT的基本蝶形（的基本蝶形（ P31 ）則有所不同，）則有所不同，這才這才是實質(zhì)的不同，是實質(zhì)的不同，DIF的復(fù)數(shù)乘法只出現(xiàn)在減法之后，的復(fù)數(shù)乘法只出現(xiàn)在減法之后，DIT則是先則是先作復(fù)

20、乘后再作加減法。作復(fù)乘后再作加減法。 47 但是，但是，DIF與與DIT就運(yùn)算量來說則是相同的，即都有就運(yùn)算量來說則是相同的，即都有M級級（列）運(yùn)算，（列）運(yùn)算， 每級運(yùn)算需每級運(yùn)算需N/2個蝶形運(yùn)算來完成，總共需要個蝶形運(yùn)算來完成，總共需要 次次復(fù)乘與復(fù)乘與 次復(fù)加次復(fù)加 ，DIF法與法與DIT法都可進(jìn)行原位運(yùn)法都可進(jìn)行原位運(yùn)算。算。 頻率抽取頻率抽取FFT算法的輸入是自然順序，輸出是倒位序的。算法的輸入是自然順序，輸出是倒位序的。 因此因此運(yùn)算完畢后，運(yùn)算完畢后， 要通過變址計算將倒位序轉(zhuǎn)換成自要通過變址計算將倒位序轉(zhuǎn)換成自然位序，然后再輸然位序，然后再輸出。出。 轉(zhuǎn)換方法與時間抽取法相同

21、。轉(zhuǎn)換方法與時間抽取法相同。 NN2log2NNLN2log48 如果將如果將DIT的基本蝶形加以轉(zhuǎn)置，就得到的基本蝶形加以轉(zhuǎn)置，就得到DIF的基本蝶形的基本蝶形; 反過來反過來, 將將DIF的基本蝶形加以轉(zhuǎn)置，的基本蝶形加以轉(zhuǎn)置， 就得到就得到DIT的基本蝶形，的基本蝶形， 因而因而DIT法與法與DIF法的基本蝶形是互為轉(zhuǎn)置的。法的基本蝶形是互為轉(zhuǎn)置的。 按照轉(zhuǎn)置定按照轉(zhuǎn)置定理，理， 兩個流圖的輸入兩個流圖的輸入-輸出特性必然相同。轉(zhuǎn)置就是將流圖輸出特性必然相同。轉(zhuǎn)置就是將流圖的所有支路方向都反向，并且交換輸入與輸出，但節(jié)點(diǎn)變量的所有支路方向都反向，并且交換輸入與輸出，但節(jié)點(diǎn)變量值不交換。值

22、不交換。 因此可以說，有多少種按時間抽取的因此可以說，有多少種按時間抽取的FFT流圖流圖就存在多就存在多少種按頻率抽取的少種按頻率抽取的FFT流圖。頻率抽取法與時間抽流圖。頻率抽取法與時間抽取法是兩種等價的取法是兩種等價的FFT運(yùn)算。運(yùn)算。 1 10NW2NWx(0)x(1)x(2)x(3) 1 1x(4)x(5)x(6)x(7)0NW1NW2NW3NWX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW 1 1 1 1 1 1 1 10NW0NW0NW0NW49n頻率抽取法輸入是自然順序，輸出是倒位序的；時間抽取法正好相反。n頻率抽取法的基本蝶形與時間抽取法的基本蝶形

23、有所不同。n頻率抽取法運(yùn)算量與時間抽取法相同。n頻率抽取法與時間抽取法的基本蝶形是互為轉(zhuǎn)置的。 504.5 N為復(fù)合數(shù)的為復(fù)合數(shù)的FFT算法算法 上面討論的是序列的點(diǎn)數(shù)上面討論的是序列的點(diǎn)數(shù)N為為2的冪次的冪次(即即N=2M)情況下，按時間抽情況下，按時間抽取和按頻率抽取的基取和按頻率抽取的基 -2FFT算法的基本原理。這種基算法的基本原理。這種基-2FFT算法在實算法在實際中使用得最多，因為它的程序簡單，效率高，使用方便。但實際上際中使用得最多，因為它的程序簡單，效率高，使用方便。但實際上無法保證總是處理長度為無法保證總是處理長度為2的整數(shù)冪次的序列。若不滿足的整數(shù)冪次的序列。若不滿足N=2

24、M，可將，可將x(n)增補(bǔ)一些零值點(diǎn)，以使增補(bǔ)一些零值點(diǎn)，以使N增長到最鄰近的一個增長到最鄰近的一個2M數(shù)值。數(shù)值。 有限長序有限長序列補(bǔ)零之后，并不影響其頻譜列補(bǔ)零之后，并不影響其頻譜X(ej)，只不過其頻譜的抽樣點(diǎn)數(shù)增加，只不過其頻譜的抽樣點(diǎn)數(shù)增加了，所造成的結(jié)果是增加了計算量而已。了，所造成的結(jié)果是增加了計算量而已。 但是，有時計算量增加太多，但是，有時計算量增加太多，浪浪費(fèi)較大。費(fèi)較大。 51 例如例如, x(n)的點(diǎn)數(shù)的點(diǎn)數(shù)N=300，則須補(bǔ)到，則須補(bǔ)到N=29=512，要補(bǔ)，要補(bǔ)212個零值點(diǎn)，個零值點(diǎn)，因而人們才研究因而人們才研究N2M時的時的FFT算法。算法。 若若N是一個復(fù)合

25、數(shù)，即它可以分解成一些因子的乘積，則可以用是一個復(fù)合數(shù)，即它可以分解成一些因子的乘積，則可以用FFT的一般算法，即混合基的一般算法，即混合基FFT算法，如庫利算法，如庫利-圖基（圖基（CooleyTukey）算法，而基）算法，而基 -2 算法只是這種一般算法的特例。算法只是這種一般算法的特例。 這里，這里， 我們我們不作詳細(xì)介紹，不作詳細(xì)介紹， 感興趣的話可參考文獻(xiàn)感興趣的話可參考文獻(xiàn)1。 總之，不管采用什么方法，總之，不管采用什么方法，計算計算DFT的高效算法是把計算長度為的高效算法是把計算長度為N的序列的的序列的DFT逐逐次分解成計算長度較短序列的次分解成計算長度較短序列的DFT。這是很多

26、高效。這是很多高效算法的標(biāo)準(zhǔn)方法和基本原理。算法的標(biāo)準(zhǔn)方法和基本原理。 52 例例4-5 已知信號已知信號x(t)=0.15 sin(2f1t)+sin(2f2t)-0.1sin(2f3t)，其中，其中, f1=1Hz， f2=2 Hz，f3=3Hz。 x(t)包含三個正弦波，但從時域波形圖包含三個正弦波，但從時域波形圖3-26（a）來看，似乎是一個正弦信號，很難看到小信號的存在，）來看，似乎是一個正弦信號，很難看到小信號的存在， 因為因為它被大信號所掩蓋。它被大信號所掩蓋。 取取fs=32 Hz作頻譜分析。作頻譜分析。 解解 因因fs=32 Hz，故，故 nnnnTxnx326sin1 .

27、0324sin322sin15. 0)()(53圖圖 已知信號的時域圖和幅頻圖已知信號的時域圖和幅頻圖 2 10120102030(a)(b)0051015816nkx(n)|X(k)|54 該信號為周期信號，其周期為該信號為周期信號，其周期為N=32?，F(xiàn)對?，F(xiàn)對x(n)作作32點(diǎn)的離散傅點(diǎn)的離散傅里葉變換里葉變換(DFT)，其幅度特性，其幅度特性|X(k)|如如下頁下頁圖所示。圖中僅給出了圖所示。圖中僅給出了k=0, 1, , 15的結(jié)果。的結(jié)果。k=16, 17, , 32的結(jié)果可由的結(jié)果可由|X(N-k)|=|X(k)|得出。因得出。因N=32, 故頻率分辨率故頻率分辨率F=fs/N=1

28、Hz。 由由下頁下頁圖圖 (b) 可知可知, k=1, 2, 3所對應(yīng)所對應(yīng)的頻譜即為頻率的頻譜即為頻率f1=1 Hz，f2=2 Hz, f3=3Hz的正弦波所對應(yīng)的頻譜。的正弦波所對應(yīng)的頻譜。幅頻圖如幅頻圖如下頁下頁圖圖(b), 該圖中小信號成該圖中小信號成分清楚顯示出來。分清楚顯示出來。 可見小信號可見小信號成分在時域中很難辨識而在頻域中容易識別。成分在時域中很難辨識而在頻域中容易識別。 幅度：幅度：2Hz，16*(2/fs)=1 1Hz，2.4*(2/fs)=0.15 3Hz，1.6*(2/fs)=0.155 如果我們事先不知道信號的最高頻率，可以根據(jù)信號的時域波形如果我們事先不知道信號

29、的最高頻率，可以根據(jù)信號的時域波形圖來估計它。例如，圖來估計它。例如， 某信號的波形如圖某信號的波形如圖 3-23 所示。所示。 先找出相鄰的波先找出相鄰的波峰與波谷之間的距離，如圖中峰與波谷之間的距離，如圖中t1，t2，t3，t4。 然后，選出其中最小的然后，選出其中最小的一個如一個如t4。這里。這里, t4可能就是由信號的最高頻率分量形成的?？赡芫褪怯尚盘柕淖罡哳l率分量形成的。 峰與峰與谷之谷之間的距離就是周期的一半。間的距離就是周期的一半。 因此，最高頻率為因此，最高頻率為 )(214Hztfh知道道fh后就能確定采樣頻率后就能確定采樣頻率 hsff256估算信號最高頻率估算信號最高頻率

30、fh tt3t4ot1t2x(t)57MN 2IDFT公式公式 10)(1NknkNWkXNkXIDFTnxDFT公式公式 nkNNnWnxnxDFTkX10)()()(比較可以看出，比較可以看出， nkNWnkNWMN211IDFT多出多出M個個1/2可分解到可分解到M級蝶形運(yùn)算中。級蝶形運(yùn)算中。585910)(1)(NknkNWkXNkXIDFT10)(1NknkNWkXN10)(1NknkNWkXN1( )( )IFFT X kFFT XkN( )X k 求共軛( )XkFFT 求( )FFT Xk( )FFT XkN 除以( )x n 求共軛60n用用FFT進(jìn)行譜分析的進(jìn)行譜分析的Ma

31、tlab實現(xiàn)實現(xiàn)n用用CZT進(jìn)行譜分析的進(jìn)行譜分析的Matlab實現(xiàn)實現(xiàn)n在在Matlab中使用的線性調(diào)頻中使用的線性調(diào)頻z變換函數(shù)為變換函數(shù)為czt，其，其調(diào)用格式為調(diào)用格式為nX= czt(x, M, W, A)n其中，其中，x是待變換的時域信號是待變換的時域信號x(n)，其長度為，其長度為N，M是變是變換的長度，換的長度，W確定變換的步長，確定變換的步長，A確定變換的起點(diǎn)。若確定變換的起點(diǎn)。若M= N，A= 1，則，則CZT變成變成DFT。61例例4.6 設(shè)模擬信號設(shè)模擬信號 ，以，以 t= 0.01n (n=0: N-1) 進(jìn)行取樣，試用進(jìn)行取樣，試用fft函數(shù)對其做頻譜分析。函數(shù)對其

32、做頻譜分析。N分別為：分別為：(1) N=45；(2) N=50；(3) N=55；(2) N=60。 ( )2sin(4)5cos(8)x ttt程序清單如下程序清單如下 %計算計算N=45的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=45;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y)title(FFT N=45)62%計算計算N=50的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=50;n=0:N-1;t=0.01*n;q=n*

33、2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y)title(FFT N=50)63%計算計算N=55的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=55;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,3)plot(q,abs(y)title(FFT N=55)64%計算計算N=60的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=60;n=0:N-1

34、;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,4)plot(q,abs(y)title(FFT N=60)6502468050100150FFT N=4502468050100150FFT N=5002468050100150FFT N=5502468050100150FFT N=60從圖中可以看出，這幾種情況下均有較好的精度。從圖中可以看出，這幾種情況下均有較好的精度。 66分析：由分析：由t=0.01n進(jìn)行取樣可得，采樣頻率進(jìn)行取樣可得，采樣頻率fs=100Hz。而連續(xù)。

35、而連續(xù)信號的最高模擬角頻率為信號的最高模擬角頻率為8 ，由，由2 f可得，最高頻可得，最高頻率為率為8 /2 =4Hz。因此，滿足采樣定理的要求。因此，滿足采樣定理的要求。 采樣序列為采樣序列為( )2cos(4)5cos(8)x nTnTn48( )2cos5cos100100 x nnn即即為周期序列，周期為周期序列，周期N=50。將程序中將程序中plot改為改為stem函數(shù)，則可以更清楚地看出頻譜。函數(shù)，則可以更清楚地看出頻譜。6702468050100150FFT N=4502468050100150FFT N=5002468050100150FFT N=550246805010015

36、0FFT N=6068 信號消噪信號消噪 假設(shè)信號在傳輸過程中，受到噪聲的干擾，則在接收假設(shè)信號在傳輸過程中，受到噪聲的干擾，則在接收端得到的信號由于受到噪聲的干擾，信號將難以辨識。用端得到的信號由于受到噪聲的干擾，信號將難以辨識。用FFT方法消噪就是對含噪信號的頻譜進(jìn)行處理，將噪聲所方法消噪就是對含噪信號的頻譜進(jìn)行處理，將噪聲所在頻段的在頻段的X(k)值全部置零后，再對處理后的值全部置零后，再對處理后的X(k)進(jìn)行離散傅進(jìn)行離散傅里葉反變換（里葉反變換（IFFT）可得原信號的近似結(jié)果。）可得原信號的近似結(jié)果。這種方法要這種方法要求噪聲與信號的頻譜不在同一頻段，求噪聲與信號的頻譜不在同一頻段，

37、 否則，否則， 將很難處理。將很難處理。 69 例例 將上述消噪原理用于語音消噪。圖將上述消噪原理用于語音消噪。圖4-30是一個實際例子，語音是一個實際例子，語音信號受到了強(qiáng)烈的嘯叫噪聲干擾，無法聽清語音意思，信號受到了強(qiáng)烈的嘯叫噪聲干擾，無法聽清語音意思， 如圖如圖4-30(a)， 信號淹沒在噪聲中（信噪比只有信號淹沒在噪聲中（信噪比只有-10dB）。試用）。試用FFT方法消噪。方法消噪。 先作先作FFT分析，得到其功率譜如圖分析，得到其功率譜如圖4-30(b)，可見在頻率，可見在頻率2.5 kHz附近有一極強(qiáng)分量。附近有一極強(qiáng)分量。 這就是嘯叫噪聲干擾。圖這就是嘯叫噪聲干擾。圖4-30(b

38、)中頻率在中頻率在30800Hz范圍是語音信號。對頻譜進(jìn)行修正，去除噪聲頻段，即將大于范圍是語音信號。對頻譜進(jìn)行修正，去除噪聲頻段，即將大于2.5 kHz部分的部分的X(k)值全部置為零，值全部置為零， 圖圖4-30(c)是去噪后的功率譜。再是去噪后的功率譜。再由反變換（由反變換（IFFT）重構(gòu)信號得到原語音信號如圖）重構(gòu)信號得到原語音信號如圖4-30(d)。這時信噪。這時信噪比為比為14dB，提高了，提高了24dB。 這就是早期的數(shù)字式錄音音樂中所采用的這就是早期的數(shù)字式錄音音樂中所采用的消噪方消噪方法。法。 70圖 4-30 語音信號消噪過程(a) 信號淹沒在嘯叫噪聲中； (b) 信號與噪

39、聲的功率譜； (c) 去噪后的功率譜； (d) 重構(gòu)原語音信號71n離散哈特萊變換定義離散哈特萊變換定義 n 設(shè)設(shè)x(n),n=0,1,N-1，為一實序列，其，為一實序列，其DHT定義為定義為1, 2 , 1 , 0)2()()()(10NkknNcasnxnxDHTkXNnH式中，cas()=cos+sin。其逆變換(IDHT)為1, 2 , 1 , 0)2()(1)()(10NnknNcaskXNkXIDHTnxNkHH72已知已知x(n)的的DHT，則，則DFT可用下式求得可用下式求得： 11( )( )()( )()22HHHHX kXkXNkj XkXNk DHT與與DFT之間的關(guān)系

40、之間的關(guān)系 731）實數(shù)序列的）實數(shù)序列的FFT以上討論的FFT算法都是復(fù)數(shù)運(yùn)算,包括序列x(n)也認(rèn)為是復(fù)數(shù),但大多數(shù)場合,信號是實數(shù)序列,任何實數(shù)都可看成虛部為零的復(fù)數(shù),例如,求某實信號x(n)的復(fù)譜,可認(rèn)為是將實信號加上數(shù)值為零的虛部變成復(fù)信號(x(n)+j0),再用FFT求其離散傅里葉變換。這種作法很不經(jīng)濟(jì),因為把實序列變成復(fù)序列,存儲器要增加一倍,且計算機(jī)運(yùn)行時,即使虛部為零,也要進(jìn)行涉及虛部的運(yùn)算,浪費(fèi)了運(yùn)算量。合理的解決方法是利用復(fù)數(shù)據(jù)FFT對實數(shù)據(jù)進(jìn)行有效計算,下面介紹兩種方法。74 （1）用 一個N點(diǎn)FFT同時計算兩個N點(diǎn)實序列的DFT 設(shè)x (n)、y (n)是彼此獨(dú)立的兩

41、個N點(diǎn)實序列,且 X (k)=DFTx (n) ， Y (k)=DFTy(n) 則X (k)、Y(k)可通過一次FFT運(yùn)算同時獲得。首先將x (n)、y(n)分別當(dāng)作一復(fù)序列的實部及虛部,令 g(n)=x (n)+jy(n)75 kjGkGkjYkjXkYkXkYkjYkjXkXkjYkXkGirriiririr kNGkGjkNGkGkNGkGkXngiirr212121ReDFT*通過通過g(n)的的FFT運(yùn)算結(jié)果運(yùn)算結(jié)果G(k),由上式可得到由上式可得到X (k) 的值。的值。 通過通過FFT運(yùn)算可獲得運(yùn)算可獲得g(n)的的DFT值值利用離散傅里葉變換的共軛對稱性利用離散傅里葉變換的共軛

42、對稱性76 kNGkGjkNGkGkNGkGkjYngjiirr212121ImDFT* kNGkGjkNGkGkYrrii2121作一次作一次點(diǎn)復(fù)序列的點(diǎn)復(fù)序列的FFT，再通過加、減法運(yùn)算就可以，再通過加、減法運(yùn)算就可以將將X(k)與與Y(k)分離出來。顯然，這將使運(yùn)算效率提高一倍。分離出來。顯然，這將使運(yùn)算效率提高一倍。77 設(shè)設(shè)x(n)是是2N點(diǎn)的實序列點(diǎn)的實序列,現(xiàn)人為地將現(xiàn)人為地將x(n)分為偶數(shù)組分為偶數(shù)組x1(n)和奇數(shù)組和奇數(shù)組x2(n) x1(n)=x(2n) n=0,1,N-1 x2(n)=x(2n+1) n=0,1,N-1然后將然后將x1(n)及及x2(n)組成一個復(fù)序列

43、：組成一個復(fù)序列： y(n)=x1(n)+jx2(n)（2） 用一個用一個N點(diǎn)的點(diǎn)的FFT運(yùn)算獲得一個運(yùn)算獲得一個2N點(diǎn)實序列的點(diǎn)實序列的DFT78通過通過N點(diǎn)點(diǎn)FFT運(yùn)算可得到：運(yùn)算可得到： Y(k)=X1(k)+jX2(k) ，N點(diǎn)點(diǎn)根據(jù)前面的討論，得到根據(jù)前面的討論，得到)()(2)()()(21)(*2*1kNYkYjkXkNYkYkX79 為求為求 2N 點(diǎn)點(diǎn) x（n）所對應(yīng)）所對應(yīng) X（k），需求出），需求出 X（k）與）與 X1（k）、）、X2（k）的關(guān)系）的關(guān)系 ： 10)12(21022120) 12()2()()(NnknNNnnkNNnnkWnxWnxWnxkX10210

44、) 12()2(NnnkNkNNnnkNWnxWWnx 1011022101011)12()()()2()()(NnnkNNnnkNNnNnnkNnkNWnxWnxkXWnxWnxkX801）由由x1(n)及及x2(n)組成復(fù)序列，經(jīng)組成復(fù)序列，經(jīng)FFT運(yùn)算求得運(yùn)算求得 Y(k)，2）利用共軛對稱性求出）利用共軛對稱性求出 X1(k)、X2(k)，3）最后利用上式求出）最后利用上式求出 X（k）, 達(dá)到用一個達(dá)到用一個N點(diǎn)的點(diǎn)的FFT計算計算一個一個2N點(diǎn)的實序列的點(diǎn)的實序列的DFT的目的。的目的。 X(k)=X1(k)+W2Nk X2(k)所以所以81 2) 線性卷積的線性卷積的FFT算法算

45、法 線性卷積是求離散系統(tǒng)響應(yīng)的主要方法之一線性卷積是求離散系統(tǒng)響應(yīng)的主要方法之一,許多重要應(yīng)許多重要應(yīng)用都建立在這一理論基礎(chǔ)上用都建立在這一理論基礎(chǔ)上,如卷積濾波等。如卷積濾波等。 以前曾討論了用圓周卷積計算線性卷積的方法歸納如下以前曾討論了用圓周卷積計算線性卷積的方法歸納如下: 將長為將長為N2的序列的序列x(n)延長到延長到L,補(bǔ)補(bǔ)L-N2個零，個零， 將長為將長為N1的序列的序列h(n)延長到延長到L,補(bǔ)補(bǔ)L-N1個零，個零， 如果如果LN1+N2-1,則圓周卷積與線性卷積相等則圓周卷積與線性卷積相等,此時此時,可用可用FFT計算線性卷積，方法如下計算線性卷積，方法如下: 82a.計算計

46、算X(k)=FFTx(n)b. 求求H(k)=FFTh(n)c. 求求Y(k)=H(k)X(k) k=0L-1d. 求求y(n)=IFFTY(k) n=0L-1 可見，只要進(jìn)行二次可見，只要進(jìn)行二次FFT,一次一次IFFT就可就可完成線性卷積計算。完成線性卷積計算。 計算表明計算表明,L32時時,上述計算線性卷積的上述計算線性卷積的方法比直接計算線卷積有明顯的優(yōu)越性方法比直接計算線卷積有明顯的優(yōu)越性,因此因此,也稱上述循環(huán)卷積方法為快速卷積法。也稱上述循環(huán)卷積方法為快速卷積法。83 上述結(jié)論適用于上述結(jié)論適用于 x(n)、h(n) 兩序列長度比較接近兩序列長度比較接近或相等的情況或相等的情況,

47、如果如果x(n)、h(n) 長度相差較多長度相差較多,例如例如, h(n) 為某濾波器的單位脈沖響應(yīng)為某濾波器的單位脈沖響應(yīng),長度有限長度有限,用來處理用來處理一個很長的輸入信號一個很長的輸入信號 x(n),或者處理一個連續(xù)不斷的或者處理一個連續(xù)不斷的信號，按上述方法信號，按上述方法, h(n) 要補(bǔ)許多零再進(jìn)行計算要補(bǔ)許多零再進(jìn)行計算,計算計算量有很大的浪費(fèi)，或者根本不能實現(xiàn)。量有很大的浪費(fèi)，或者根本不能實現(xiàn)。 為了保持快速卷積法的優(yōu)越性為了保持快速卷積法的優(yōu)越性,可將可將 x(n) 分為許分為許多段多段,每段的長度與每段的長度與 h(n) 接近接近 , 處理方法有兩種：處理方法有兩種：84

48、h(n)x(n)85 假定假定 xi(n) 表示表示 x（n）序列的第）序列的第i段段 ：01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnhnxny)()(*)()(*)()()(*)()(nhnxnyii其中其中 于是輸出可分解為：于是輸出可分解為： 則輸入序列可表為：則輸入序列可表為： 86 1）先對）先對 h（n）及）及 xi（n）補(bǔ)零，補(bǔ)到具有）補(bǔ)零，補(bǔ)到具有N點(diǎn)長度，點(diǎn)長度，N=N1+N2-1。 一般選一般選 N=2M。 2）用基）用基2 FFT計算計算 yi（n）=xi（n）*h（n）。）。 3）重疊部分相加構(gòu)成最后的輸出序列。）重疊部分相加構(gòu)

49、成最后的輸出序列。 iinyny)()(由于由于 yi（n）的長度為）的長度為N，而，而xi（n）的長度為）的長度為N2，因此相鄰，因此相鄰兩段兩段 yi（n）序列必然有）序列必然有N-N2=N1-1點(diǎn)發(fā)生重疊。點(diǎn)發(fā)生重疊。87計算步驟：計算步驟：a. 事先準(zhǔn)備好濾波器參數(shù)事先準(zhǔn)備好濾波器參數(shù) H（k）=DFTh(n)，N點(diǎn)點(diǎn)b. 用用N點(diǎn)點(diǎn)FFT計算計算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用用N點(diǎn)點(diǎn)IFFT求求 yi（n）=IDFTYi(k)e. 將重疊部分相加將重疊部分相加8889 這種方法和第一種方法稍有不同，即將上面分段這種方法和第一種方法稍有不同，即

50、將上面分段序列中補(bǔ)零的部分不是補(bǔ)零，而是保留原來的輸入序列中補(bǔ)零的部分不是補(bǔ)零，而是保留原來的輸入序列值，這時，如利用序列值，這時，如利用DFT實現(xiàn)實現(xiàn)h(n)和和xi(n)的循環(huán)的循環(huán)卷積，則每段卷積結(jié)果中有卷積，則每段卷積結(jié)果中有N1-1個點(diǎn)不等于線性卷個點(diǎn)不等于線性卷積值需舍去。積值需舍去。 重疊保留法與重疊相加法的計算量差不多，但省重疊保留法與重疊相加法的計算量差不多，但省去了重疊相加法最后的相加運(yùn)算。去了重疊相加法最后的相加運(yùn)算。 9091 4）用）用FFT計算相關(guān)函數(shù)計算相關(guān)函數(shù) 相關(guān)的概念很重要，互相關(guān)運(yùn)算廣泛應(yīng)用于信號相關(guān)的概念很重要，互相關(guān)運(yùn)算廣泛應(yīng)用于信號分析與統(tǒng)計分析，如通過相關(guān)函數(shù)峰值的檢測測量兩個分析與統(tǒng)計分析，如通過相關(guān)函數(shù)峰值的檢測測量兩個信號的時延差等。信號的時延差等。 兩個長為兩個長為N的實離散時間序列的實離散時間序列 x（n）與）與y（n）的）的互相關(guān)函數(shù)定義為互相關(guān)函數(shù)定義為 1