時域和頻域特征提取Matlab編程實例

1、第一章 緒論1.1 概述機械信號是指機械系統(tǒng)在運行過程中各種隨時間變化的動態(tài)信息，經(jīng)各種測試儀器拾取并記錄和存儲下來的數(shù)據(jù)或圖像。機械設(shè)備是工業(yè)生產(chǎn)的基礎(chǔ)，而機械信號處理與分析技術(shù)則是工業(yè)發(fā)展的一個重要基礎(chǔ)技術(shù)。隨著各行各業(yè)的快速發(fā)展和各種各樣的應(yīng)用需求，信號分析和處理技術(shù)在信號處理速度、分辨能力、功能范圍以及特殊處理等方面將會不斷進(jìn)步，新的處理激素將會不斷涌現(xiàn)。當(dāng)前信號處理的發(fā)展主要表現(xiàn)在：1.新技術(shù)、新方法的出現(xiàn)；2.實時能力的進(jìn)一步提高；3.高分辨率頻譜分析方法的研究三方面。信號處理的發(fā)展與應(yīng)用是相輔相成的，工業(yè)方面應(yīng)用的需求是信號處理發(fā)展的動力，而信號處理的發(fā)展反過來又拓展了它的應(yīng)用領(lǐng)

2、域。機械信號的分析與處理方法從早期模擬系統(tǒng)向著數(shù)字化方向發(fā)展。在幾乎所有的機械工程領(lǐng)域中，它一直是一個重要的研究課題。機械信號分析與處理技術(shù)正在不斷發(fā)展，它已有可能幫助從事故障診斷和監(jiān)測的專業(yè)技術(shù)人員從機器運行記錄中提取和歸納機器運行的基本規(guī)律，并且充分利用當(dāng)前的運行狀態(tài)和對未來條件的了解與研究，綜合分析和處理各種干擾因素可能造成的影響，預(yù)測機器在未來運行期間的狀態(tài)和動態(tài)特性，為發(fā)展預(yù)知維修制度、延長大修期及科學(xué)地制定設(shè)備的更新和維護計劃提供依據(jù)，從而更為有效地保證機器的穩(wěn)定可靠運行，提高大型關(guān)鍵設(shè)備的利用率和效率。機械信號處理是通過對測量信號進(jìn)行某種加工變換，削弱機械信號中的無用的冗余信號，

3、濾除混雜的噪聲干擾，或者將信號變成便于識別的形式以便提取它的特征值等。機械信號處理的基本流程圖如圖1.1所示。 測量信號信號處理 特征提取獲得設(shè)備結(jié)構(gòu)特征或運行狀態(tài)機械設(shè)備圖1.1 機械信號處理的基本流程 本文主要就第三、第四步驟展開討論。第2章 機械信號的時域處理及其分析方法2.1 時域統(tǒng)計特征參數(shù)處理通過時域波形可以得到的一些特征參數(shù)，它們常用于對機械進(jìn)行快速評價和簡易診斷。 有量綱的幅值參數(shù)有量綱的幅值參數(shù)包括方根幅值、平均幅值、均方幅值和峰值等。若隨機過程x(t)符合平穩(wěn)、各態(tài)歷經(jīng)條件且均值為零，設(shè)x為幅值，p(x)為概率密度函數(shù)，有量綱型幅值參數(shù)可定義為xd=式中：xr為方根均值，為

4、均值，為均方值，為峰值。由于有量綱型幅值參數(shù)來描述機械狀態(tài)，不但與及其的狀態(tài)有關(guān)，而且與機器的運動參數(shù)（如轉(zhuǎn)速、載荷等）有關(guān)，因此直接用它們評價不同工況的機械無法得出統(tǒng)一的結(jié)論。 無量綱型參數(shù)無量綱型參數(shù)具有對機械工況變化不敏感的特點，這就意味這，理論上它們與機械的運動條件無關(guān)，它們只依賴于概率密度函數(shù)p(x)的形狀，所以無量綱型參數(shù)是一種較好的評價參數(shù)。一般它可定義為，由此公式，可得到如下的一些指標(biāo)波形指標(biāo)l=2,m=1 K=峰值指標(biāo)l,m=2 C=脈沖指標(biāo)l,m=1 I=裕度指標(biāo)l,m=1/2 L=峭度指標(biāo) K=式中為信號標(biāo)準(zhǔn)差=2.2 相關(guān)分析方法以及應(yīng)用所謂相關(guān)，就是指變量之間的線性關(guān)

