1、MEASUREMENTINFORMATION SIGNAL ANALYSIS IN MECHANICAL ENGINEERING 機械工程測試機械工程測試信息信息信號分析信號分析 機械科學與工程學院機械科學與工程學院 機械電子信息工程系機械電子信息工程系課件資料下載：課件資料下載：郵箱地址：郵箱地址：密碼：密碼：注意下載時不要刪除原始文件注意下載時不要刪除原始文件 第六章第六章 數字信號分析數字信號分析(I)DFT與與FFT46-5 現代譜分析方法現代譜分析方法-最大熵譜估計最大熵譜估計6-3 FFT6-4 譜分析與譜估計譜分析與譜估計6-2 DFT6-1 模擬信號離散化模擬信號離散化目錄目錄

2、56-3 FFT算法6.3.1、DFT的計算工作量 1, 1 , 0 ,)()(10NkWnxkXNnnkN101, 1 , 0 ,)(1)(NknkNNnWkXNnxqFFT背景背景2jNNWe兩者的差別僅在指數的符號和因子兩者的差別僅在指數的符號和因子1/N.6q通常通常x(n)和和WNnk都是復數，所以計算一個都是復數，所以計算一個x(k)的值的值需要需要N次次復數乘法運算，和復數乘法運算，和N-1次次復數加法運算。復數加法運算。那么，所有的那么，所有的X(k)就要就要N2次復數乘法運算，次復數乘法運算，N(N-1)次復數加法運算。當次復數加法運算。當N很大時，運算量將是驚人的很大時，運

3、算量將是驚人的,如如N=1024，則要完成，則要完成1048576次次(一百多萬次一百多萬次)運算。運算。難以做到實時處理。難以做到實時處理。q計算一個計算一個X(k)的值的計算量，如的值的計算量，如X(1)，k=11210) 1()2() 1 ()0() 1 (NNNNNWNxWxWxWxX76.3.2、改進的途徑、改進的途徑 2. WNnk的對稱性和周期性的對稱性和周期性()2()(),;NnknkNNnkn N kn k NNNNWWWWW .),1( 1),1()2/(2/)(2)()(2kNNkNjNNNnkNnNNkNnkNknNNkNnNWWeWWeWWWWWN得得: :對稱性對

4、稱性: :周期性周期性: :1. WN0=1； WNN/2=e-j2 /NN/2 = -1 不必運算不必運算8q利用上述特性，可以將有些項合并，并將DFT分解為短序列，從而降低運算次數，提高運算速度。1965年，庫利(cooley)和圖基(Tukey)首先提出FFT算法。對于N點DFT，僅需(N/2)log2N次復數乘法運算。例如N=1024=210時，需要(1024/2)log2210 =512*10= 5120 次。5120/1048576=4.88%，速度提高20倍。96.3.3、庫利、庫利-圖基算法圖基算法一、算法原理一、算法原理(基基2FFT)(一) N/2點DFT1. 先將先將x(

5、n)按按n的奇偶分為兩組作的奇偶分為兩組作DFT，設設N=2L ，不足時，可，不足時，可補些零。有補些零。有: n為偶數時為偶數時： n為奇數時：為奇數時：1, 1 , 0 ),() 12(1, 1 , 0 ),()2(2221NNrrxrxrrxrx10( ) ( )( )NnkNnX kDFT x nx n W因此，因此，按時間抽取按時間抽取(DIT)的的FFT算法算法 庫利庫利-圖基算法圖基算法10庫利庫利-圖基算法圖基算法所以，上式可表示為：222)/(222NNNWeeWjjN由于：)()()()()(211021012222kXWkXWrxWWrxkXkNrrkkNrrkNNNN1

6、022102110)12(10210102222)()() 12()2()()()(NNNNrrkNkNrrkNrkrNrrkNNnNnnkNnkNWrxWWrxWrxWrxWnxWnxkX(n為偶數)(n為奇數)11庫利庫利-圖基算法圖基算法2. 兩點結論兩點結論: (1) X1(k)，X2(k)均為均為N/2點的點的DFT。 (2) X(k)=X1(k)+WNkX2(k)只能確定出只能確定出X(k)的的k=N/2-1個、即個、即前一半前一半的結果。的結果。2(0), (1), (1)Nxxx10102210101122222222) 12()()()2()()(NNNNNNNNrrkrrk

7、rrkrrkWrxWrxkXWrxWrxkX其中其中,12庫利庫利-圖基算法圖基算法q同理，同理，X2(N/2+k)=X2(k)，即即X1(k)，X2(k)的后一的后一半，分別等于其前一半的值。半，分別等于其前一半的值。)()()()2(1101)(101122222kXWrxWrxkNXNNNNNrrkkrr由于WN/2r(k+N/2)=WN/2rk (周期性）周期性）, ,所以:(3) X(k)的后一半的確定的后一半的確定13庫利庫利-圖基算法圖基算法q可見，見，X(k)的后一半，也完全由的后一半，也完全由X1(k)，X2 (k)的前一的前一半所確定。半所確定。q*N點的點的DFT可由兩個

8、可由兩個N/2點的點的DFT來計算。來計算。1, 1 , 0 ),()()2()2()2(221212NkNkNkkXWkXNkXWNkXNkXN又由于WN(N/2+k)=WNN/2WNk= -WNk，所以144. 蝶形運算蝶形運算實現上式運算的流圖稱作蝶形運算實現上式運算的流圖稱作蝶形運算前一半前一半后一半后一半（N/2個蝶形個蝶形)( (前一半前一半) )( (后一半后一半) )1 1 11-1-1)()()()()()(2121kXWkXkXkXWkXkXkNkN)1,()1, 1 , 0(22 N Nk kk kN NN N)()()(21kXWkXkXkN)()()2(21kXWkX

9、kNXkN)(1kX)(2kXkNW由X1(k)、X 2(k)表示X(k)的運算是一種特殊的運算-碟形運算碟形運算155.計算工作量分析計算工作量分析(1)N/2點的點的DFT運算量：復乘次數運算量：復乘次數:(N/2)2=N2/4 復加次數復加次數:N/2(N/2-1)(2)兩個兩個N/2點的點的DFT運算量：復乘次數運算量：復乘次數:N2/2 復加次數復加次數: N/2(N/2-1)(3)N/2個蝶形運算的運算量：復乘次數個蝶形運算的運算量：復乘次數:N/2 復加次數復加次數:2.N/2=N復乘:復加:2)12(2NNNN22222NNN總共運算量總共運算量: :按奇、偶分組后的計算量：按

10、奇、偶分組后的計算量：p但是，但是，N點點DFT的復乘為的復乘為N2；復加為；復加為N(N-1)；與分解后相比；與分解后相比可知，計算工作點差不多減少可知，計算工作點差不多減少 一半。一半。16q例如：例如：N=8 時的時的DFT,可以分解為兩個可以分解為兩個N/2 = 4點的點的DFT。具體方法如下。具體方法如下： (1) n為偶數時，即為偶數時，即x(0)，x(1)，x(2)，x(3)； 分別記作分別記作：33114400( )( )(2 ),0,1,2,3rkrkrrX kx r Wxr Wk);6()3(),4()2(),2()1(),0()0(1111xxxxxxxx進行進行N/2=

11、4點的點的DFT，得，得X1(k)17 (2) n為奇數時為奇數時，分別記作，分別記作:2222(0)(1),(1)(3),(2)(5),(3)(7);xxxxxxxx33224400( )( )(21),0,1,2,3rkrkrrXkx r WxrWk1212( )( )( )(4)( )( ),0,1,2,3kNkNX kX kW XkX kX kW Xkk進行進行N/2=4點的點的DFT，得，得X2(k)18x1(0)= x(0) x1(1)=x(2)x1(2)=x(4) x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) x2(2)=x(5) x2(3)=x(7) X X

