MATLAB在機械振動信號中的應(yīng)用

1、MATLAB在機械振動信號中的應(yīng)用中振（山東理工大學(xué)交通與車輛工程學(xué)院）摘要：綜述了現(xiàn)代信號分析處理理論、方法如時域分析（包括時域參數(shù)識別、相關(guān)分析等）、頻域分析（包括傅立葉變換、功率譜分解等），并結(jié)合MATLAB中的相關(guān)函數(shù)來對所擬合的振動信號進(jìn)行時域分析和頻域分析，并對繪出的頻譜圖進(jìn)行說明。關(guān)鍵詞：時域分析頻域分析MATLAB信號是信息的載體，采用合適的信號分析處理方法以獲取隱藏于傳感觀測信號中的 重要信息（包括時域與頻域信息等），對于許多工程應(yīng)用領(lǐng)域均具有重要意義。對獲取振 動噪聲信號的分析處理，是進(jìn)行狀態(tài)監(jiān)測、故障診斷、質(zhì)量檢查、源識別、機器產(chǎn)品的 動態(tài)性能測試與優(yōu)化設(shè)計等工作的重要環(huán)

2、節(jié)，它可以預(yù)先發(fā)現(xiàn)機械部件的磨損和缺陷等 故障，從而可以提高產(chǎn)品的質(zhì)量，降低維護(hù)費用。隨著測試技術(shù)的迅速發(fā)展，各種信號 分析方法也隨之涌現(xiàn)，并廣泛應(yīng)用在各個領(lǐng)域10時域描述簡單直觀，只能反映信號的幅值隨時間的變化，而不能明確的揭示信號隨時間的變化關(guān)系。為了研究信號的頻率組成和各頻率成分的幅值大小、相位關(guān)系，應(yīng)對 信號進(jìn)行頻譜分析，即把時域信號通過適當(dāng)?shù)臄?shù)學(xué)方法處理變成頻率f （或角頻率 ）為獨立變量， 相應(yīng)的幅值或相位為因變量的頻域描述。頻域分析法將時域分析法中的 微分或差分方程轉(zhuǎn)換為代數(shù)方程，有利于問題的分析 2 0MATLAB是MathWorks 公司于1982年推出的一種功能強大、 效率

3、高、交互性好 的數(shù)值計算和可視化計算機高級語力，它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機地融合為一體，形成了一個極其方便、用戶界面良好的操作環(huán)境。隨著其自身版本的不斷提高，MATLAB的功能越來越強大，應(yīng)用范圍也越來越廣，如廣泛應(yīng)用于信號處理、數(shù)字圖像處理、仿真、自動化控制、小波分析及神經(jīng)網(wǎng)絡(luò)等領(lǐng)域30本文主要運用了 MATLAB R2014a對機械振動信號進(jìn)行分析。分析過程包括時域分析和頻域分析兩大部分，時域分析的指標(biāo)包括隨機信號的均值、方差以及均方值。頻 域分析的性能指標(biāo)包括對功率譜分析、倒頻譜分析。在進(jìn)行上述分析之前先要對振動信 號進(jìn)行擬合。機械振動分為確定性振動和隨機振動，確定性

