1、2.2.7.3 相關分析一、相關 二、互相關函數與自相關函數 三、相關函數的工程意義及應用 一、相關(correlation)相關：用來描述一個隨機過程自身在不同時刻的狀態間，或者兩個隨機過程在某個時刻狀態間線性依從關系的數字特征。 圖2.52 變量x和y的相關性 （a）精確相關 （b）中等程度相關 （c）不相關 評價變量x和y間線性相關程度的經典方法：n協方差xy： 式中，E表示數學期望值； x=Ex為隨機變量x的均值； y=Ey為隨機變量y的均值； n相關函數xy：式中x、y分別為x、y的標準偏差，而x和y的方差x2和y2則分別為 NiyixiNyxxyyxNyxE11lim（2.142）

2、11xyyxxyxy（2.143）22xxxE22yyxE（2.144）（2.145）利用柯西許瓦茲不等式(Cauchy-Schwarz inequality) 可知|xy|1。當xy=1時，所有數據點均落在y-y=m(x- x)的直線上，因此x,y兩變量是理想的線性相關。當xy=0時，(xi-x)與(yi-y)的正積之和等于其負積之和，因而其平均積xy為0，表示x,y之間完全不相關。 222yxyxyExEyxE（2.146）二、互相關函數與自相關函數 對于各態歷經過程，可定義時間變量x(t)和y(t)的互協方差(cross-covariance)函數為 式中稱x(t)與y(t)的互相關(c

3、ross-correlation)函數，自變量稱為時移。 yxxyTxxTyxxyRdttytxTtytxEC01lim（2.147） TTxydttytxTR01lim（2.148）當y(t) x(t)時，得自協方差(auto-covariance)函數 其中 稱為x(t)的自相關(auto-correlation)函數。周期函數的自相關函數仍為周期函數，且兩者的頻率相同，但丟掉了相角信息。 同頻相關，不同頻不相關。 201limxxTxxTxRdttxtxTC（2.149） TTxdttxtxTR01lim（2.150）圖2.53 典型的自相關函數和互相關函數曲線(a)自相關函數 （b）互