12、1(0)(0)X X1(1)(1)X X1(2)(2)X X1(3)(3)X X2(0)(0)X X2(1)(1)X X2(2)(2)X X2(3)(3)-1-1-1-1X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)(3)對對X1(k)和和X2(k)進行蝶形運算，前半部為進行蝶形運算，前半部為X(0) X(3)，后半部分為后半部分為X(4) X(7)，整個過程如下圖所示整個過程如下圖所示:N/2點DFTN/2點DFT0NW1NW2NW3NW19(二二) N/4點點DFTq由于由于N=2L，所以所以N/2仍為偶數，

13、可以進一步把每個仍為偶數，可以進一步把每個N/2點的序列再按其奇偶部分分解為兩個點的序列再按其奇偶部分分解為兩個N/4的子序的子序列。例如，列。例如，n為偶數時的為偶數時的 N/2點分解為點分解為:134144(2 )( ), 0,1,1(21)( ), 0,1,1NNxlx lxlx l進行N/4N/4點的DFTDFT，得到klNllkNllkNllkNlWlxWlxkXWlxWlxkXNNNN) 12(2/1014/104422/1014/1033) 12()()()2()()(4444( (偶中偶) )(偶中奇)20)()()(4312kXWkXkXkN1, 1 , 04Nk從而可得到前

14、N/4點的X1(k)()()4(4312kXWkXkNXkN后N/4點的X1(k)為：1, 1 , 04Nk21(奇中偶奇中偶)104/64/1026104/5104/254444)() 12()()()2()(NNNNllkNlkNlllkNllkNWlxWlxkXWlxWlxkX(奇中奇奇中奇)q同樣對同樣對n為奇數時，為奇數時，N/2點分為兩個點分為兩個N/4點的序列進行點的序列進行DFT，則有：，則有：k25N/26k25N/26X (k)X (k)WX (k) ;k0,1,14NX (k)X (k)WX (k) ;k0,1,144NN進行碟形運算，得：、由)()(65kXkX22q例

15、如，例如，N=8時的時的DFT可分解為四個可分解為四個N/4的的DFT，具體，具體步驟如下：步驟如下：)4()2() 1 ()0()0()0()()()(131313xxxxxxnxrxlx(1) 將原序列x(n)的“偶中偶”部分：構成N/4點DFT，從而得到X3(0)， X3(1)。23)6()3()1()2()1()0()()()(141414xxxxxxnxrxlx(2) 將原序列x(n)的“偶中奇”部分:構成N/4點DFT，從而得到X4(0), X4(1)。(3) 將原序列x(n)的“奇中偶”部分:)5()2()1()1()0()0()()()(252525xxxxxxnxrxlx構成

16、N/4點DFT，從而得到X5(0), X5(1)。24(4) 將原序列x(n)的“奇中奇”部分:)7()3()1()3()1()0()()()(262626xxxxxxnxrxlx構成N/4點DFT，從而得到X6(0), X6(1)。（5）由 X3(0), X3(1), X4(0), X4(1)進行碟形運算, 得到X1(0), X1(1), X1(2), X1(3)。（6）由 X5(0), X5(1), X6(0), X6(1)進行碟形運算, 得到X2(0), X2(1), X2(2), X2(3)。25-1-1-1-2-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)（

17、7）由X1(0), 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) 。x3(0)=x1(0)=x(0)x3(1)=x1(2)=x(4)x4(0)=x1(1)=x(2)x4(1)=x1(3)=x(6)x5(0)=x2(0)=x(1) x5(1)=x2(2)=x(5) x6(0)=x2(1)=x(3)x6(1)=x2(3)=x(7)N/4DFTN/4DFTN/4DFTN/4DFT0NW1NW2NW3NW1(0)X1(0)X1(0)X1(0)