4、振動又分為周期振動和非周 期振動，周期振動又進(jìn)一步分為簡諧振動和復(fù)雜的周期振動。所以可以根據(jù)上述的分類來擬合振動信號2。在設(shè)計信號的處理程序時，運用 MATLAB中的相關(guān)函數(shù)來對所擬合的振動信號進(jìn)行時域分析和頻域分析，并對繪出的頻譜圖進(jìn)行說明1時域分析1.1均值 對于一個各態(tài)歷經(jīng)隨機隨機信號x(t),其均值x為(D1 TTim T 0x(t)dt式中Mt)樣本函數(shù);T 觀測時間;常值分量。1.2方差2是描述隨機信號的波動分量,定義為1 TTim T ox(t)2x出(1)稱為標(biāo)準(zhǔn)差它表示信號x(t)偏離其均值x平方的均值，方差的正平方根1.3均方值X2是隨機信號X(t)平方的平均值，定義為(3

5、)21 T 2X Tim T 0 X (t)dt它描述信號的能量或強度，是X(t)平方的均值。均方值的正平方根值稱為均方根值X rms參數(shù)X、2、 X2之間的關(guān)系為(4)1.4時域統(tǒng)計分析 概率密度分析是以幅值大小為橫坐標(biāo)，以每個幅值間隔內(nèi)出現(xiàn)的概率 為縱坐標(biāo)進(jìn)行統(tǒng)計分析的方法。它反映了信號落在不同幅值強度區(qū)域內(nèi)的概率情況。計算方法如下:P(X)lim PX X(t) XXX 0limIlim0(Tti)(5)概率密度函數(shù)p(X)給出了信號取不同幅值大小的概率，是隨機信號的主要特征參數(shù)之 一。不同的隨機信號有不同的概率密度函數(shù)圖形，可以借此來識別信號的性質(zhì)，如正弦 信號加隨機噪聲、窄帶隨機信號

6、及寬帶隨機信號等。概率分布函數(shù)是信號幅值小于或等 于某值R的概率，定義為：(6)F (x) p(X)dX概率分布函數(shù)又稱為累積概率函數(shù)，表示了信號幅值落在某一區(qū)間的概率4 O2頻域分析2.1 傅里葉變換任何周期函數(shù)，均可展開成正交函數(shù)線性組合的無窮級數(shù)，如三角函數(shù)集的傅里葉 級數(shù)葉級數(shù)的表達(dá)形式如下:(7)(8)(9)，、 .， 、x(t) AnSin(n 0tn)2 n 1A0a0Andan b：(n 1,2,3,L)ann arctan 一 bn對于非周期信號或瞬變信號，利用如下的傅立葉變換進(jìn)行頻譜分析:X(f) x(t)ej2ftdtx(t) X(f)ej2 ftdf2.2 功率譜分析2

7、.2.1 經(jīng)典功率譜估計方法若x(t)為平穩(wěn)隨機信號，當(dāng)自相關(guān)函數(shù)為絕對可積時，自相關(guān)函數(shù)Rxx()和功率譜密度Sx()為一個傅里葉變換對，即Sx( )Rxx( )e j d1(10)Rxx( ) 2-Sx( )ej d同理，在頻域描述兩個隨機信號x(t)和y(t)相互關(guān)聯(lián)程度的數(shù)字特征，可以定義為 互譜功率密度簡稱互譜密度。而且，互相關(guān)函數(shù)與互譜密度是一個傅里葉變換對。Sxy( )Rxy( )e j d1(11)Rxy( ) Sxy( )ej d2.2.2 改進(jìn)的直接估計法直接法和間接法的方差性能很差，而且當(dāng)數(shù)據(jù)長度太大時，譜曲線起伏加?。蝗魯?shù) 據(jù)長度太小，則譜的分辨率又不好，所以需要改進(jìn)3

8、。提高的周期圖法估計的另一種方式就是采用對采樣數(shù)據(jù)分段使用非矩形窗，即Welch法。由于非矩形窗在邊沿趨近于零，從而減少了分段對重疊的依賴。選擇合適的 窗函數(shù)，采用每段一半的重疊率能大大降低譜估計的方差。這種方法中，記錄數(shù)據(jù)仍分 成K N段，即Kx(i)(n) x(n iM N) 0 n M 1,1 i K(12)每段M個取樣。窗函數(shù)w(n)在計算周期圖之前就與數(shù)據(jù)段相乘，于是可定義K個修正周期圖M M 12JM)( )x(i)(n)w(n)e j n i 1,2,L ,K(13)MU n 0U是窗口序列函數(shù)的平均能量1M 1 2U w2(n)(14)M n 0則定義譜估計為、4，1 K小Bx

