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

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用申振（山東理工大學(xué)交通與車輛工程學(xué)院）摘要：綜述現(xiàn)代信號(hào)分析處理理論方法如時(shí)域分（包括時(shí)域參數(shù)識(shí)別、相分等域分析（包括傅立葉變換功率譜分解)并合MATLAB中的相關(guān)函數(shù)來(lái)對(duì)所擬合的振動(dòng)信號(hào)進(jìn)時(shí)分析和頻域分析并繪的頻譜圖進(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é)，它可以預(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é)方法處理變成頻率

f

(或頻率）為獨(dú)立變量，相應(yīng)的幅值或相位為因變量的頻域描述。頻域分析法將時(shí)域分析法中的微分或差分方程轉(zhuǎn)換為代數(shù)方程，有利于問(wèn)題的分析[2

.MATLAB是MathWorks公司于1982年出的一種功能強(qiáng)效高交互性好的數(shù)值計(jì)算和可視化計(jì)算機(jī)高級(jí)語(yǔ)言,它將數(shù)值分析、矩陣運(yùn)算、信號(hào)處理和圖形顯示有機(jī)地融合為一體形成了一個(gè)極其方便、用戶界面良好的操作環(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對(duì)機(jī)械振信號(hào)進(jìn)行分析分析過(guò)程包括時(shí)域分析和頻域分析兩大部分，時(shí)域分析的指標(biāo)包括隨機(jī)信號(hào)的均值、方差以及均方值。頻域分析的性能指標(biāo)包1

T12Txx(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用T12Txx括對(duì)功率譜分析、倒頻譜分析。在進(jìn)行上述分析之前先要對(duì)振動(dòng)信號(hào)進(jìn)行擬合.機(jī)械振動(dòng)分為確定性振動(dòng)和隨機(jī)振動(dòng)確定性振動(dòng)又分為周期振動(dòng)和非周期振動(dòng)，周期振動(dòng)又進(jìn)一步分為簡(jiǎn)諧振動(dòng)和復(fù)雜的周期振動(dòng)。所以可以根據(jù)上述的分類來(lái)擬合振動(dòng)信號(hào)[

.在計(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)(t)

，其均值為x

x

T

1

0

x(t

（1）式中xt)

——樣本函數(shù)T

-—測(cè)時(shí)間；—-常值分量。x1.2方差是述隨機(jī)信號(hào)的波動(dòng)分量定義為xlimt)]T0

（1）它表示信號(hào)(t)

偏離其均值平方的均值，方差的正平方根稱為標(biāo)準(zhǔn)差。xx1均值

x

是隨機(jī)信號(hào)(t)

平方的平均值，定義為

x

2

limT

1

T0

()

(3)它描述信號(hào)的能量或強(qiáng)度是(t)

平方的均值方值的正平方根值稱為均根值

rms

參數(shù)、、之間的關(guān)系為xxx

2x

-

x

2

（4）1.4

時(shí)域統(tǒng)計(jì)分析

概率密度分析是以幅值大小為橫坐標(biāo)，以每個(gè)幅值間隔內(nèi)出現(xiàn)的概率為縱坐標(biāo)進(jìn)行統(tǒng)計(jì)分析的方法。它反映了信號(hào)落在不同幅值強(qiáng)度區(qū)域內(nèi)的概率情況。計(jì)算方法如下：(0

P[xx(

lim

/T

(0

i

)i

（5)2

n2nnn2nn概率密度函數(shù)()

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用給出了信號(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)幅值小于或等于某值R的概率定義為：F)

)dx

(6）概率分布函數(shù)又稱為累積概率函數(shù),表示了信號(hào)幅值落某一區(qū)間的概率［］

.22.1

頻域分析傅里葉變換任何周期函數(shù)，均可展開(kāi)成正交函數(shù)線性組合的無(wú)窮級(jí)數(shù)，如三角函數(shù)集的傅里葉級(jí)數(shù)。葉級(jí)數(shù)的表達(dá)形式如下：a(t)Asin(nt0n

（7)A022nnarctann

（n

（8）對(duì)于非周期信號(hào)或瞬變信號(hào)，利用如下的傅立葉變換進(jìn)行頻譜分析：

)t)ejdtxt)X)ej

2.2

功率譜分析2

經(jīng)典功率譜估方法若xt)

為平穩(wěn)隨機(jī)信號(hào)自相關(guān)函數(shù)為絕對(duì)可時(shí)相關(guān)函數(shù)R)和功率譜密度S()xx為一個(gè)傅里葉變換對(duì)即3

