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

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

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

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

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

5、析 概率密度分析是以幅值大小為橫坐標(biāo)，以每個幅值間隔內(nèi)出現(xiàn)的概率為縱坐標(biāo)進行統(tǒng)計分析的方法。它反映了信號落在不同幅值強度區(qū)域內(nèi)的概率情況。計算方法如下： （5） 概率密度函數(shù)給出了信號取不同幅值大小的概率，是隨機信號的主要特征參數(shù)之一。不同的隨機信號有不同的概率密度函數(shù)圖形，可以借此來識別信號的性質(zhì)，如正弦信號加隨機噪聲、窄帶隨機信號及寬帶隨機信號等。概率分布函數(shù)是信號幅值小于或等于某值的概率，定義為： （6）概率分布函數(shù)又稱為累積概率函數(shù)，表示了信號幅值落在某一區(qū)間的概率4。2 頻域分析2.1 傅里葉變換任何周期函數(shù)，均可展開成正交函數(shù)線性組合的無窮級數(shù)，如三角函數(shù)集的傅里葉級數(shù)。 葉級數(shù)的

6、表達形式如下： （7） （8）對于非周期信號或瞬變信號，利用如下的傅立葉變換進行頻譜分析： （9）2.2 功率譜分析2.2.1 經(jīng)典功率譜估計方法若為平穩(wěn)隨機信號，當(dāng)自相關(guān)函數(shù)為絕對可積時，自相關(guān)函數(shù)和功率譜密度為一個傅里葉變換對，即 （10）同理，在頻域描述兩個隨機信號和相互關(guān)聯(lián)程度的數(shù)字特征，可以定義為互譜功率密度簡稱互譜密度。而且，互相關(guān)函數(shù)與互譜密度是一個傅里葉變換對。 （11）2.2.2 改進的直接估計法直接法和間接法的方差性能很差，而且當(dāng)數(shù)據(jù)長度太大時，譜曲線起伏加??；若數(shù)據(jù)長度太小，則譜的分辨率又不好，所以需要改進3。提高的周期圖法估計的另一種方式就是采用對采樣數(shù)據(jù)分段使用非矩形

7、窗，即Welch法。由于非矩形窗在邊沿趨近于零，從而減少了分段對重疊的依賴。選擇合適的窗函數(shù)，采用每段一半的重疊率能大大降低譜估計的方差。這種方法中，記錄數(shù)據(jù)仍分成段，即 （12）每段個取樣。窗函數(shù)在計算周期圖之前就與數(shù)據(jù)段相乘，于是可定義個修正周期圖 （13）是窗口序列函數(shù)的平均能量 （14）則定義譜估計為 （15）2.2.3 AR模型功率譜估計法傳統(tǒng)的功率譜估計方法是利用加窗的數(shù)據(jù)或加窗的相關(guān)函數(shù)估計值的傅里葉變換來計算的，具有一定的優(yōu)勢，如計算效率高，估計值正比于正弦波信號的功率等。但是同時也存在許多缺點，主要缺點就是方差性能較差、譜分辨率低。而參數(shù)模型法可大大提高功率譜估計的分辨率，是

8、現(xiàn)代譜估計的主要研究內(nèi)容，在語音分析、數(shù)據(jù)壓縮以及通信等領(lǐng)域有著廣泛的應(yīng)用3。按照模型化進行功率譜估計，其主要思想如下：（1） 選擇模型；（2） 從給出的數(shù)據(jù)樣本估計假設(shè)的模型；（3） 將估計的模型參數(shù)打入模型的理論功率譜公式中得出一個較好的譜估計值。假設(shè)產(chǎn)生隨機序列的系統(tǒng)模型為一個線性差分方程，即 （16）式中，表示白噪聲序列，對上式進行Z變換，可得 （17）所以系統(tǒng)的傳遞函數(shù)為 （18） 式中， （19） （20）假定輸入白噪聲功率譜密度為，那么輸出功率譜密度為 （21）又根據(jù)，所以得 （22）這樣，當(dāng)確定了系數(shù)、和后，就可以求解得到隨機信號的功率譜密度了。通過上式可知，如果，時，則系統(tǒng)的