5、系，它是一個非常重要的概念。對于確定性信號，兩個變量之間可以用函數(shù)關(guān)系來描述，兩者一一對應(yīng)并為確定的數(shù)值。而兩個隨即變量之間不具有確定的關(guān)系。但是，如果這兩個變量之間存在著某種不確定但卻有著表征其特性的近似關(guān)系，這兩個變量之間會有一定的線性關(guān)系。這時，對于一個隨機機械信號，可以采用相關(guān)性函數(shù)來描述其在不同時間的幅值變化相關(guān)程度。 自相關(guān)函數(shù)的概念和性質(zhì) x(t)是各態(tài)歷經(jīng)隨機過程的一個樣本函數(shù)，x(t+t)是x(t)時移t后的樣本(圖2.6)，把相關(guān)系數(shù)rx(t)x(t+t)簡寫為rx(t)，那么就有:   圖2.6 波形圖若用Rx(t)表示自相關(guān)函數(shù)，其定義為: 信號的性質(zhì)不同，自

6、相關(guān)函數(shù)有不同的表達(dá)形式。如對周期信號（功率信號）：非周期信號（能量信號）： 圖2.7給出了自相關(guān)函數(shù)具有的性質(zhì)。正弦函數(shù)的自相關(guān)函數(shù)是一個余弦函數(shù)，在=0時具有最大值。它保留了幅值信息和頻率信息,但丟失了原正弦函數(shù)中的初始相位信息。2.3 Matlab編程實驗結(jié)果 構(gòu)造加噪周期信號，時域特征分析，自相關(guān)函數(shù)特性的驗證，（程序1）圖2.8 噪聲-自相關(guān).jpg 如圖所示：自相關(guān)函數(shù)消除了大量的噪聲，周期成分變得非常明顯。原始信號的時域處理結(jié)果：平均值：0.0184極小值：-2.8138極大值：2.8557標(biāo)準(zhǔn)差：1.0103 方差： 1.0207峰峰值：5.6695第3章 機械信號的頻域處理方

7、法及其應(yīng)用信號處理中，傅立葉變換把一個隨機信號解析成不同頻率的正弦波，使信號的頻域分析稱為可能。由于計算機技術(shù)的發(fā)展，在微機上直接使用離散傅立葉變換變得非常方便，這使得頻域分析稱為常用的處理方法。常用的頻域分析方法包括自譜、功率譜、倒譜等。3.1 頻譜的分析方法DFT和FFT 離散傅立葉變換DFT傅立葉變換及其逆變換都不適合用數(shù)字計算機計算。要進(jìn)行數(shù)字計算和處理，必須將連續(xù)信號離散化，無限數(shù)據(jù)有限化。這種對有限個離散數(shù)據(jù)的傅立葉變換，稱為有限離散傅立葉變換，簡稱DFT（Discrete Fourier Trasform）。3.1.2 快速傅立葉變換FFT1965年J.W.Cooley 和J.W

8、.Tukey研究一種DFT的快速算法，稱為快速傅立葉變換，簡稱FFT(FastFourier Transform)。FFT的迅速發(fā)展，使數(shù)字頻譜分析取得了突破性的進(jìn)展。根據(jù)FFT快速變換的指導(dǎo)思想，就可以編制FFT的計算程序。時間序列從時域到頻域要用FFT變換，從頻域到時域要用逆變換IFFT，F(xiàn)FT和IFFT的公式可以統(tǒng)一。 功率譜密度函數(shù)的物理意義Sx(f)和Sxy(f)是隨機信號的頻域描述函數(shù)。Sx(f)表示信號的功率密度沿頻率軸的分布，故又稱Sx(f)為功率譜密度函數(shù)。3.2 功率譜方法以及應(yīng)用功率譜的定義式為若X（）=DFTx(m)，x(n)為N點序列。則X（） =DFTx(-m)從而

