時(shí)域和頻域特征提取Matlab編程實(shí)例

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

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

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

4、均值，為均方值，為峰值。由于有量綱型幅值參數(shù)來描述機(jī)械狀態(tài)，不但與及其的狀態(tài)有關(guān)，而且與機(jī)器的運(yùn)動(dòng)參數(shù)（如轉(zhuǎn)速、載荷等）有關(guān)，因此直接用它們?cè)u(píng)價(jià)不同工況的機(jī)械無法得出統(tǒng)一的結(jié)論。 無量綱型參數(shù)無量綱型參數(shù)具有對(duì)機(jī)械工況變化不敏感的特點(diǎn)，這就意味這，理論上它們與機(jī)械的運(yùn)動(dòng)條件無關(guān)，它們只依賴于概率密度函數(shù)p(x)的形狀，所以無量綱型參數(shù)是一種較好的評(píng)價(jià)參數(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=式中為信號(hào)標(biāo)準(zhǔn)差=2.2 相關(guān)分析方法以及應(yīng)用所謂相關(guān)，就是指變量之間的線性關(guān)

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

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

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

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

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

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

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

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

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

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

15、3-10 原始信號(hào)FFT圖3-11 ZOOM-FFT程序清單程序1.構(gòu)造加噪周期信號(hào)，時(shí)域特征分析，自相關(guān)函數(shù)特性的驗(yàn)證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); %計(jì)算平均值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); %帶白噪聲的余弦信號(hào)自相關(guān)subplot(2,2,3)plot(t,z)xlabel(

17、9;t');ylabel('z(t)');title('原始信號(hào)');subplot(2,2,4)plot(maxlags/fs,c)xlabel('t');ylabel('rz(t)'); title('自相關(guān)'); 程序2: 產(chǎn)生一組由60HZ和150HZ的正弦信號(hào)和隨機(jī)噪聲組成的信號(hào)，對(duì)其做頻譜分析、倒譜分析以及幾種種功率譜估計(jì)方法的比較。%1.(頻譜分析)產(chǎn)生一組由60HZ和150HZ的正弦信號(hào)和隨機(jī)噪聲組成的信號(hào)，觀察其時(shí)域波形與頻譜。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); %實(shí)倒譜figure(2)subplot(2,1,1)plot(t,D) E=cceps(x); %復(fù)倒譜subplot(2,1,2)plot(t,E)%3.1功率譜估計(jì)（Welch法）Pxx=abs(fft(x,N).2/N

19、; %采樣點(diǎn)數(shù)為1024Pxx_short=abs(fft(x,256).2/256; %采樣點(diǎn)數(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將信號(hào)分段估計(jì)并進(jìn)行平均來減少功率譜估計(jì)的協(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將信號(hào)分為六段作功率譜估計(jì)再平均。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對(duì)數(shù)據(jù)分段家非矩形創(chuàng)形成修正的周期突法。窗在其邊沿為零，這減少了分段對(duì)混疊的依賴效果。用合適的窗函數(shù)（如海明窗，漢寧窗），采用分段長(zhǎng)度一半的混疊率能%極大地降低估計(jì)的協(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長(zhǎng)度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),