MATLAB在機械振動信號中的應用

(完整版)MATLAB在機械振動信號中的應用申振（山東理工大學交通與車輛工程學院）摘要：綜述現(xiàn)代信號分析處理理論方法如時域分（包括時域參數(shù)識別、相分等域分析（包括傅立葉變換功率譜分解)并合MATLAB中的相關(guān)函數(shù)來對所擬合的振動信號進時分析和頻域分析并繪的頻譜圖進行說明關(guān)鍵詞時域分析頻域分析MATLAB信號是信息的載體，采用合適的信號分析處理方法以獲取隱藏于傳感觀測信號中的重要信息括時域與頻域信息等對于許多工程應用領域均具有重要意義。對獲取振動噪聲信號的分析處理，是進行狀態(tài)監(jiān)測、故障診斷、質(zhì)量檢查、源識別、機器產(chǎn)品的動態(tài)性能測試與優(yōu)化設計等工作的重要環(huán)節(jié)，它可以預先發(fā)現(xiàn)機械部件的磨損和缺陷等故障，從而可以提高產(chǎn)品的質(zhì)量，降低維護費用。隨著測試技術(shù)的迅速發(fā)展,種信號分析方法也隨之涌現(xiàn),并廣泛應用在各個領域［1]。時域描述簡單直觀，只能反映信號的幅值隨時間的變化,而不能明確的揭示信號隨時間的變化關(guān)系。為了研究信號的頻率組成和各頻率成分的幅值大小、相位關(guān)系，應對信號進行頻譜分析把時域信號通過適當?shù)臄?shù)學方法處理變成頻率

f

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

.MATLAB是MathWorks公司于1982年出的一種功能強效高交互性好的數(shù)值計算和可視化計算機高級語言,它將數(shù)值分析、矩陣運算、信號處理和圖形顯示有機地融合為一體形成了一個極其方便、用戶界面良好的操作環(huán)境。隨著其自身版本的不斷提高，MATLAB功能越來越強大應用范圍也越來越廣，如廣泛應用于信號處理、數(shù)字圖像處理、仿真、自動化控制、小波分析及神經(jīng)網(wǎng)絡等領域[3]

.本文主要運用了MATLAB對機械振信號進行分析分析過程包括時域分析和頻域分析兩大部分，時域分析的指標包括隨機信號的均值、方差以及均方值。頻域分析的性能指標包1

T12Txx(完整版)MATLAB在機械振動信號中的應用T12Txx括對功率譜分析、倒頻譜分析。在進行上述分析之前先要對振動信號進行擬合.機械振動分為確定性振動和隨機振動確定性振動又分為周期振動和非周期振動，周期振動又進一步分為簡諧振動和復雜的周期振動。所以可以根據(jù)上述的分類來擬合振動信號[

.在計信號的處理程序時，運用MATLAB中的相關(guān)函數(shù)來對所擬合的振動信號進行時域分析和頻域分析，并對繪出的頻譜圖進行說明。1

時域分析1.1均值

對于一個各態(tài)歷經(jīng)隨機隨機信號(t)

，其均值為x

x

T

1

0

x(t

（1）式中xt)

——樣本函數(shù)T

-—測時間；—-常值分量。x1.2方差是述隨機信號的波動分量定義為xlimt)]T0

（1）它表示信號(t)

偏離其均值平方的均值，方差的正平方根稱為標準差。xx1均值

x

是隨機信號(t)

平方的平均值，定義為

x

2

limT

1

T0

()

(3)它描述信號的能量或強度是(t)

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

rms

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

2x

-

x

2

（4）1.4

時域統(tǒng)計分析

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

P[xx(

lim

/T

(0

i

)i

（5)2

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

(完整版)MATLAB在機械振動信號中的應用給出了信號取不同幅值大小的概率是隨機信號的主要特征參數(shù)之一不同的隨機信號有不同的概率密度函數(shù)圖形可以此來識信號的性質(zhì)，如正弦信號加隨機噪聲、窄帶隨機信號及寬帶隨機信號等。概率分布函數(shù)是信號幅值小于或等于某值R的概率定義為：F)

)dx

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

.22.1

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

（7)A022nnarctann

（n

（8）對于非周期信號或瞬變信號，利用如下的傅立葉變換進行頻譜分析：

)t)ejdtxt)X)ej

2.2

功率譜分析2

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

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

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

（10)同理，在頻域描述兩個隨機信號xt)

和(t)

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

（11）2.2。2

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

K

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

0K

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

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

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

1

Mn

x(i)(n)(e

2

iK

（13)U

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

1

Mn

w2(n)

(14）則定義譜估計為w

1K

i

(iM

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

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

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

（16式中，(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在機械振動信號中的應用Px

2

(ej(j)

2

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

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

（23）H(

()z)Az)

zj

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

(j

)

qaj

j

2

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

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

t

ft1

e

t

ft2

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

疊加而成。f100Hzf3001

。其時域波形如圖1所（程序詳見附錄16

(完整版)MATLAB在機械振動信號中的應用圖1時域波形時域分析結(jié)果：序列的平均值為0。5050序列的最小值為-10.7448序列的最大值為12。0222序列的標準差為2。9153序列的方差為8。4992序列的均方值為2.9580圖2

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

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

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

結(jié)論信號是信息的載體，因此采用合適的信號分析處理方法以獲取隱藏于傳感觀測信號中的重要信息包括時域與頻域信息等于許工程應用領域均具有重要意義本在研究現(xiàn)代信號分析處理理論方法如時域分析包括時域參數(shù)識別相關(guān)析以及統(tǒng)計分析等域分析(包括傅立葉變換、功率譜分解等）的基礎上，結(jié)合仿真數(shù)據(jù)對機械振動信號分析處理，具有一定9

(完整版)MATLAB在機械振動信號中的應用的參考價值。參考文獻[1]凱工測試技術(shù)［M].西安:北業(yè)大學出版社，［2］許同樂。機械工程測試技術(shù)。北京：機械工業(yè)出版社,2010.［3]薛年喜MATLAB在數(shù)字信號處理中的應用（第二版）北京：清華大學出版社,2008［4]焦衛(wèi)東。旋轉(zhuǎn)機械振動信號分析。浙江.嘉興學院學報.2007.附錄附錄一時域分析、頻域分析程序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（’振動信號的波形’xlabel('時間/秒')；ylabel（’幅度grid;on；％時域分析mi=min(Y);disp(mi）;%最小值mx=max（Y)；disp（mx)；％最大值st=std(Y)；disp標準差m=meandisp（m);％均值vr=var（Ydisp；％方差rm=rmsdisp(rm);％均方值10

(完整版)MATLAB在機械振動信號中的應用％頻域分析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方法進行估計程序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