|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 longyi2003 于 2013-5-8 21:17 编辑
这几个图的效果都好差啊,,各位大牛们,,能告诉我这是什么原因造成的吗?
程序如下:- clear;
- N=3200;
- t=linspace(0,1,N);
- delta=t(2)-t(1); % 采样周期
- fs=1/delta; % 采样频率
- x=50*sin(2*pi*50*t);
- y=15*sin(2*pi*150*t);
- y1=8*sin(2*pi*250*t);
- y2=5*sin(2*pi*350*t);
- y3=2*sin(2*pi*450*t);
- y4=1.8*sin(2*pi*550*t);
- z=x+y+y1+y2+y3+y4;
- % plot(t,z);
- imf=emd1(z);
- cemd_visu(z,1:length(z),imf);
- [B,f,tt]=hhspectrum(imf);
- [im,tt,Cenf]=toimage(B,f);
- disp_hhs(im,[],fs);
- colormap(flipud(gray)) %绘出Hilbert谱
- %画出边际谱
- %N=length(Cenf);%设置频率点数 %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
- for k=1:size(im,1)
- bjp(k)=sum(im(k,:))*1/fs;
- end
- figure;
- plot(Cenf(1,:)*fs,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
- xlabel('频率 / Hz');
- ylabel('幅值');
- s = size(imf);
- k = s(1);
- % figure
- %imf分量幅度
- for j=1:k
- %imf1的Hilbert变换
- xn(j,:)=hilbert(imf(j,:));
- xr(j,:)=real(xn(j,:));
- xi(j,:)=imag(xn(j,:));
- A(j,:)=sqrt(xr(j,:).^2+xi(j,:).^2);
- % subplot(k,1,j);plot(t,A(j,:));
- % plot(t,xn(j,:));
- % title(['imf',int2str(j)]);
- % hold on;
- end
- % hold off;
- %imf1的瞬时相位
- % figure,
- for j=1:k
- P(j,:)=atan2(xi(j,:),xr(j,:));
- % subplot(k,1,j);plot(t,P(j,:));
- % title('imf瞬时相位')
- end
- %imf1瞬时频率
- figure,
- for j=1:k;
- xh(j,:)=unwrap(P(j,:));
- % fs=6400;
- xhd(j,:)=fs*diff(xh(j,:))/(2*pi);
- % subplot(k,1,j);
- plot(t(1:N-1),xhd(j,:));title('imf瞬时频率')
- hold on;
- end
- hold off
复制代码 |
|