将式% T=0.25*2中的*2去掉即可。
因为fs=2*pi/T;这句话是为了将频率转化为角频率,从楼主给的式t=0:0.25:50中
可以看出采样频率,1/fs=0.25,即fs=1/0.25;
而与楼主用以下两句所得的采样频率不一致:
% T=0.25*2; % 采样周期 -采样间隔0.25s
% fs=2*pi/T; % 对应的采样频率(角频率)
效果见附图:从谱图上看,角频率2和4处各有一条线且比较明显。
为了方便理解,我把楼主的程序编辑修改了一下:
clc;clear all;
fs=4;
dt=1/fs;
N=200;
n=0:N-1;
t=n*dt;%以上语句与式t=0:0.25:50等效;
y=sin(2*t)+3*cos(4*t);
z=y';
imf=emd(z); % 求各IMF
emd_visu(z,1:length(z),imf) % 画各IMF时间序列 和 fine2course ,course2fine 的信号重建(3 figures)
[A,f,tt]=hhspectrum(imf); % 对每一IMF求局部振幅和局部频率
[im,tt]=toimage(A,f); % 求Hilbert-Huang谱(转换成image)
w=2*pi*fs; % 对应的采样频率(角频率)
disp_hhs(im,[],w); % 画谱图
colormap(flipud(gray))
[ 本帖最后由 wains 于 2009-5-27 20:21 编辑 ] |