18、X2(0)X2(1)X2(2)X2(3)X0NW2NW0NW2NW3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)X26q這樣這樣, ,又一次分解又一次分解, ,得到四個得到四個N/4點點DFT，兩級兩級蝶蝶形運算，其運算量有大約減少一半即為形運算，其運算量有大約減少一半即為N點點DFT的的1/4。q對于對于N=8時時DFT，N/4點即為兩點點即為兩點DFT，因此因此0,1k ,)()(0,1k ,)()(0,1k ,)()(0,1k ,)()(4/1066104/55104/44104/33lkNlllkNllkNllkNWlxkXWlxkXWlxkXWlxkX

19、27亦即亦即, , )4()0() 1 ()0() 1 ()4()0() 1 ()0()0(031233030233xWxxWxXxWxxWxXNN)6()2() 1 ()0() 1 ()6()2() 1 ()0()0(041244040244xWxxWxXxWxxWxXNN)5() 1 () 1 ()0() 1 ()5() 1 () 1 ()0()0(051255050255xWxxWxXxWxxWxXNN)7() 3() 1 ()0() 1 ()7() 3() 1 ()0()0(061266060266xWxxWxXxWxxWxXNN28-1-1-1-1-1-1-1-1-1-1-1-1X(

20、0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)xxx因此因此,8 8點點DFTDFT的的FFTFFT的運算流圖如下的運算流圖如下0NW0NW0NW0NWx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)3NW2NW1NW0NW0NW2NW0NW2NWx3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)29q此此FFT算法，是在時間上對輸入序列的次序算法，是在時間上對輸入序列的次序是屬于偶數還是屬于奇數來進行分解的，是屬于偶數還是屬于奇數來進行分解的，所以稱作按

