華科-機(jī)械工程測(cè)試信息信號(hào)分析-課件-ch6-03-數(shù)字信號(hào)分析

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

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

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

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

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

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

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

8、N-1)；與分解后相比可知，計(jì)算工作點(diǎn)差不多減少 一半。,例如：N=8 時(shí)的DFT,可以分解為兩個(gè)N/2 = 4點(diǎn)的DFT。具體方法如下：(1) n為偶數(shù)時(shí)，即x(0)，x(1)，x(2)，x(3)； 分別記作：,進(jìn)行N/2=4點(diǎn)的DFT，得X1(k),(2) n為奇數(shù)時(shí)，分別記作:,進(jìn)行N/2=4點(diǎn)的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)對(duì)X1(k)和X2(k)進(jìn)行蝶形運(yùn)算，前半部為X(0) ? X(3)，后半部分為X(4) ?X(7)，整個(gè)過(guò)程如下圖所示:,N/2點(diǎn)DFT,N/2點(diǎn)DFT,(二) N/4點(diǎn)DFT由于N=2L，所以N/2仍為偶數(shù)，可以進(jìn)一步把每個(gè)N/2點(diǎn)的序列再按其奇偶部分分解為兩個(gè)N/4的子序列。例如，n為偶數(shù)時(shí)

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

12、:,構(gòu)成N/4點(diǎn)DFT，從而得到X4(0), X4(1)。,(3) 將原序列x(n)的“奇中偶”部分:,構(gòu)成N/4點(diǎn)DFT，從而得到X5(0), X5(1)。,(4) 將原序列x(n)的“奇中奇”部分:,構(gòu)成N/4點(diǎn)DFT，從而得到X6(0), X6(1)。,（5）由 X3(0), X3(1), X4(0), X4(1)進(jìn)行碟形運(yùn)算, 得到X1(0), X1(1), X1(2), X1(3)。,（6）由 X5(0), X

13、5(1), X6(0), X6(1)進(jìn)行碟形運(yùn)算, 得到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)進(jìn)行碟形運(yùn)算，得到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、次分解,得到四個(gè)N/4點(diǎn)DFT，兩級(jí)蝶形運(yùn)算，其運(yùn)算量有大約減少一半即為N點(diǎn)DFT的1/4。對(duì)于N=8時(shí)DFT，N/4點(diǎn)即為兩點(diǎn)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點(diǎn)DFT的FFT的運(yùn)算流圖如下,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í)間上對(duì)輸入序列的次序是屬于偶數(shù)還是屬于奇數(shù)來(lái)進(jìn)行分解的，所以稱(chēng)作按時(shí)間抽取的算法(DIT)。二、運(yùn)算量

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

18、,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,-1,-1,-1,-1,-1,-1,-1,-1,...,..,1.原位運(yùn)算 輸入數(shù)據(jù)、中間運(yùn)算結(jié)果和最后輸出均用同一存儲(chǔ)器。,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),由運(yùn)算流圖可知，一共有N個(gè)輸入/出行，一共有l(wèi)og2N=L列(級(jí))蝶形運(yùn)算(基本迭代運(yùn)算

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

21、。 這樣，蝶形運(yùn)算的兩個(gè)輸出值仍可放回蝶形運(yùn)算的兩個(gè)輸入所在的存儲(chǔ)器中，即實(shí)現(xiàn)所謂原位運(yùn)算。每一組(列)有N/2個(gè)蝶形運(yùn)算，所以只需N個(gè)存儲(chǔ)單元，可以節(jié) 省存儲(chǔ)單元。,2. 倒位序規(guī)律由圖可知，輸出X(k)按正常順序排列在存儲(chǔ)單元，而輸入是按順序：,這種順序稱(chēng)作倒位序，即二進(jìn)制數(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為例說(shuō)明如下：,n0=1，(偶),3.倒位序?qū)崿F(xiàn)輸入序列先按自然順序存入存儲(chǔ)單元，然后經(jīng)變址運(yùn)算來(lái)實(shí)現(xiàn)倒位序排列設(shè)輸入序列的序號(hào)為n，二進(jìn)制為(n2 n1 n0 )2，倒位序順序用 表示，其倒位序

23、二進(jìn)制為(n0 n1 n2 )2 ，例如，N=8時(shí)如下表：,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 二進(jìn)制n2 n1 n0 倒位序二進(jìn)制n0 n1 n2 倒位順序,變址處理方法,存儲(chǔ)單元,自然順序,變址,倒位序,4

25、. 蝶形運(yùn)算兩節(jié)點(diǎn)的距離：2m-1其中，m表示第m列，且m =1,… ,L例如N=8=23，第一級(jí)(列)距離為21-1=1， 第二級(jí)(列)距離為22-1=2， 第三級(jí)(列)距離為23-1=4。,5. WNr 的確定(僅給出方法)考慮蝶形運(yùn)算兩節(jié)點(diǎn)的距離為2m-1，蝶形運(yùn)算可表為： 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)位，右邊位置補(bǔ)零，就可得到(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.存儲(chǔ)單元存輸入序列 (n)，n=0, 1,…, N-1，計(jì)N個(gè)單元；存放系數(shù)WNr，r=0, 1, … , (N/2)-1，需N/2個(gè)存儲(chǔ)單元；共計(jì)(N+N/2)個(gè)存儲(chǔ)單元。,6.3.4 IFFT算法,一. 稍微變動(dòng)

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

29、；然后再取一次共軛；最后再乘1/N，即得 x(n)。所以FFT，IFFT可用一個(gè)子程序。,6.3.5 線性卷積的FFT算法,一、線性卷積的長(zhǎng)度設(shè)一離散線性移不變系統(tǒng)的沖激響應(yīng)為h(n)，其輸入信號(hào)為x(n)，其輸出為y(n)，并且x(n)的長(zhǎng)度為L(zhǎng)點(diǎn)，h(n)的長(zhǎng)度為M點(diǎn)，則：,以實(shí)例說(shuō)明：,,,,,,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,,,。,可見(jiàn)，y(n)的長(zhǎng)度為L(zhǎng)+M-1,二、用FFT算y(n)的步驟,1、將x(n)，y(n)補(bǔ)零點(diǎn)，至少為N=M+L-1點(diǎn)；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),流程圖,三、幾點(diǎn)說(shuō)明,1、當(dāng) M=L 時(shí)，用圓周卷積計(jì)算線性卷積的速度快，點(diǎn)數(shù)越多，速度越快，N≥64時(shí)，速度增加明顯。2、L>>M 時(shí)，速度增加不太明顯，為了增加速度，采用 (1)重疊相加法 (2)重疊保留法(從略