9、( ) - JM彳 )(15)K i 12.2.3 AR模型功率譜估計法傳統(tǒng)的功率譜估計方法是利用加窗的數(shù)據(jù)或加窗的相關(guān)函數(shù)估計值的傅里葉變換 來計算的，具有一定的優(yōu)勢，如計算效率高，估計值正比于正弦波信號的功率等。但是 同時也存在許多缺點，主要缺點就是方差性能較差、譜分辨率低。而參數(shù)模型法可大大提高功率譜估計的分辨率，是現(xiàn)代譜估計的主要研究內(nèi)容，在語音分析、數(shù)據(jù)壓縮以及 通信等領(lǐng)域有著廣泛的應(yīng)用網(wǎng)按照模型化進(jìn)行功率譜估計，其主要思想如下：(1)選擇模型；(2)從給出的數(shù)據(jù)樣本估計假設(shè)的模型；(3)將估計的模型參數(shù)打入模型的理論功率譜公式中得出一個較好的譜估計值。假設(shè)產(chǎn)生隨機序列x(n)的系統(tǒng)

10、模型為一個線性差分方程，即x(n)biw(n i) ajx(n j)(16)式中，w(n)表示白噪聲序列，對上式進(jìn)行Z變換，可得所以系統(tǒng)的傳遞函數(shù)為q kajX(z)zj 0qbW (z)z ii 0(17)式中,H(z)器B(z)A(z)(18)假定輸入白噪聲功率譜密度為Pw (z)qA(z)ajz j j oqB(z)bz ii 0(19)(20)。，那么輸出功率譜密度為_1(21)Px(Z)W*旦()又根據(jù)z ej ,所以得(22)這樣，當(dāng)確定了系數(shù)aj> b和W后，就可以求解得到隨機信號的功率譜密度 px( )了通過上式可知，如果i 1 , bi0時，則系統(tǒng)的差分方程變?yōu)閝(23

11、)x(n) ajx(n j) w(n)j i上式即為自回歸模型，簡稱為 AR (Auto-Regressive )模型，再將該式進(jìn)行Z變換,H(z)X(z) 1W(z) A(z)(24)所以，AR模型又稱為全極點模型。AR模型的輸出功率譜為R()2wA(ej )qj k1 ajej 1(25)顯然，計算出w和aj后，就可以求解得到隨機信號的功率譜密度 Px()本文采用AR模型的一種Burg法進(jìn)行功率譜估計。3仿真研究仿真帶噪聲信號如下：x(t) 6e 1.0t sin(2 f1t) 8e2.0t sin(2 fzt) randn(t)該仿真帶噪聲信號由兩個正弦信號6e1.0tsin(2m)、8