9、有 DFTR(M)= DFTx(m) DFTx(-m)即 （）= X（）X（）=| X（）|2綜上所述，先用FFT求出隨機離散序列的DFT，再計算幅頻特性的平方，再除以N，即得到該隨機信號的功率譜估計。3.3 倒頻譜分析方法倒頻譜實際上是頻域信號取對數(shù)的傅立葉變換再處理，或稱為“頻域信號的傅立葉變換再變換”。對功率譜密度函數(shù)取對數(shù)的目的是使再變換以后信號的能量更加集中。倒頻譜可以分析復(fù)雜頻譜上的周期成分，分離和提取在密集泛頻信號中的成分。對于具有同族諧頻和異族諧頻等復(fù)雜信號的分析，效果很好。倒頻譜用于對語音分析中的語言音調(diào)的測定和檢測、機械振動譜圖中的諧波分量作故障檢測和診斷以及排除回波等方面

10、是很有效的。 倒頻譜的數(shù)學(xué)描述倒頻譜函數(shù)CF(q)(power cepstrum)其數(shù)學(xué)表達(dá)式為： CF(q)又叫功率倒頻譜，或叫對數(shù)功率譜的功率譜。工程上常用的是式(2.67)的開方形式，即： C0(q)稱為幅值倒頻譜，有時簡稱倒頻譜。倒頻譜自變量q的物理意義為了使其定義更加明確，還可以定義：即倒頻譜定義為信號的雙邊功率譜對數(shù)加權(quán)，再取其傅里葉逆變換，聯(lián)系一下信號的自相關(guān)函數(shù)：看出，這種定義方法與自相關(guān)函數(shù)很相近，變量q與在量綱上完全相同。 為了反映出相位信息，分離后能恢復(fù)原信號，又提出一種復(fù)倒頻譜的運算方法。若信號x(t)的傅里葉變換為X(f): x(t)的倒頻譜記為：

11、顯而易見，它保留了相位的信息。 倒頻譜與相關(guān)函數(shù)不同的只差對數(shù)加權(quán)，目的是使再變換以后的信號能量集中，擴大動態(tài)分析的頻譜范圍和提高再變換的精度。還可以解卷積(褶積)成分，易于對原信號的分離和識別。3.4 細(xì)化譜分析方法細(xì)化譜分析法是增加頻譜中某些部分分辨能力的方法，即“局部放大”的方法。所謂細(xì)化分析室只對固定某窄帶部分進(jìn)行放大，像照相機將照片的個別部分放大一樣，使其動態(tài)范圍和分辨率都提高。細(xì)化的分析過程中，首先像通常的FFT做法那樣，選用采樣頻率fs=1/h進(jìn)行采樣，可得到N點離散序列xn.假設(shè)我們感興趣的譜中心頻率為fk的一個窄帶f,然后用一個復(fù)正弦序列（單位旋轉(zhuǎn)矢量）exp-j2fknh乘

12、以xn的yn新的N點離散序列。根據(jù)頻移定理，即將頻率原點有效地移至頻率fk（即復(fù)調(diào)制）。fk成為新的頻率坐標(biāo)原點。正、負(fù)采樣頻率fs也同樣移動了一個量fk。低通濾波后得到gm序列所保留下來的窄頻帶，若濾波后的總帶寬小于采樣頻率的1/D倍，就有可能把采樣頻率降低到1/D，而不會再新的乃奎斯特頻率附近產(chǎn)生混疊。然后再重新采樣，用fs2= fs/D的頻率來采樣，即降低了采樣頻率。由采樣定理可知，降低采樣頻率而又保持同樣的采樣點數(shù)N時，就相當(dāng)于總的時間窗增長D倍，那么，頻率分辨率也提高了D倍。所以，對經(jīng)過重新采樣后獲得的新的離散序列rm進(jìn)行復(fù)數(shù)FFT計算，即可得到細(xì)化后的譜線，這些譜線就代表中心頻率為