4、相關函數 例1 求正弦函數x(t)=Asin(t+)的自相關函數。解：正弦函數x(t)是一個均值為零的各態歷經隨機過程，其各種平均值可用一個周期內的平均值來表示。令t+=,則dt=d/,由此得正弦函數的自相關函數是一個與原函數具有相同頻率的余弦函數，它保留了原信號的幅值和頻率信息，但失去了原信號的相位信息。 自相關函數可用來檢測淹沒在隨機信號中的周期分量。 0000)(sin)sin(1)()(1lim)(TTTxdtttATdttxtxTRcos2)sin(sin2)(2202AdARx自相關和互相關函數的估計 和 具有限個數據點N的相關函數估計的數字處理表達式則為： )(xR)(xyR T

5、xdttxtxTR01 TxydttytxTR01（2.160）（2.161） 11NonxrnxnxNrR 11NonxyrnynxNrRNrr, 2 , 1 , 0（2.162）（2.163）三、相關函數的工程意義及應用 不同類別信號的辨識 圖2.55 典型信號的自相關函數 相關濾波(filtering by correlation) 圖4.79 相關濾波頻譜分析儀原理框圖 相關測速和測距 圖2.56 相關法測量聲傳播距離 圖2.57 帶鋼測速系統 測量流速和流量 圖2.58 相在法測定流量 2.2.7.4 功率譜分析 2.2.7.4 功率譜(power spectrum)分析 一、自功率

6、譜密度函數 二、巴塞伐爾（Parseval）定理 三、互功率譜密度函數 四、自譜和互譜的估計五、工程應用 一、自功率譜密度函數 設x(t)為一零均值的隨機過程，且x(t)中無周期性分量，則其自相關函數Rx()在當時有 該自相關函數Rx()滿足傅里葉變換的條件 。對作傅里葉變換可得 其逆變換為 0)(xRdRx)( deRfSfjxx2（2.167） dfefSRfjxx2（2.168）Sx(f)為x(t)的自功率譜密度函數(auto power spectrum)，簡稱自譜或功率譜。功率譜Sx(f)與自相關函數Rx()之間是傅里葉變換對的關系，亦即式（2.167）和（2.168）稱為維納辛欽（

7、Wiener-Khintchine）公式。由于Rx()為實偶函數，因此亦為Sx(f)實偶函數。 )()(fSRxFTIFTx圖2.59 單邊功率譜和雙邊功率譜 當=0時，根據自相關函數Rx()和自功率譜密度函數Sx(f)的定義，可得Sx(f)曲線下面和頻率軸所包圍的面積即為信號的平均功率；Sx(f)就是信號的功率譜密度沿頻率軸的分布，故也稱為功率譜。 dffSdttxTRxTTTx2221lim0（2.169）二、巴塞伐爾（Parseval）定理設有變換對： 按頻域卷積定理有 令k=0，有 又令h(t)=x(t)，得 )()(fXtx)()(fHth)(*)()()(fHfXthtxdffkH

8、fXdtethtxktj)()()()(2dffHfXdtthtx)()()()(dffXfXdttx)()()(2x(t)為實函數，故X(-f)=X*(f)，于是有v巴塞伐爾定理：信號在時域中計算的總巴塞伐爾定理：信號在時域中計算的總能量等于它在頻域中計算的總能量。能量等于它在頻域中計算的總能量。式（2.170）又稱信號能量等式。|X(f)|2稱能量譜，它是沿頻率軸的能量分布密度。在整個時間軸上信號的平均功率可計算為自譜密度函數與幅值譜之間的關系為 dffXdffXfXdttx2*2（2.170） dffxTdttxTPTTTT22221lim1lim（2.171） 21limfxTfSTx

9、（2.172） 對于單邊功率譜G(f) 也應滿足巴塞伐爾定理，故有由此規定 Gx(f)的圖形如圖2.59中所示。 dffGdffSPxx（2.173）0)(2)(ffSfGxx圖2.59 單邊功率譜和雙邊功率譜 根據信號功率（或能量）在頻域中的分布情況，將隨機過程區分為窄帶隨機、寬帶隨機和白噪聲等幾種類型。窄帶過程的功率譜（或能量）集中于某一中心頻率附近，寬帶過程的能量則分布在較寬的頻率上，而白噪聲過程的能量在所分析的頻域內呈均勻分布狀態。 三、互功率譜密度函數 若互相關函數Rxy()滿足傅里葉變換的條件 ，則定義Rxy()的傅里葉變換 為信號x(t)和y(t)的互功率譜密度函數，簡稱互譜密度

10、函數(cross power spectrum)或互譜。 根據維納辛欽關系，互譜與互相關函數也是一個傅里葉變換對，即 因此Sxy(f)的傅里葉逆變換為： dRxy)( dfeRSftjxyxy2（2.175）)()(fSRxyFTIFTxy dfefSRfjxyxy2（2.176） 定義信號x(t)和y(t)的互功率為因此互譜和幅值譜的關系為 正如Ryx()Rxy()一樣，當x和y的順序調換時，Syx()Sxy() 。但根據Rxy(-)=Ryx() 及維納辛欽關系式，不難證明：其中 dffXfYTdttytxTPTTTT*221lim1lim（2.177） fXfYTfSTxy*1lim（2.

11、178） fSfSfSyxxyxy*（2.179） )()(1lim*fYfXTfStxy Sxy(f)也是含正、負頻率的雙邊互譜，實用中也常取只含非負頻率的單邊互譜Gxy(f)，由此規定 自譜是f的實函數，而互譜則為f的復函數，實部Cxy(f)稱為共譜(cospectrum)，虛部Qxy(f)稱為重譜(quad spectrum)，即寫為幅頻和相頻的形式： 02ffSfGxyxy（2.180） fjQfCfGxyxyxy（2.181） fCfQarctgffQfCfGefGfGxyxyxyxyxyxyfjxyxyxy22（2.182）四、自譜和互譜的估計定義功率譜亦即自譜的估計值 互譜的估計

12、為 21fXTfSx（2.183） fYfXTfSxy*1（2.184） fXfYTfSyx*1（2.185）五、工程應用 1.求取系統的頻響(frequency response)函數 線性系統的傳遞函數H(s)或頻響函數H(j)十分重要，在機器故障診斷等多個領域常要用到它。w例1：機器由于其軸承的缺陷而在機器運行中會造成沖擊脈沖信號，此時若用安裝在機殼外部的加速度傳感器來接收時，必須考慮機殼的傳遞函數。w例2：當信號經過一個復雜系統被傳輸時，系統各環節的傳遞函數便必須要加以考慮。 一個線性系統的輸出y(t)等于其輸入x(t)和系統的脈沖響應h(t)的卷積，即 根據卷積定理，上式在頻域中化為

13、式中H(f)即為系統的頻響函數。 thtxty*（2.192） fXfHfY（2.193）通過自譜和互譜來求取H(f)： 對式（2.193）兩端乘以各自的復共軛并取期望值有 上式反映出輸入與輸出的功率譜密度和頻響函數間的關系；式中沒有頻響函數的相位信息，因此不可能得到系統的相頻特性。 fSfHfSxy2（2.194） 如果在式（2.193）兩端乘以x(f)的復共軛并取期望值，則有 由于Sx(f)為實偶函數，因此頻響函數的相位變化完全取決于互譜密度函數的相位變化。式（2.195）將輸入、輸出的相位關系完全保留了下來，且在這里輸入的形式并不一定限制為確定性信號，也可以是隨機信號。 fXfXfHfX

14、fY* fSfHfSxxy（2.195）通常一個測試系統往往受到內部和外部噪聲的干擾。從而輸出也會帶入干擾。輸入信號與噪聲是獨立無關的，因此它們的互相關為零。結論：在用結論：在用和和求取系統頻響函求取系統頻響函數時不會受到系統干擾的影響。數時不會受到系統干擾的影響。 2.旋轉機械振動特性檢測 n旋轉機械的轉軸部件從起動、升速到額定轉速的過程共經歷了全部轉速的變化，因此在各個轉速下的振動狀態可用來對機器的臨界轉速、固有頻率和阻尼比等各參數進行辨識。n起動和停車過程則包含了豐富的信息。是常規運行狀態下所無法獲得的。n“瀑布圖(waterfall plot)法”：在機械振動或停車過程中將不同轉速下振

15、動的功率譜圖迭加而形成的一種圖。圖2.61 旋轉機械的瀑布圖 由圖可見機器的回轉頻率n(r/min)及其各次諧波下譜峰高度，由此來得出機器的臨界轉速、固有頻率及阻尼比等數據。從圖可見，機器臨界轉速約為4000r/min，機器振動的高次諧波分量很小，主要是回轉頻率處的譜峰，因此可判斷轉子存在有較嚴重的失衡。此外還可看到圖中頻率60HZ處有一譜峰值，它不隨轉速升高而改變，判斷為電源的脈動干擾。 2.3 數字信號處理數字信號處理(digital signal processing)：利用計算機或專用信號處理設備，以數值計算的方法對信號作采集、變換、綜合、估值與識別等處理。一、離散傅里葉變換（DFT）

16、二、離散傅里葉變換的性質 三、采樣定理 四、泄漏與加窗處理 五、柵欄效應六、快速傅里葉變換（FFT） 一、離散傅里葉變換（DFT）對于一個非周期的連續時間信號x(t)來說，它的傅里葉變換應該是一個連續的頻譜X(f)，其運算公式根據第二章的內容有 dtetxfXFTftj2: dfefXtxIFTftj2:（2.199）（2.200）圖2.63 傅里葉變換的幾種類型 對于無限連續信號的傅里葉變換共有四種情況 ：n對于非周期連續信號X(t)，頻譜X(f)是連續譜；n對于周期連續信號，傅里葉變換轉變為傅里葉級數，因而其頻譜是離散的；n對于非周期離散信號，其傅里葉變換是一個周期性的連續頻譜；n對于周期

17、離散的時間序列，其頻譜也是周期離散的 。結論：若x(t)是周期的，頻域中X(f)必然是離散的，反之亦然。若x(t)是非周期的，則X(f)一定是連續的，反之亦然。第四種亦即時域和頻域都是離散的信號，且都是周期的，給我們利用計算機實施頻譜分析提供了一種可能性。對這種信號的傅里葉變換，我們只需取其時域上一個周期（N個采樣點）和頻域一個周期（同樣為N個采樣點）進行分析，便可了解該信號的全部過程。 DFT的定義：對有限長度的離散時域或頻域信號序列進行傅里葉變換或逆變換，得到同樣為有限長度的離散頻域或時域信號序列的方法，便稱為離散傅里葉變換（digital Fourier transform, DFT）或

18、其逆變換（IDFT）。離散傅里葉變換的公式：式中x(n)和X(k)分別為 和 的一個周期，此處將t和f0均歸一化為1。 nkNNonNonnkNjWnxenxkX112 11211NoKNoKnkNnkNjWkXNekXNnx（2.205）（2.206）NjNeW2)( tnx)(0kfX離散傅里葉變換意義：可以對任意連續的時域信號進行采樣和截斷并對其作離散傅里葉變換的運算，得到離散的頻譜，該頻譜的包絡即是對原連續信號真正頻譜的估計。 離散傅里葉變換的過程：n時域采樣(sampling in t-domain)；n時域截斷(truncation in t-doman)；n頻域采樣(sampli

19、ng in f-domain)。 圖2.64 離散傅里葉變換的圖解過程 （一）圖2.64 離散傅里葉變換的圖解過程（二） 圖2.64 離散傅里葉變換的圖解過程（三） 三、采樣定理(sampling theorem) 混疊（aliasing）：若采樣率過低即采樣間隔大，則系列的離散時間序列可能不能真正反映原始信號的波形特征，在頻域處理時會出現頻率混淆。 圖2.65 不同采樣率對采樣信號產生的影響（一） 圖2.65 不同采樣率對采樣信號產生的影響（二） 采樣定理：為避免混疊產生，要求的采樣頻率fs必須高于信號頻率成分中最高頻率fmax的兩倍，即 乃奎斯特（Nyquist）頻率：在給定的采樣頻率fs

20、條件下，信號中能被分辨的最高頻率。只有低于乃奎斯特頻率的頻率成分才能被精確地采樣，亦即為避免頻率混淆，應使被分析信號的最高頻率fmax低于乃奎斯特頻率。 max2 ffs（2.222）2sNyqff（2.223）圖2.66 混疊產生的條件 四、泄漏(leakage)與加窗(windowing) 圖2.67 余弦信號加窗截斷造成的泄漏現象 抑制或減小泄漏效應的方法：選擇性能更好的特殊窗來替代矩形窗，亦即加窗處理。評價窗函數的性能指標：1.3dB帶寬B：它是主瓣歸一化的幅值下降至-3dB時的帶寬。歸一化|W(f)|=20lg|W(f)/W(0)|，帶寬B的單位為或f。 2.旁瓣幅度A(dB)，表示

21、為最大旁瓣峰值Asmax與主瓣峰值Am之比，即20lg(Asmax/Am)。3.旁瓣峰值衰減率D（dB/decade），表示為最大旁瓣峰值與相距十倍頻處的旁瓣峰值之比，也是以分貝表示。 理想的窗函數應具有最小的B和A以及最大的D。 圖2.69 常用窗函數的時域圖像 圖2.70 常用窗函數的頻譜 五、柵欄效應(picket fence effect) 柵欄效應：若信號中某頻率成分的頻率fi等于，k/T即它與輸出的頻率采樣點相重合，那么該譜線便可被精確地顯示出來；反之若fi與頻率采樣點不重合，便得不到顯示，所得的頻譜便會產生誤差。頻率分辨率f：兩條譜線間的距離。當被分析的時域信號長度T（即窗寬T=

22、NTs）和采樣頻率fs被確定之后，則頻率分辨f也被確定： TNffs1（2.231）例：對余弦信號cos2f0t作DFT。圖2.72 周期信號作整周期截取的DFT（一） 圖2.72 周期信號作整周期截取的DFT （二）圖2.73 周期函數作非整周期截取的DFT 結論：對周期信號作整周期截取是獲取正確頻譜的先決條件。 六、快速傅里葉變換（FFT） 離散傅里葉變換的計算公式為：式中 N個點的X(k)需做N2次復數乘法和N(N-1)次復數加法。而做一次復數乘法需要做四次實數相乘和兩次實數相加，做一次復數加法需要做兩次實數相加。 例：N=1024時，則需要總共1,048,576次復數乘，即4,194,

23、304次實數乘法。 112:NonnkNNonnkNjWnxenxkXDFT 11211:NoknkNNoknkNjWkXNekXNnxIDFT1, 1 ,0Nk1, 1 ,0Nn（2.232）（2.233）NjNeW2快速傅里葉變換(FFT，Fast Fourier Transform)算法的本質:充分利用因子WN的周期性和對稱性。n對稱性： n周期性：FFT算法的基本思想：避免運算中的重復運算，將長序列的DFT分割為短序列的DFT的線性組合，從而達到整體降低運算量的目的。效果：使原來的N點DFT的乘法計算量從N2次降至為N/2log2N次，如N=1024，則計算量現在為5120次，僅為原計

24、算量的4.88% 。nkNNnkNWW2nkNnkNNWW（2.234）（2.235）時間抽取(decimation-in-time)基2算法 對式(2.232)，令N=2M,將x(n)序列分割成長度各為N/2的奇序列和偶序列，即令n=2r和n=2r+1,，r=0,1, ,N/2-1則式（2.232）重寫為式中 這是因為 1221221212122122122NorrkNNorkNrkNNorNorrNrkNWrxWWrxkWrxWrxkX（2.236）NjNjNeeW4222222222NNjNjNWeeW令則式（2.236）可改寫為而 因此將式（2.239）完整地寫成 12212, 1 ,

25、 02NorrkNNkWrxkA 12212, 1 ,012NorrkNNkWrxkB（2.237）（2.238） 1, 1 ,0NkkBWkAkXkN（2.239） kBkNBkAkNA22 12,2, 1 ,022NkkBWkANkXkBWkAkXNkNkN（2.240）又因為 ，因此最終可得 12NNW 12,2, 1 ,02NkkBWkANkXkBWkAkXkNkN（2.241）圖2.74 分割一次后的A(k)、B(k)及X(k)之間的關系（N=8） 按照上述思路繼續對A(k)和B(k)作奇偶序列分解。令r=2l，r=2l+1，l=0,1, ,N/4-1，則有：令 則 14042120414012214022244244NllkNkNNllkNNlklNNllkNWlxWWlxWlxWlxkA 14, 1 , 04414NkWxkCkNNo（2.242） 14, 1 ,02