9、差分方程變?yōu)?（23）上式即為自回歸模型，簡稱為AR（Auto-Regressive）模型，再將該式進行Z變換，得 （24）所以，AR模型又稱為全極點模型。AR模型的輸出功率譜為 （25）顯然，計算出和后，就可以求解得到隨機信號的功率譜密度。本文采用AR模型的一種Burg法進行功率譜估計。3 仿真研究仿真帶噪聲信號如下：該仿真帶噪聲信號由兩個正弦信號、和一個服從正態(tài)分布的高斯白噪聲信號疊加而成。其時域波形如圖1所示（程序詳見附錄1）。圖1 時域波形圖時域分析結(jié)果：序列的平均值為 0.5050 序列的最小值為 -10.7448序列的最大值為 12.0222序列的標(biāo)準(zhǔn)差為 2.9153序列的方差為

10、 8.4992序列的均方值為 2.9580圖2 經(jīng)典功率譜估算圖在功率譜中可以很明顯的看到振動信號中有100Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分。如圖2圖3 FFT頻譜圖上圖3為FFT頻譜圖，從該頻譜中可以看到有三個主要高峰值，即在0Hz，100Hz，300Hz處。用Burg法進行PSD估計功率譜圖如圖4，從中可以很明顯的看到振動信號中有100Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分（程序詳見附錄3）：圖4 Burg法進行PSD估計功率譜圖在Welch法進行PSD功率譜估計，當(dāng)采用不同窗函數(shù)時的結(jié)果。從中可以很明顯的看到振動信號中有100

11、Hz和300Hz兩個主要的頻率。表明信號中含有這兩個頻率的周期成分。且海寧窗和布萊克曼窗較為明顯（程序詳見附錄2）。圖5 Welch法進行PSD功率譜估計功率譜圖圖5 倒譜圖理論上，傅立葉變換用于頻譜分析，可以找出受噪聲干擾的信號的頻率成分，而這用時域分析是不能分辨的。對傅立葉變換做復(fù)共軛運算，即可得到信號的功率譜密度函數(shù)，以顯示各頻率分量的能量分布。仿真帶噪信號的傅立葉變換與功率譜分解結(jié)果如圖3和圖4、5、6所示。從圖3和圖4、5、6可以清楚看到，約在頻率為100Hz、300Hz（即振動信號頻率的倍頻)處頻譜幅值和能量出現(xiàn)局部極大值，對應(yīng)機械振動的主振動源所在。4 結(jié)論信號是信息的載體，因此

12、采用合適的信號分析處理方法以獲取隱藏于傳感觀測信號中的重要信息(包括時域與頻域信息等)，對于許多工程應(yīng)用領(lǐng)域均具有重要意義。本文在研究現(xiàn)代信號分析處理理論、方法如時域分析 (包括時域參數(shù)識別、相關(guān)分析以及統(tǒng)計分析等)、頻域分析(包括傅立葉變換、功率譜分解等)的基礎(chǔ)上，結(jié)合仿真數(shù)據(jù)對機械振動信號分析處理，具有一定的參考價值。參考文獻1 馮凱.工程測試技術(shù)M.西安：西北工業(yè)大學(xué)出版社，2003.2 許同樂.機械工程測試技術(shù).北京：機械工業(yè)出版社，2010.3 薛年喜.MATLAB在數(shù)字信號處理中的應(yīng)用（第二版）.北京：清華大學(xué)出版社，20084 焦衛(wèi)東.旋轉(zhuǎn)機械振動信號分析.浙江.嘉興學(xué)院學(xué)報.2

13、007.附錄附錄一：時域分析、頻域分析程序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).*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);%最小值m

14、x=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=length(Y);r=fft(Y)/l;r=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

15、); 附錄二：Welch方法進行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=A2*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

16、,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('布萊克曼窗'); 附錄三：Burg方法進行PSD估計程序A1=6;A2=8;f1=100;