13、fk的一窄帶f間的細(xì)化譜。3.5 Matlab編程實驗結(jié)果 產(chǎn)生一組由60HZ和150HZ的正弦信號和隨機噪聲組成的信號，觀察其時域波形與頻譜。（程序2） 圖3-1 原始信號的時域波形圖 圖3-2 原始信號的頻譜圖圖3-1看不出信號 的周期成分；圖3-2可以清除看到，在頻率為60HZ和150HZ處有兩個尖峰，即為信號的兩個頻率分量。 功率譜估計（周期圖法）：1.利用上圖的帶噪原始信號的傅里葉變換后結(jié)果幅值，將幅值平方，即可得功率譜的估計值（Welch法） 圖3-3 采樣點數(shù)為1024時的估計功率譜 圖3-4 采樣點數(shù)為256時的估計功率譜由圖3-3與3-4可看出：2.為提高周期圖的平滑性，將信

14、號分段估計并進(jìn)行平均來減少功率譜估計的協(xié)方差，得到平均周期圖。 圖3-5 三段平均的估計功率譜圖3-6 六段平均的估計功率譜由圖3-5與3-6看出：分段平均法提高了功率譜圖的平滑性，分段數(shù)越多，平滑效果越好，信號細(xì)節(jié)更易丟失。3.對數(shù)據(jù)分段加非矩形創(chuàng)形成修正的功率譜估計法：圖3-7 加漢寧窗的估計功率譜由于窗在其邊沿為零，這減少了分段對混疊的依賴效果。用合適的窗函數(shù)，采用分段長度一半的混疊率能極大地降低估計的協(xié)方差。 倒頻譜分析：圖3-8 實倒譜圖3-9 復(fù)倒譜 正弦信號，其第一個功率譜變換為一脈沖，經(jīng)濾波后進(jìn)入第二次功率譜變換，其輸出為幅度很低的三角波輸出，因而檢測不到其存在。細(xì)化譜分析：圖

15、3-10 原始信號FFT圖3-11 ZOOM-FFT程序清單程序1.構(gòu)造加噪周期信號，時域特征分析，自相關(guān)函數(shù)特性的驗證fs=1000;t=0:1/fs:(1-1/fs);maxlag=100;x=randn(1,fs);c,maxlags=xcorr(x,maxlag); %白噪聲的自相關(guān)性z=cos(2*pi*20*t)+0.7*randn(1,1000);%加白噪聲m= mean(z); disp (m); %計算平均值mi = min(z); disp (mi); %極小值mx = max(z); disp (mx); %極大值st = std(z); disp (st); %標(biāo)準(zhǔn)差f