21、所以稱作按時間時間抽取的算法抽取的算法(DIT)。二、運算量二、運算量 由上述分析可知，由上述分析可知，N=8需三級蝶形運算需三級蝶形運算N=2 =8，由此可知，由此可知，N=2L 共需共需L級蝶形運算級蝶形運算,而且每級都而且每級都由由N/2個蝶形運算組成，每個蝶形運算有一次個蝶形運算組成，每個蝶形運算有一次復乘，兩次復加。復乘，兩次復加。3 330q因此，因此，N點的點的FFT的運算量為：的運算量為： 復乘復乘: mF=(N/2)L=(N/2)log2N 復加復加: aF=N L=Nlog2Nq由于計算機的乘法運算比加法運算所需由于計算機的乘法運算比加法運算所需的時間多得多，故以乘法運算作

22、為比較的時間多得多，故以乘法運算作為比較基準。基準。31三、三、DIT的的FFT算法的特點算法的特點-1-1-1-1-1-1-1-1.1.1.原位運算原位運算 輸入數據、中間運算結果和最后輸出均用同一存儲器。輸入數據、中間運算結果和最后輸出均用同一存儲器。x(0)=X0(0)x(4)=X0(1) x(2)=X0(2) x(6)=X0(3)x(1)=X0(4)x(5)=X0(5)x(3)=X0(6)x(7)=X0(7)0NW0NW0NW0NW0NW2NW0NW2NW0NW1NW2NW3NWX2(0)X2(1) X2(2) X2(3)X2(4)X2(5)X2(6)X2(7)X3(0)=X(0)X3

23、(1)=X(1) X3(2)=X(2)X3(3)=X(3)X3(4)=X(4)X3(5)=X(5)X3(6)=X(6)X3(7)=X(7)X1(0)X1(1) X1(2) X1(3)X1(4)X1(5)X1(6)X1(7)32rNmmmrNmmmWjXkXjXWjXkXkX)()()()()()(1111),3()6(),2()2(),1 ()4(),0()0(0000XxXxXxXx).7()7(),6() 3(),5()5(),4() 1 (0000XxXxXxXxq由運算流圖可知，一共有由運算流圖可知，一共有N個輸入個輸入/出行，一共有出行，一共有log2N=L列列(級級)蝶形運算蝶形運

24、算(基本迭代運算基本迭代運算)。 q設用設用m(m=1,2, ,L)表示第表示第m列；用列；用k,j表示蝶形輸入數據所在表示蝶形輸入數據所在的的(上上/下下)行數行數(0,1,2, ,N-1)；這時任何一個蝶形運算可用下；這時任何一個蝶形運算可用下面通用式表示，即：面通用式表示，即：33所以，當所以，當m=1時，則有（前兩個蝶形）則有（前兩個蝶形）00010001)1()0()1()1()0()0(NNWXXXWXXX00010001)3()2()3()3()2()2(NNWXXXWXXX34當當m=2時時，則有，則有(前兩個蝶形前兩個蝶形)2112211201120112)3()1()3()

