MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用

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

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

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

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

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

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

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

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

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

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

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

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

13、007.附錄附錄一：時(shí)域分析、頻域分析程序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('振動(dòng)信號(hào)的波形');xlabel('時(shí)間/秒');ylabel('幅度');grid; hold on;%時(shí)域分析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方法進(jìn)行PSD估計(jì)程序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方法進(jìn)行PSD估計(jì)程序A1=6;A2=8;f1=100;