www.pudn.com > HHT_power-system_power-quality_disturbances-detect > IMFeg.m, change:2010-06-18,size:1692b

```
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=100*2^0.5*cos(500*pi*t);

q=heaviside(t-0.1)-heaviside(t-0.7);
y3=y2.*q;

y4=-140*2^0.5*cos(100*pi*t);
p=heaviside(t-0.35)-heaviside(t-0.55);
y5=y4.*p;

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模拟给出）');

subplot(4,1,3);
plot(t,y5);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分解');xlabel('时间(ms)')

subplot(5,1,2);
plot(imf(1,:));
ylabel('IMF1');
set(gca,'box','off');xlabel('时间(ms)');
subplot(5,1,3);
plot(imf(2,:));
ylabel('IMF2');
set(gca,'box','off');xlabel('时间(ms)');
subplot(5,1,4);
plot(imf(3,:));
ylabel('IMF3');
set(gca,'box','off');xlabel('时间(ms)');
subplot(5,1,5);
plot(imf(4,:));
ylabel('R');
set(gca,'box','off');xlabel('时间(ms)');

[A,fa,tt]=hhspectrum(imf(1:end-1,:));
figure;
plot(fa(1,:)*fs);ylabel('频率(Hz)');xlabel('时间(ms)');title('时间-频率特性曲线');
p1=(fa(1,:)*fs)';
% figure;
% plot(fa(2,:)*fs);
% figure;
% plot(fa(3,:)*fs);

figure;
plot(A(1,:)');ylabel('幅值(v)');xlabel('时间(ms)');title('时间-幅值特性曲线');

```