華科-機械工程測試信息信號分析-課件-ch6-03-數(shù)字信號分析

1、MEASUREMENTINFORMATION SIGNAL ANALYSIS IN MECHANICAL ENGINEERING,機械工程測試?信息?信號分析,機械科學(xué)與工程學(xué)院 機械電子信息工程系,課件資料下載：郵箱地址：密碼：注意下載時不要刪除原始文件,第六章 數(shù)字信號分析(I)DFT與FFT,§6-5 現(xiàn)代譜分析方法-最大熵譜估計,§6-3 FFT,§6-4 譜分析與

2、譜估計,§6-2 DFT,§6-1 模擬信號離散化,目錄,6-3 FFT算法,6.3.1、DFT的計算工作量,,FFT背景,兩者的差別僅在指數(shù)的符號和因子1/N.,通常x(n)和WNnk都是復(fù)數(shù)，所以計算一個x(k)的值需要N次復(fù)數(shù)乘法運算，和N-1次復(fù)數(shù)加法運算。那么，所有的X(k)就要N2次復(fù)數(shù)乘法運算，N(N-1)次復(fù)數(shù)加法運算。當(dāng)N很大時，運算量將是驚人的,如N=1024，則要完成1048576次(一百

3、多萬次)運算。難以做到實時處理。,計算一個X(k)的值的計算量，如X(1)，k=1,6.3.2、改進的途徑,2. WNnk的對稱性和周期性,得:,對稱性:,周期性:,1. WN0=1； WNN/2=[e-j2?/N]N/2 = -1 不必運算,利用上述特性，可以將有些項合并，并將DFT分解為短序列，從而降低運算次數(shù)，提高運算速度。1965年，庫利(cooley)和圖基(Tukey)首先提出FFT算法。對于N點DFT，僅需(N/2)lo

4、g2N次復(fù)數(shù)乘法運算。例如N=1024=210時，需要(1024/2)log2210 =512*10= 5120 次。5120/1048576=4.88%，速度提高20倍。,6.3.3、庫利-圖基算法,一、算法原理(基2FFT)(一) N/2點DFT1. 先將x(n)按n的奇偶分為兩組作DFT，設(shè)N=2L ，不足時，可補些零。有: n為偶數(shù)時： n為奇數(shù)時：,因此，,按時間抽取(DIT)的FFT算法 —庫利-圖基算法,庫利

5、-圖基算法,所以，上式可表示為：,庫利-圖基算法,2. 兩點結(jié)論: (1) X1(k)，X2(k)均為N/2點的DFT。 (2) X(k)=X1(k)+WNkX2(k)只能確定出X(k)的k=N/2-1個、即前一半的結(jié)果。,庫利-圖基算法,同理，X2(N/2+k)=X2(k)，即X1(k)，X2(k)的后一半，分別等于其前一半的值。,由于WN/2r(k+N/2)=WN/2rk (周期性）,所以:,(3) X(k)的后一半的確定,庫利

6、-圖基算法,可見，X(k)的后一半，也完全由X1(k)，X2 (k)的前一半所確定。*N點的DFT可由兩個N/2點的DFT來計算。,又由于WN(N/2+k)=WNN/2WNk= -WNk，所以,4. 蝶形運算,實現(xiàn)上式運算的流圖稱作蝶形運算,,前一半,后一半,（N/2個蝶形),,,,,,,,,,,(前一半),(后一半),1 1,1,1,-1,由X1(k)、X 2(k)表示X(k)的運算是一種特殊的運算-碟形運算

7、,5.計算工作量分析,(1)N/2點的DFT運算量：復(fù)乘次數(shù):(N/2)2=N2/4 復(fù)加次數(shù):N/2(N/2-1)(2)兩個N/2點的DFT運算量：復(fù)乘次數(shù):N2/2 復(fù)加次數(shù): N/2(N/2-1)(3)N/2個蝶形運算的運算量：復(fù)乘次數(shù):N/2 復(fù)加次數(shù):2.N/2=N,,復(fù)乘:,復(fù)加:,總共運算量:,按奇、偶分組后的計算量：,但是，N點DFT的復(fù)乘為N2；復(fù)加為N(

8、N-1)；與分解后相比可知，計算工作點差不多減少 一半。,例如：N=8 時的DFT,可以分解為兩個N/2 = 4點的DFT。具體方法如下：(1) n為偶數(shù)時，即x(0)，x(1)，x(2)，x(3)； 分別記作：,進行N/2=4點的DFT，得X1(k),(2) n為奇數(shù)時，分別記作:,進行N/2=4點的DFT，得X2(k),x1(0)= x(0)