K(完整版)MATLAB在機(jī)KS(1((e

（10)同理，在頻域描述兩個(gè)隨機(jī)信號(hào)xt)

和(t)

相互關(guān)聯(lián)程度的數(shù)字特征，可以定義為互譜功率密度簡(jiǎn)稱互譜密度。而且互相關(guān)函數(shù)與互譜密度是一個(gè)傅里葉變換對(duì).S(1((2

（11）2.2。2

改進(jìn)的直接估法直接法和間接法的方差性能很差而且當(dāng)數(shù)據(jù)長(zhǎng)度太大時(shí),曲線起伏加?。蝗魯?shù)據(jù)長(zhǎng)度太小，則譜的分辨率又不好所以需要改進(jìn)[提高的周期圖法估計(jì)的另一種方式就是采用對(duì)采樣數(shù)據(jù)分段使用非矩形窗,即Welch法。由于非矩形窗在邊沿趨近于零，從而減少了分段對(duì)重疊的依賴。選擇合適的窗函數(shù)，采用每段一半的重疊率能大大降低譜估計(jì)的方差。這種方法中，記錄數(shù)據(jù)仍分成

K

段，即x(i)()x(N)

0K

（12每段M個(gè)取。窗函數(shù)w()

在計(jì)算周期圖之前就與數(shù)據(jù)段相乘，于是可定義K

個(gè)修正周期圖J(i))M

1

Mn

x(i)(n)(e

2

iK

（13)U

是窗口序列函數(shù)的平均能量U

1

Mn

w2(n)

(14）則定義譜估計(jì)為w

1K

i

(iM

（15）2.2.3AR模型功率譜估計(jì)法傳統(tǒng)的功率譜估計(jì)方法是利用加窗的數(shù)據(jù)或加窗的相關(guān)函數(shù)估計(jì)值的傅里葉變換來(lái)計(jì)算4

q(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)q的，具有一定的優(yōu)勢(shì)，如計(jì)算效率高，估計(jì)值正比于正弦波信號(hào)的功率等。但是同時(shí)也存在許多缺點(diǎn)，主要缺點(diǎn)就是方差性能較差、譜分辨率低。而參數(shù)模型法可大大提高功率譜估計(jì)的分辨率是現(xiàn)代譜估計(jì)的主要研究?jī)?nèi)容語(yǔ)音分析據(jù)壓縮以及通信等領(lǐng)域有著廣泛的應(yīng)用[按照模型化進(jìn)行功率譜估計(jì)主要思想如下：（1）選擇模型；（2）從給出的數(shù)據(jù)樣本估計(jì)假設(shè)的模型;（3）將估計(jì)的模型參數(shù)打入模型的理論功率譜公式中得出一個(gè)較好的譜估計(jì)值。假設(shè)產(chǎn)生隨機(jī)序列x(n

的系統(tǒng)模型為一個(gè)線性差分方程即x(n((j)ijij

（16式中，(n)

表示白噪聲序列，對(duì)上式進(jìn)行Z變換可得

a()zj

bW(z)i

（17)j

i所以系統(tǒng)的傳遞函數(shù)為H(z)

X(z)W(z)

(z)(

（18式中，q(z)zjj

j

(19）q(z)bzii

（20）假定輸入白噪聲功率譜密度為(z),那么輸出功率譜密度為ww(

(zB((z)A

))

(21）又根據(jù)zj，所以得5

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

2

(ej(j)

2

（22這樣,當(dāng)確定了系數(shù)a、和2后就可以求解得到隨機(jī)信號(hào)的功率譜密度了通過(guò)上可jiw知，如果ii

時(shí)，則系統(tǒng)的差分方程變?yōu)閝x(nx(nj)(njj上式即為自回歸模型簡(jiǎn)稱為AR（Auto—Regressive）模型，再將該式進(jìn)行變換得

（23）H(

()z)Az)

zj

（24)j所以，AR型稱為全極點(diǎn)模型AR型輸出功率譜為

(j

)

qaj

j

2

（25)j顯然，計(jì)算出2和后,就可以求解得到隨機(jī)信號(hào)的功率譜密度p。wj本文采用型的一種Burg法進(jìn)行功率譜估計(jì)。3

仿真研究仿真帶噪聲信號(hào)如下：(t)tt)ett)randn(該仿真帶噪聲信號(hào)由兩個(gè)正弦信號(hào)e

t

ft1

e

t

ft2

和一個(gè)服從正態(tài)分布的高斯白噪聲信號(hào)(t

疊加而成。f100Hzf3001

。其時(shí)域波形如圖1所（程序詳見(jiàn)附錄16

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用圖1時(shí)域波形時(shí)域分析結(jié)果：序列的平均值為0。5050序列的最小值為-10.7448序列的最大值為12。0222序列的標(biāo)準(zhǔn)差為2。9153序列的方差為8。4992序列的均方值為2.9580圖2

經(jīng)典功率譜估圖在功率譜中可以很明顯的看到振動(dòng)信號(hào)中有100Hz和兩個(gè)主要的頻率。表明信號(hào)中含有這兩個(gè)頻率的周期成分。如圖27

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用圖3FFT頻圖上圖3為頻譜圖從該譜中可以看到有三個(gè)主要峰值,即在0Hz，300Hz處。用Burg法進(jìn)行PSD估計(jì)率譜圖如圖4，從中可以很明顯的看到振動(dòng)信號(hào)中有100Hz和300Hz個(gè)主要的頻率.表明號(hào)中含有這兩個(gè)頻率的周期成分（程序詳見(jiàn)附錄圖4Burg法進(jìn)行PSD估計(jì)功率譜圖在Welch法進(jìn)行PSD功率譜估計(jì)當(dāng)采用不同窗函數(shù)時(shí)的結(jié)果。從中可以很明顯的看到振動(dòng)信號(hào)中有100Hz和300Hz兩個(gè)主要的頻率。表明信號(hào)中含有這兩個(gè)頻率的周期成分。且海寧窗和布萊克曼窗較為明顯（程序詳見(jiàn)附錄8

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用圖5Welch法行PSD功率譜估計(jì)功率譜圖圖5倒譜圖理論上，傅立葉變換用于頻譜分析，可以找出受噪聲干擾的信號(hào)的頻率成分，而這用時(shí)域分析是不能分辨的。對(duì)傅立葉變換做復(fù)共軛運(yùn)算,可得到信號(hào)的功譜密度函數(shù)以顯示各頻率分量的能量分布。仿真帶噪信號(hào)的傅立葉變換與功率譜分解結(jié)果如圖和圖、5、6所示。從圖3和圖4、6可以清楚看到，約在頻率為100Hz、300Hz(即振動(dòng)信號(hào)頻率的倍頻）處頻譜幅值和能量出現(xiàn)局部極大值，對(duì)應(yīng)機(jī)械振動(dòng)的主振動(dòng)源所在。4

結(jié)論信號(hào)是信息的載體，因此采用合適的信號(hào)分析處理方法以獲取隱藏于傳感觀測(cè)信號(hào)中的重要信息包括時(shí)域與頻域信息等于許工程應(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)分析處理，具有一定9

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用的參考價(jià)值。參考文獻(xiàn)[1]凱工測(cè)試技術(shù)［M].西安:北業(yè)大學(xué)出版社，［2］許同樂(lè)。機(jī)械工程測(cè)試技術(shù)。北京：機(jī)械工業(yè)出版社,2010.［3]薛年喜MATLAB在數(shù)字信號(hào)處理中的應(yīng)用（第二版）北京：清華大學(xué)出版社,2008［4]焦衛(wèi)東。旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)分析。浙江.嘉興學(xué)院學(xué)報(bào).2007.附錄附錄一時(shí)域分析、頻域分析程序A1=6；A2=8f1=100;f2=300；fs=1000t=0:1/fs:2N=length(t）;X1=A1*exp（-1.0*t).＊sin＊pi*f1X2=A2*exp（—2。0＊t)＊f2＊t)；R=randY=X1+X2+R；figure(1)；plot(t，Y);title（’振動(dòng)信號(hào)的波形’xlabel('時(shí)間/秒')；ylabel（’幅度grid;on；％時(shí)域分析mi=min(Y);disp(mi）;%最小值mx=max（Y)；disp（mx)；％最大值st=std(Y)；disp標(biāo)準(zhǔn)差m=meandisp（m);％均值vr=var（Ydisp；％方差rm=rmsdisp(rm);％均方值10

(完整版)MATLAB在機(jī)械振動(dòng)信號(hào)中的應(yīng)用％頻域分析l=length（Yr=fft(Y)/l;r=fftshift(rf=linspace(—fs/2,fs/2,l)figure（2，abs(r)）;gridholdon;figure(3）;psd（Y,2048，kaiser(512，5),0,0.95figure(4)；yc=rcepsplot（yc附錄二Welch方法進(jìn)行估計(jì)程序A1=6；A2=8f1=100；f2=300；fs=1000；nfft=1024；t=0：1/fs；N=length(t)X1=A1＊exp（—1.0。*sin(2*pi＊f1＊tX2=A2*exp（—1。5*t)＊sin（2＊piR=rand（1，NY=X1+X2+R；window1=boxcar(100)；window2=hamming（10