直接贴出m文件吧。
上面的帖子没法修改。
clc;
clear all;
fa=10;
fb=20;
a=1;
b=1;
dt=0.002;
fs=1/dt;
t=0:dt:5;
tspan=dt*(length(t)-1);
pi=3.1415926;
x=a*sin(2*pi*fa*t)+b*sin(2*pi*fb*t);
imf=emd(x);
idx=size(imf);
han=idx(1)+1;
subplot(han,1,1)
plot(t,x)
ylabel('original')
%++++++++++++++++++++++++++++++++
for j=1:han-1
subplot(han,1,j+1)
plot(t,imf(j,:))
ylabel('imf')
end
[A,f,tt]=hhspectrum(imf); %HHT 时频谱计算
[im,tt]=toimage(A,f,tt,length(tt));
disp_hhs(im,[],fs);
ylim([0,30])
NN=size(im,1);
for k=1:NN
bjp(k)=sum(im(k,:))*1/fs;
end
fbjp=(0:NN-1)/NN*(fs/2);
%----------------------边际谱--------------------------
figure(3)
%bjp=bjp/max(bjp);
plot(fbjp,bjp)
xlim([0,30])
xlabel('频率/Hz')
ylabel('幅值')
title('边际谱')
|