9、 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),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,X1(0),X1(1),X1(2),X1(3),X2(0),X2(1),

10、X2(2),X2(3),-1,-1,-1,-1,X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),(3)對X1(k)和X2(k)進行蝶形運算，前半部為X(0) ? X(3)，后半部分為X(4) ?X(7)，整個過程如下圖所示:,N/2點DFT,N/2點DFT,(二) N/4點DFT由于N=2L，所以N/2仍為偶數(shù)，可以進一步把每個N/2點的序列再按其奇偶部分分解為兩個N/4的子序列。例如，n為偶數(shù)時

11、的 N/2點分解為:,進行N/4點的DFT，得到,,(偶中偶),(偶中奇),從而可得到前N/4點的X1(k),后N/4點的X1(k)為：,,(奇中偶),(奇中奇),同樣對n為奇數(shù)時，N/2點分為兩個N/4點的序列進行DFT，則有：,,例如，N=8時的DFT可分解為四個N/4的DFT，具體步驟如下：,(1) 將原序列x(n)的“偶中偶”部分：,構(gòu)成N/4點DFT，從而得到X3(0)， X3(1)。,(2) 將原序列x(n)的“偶中奇”部分

12、:,構(gòu)成N/4點DFT，從而得到X4(0), X4(1)。,(3) 將原序列x(n)的“奇中偶”部分:,構(gòu)成N/4點DFT，從而得到X5(0), X5(1)。,(4) 將原序列x(n)的“奇中奇”部分:,構(gòu)成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), X

13、5(1), X6(0), X6(1)進行碟形運算, 得到X2(0), X2(1), X2(2), X2(3)。,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,-1,-1,-1,-2,-1,-1,X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),（7）由X1(0), X1(1), X1(2)，X1(3)，X2(0)， X2(1)，X2(2)，X2(3)進行碟形運算，得到X(0

14、), 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/4DFT,N/4DFT,N/4DFT,N/4DFT,這樣,又一

15、次分解,得到四個N/4點DFT，兩級蝶形運算，其運算量有大約減少一半即為N點DFT的1/4。對于N=8時DFT，N/4點即為兩點DFT，因此,亦即,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,-1,-1,-1,-1,,-1,-1,-1,-1,-1,-1,-1,-1,X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7)

16、,因此,8點DFT的FFT的運算流圖如下,x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7),x3(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),此FFT算法，是在時間上對輸入序列的次序是屬于偶數(shù)還是屬于奇數(shù)來進行分解的，所以稱作按時間抽取的算法(DIT)。二、運算量

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

18、,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,-1,-1,-1,-1,-1,-1,-1,-1,...,..,1.原位運算 輸入數(shù)據(jù)、中間運算結(jié)果和最后輸出均用同一存儲器。,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),X2(0)X2(1)

19、X2(2) X2(3)X2(4)X2(5)X2(6)X2(7),X3(0)=X(0)X3(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),由運算流圖可知，一共有N個輸入/出行，一共有l(wèi)og2N=L列(級)蝶形運算(基本迭代運算

20、)。 設(shè)用m(m=1,2, … ,L)表示第m列；用k,j表示蝶形輸入數(shù)據(jù)所在的(上/下)行數(shù)(0,1,2,… ,N-1)；這時任何一個蝶形運算可用下面通用式表示，即：,所以，當(dāng)m=1時，則有（前兩個蝶形）,,,當(dāng)m=2時，則有(前兩個蝶形),當(dāng)m=3時，則有(前兩個蝶形),可見，在某列進行蝶形運算的任意兩個節(jié)點(行)k和j的節(jié)點變量Xm-1(k)，Xm-1(j)就完全可以確定蝶形運算的結(jié)果Xm(k)，Xm(j)，與其它行(節(jié)點)無關(guān)