12、e 2.0t sin(2 fR和一個服從正態(tài)分布的高斯白噪聲信號randn(t)疊加而成。G 100Hz, f2 300Hz。其時域波形如圖1所示（程序詳見附錄1）圖1時域波形圖時域分析結(jié)果：序列的平均值為0.5050序列的最小值為-10.7448序列的最大值為12.0222序列的標(biāo)準(zhǔn)差為2.9153序列的方差為8.4992序列的均方值為2.9580*圖2經(jīng)典功率譜估算圖陋幡N露我在功率譜中可以很明顯的看到振動信號中有 100Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分。如圖 2圖3 FFT頻譜圖上圖3為FFT頻譜圖，從該頻譜中可以看到有三個主要高峰值，即在0Hz , 1

13、00Hz , 300Hz 處。用Burg法進(jìn)行PSD估計功率譜圖如圖4,從中可以很明顯的看到振動信號中有100Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分(程序詳見附錄3):圖4 Burg法進(jìn)行PSD估計功率譜圖在Welch法進(jìn)行PSD功率譜估計，當(dāng)采用不同窗函數(shù)時的結(jié)果。從中可以很明顯的看到振動信號中有100Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分。且海寧窗和布萊克曼窗較為明顯（程序詳見附錄2）。圖5 Welch法進(jìn)行PSD功率譜估計功率譜圖CD1sc圖5倒譜圖理論上，傅立葉變換用于頻譜分析，可以找出受噪聲干擾的信號的頻率成分，而這 用時域分析

14、是不能分辨的。對傅立葉變換做復(fù)共腕運算，即可得到信號的功率譜密度函 數(shù)，以顯示各頻率分量的能量分布。仿真帶噪信號的傅立葉變換與功率譜分解結(jié)果如圖 3和圖4、5、6所示。從圖3和圖4、5、6可以清楚看到，約在頻率為100Hz、300Hz （即振動信號頻率的倍頻）處頻譜幅值和能量出現(xiàn)局部極大值，對應(yīng)機械振動的主振動源所在。4結(jié)論信號是信息的載體，因此采用合適的信號分析處理方法以獲取隱藏于傳感觀測信號 中的重要信息（包括時域與頻域信息等），對于許多工程應(yīng)用領(lǐng)域均具有重要意義。本文 在研究現(xiàn)代信號分析處理理論、方法如時域分析（包括時域參數(shù)識別、相關(guān)分析以及統(tǒng)計分析等）、頻域分析（包括傅立葉變換、功率譜

15、分解等）的基礎(chǔ)上，結(jié)合仿真數(shù)據(jù)對機械 振動信號分析處理，具有一定的參考價值。參考文獻(xiàn)1馮凱.工程測試技術(shù)M.西安：西北工業(yè)大學(xué)出版社，2003.2許同樂.機械工程測試技術(shù).北京：機械工業(yè)出版社，2010.3薛年喜.MATLAB在數(shù)字信號處理中的應(yīng)用(第二版).北京：清華大學(xué)出版社，20084焦衛(wèi)東.旋轉(zhuǎn)機械振動信號分析.浙江.嘉興學(xué)院學(xué)報.2007.附錄附錄一：時域分析、頻域分析程序A1=6;A2=8;f1=100;f2=300;fs=1000;t=0:1/fs:2;N=length(t);X1=A1*exp(-1.0*t).*sin(2*pi*f1*t);X2=A2*exp(-2.0*t).

16、*sin(2*pi*f2*t);R=rand(1,N);Y=X1+X2+R;figure(1);plot(t,Y);title('振動信號的波形');xlabel('時間/秒')；ylabel('幅度');grid; hold on;%時域分析mi=min(Y); disp(mi);% 最小值mx=max(Y); disp(mx);% 最大值st=std(Y); disp(st);% 標(biāo)準(zhǔn)差m=mean(Y); disp(m); % 均值vr=var(Y); disp(vr);% 方差rm=rms(Y); disp(rm);% 均方值%頻域分析l

17、=length(Y);r=fft(Y)/l尸fftshift(r);f=linspace(-fs/2,fs/2,l);figure(2);plot(f,abs(r); grid; hold on;figure(3);psd(Y,2048,1000,kaiser(512,5),0,0.95);figure(4);yc=rceps(Y);plot(yc);附錄二：Welch方法進(jìn)行PSD估計程序A1=6;A2=8;f1=100;f2=300;fs=1000;nfft=1024;t=0:1/fs:2;N=length(t);X1=A1*exp(-1.0*t).*sin(2*pi*f1*t);X2=A

18、2*exp(-1.5*t).*sin(2*pi*f2*t);R=rand(1,N);Y=X1+X2+R;window1=boxcar(100);window2=hamming(100);window3=blackman(100);noverlap=20;Pxx1,f1=pwelch(Y,window1,noverlap,nfft,fs);Pxx2,f2=pwelch(Y,window2,noverlap,nfft,fs);Pxx3,f3=pwelch(Y,window3,noverlap,nfft,fs);PXX1=10*log10(Pxx1);PXX2=10*log10(Pxx2);PXX3=10*log10(Pxx3);subplot(3,1,1)plot(f1,PXX1);title('矩形窗')；subplot(3,1,2)plot(f2,PXX2);title('海寧窗')；subplot(3,1,3)plot(f3,PXX3);xlabel('頻率(Hz)');ylabel('幅度(dB)');title('布萊克曼窗&#