故障诊断作业5

故障诊断作业5

1. 仿照例5.2，例5.3对下面的信号函数进行EMD与小波分解对比，并说明。

x=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t)

第一个示例（仿照例5.2）：

代码如下：

clc;  
clear all;   
close all;       
  
%% 构造一个x=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t)信号  
nfft=256;  
fs=2048;  
N=1024;   
t=0:1/fs:(N-1)/fs;  
x=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t);  
imf=emd(x);  
[m,n]=size(imf);  
figure(1)  
plot(x);  
xlim([1 1024]);  
xlabel('时间（点数）','FontSize',12);  
ylabel('初始信号','FontSize',12);  
  
%% 求Hilbert谱  
figure(2)  
dt=1/fs;  
h=nspabz(imf',500,0,500,0,N/fs);  
subplot(211);  
surf(h(1:200,50:end-20));  
shading interp; %通过在每个线条或面中对颜色图索引或真彩色值进行插值来改变该线条或面中的颜色。  
xlabel('时间（点数）','FontSize',12);  
ylabel('频率/Hz','FontSize',12);  
zlabel('幅值','FontSize',12);  
title('(a) Hilbert时频图','FontSize',12);  
view([-75,25]);  
yt=subplot(223);  
% imagesc(h(1:200,:));  
contour((1:1024)/10,(1:200)/200,h(1:200,:)); %Hilbert谱  
xlabel('时间（点数）','FontSize',12);  
ylabel('频率/Hz','FontSize',12);  
set(yt,'ydir','nor')  
title('(b) Hilbert时频谱（等高线图）','FontSize',12)  
  
%% 求Hilbert边际谱  
ms=mspc(h); %Hilbert边际谱  
subplot(224);   
plot((2/500:2/500:length(ms)*2/500),ms);  
shading flat; %每个网格线段和面具有恒定颜色，该颜色由该线段的端点或该面的角边具有最小索引的颜色值确定。  
xlabel('频率/Hz','FontSize',12);  
ylabel('幅值','FontSize',12);  
title('(c) Hilbert 边际谱','FontSize',12);  
% sgtitle('仿真信号的Hilbert谱和 Hilbert时频谱');  
  
%% 小波分解方法处理  
T=wpdec(x,4,'db8');%采用DB8进行4层小波包分解  
%分解系数重构及小波包节点排序  
y(:,1)=wprcoef(T,[4,0]);  
y(:,2)=wprcoef(T,[4,1]);  
y(:,3)=wprcoef(T,[4,2]);  
y(:,4)=wprcoef(T,[4,3]);  
y(:,5)=wprcoef(T,[4,4]);  
y(:,6)=wprcoef(T,[4,5]);  
y(:,7)=wprcoef(T,[4,6]);  
y(:,8)=wprcoef(T,[4,7]);  
y(:,9)=wprcoef(T,[4,8]);  
y(:,10)=wprcoef(T,[4,9]);  
y(:,11)=wprcoef(T,[4,10]);  
y(:,12)=wprcoef(T,[4,11]);  
y(:,13)=wprcoef(T,[4,12]);  
y(:,14)=wprcoef(T,[4,13]);  
y(:,15)=wprcoef(T,[4,14]);  
y(:,16)=wprcoef(T,[4,15]);  
  
figure(3)  
dt=1/fs;  
h=nspabz(y,500,0,500,0,N/fs);  
subplot(211)  
surf(h(1:200,50:end-20));  
shading interp  
xlabel('时间（点数）','FontSize',12);  
ylabel('频率/Hz','FontSize',12);  
zlabel('幅值','FontSize',12);  
title('（a）小波时频图','FontSize',12);  
view([-75,25]);  
yt=subplot(223);  
contour((1:1024)/10,(1:200)/200,h(1:200,:));  
xlabel('时间（点数）','FontSize',12);  
ylabel('频率/Hz','FontSize',12);  
set(yt,'ydir','nor')  
title('（b）小波时频谱（等高线图）','FontSize',12)  
ms=mspc(h);  
subplot(224);  
plot((1:length(ms))/250,ms)  
shading flat;  
xlabel('频率/Hz','FontSize',12);  
ylabel('幅值','FontSize',12);  
title('（c）小波边际谱','FontSize',12)  
% sgtitle('仿真信号的小波谱');  

打印输出结果：

图1-1 仿真信号时域波形图

图1-2 仿真信号的Hilbert谱和Hilbert时频谱

图1-3 仿真信号的小波谱和小波时频谱

总结说明：

图1-1所示为x(t)的时域波形图,对其进行EMD并作出Hilbert时频谱和Hilbert 边际谱,如图1-2所示。从图1-2(a)和图1-2(b)中可以清晰分辨原始信号在不同时间和不同频段上的分布,与原始信号基本吻合;图1-2(c)所示的Hilbert边际谱中,4个谱峰对应的频率正好是原信号中4个时频分量的瞬时频率0.04 Hz、0.1 Hz、0.2 Hz和0.4 Hz。图1-2给出了仿真信号的完整时频解释。

将仿真信号进行小波包分解,得到的小波谱如图1-3所示。由图1-3(a)和图1-3(b)可以看出,信号中的在时频面上的分布分散,出现了许多并不存在的虚假频率成分,而在真实频率0.2 Hz外的聚集性变差,越到后面越差。图1-3(c)所示时频分布在时间轴上的积分,显示出原始信号频率成分的大致分布,第1个峰值时段分量的固有频率 0.15 Hz ,第2个峰值时段分量的固有频率 0.4 Hz，剩下的频率区域出现频率模糊的现象，难以辨别信号。该仿真实例表明了Hilbert谱良好的局部化特性。

Hilbert-Huang变换谱所有成分的频率分辨率是致的,时间分辨率和频率分辨率相互独立,彼此没有干扰;而小波分析的时间分辨率和频率分辨率是相互影响的。因此Hilbert-Huang变换在处理强间歇性信号时具有优势。

第二个示例（仿照例5.3）：

假设s1为正常信号，s2，s3为非正常信号

s1=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t);  

