|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 汀小汀 于 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');
- [A,f,t]=hhspectrum(imf);
- [E,t,cenf]=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('边际谱')
复制代码 |
-
x(t)
-
imf
-
瞬时幅值
-
瞬时频率
-
Hilbert谱
-
边际谱
|