25、3()1()1()2()0()2()2()0()0(NNNNWXXXWXXXWXXXWXXX1223122302230223)5()1()5()5()1()1()4()0()4()4()0()0(NNNNWXXXWXXXWXXXWXXX當當m=3時時，則有，則有(前兩個蝶形前兩個蝶形) 35q可見，在某列進行蝶形運算的任意兩個節點可見，在某列進行蝶形運算的任意兩個節點(行行)k和和j的節點變量的節點變量Xm-1(k)，Xm-1(j)就完全可以確定蝶形就完全可以確定蝶形運算的結果運算的結果Xm(k)，Xm(j)，與其它行，與其它行(節點節點)無關。無關。 這樣，蝶形運算的兩個輸出值仍可放回蝶形運

26、算這樣，蝶形運算的兩個輸出值仍可放回蝶形運算的兩個輸入所在的存儲器中，即實現所謂原位運的兩個輸入所在的存儲器中，即實現所謂原位運算。每一組算。每一組(列列)有有N/2個蝶形運算，所以只需個蝶形運算，所以只需N個存個存儲單元，可以節儲單元，可以節 省存儲單元。省存儲單元。362. 倒位序規律倒位序規律 由圖可知，輸出由圖可知，輸出X(k)按正常順序排列在存儲單元，按正常順序排列在存儲單元，而輸入是按順序：而輸入是按順序：;7 (),3 (),5 (),1 (;6 (),2 (),4 (),0 (）xxxxxxxx這種順序稱作倒位序，即二進制數倒位。這種順序稱作倒位序，即二進制數倒位。37n0=0

27、，(偶)n1=0n1 =1n1 =0n1 =101010101 ),(012nnnx(n2)x(000) 0 乾x(100) 4 兌x(010) 2 離x(110) 6 震x(001) 1 巽x(101) 5 坎x(011) 3 艮x(111) 7 坤這是由奇偶分組造成的，以N=8為例說明如下：n0=1，(偶)383.3.倒位序實現倒位序實現 輸入序列先按自然順序存入存儲單元，然后經輸入序列先按自然順序存入存儲單元，然后經變址運算來實現變址運算來實現倒位序排列倒位序排列 設輸入設輸入序列的序號為序列的序號為n，二進制為，二進制為(n2 n1 n0 )2，倒位序倒位序順序用順序用 表示，其表示，

28、其倒位序倒位序二進制為二進制為(n0 n1 n2 )2 ，例如，例如，N=8時如下表：時如下表： n390 0 0 0 0 0 0 0 1 0 0 1 1 0 0 4 2 0 1 0 0 1 0 2 3 0 1 1 1 1 0 6 4 1 0 0 0 0 1 1 5 1 0 1 1 0 1 5 6 1 1 0 0 1 1 3 7 1 1 1 1 1 1 7 自然順序自然順序n 二進制二進制n2 n1 n0 倒位序二進制倒位序二進制n0 n1 n2 倒位順序倒位順序 n40A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8)x(0) x(1) x(2) x(3) x(4

29、) x(5) x(6) x(7)x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)變址處理方法變址處理方法存儲單元存儲單元自然順序自然順序變址變址倒位序倒位序4. 蝶形運算兩節點的距離：蝶形運算兩節點的距離：2m-1其中，其中，m表示第表示第m列，且列，且m =1, ,L例如例如N=8=23，第一級，第一級(列列)距離為距離為21-1=1， 第二級第二級(列列)距離為距離為22-1=2， 第三級第三級(列列)距離為距離為23-1=4。415. WNr 的確定的確定(僅給出方法僅給出方法) 考慮蝶形運算兩節點的距離為考慮蝶形運算兩節點的距離為2m-1，蝶形運算可表為：

30、蝶形運算可表為： Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WNr Xm(k+2m-1)=Xm-1(k)-Xm-1(k+2m-1)WNr 由于由于N為已知，所以將為已知，所以將r的值確定即可。的值確定即可。 為此，令為此，令k=(n2 n1 n0)2，再將，再將k= (n2 n1 n0)2 左移左移(L-m)位，位，右邊位置補零，就可得到右邊位置補零，就可得到(r)2的值，即的值，即(r)2 =(k)22L-m 。 42 例如：例如：N=8=23 (1) k=2，m=3的的r值，因值，因k=2=(010)2 左移左移L-m=3-3=0 ，故故 r=(010)2=2； (2) k=3，

31、m=3的的r值；因值；因k= 3=(011)2左移左移0位，位，r =3； (3) k=5，m=2的的r值；因值；因k=5=(101)2 左移左移L-m=1位，位，r=(010)2 =2。436.存儲單元存儲單元 存輸入序列存輸入序列 (n)，n=0, 1, N-1，計計N個單元個單元； 存放系數存放系數WNr，r=0, 1, , (N/2)-1，需需N/2個存個存儲單元；儲單元； 共計共計(N+N/2)個存儲單元個存儲單元。x446.3.4 IFFT算法算法一一. 稍微變動稍微變動FFT程序和參數可實現程序和參數可實現IFFTnkNNkNnnkNWkXNkXIDFTnxWnxnxDFTkX1

32、010)(1)()()()()(q比較兩式可知比較兩式可知,只要只要DFT的每個系數的每個系數WNnk 換成換成WN-nk，最后再，最后再乘以常數乘以常數1/N就可以得到就可以得到IDFT的快速算法的快速算法-IFFT。q另外，可以將常數另外，可以將常數1/N分配到每級運算中，分配到每級運算中，1/N =1/2L=(1/2)L，也就是每級也就是每級蝶形運算均乘以蝶形運算均乘以1/2。 45利用利用FFT程序實現程序實現IFFT二二. 不改不改(FFT)的程序直接實現的程序直接實現IFFT 1*0101( )( )1( )NnkNkNnkNkxnX k WNXk WN101( )( )1( )NnkNkx nXk WNDFTXkN,nknkNNWWA BAB因為因為所以所以因此因此q步驟為：先將步驟為：先將X(