s2=cos(2*pi*10*t)+cos(2*pi*50*t)+cos(2*pi*100*t)+cos(2*pi*200*t);  

s3=sin(2*pi*35*t)+sin(2*pi*75*t)+sin(2*pi*150*t)+sin(2*pi*250*t);  

代码如下：

clc;  
clear all;   
close all;   
  
%% 正常信号输出和非正常信号的输出（s1为正常信号，s2和s3为非正常信号）  
fs=256;  
n=1200;  
t=0:1/fs:(n-1)/fs;  
s1=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t);  
s2=cos(2*pi*10*t)+cos(2*pi*50*t)+cos(2*pi*100*t)+cos(2*pi*200*t);  
s3=sin(2*pi*35*t)+sin(2*pi*75*t)+sin(2*pi*150*t)+sin(2*pi*250*t);  
  
figure(1)  
subplot(311);  
plot(t,s1);  
xlabel('t/s');title('正常信号s1');  
xlim([0 (n-1)/fs]);ylim([-5 5]);  
subplot(312);  
plot(t,s2);  
xlabel('t/s');title('非正常信号s2');  
xlim([0 (n-1)/fs]);ylim([-5 5]);  
subplot(313);  
plot(t,s3);  
xlabel('t/s');title('非正常信号s3');  
xlim([0 (n-1)/fs]);ylim([-5 5]);  
  
%% 进行EMD分解，并对前6个IMF分量进行AR谱分析  
imf1=emd(s1);  
figure(2) %正常信号的前6个IMF分量  
title('正常信号的前6个IMF分量');  
for(i=1:6)  
    subplot(6,1,i);  
    plot(t,imf1(i,:));  
    xlabel('t/s');  
    xlim([0 (n-1)/fs]);  
end  
% sgtitle('First 6 IMF components of normal signal');  
  
%%  
imf2=emd(s2);  
imf3=emd(s3);  
%对前6个IMF分量进行AR谱分析  
N=1024;  
for (i=1:6)  
    xpsd=pburg(imf1(i,:),10,N);  
    xpsd1(i,:)=10*log10(xpsd(1:400)+0.000001);  
end  
for (i=1:6)  
    xpsd=pburg(imf2(i,:),10,N);  
    xpsd2(i,:)=10*log10(xpsd(1:400)+0.000001);  
end  
for (i=1:6)  
    xpsd=pburg(imf3(i,:),10,N);  
    xpsd3(i,:)=10*log10(xpsd(1:400)+0.000001);  
end  
ss=N/2;  
dd=(0:1/ss:1)*fs/2;  
d=dd(1:400); %频率选择在10000Hz以内  
  
figure(3)  
for(i=1:6)  
    subplot(3,2,i)  
    a={'IMF1','IMF2','IMF3','IMF4','IMF5','IMF6'};  
    plot(d,xpsd1(i,:),'b-',d,xpsd2(i,:),'g--',d,xpsd3(i,:),'r-.');  
    legend('正常信号s1','非正常信号s2','非正常信号s3');  
    xlabel('频率/Hz');ylabel('EMD-AR谱/dB');  
    title(a{i});  
end  
  
%% 将前6个IMF分量的AR谱累加  
q1=sum(xpsd1);  
q2=sum(xpsd2);  
q3=sum(xpsd3);  
figure(4)  
plot(d,q1,'b-',d,q2,'g--',d,q3,'r-.');  
legend('正常信号s1','非正常信号s2','非正常信号s3');  
xlabel('频率/Hz');ylabel('EMD-AR谱、dB');  
% title('3种信号下前6个IMF分量EMD-AR谱累加能量对比图');  

打印结果如图：

图2-1 3种信号的振动号波形图

图2-2 3种信号的振动号波形图

图2-3 3种信号前6个IMF分量EMD-AR谱对比图

图2-4 3种信号前6个IMF分量EMD-AR谱累加能量对比图

总结：

(1)从图2-3中可以看出,能量由正常信号s1依次递增的特征频带有: IMF1分量0 ~ 50 Hz，90 ~ 100 Hz, IMF2分量30 ~ 50 Hz; IMF3分量0 ~ 24.75Hz; IMF4分量0 ~ 10 Hz; IMF5 分量0 ~ 4.5 Hz; IMF6 分量0 ~ 2Hz。依次递减的特征频带有: IMF1分量53 ~ 80 Hz; IMF2分量49.75 ~ 80 Hz; IMF3分量24.75 ~ 60 Hz, 76 ~ 100 Hz; IMF4分量10 ~ 70 Hz; IMF5 分量4.5 ~ 20 Hz; IMF6 分量2 ~ 8 Hz；

EMD与AR谱的结合，很容易能分辨出正常信号与非正常信号的区别，能够有效地分析信号，克服了Hilbert分离法存在加窗效应的问题；

(2)从图2-4中可以看出，前6个IMF分量EMD-AR谱累加能量对比，能够很直观提取非正常信号的特征。