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

```
clc
close all
clear all
fs=1000;
t=1/fs:1/fs:0.5;
N=1000;

y1=220*2^0.5*cos(100*pi*t);

y2=250*cos(200*pi*t)+30*cos(300*pi*t);
% y2=50;
% q=heaviside(t-0.1)-heaviside(t-0.3);
% y3=y2.*q;

% y4=200*cos(200*pi*t);
% p=heaviside(t-0.3)-heaviside(t-0.4);
% y5=y4.*p;
z=y1+y2;
% z=y1+y3+y5;

figure;
subplot(2,1,1);
plot(t,y1);ylabel('幅值(v)');xlabel('时间(ms)');title('工频标准信号（MATLAB模拟给出）');
subplot(2,1,2)
plot(t,y3);ylabel('幅值(v)');xlabel('时间(ms)');title('加入xin1（MATLAB模拟给出）');

plot(t,y5);ylabel('幅值(v)');xlabel('时间(ms)');title('加入2（MATLAB模拟给出）');

plot(t,z);ylabel('幅值(v)');xlabel('时间(ms)');title('信号（MATLAB模拟给出）');

imf=emd(z);
figure
subplot(5,1,1)
plot(z);
ylabel('x(t)');
set(gca,'box','off')
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('IMF4');
set(gca,'box','off')

[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('频率(Hz)');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,:));

```