16、c = st.2； %方差figure(1)subplot(2,2,1) %2*2第一張圖plot(t,x) %圖片區(qū)域大小xlabel('t');ylabel('x(t)');title('白噪聲'); %加標(biāo)題subplot(2,2,2)plot(maxlags/fs,c)xlabel('t');ylabel('r(t)');title('白噪聲自相關(guān)');c,lags=xcorr(z,maxlag); %帶白噪聲的余弦信號自相關(guān)subplot(2,2,3)plot(t,z)xlabel(

17、9;t');ylabel('z(t)');title('原始信號');subplot(2,2,4)plot(maxlags/fs,c)xlabel('t');ylabel('rz(t)'); title('自相關(guān)'); 程序2: 產(chǎn)生一組由60HZ和150HZ的正弦信號和隨機噪聲組成的信號，對其做頻譜分析、倒譜分析以及幾種種功率譜估計方法的比較。%1.(頻譜分析)產(chǎn)生一組由60HZ和150HZ的正弦信號和隨機噪聲組成的信號，觀察其時域波形與頻譜。fs=1000;N=1024;t=(0:N-1)/fs;f1=

18、60;f2=150;s1=sin(2*pi*f1*t)+sin(2*pi*f2*t);s2=2*randn(size(t);x=s1+s2;figure(1)subplot(2,1,1)plot(t,x)X=abs(fft(x);f=(0:N/2-1)*fs/N;subplot(2,1,2)plot(f(1:N/2),X(1:N/2)%2.倒譜分析D=rceps(x); %實倒譜figure(2)subplot(2,1,1)plot(t,D) E=cceps(x); %復(fù)倒譜subplot(2,1,2)plot(t,E)%3.1功率譜估計（Welch法）Pxx=abs(fft(x,N).2/N

19、; %采樣點數(shù)為1024Pxx_short=abs(fft(x,256).2/256; %采樣點數(shù)為256figure(3)subplot(2,1,1)plot(0:N-1)/N*fs,10*log10(Pxx)subplot(2,1,2)plot(0:255)/256*fs,10*log10(Pxx_short)*10)%3.2將信號分段估計并進(jìn)行平均來減少功率譜估計的協(xié)方差，得到平均周期圖。Pxx=(abs(fft(x(1:256).2+abs(fft(x(257:512).2+abs(fft(x(513:768).2)/256/3;figure(4)plot(0:255)/256*fs,

20、10*log10(Pxx) %3.3將信號分為六段作功率譜估計再平均。Pxx=(abs(fft(x(1:128).2+abs(fft(x(129:256).2+abs(fft(x(257:384).2+abs(fft(x(385:512).2+abs(fft(x(513:640).2+abs(fft(x(641:768).2)/256/6;figure(4)plot(0:127)/128*fs,10*log10(Pxx)%3.4對數(shù)據(jù)分段家非矩形創(chuàng)形成修正的周期突法。窗在其邊沿為零，這減少了分段對混疊的依賴效果。用合適的窗函數(shù)（如海明窗，漢寧窗），采用分段長度一半的混疊率能%極大地降低估計的協(xié)

21、方差。漢寧法：w=hanning(256)'Pxx=(abs(fft(w.*x(1:256).2+abs(fft(w.*x(129:384).2+abs(fft(w.*x(257:512).2+abs(fft(w.*x(385:640).2+abs(fft(w.*x(513:768).2+abs(fft(w.*x(641:896).2)/(norm(w)2*6);figure(5)plot(0:255)/256*fs,10*log10(Pxx)程序3：%ZOOM-FFTfs=200;N=1024;n=0:N-1;t=n/fs;f=(0:N-1)*fs/N;f1=7;f2=7.2;f3=

22、8;s1=sin(2*pi*t*f1);s2=sin(2*pi*t*f2);s3=sin(2*pi*t*f3);x=s1+s2+s3;load zoomfftdata;fi=6;%最小細(xì)化截止頻率np=10;%放大倍數(shù)nfft=512;%fft長度nt=length(x);fa=fi+0.5*fs/np;%最大細(xì)化截止頻率nf=2nextpow2(nt);%？na=round(0.5*nf/np+1);%頻移n=0:nt-1;b=n*pi*(fi+fa)/fs;%確定旋轉(zhuǎn)因子y=x.*exp(-i*b);b=fft(y,nf);%fft變換a(1:na)=b(1:na);%正頻率帶通內(nèi)的元素賦

23、值a(nf-na+1:nf)=b(nf-na+1:nf);b=ifft(a,nf);%負(fù)頻率帶通內(nèi)的元素賦值c=b(1:np:nt);%重采樣a=fft(c,nfft)*2/nfft;%進(jìn)行ZOOM-FFT(nfft啥玩意兒？)%變換結(jié)果重新排序:y2=zeros(1,nfft/2);y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);y2(nfft/4+1:nfft/2)=a(1:nfft/4);n=0:(nfft/2-1);f2=fi+n*2*(fa-fi)/nfft;%FFT變換y1=fft(x,nfft)*2/nfft;f1=n*fs/nfft;ni=round(fi*nfft/fs+1);na=round(fa*nfft/fs+1);%輸出波形subplot(2,1,1)plot(t,x);subplot(2,1,2)nn=ni:na;plot(f1(nn),