www.pudn.com > HHT_power-system_power-quality_disturbances-detect > IMFeg.asv, change:2010-06-13,size:2299b

 
clc 
close all 
clear all 
fs=1000; 
t=1/fs:1/fs:1; 
N=1000; 
 
y1=220*2^0.5*cos(100*pi*t); 
 
 
y2=80*2^0.5*cos(500*pi*t); 
% % y2=50; 
q=heaviside(t-0.1)-heaviside(t-0.65); 
y3=y2.*q; 
y4=y1+y3; 
y5=-y4; 
 
% % y6=y1+y3; 
% y4=-110*2^0.5; 
 
p=heaviside(t-0.2)-heaviside(t-0.55); 
y6=y5.*p; 
 
z=y4+y6; 
% z=y1+y3+y5; 
 
 
 
 
figure; 
subplot(4,1,1); 
plot(t,y1);ylabel('幅值(v)');xlabel('时间(s)');title('工频标准信号（MATLAB模拟给出）'); 
subplot(4,1,2); 
plot(t,y3);ylabel('幅值(v)');xlabel('时间(s)');title('加入谐波信号（MATLAB模拟给出）'); 
% plot(t,y3);ylabel('幅值(v)');xlabel('时间(ms)');title('加入1（MATLAB模拟给出）'); 
 
subplot(4,1,3); 
plot(t,y6);ylabel('幅值(v)');xlabel('时间(s)');title('加入中断信号（MATLAB模拟给出）'); 
 
subplot(4,1,4); 
plot(t,z);ylabel('幅值(v)');xlabel('时间(s)');title('目标信号（MATLAB模拟给出）'); 
 
 
 
imf=emd(z); 
figure; 
 
 
  subplot(5,1,1) 
     plot(z); 
     ylabel('x(t)'); 
     set(gca,'box','off');title('对目标信号进行EMD分解'); 
 
  subplot(5,1,2); 
     plot(imf(1,:)); 
     ylabel('IMF1'); 
     set(gca,'box','off') 
  subplot(5,1,3); 
     plot(imf(2,:)); 
     ylabel('IMF2'); 
     set(gca,'box','off') 
  subplot(5,1,4); 
     plot(imf(3,:)); 
     ylabel('IMF3'); 
     set(gca,'box','off') 
  subplot(5,1,5); 
     plot(imf(4,:)); 
     ylabel('R'); 
     set(gca,'box','off');xlabel('时间(ms)') 
 
 
  
  
  
  
 [A,fa,tt]=hhspectrum(imf(1:end-1,:));          % 求每一层imf的HHT谱 
figure; 
plot(fa(1,:)*fs);ylabel('频率(Hz)');xlabel('时间(ms)');title('时间-频率特性曲线'); 
 
% figure; 
% plot(fa(2,:)*fs); 
% figure; 
% plot(fa(3,:)*fs); 
 
 
 
 
 
figure; 
plot(A(1,:)');ylabel('幅值(v)');xlabel('时间(ms)');title('时间-幅值特性曲线'); 
 
 
 
% % 求边际谱用的是矩形积分公式，输入的应是调用了toimage后的结果，而不是调用了hhspectrum后的结果 
% [E,tt1]=toimage(A,fa,tt,length(tt)); 
% E=flipud(E); 
% for k=1:size(E,1)% 
% [E,tt1]=toimage(A,fa,tt,length(tt)); 
%  
% for k=1:size(E,1) 
%     bjp(k)=sum(E(k,:))*1/fs;   
% end 
% f=(0:N-3)/N*(fs/2); 
%  figure 
% plot(f,bjp); 
% xlabel('频率 / Hz'); 
% ylabel('幅值(V)'); 
%  
%     bjp(k)=sum(E(k,:))*1/fs;   
% end 
% f=(0:N-3)/N*(fs/2); 
%  figure 
% plot(f,bjp); 
% xlabel('频率 / Hz'); 
% ylabel('幅值(V)'); 
 
% [A,fa,tt]=hhspectrum(imf(1,:));