一个EMD的程序，大家帮我看看结果对不对？
本帖最后由 汀小汀 于 2011-5-10 12:19 编辑本人初学EMD
程序是根据论坛里面的一个程序改动的
仿真信号为x(t)=2sin(30πt)+4sin(20πt)sin(2πt/10)+sin(10πt)
请大家指教下，先谢了！
fs=1000;
l=1;
t=1/fs:1/fs:l;
N=length(t);
x=2*sin(2*pi*15*t)+4*sin(2*pi*10*t).*sin(2*pi*t/10)+sin(2*pi*5*t);
plot(t,x);
xlabel('t/s');ylabel('x(t)');
imf=emd(x);
cemd_visu(x,1:length(x),imf);
xn1=hilbert(imf(1,:));
xr1=real(xn1);
xi1=imag(xn1);
A1=sqrt(xr1.^2+xi1.^2);
figure,subplot(3,1,1);plot(t,A1);
xlabel('t/s');ylabel('瞬时振幅');title('imf1');
xn2=hilbert(imf(2,:));
xr2=real(xn2);
xi2=imag(xn2);
A2=sqrt(xr2.^2+xi2.^2);
subplot(3,1,2);plot(t,A2);
xlabel('t/s');ylabel('瞬时振幅');title('imf2');
xn3=hilbert(imf(3,:));
xr3=real(xn3);
xi3=imag(xn3);
A3=sqrt(xr3.^2+xi3.^2);
subplot(3,1,3);plot(t,A3);
xlabel('t/s');ylabel('瞬时振幅');title('imf3');
P1=atan2(xi1,xr1);
xh1=unwrap(P1);
fs=1000;
xhd1=fs*diff(xh1)/(2*pi);
figure,subplot(3,1,1);plot(t(1:999),xhd1);
xlabel('t/s');ylabel('瞬时频率');title('imf1');
P2=atan2(xi2,xr2);
xh2=unwrap(P2);
fs=1000;
xhd2=fs*diff(xh2)/(2*pi);
subplot(3,1,2);plot(t(1:999),xhd2);
xlabel('t/s');ylabel('瞬时频率');title('imf2');
P3=atan2(xi3,xr3);
xh3=unwrap(P3);
fs=1000;
xhd3=fs*diff(xh3)/(2*pi);
subplot(3,1,3);plot(t(1:999),xhd3);
xlabel('t/s');ylabel('瞬时频率');title('imf3');
=hhspectrum(imf);
=toimage(A,f);
disp_hhs(E,[],1000);
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/l*1/fs;
end
H=size(E,1);
f=(0:H-1)/H*(fs/2);
figure();
plot(f,bjp);
xlabel('频率/HZ');
ylabel('幅值');title('边际谱')
哥们，你首先把刻度值标上啊！
回复 2 # 杨德昌 的帖子
不好意思啊，马上改下图 回复 2 # 杨德昌 的帖子
页:
[1]