清華大學(xué)測試技術(shù)第三章

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

24、算量的4.88% 。nkNNnkNWW2nkNnkNNWW（2.234）（2.235）時(shí)間抽取(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）重寫為式中 這是因?yàn)?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）又因?yàn)?，因此最終可得 12NNW 12,2, 1 ,02NkkBWkANkXkBWkAkXkNkN（2.241）圖2.74 分割一次后的A(k)、B(k)及X(k)之間的關(guān)系（N=8） 按照上述思路繼續(xù)對A(k)和B(k)作奇偶序列分解。令r=2l，r=2l+1，l=0,1, ,N/4-1，則有：令 則 14042120414012214022244244NllkNkNNllkNNlklNNllkNWlxWWlxWlxWlxkA 14, 1 , 04414NkWxkCkNNo（2.242） 14, 1 ,02