21、。 這樣，蝶形運算的兩個輸出值仍可放回蝶形運算的兩個輸入所在的存儲器中，即實現(xiàn)所謂原位運算。每一組(列)有N/2個蝶形運算，所以只需N個存儲單元，可以節(jié) 省存儲單元。,2. 倒位序規(guī)律由圖可知，輸出X(k)按正常順序排列在存儲單元，而輸入是按順序：,這種順序稱作倒位序，即二進制數(shù)倒位。,,,,,,,,,,,,,,,,,,,,,,,,n0=0，(偶),n1=0,n1 =1,n1 =0,n1 =1,0,1,0,1,0,1,0,1,,,(

22、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，(偶),3.倒位序?qū)崿F(xiàn)輸入序列先按自然順序存入存儲單元，然后經(jīng)變址運算來實現(xiàn)倒位序排列設(shè)輸入序列的序號為n，二進制為(n2 n1 n0 )2，倒位序順序用 表示，其倒位序

23、二進制為(n0 n1 n2 )2 ，例如，N=8時如下表：,0 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

24、 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 倒位順序,變址處理方法,存儲單元,自然順序,變址,倒位序,4

25、. 蝶形運算兩節(jié)點的距離：2m-1其中，m表示第m列，且m =1,… ,L例如N=8=23，第一級(列)距離為21-1=1， 第二級(列)距離為22-1=2， 第三級(列)距離為23-1=4。,5. WNr 的確定(僅給出方法)考慮蝶形運算兩節(jié)點的距離為2m-1，蝶形運算可表為： Xm(k)=X

26、m-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 。,,例如：N=8=23 (1) k=2，m=3的r值，因k=2=(010)2 左移L-m=3-3=

27、0 ，故 r=(010)2=2；(2) k=3，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。,6.存儲單元存輸入序列 (n)，n=0, 1,…, N-1，計N個單元；存放系數(shù)WNr，r=0, 1, … , (N/2)-1，需N/2個存儲單元；共計(N+N/2)個存儲單元。,6.3.4 IFFT算法,一. 稍微變動

28、FFT程序和參數(shù)可實現(xiàn)IFFT,比較兩式可知,只要DFT的每個系數(shù)WNnk 換成WN-nk，最后再乘以常數(shù)1/N就可以得到IDFT的快速算法--IFFT。另外，可以將常數(shù)1/N分配到每級運算中，1/N =1/2L=(1/2)L，也就是每級蝶形運算均乘以1/2。,利用FFT程序?qū)崿F(xiàn)IFFT,二. 不改(FFT)的程序直接實現(xiàn)IFFT,,因為,所以,因此,步驟為：先將X(k)取共軛，即將X(k)的虛部乘-1，直接利用FFT程序計算DFT

29、；然后再取一次共軛；最后再乘1/N，即得 x(n)。所以FFT，IFFT可用一個子程序。,6.3.5 線性卷積的FFT算法,一、線性卷積的長度設(shè)一離散線性移不變系統(tǒng)的沖激響應(yīng)為h(n)，其輸入信號為x(n)，其輸出為y(n)，并且x(n)的長度為L點，h(n)的長度為M點，則：,以實例說明：,,,,,,0 1,1,2,3,2,3,,,,,,,,,,,,,,,,,,,,,,,,0 1 2 3,。,,,,,,0 1,1,

30、2,3,2,3,。,,,,,,,0 1 2 3 4,,,,,,0 1,1,2,3,2,3,,,,,,,,0 1 2 3 4 5,。,。,,,,,,,,0 1 2 3 4 5 6,1,3,3,5,6,,,。,可見，y(n)的長度為L+M-1,二、用FFT算y(n)的步驟,1、將x(n)，y(n)補零點，至少為N=M+L-1點；2、求H(k)=FFT[h(n)]；3、求X(k

31、)=FFT[x(n)]；4、求Y(k)=X(k)H(k)；5、求y(n)=IFFT[Y(k)]；,FFT,FFT,IFFT,,,,x,,,,x(n),h(n),,,,y(n),X(k),H(k),Y(k),流程圖,三、幾點說明,1、當(dāng) M=L 時，用圓周卷積計算線性卷積的速度快，點數(shù)越多，速度越快，N≥64時，速度增加明顯。2、L>>M 時，速度增加不太明顯，為了增加速度，采用 (1)重疊相加法 (2)重疊保留法(從略