HHT分析实例；请高人指导
对一简单的信号做HHT分析，该信号含有两个频率成分10Hz和35Hz，得到的希尔伯特谱（二维）纵轴的频率怎么不是10和35；程序和图像如下：N=2048;
t=linspace(0,1,N);
delta=t(2)-t(1); % 采样周期
fs=1/delta; % 采样频率
s=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
imf=emd(s);
=hhspectrum(imf);
=toimage(A,fa);
figure;
myimage(E,'2D'); % 画出二维希尔伯特谱
figure;
myimage(E,'3D'); %三维
图1 原始信号
图2 经验模态分解
图3 二维HHT频谱
应该可以分开吧，IMF余量能量不大 可以分出来，imf1的频率为35Hz，imf2的频率为10Hz；但是在第三个图中，红色线条的纵坐标应该有很明显的两个：一个是35，一个是10。怎么做呢？ 频率太接近了，换一个差距大一点的会好一些 回复 4 # zhangnan3509 的帖子
假如是s=5*sin(2*pi*50*t)+5*sin(2*pi*200*t);的话，出来的怎么一个是对应20，一个对应80呢？ 回复 1 # zhangshun5233 的帖子
建议楼主把myimage的代码发出来看看，这样才能更好的解决问题，呵呵！~ 需要myimage的函数，楼主发我一个，zhangjianjohnny@gmail.com 本帖最后由 zhangshun5233 于 2011-12-7 11:20 编辑
“myimage”，我曾经在这个论坛里面找的，顺便谢谢那个楼主哦，忘记是谁了哈~
function h=myimage(dummyCoefs,plotmode,cmap)
% MYIMAGE 3-D visualization of coefficients.
% DUMMYCOEFS is the coefficients matrix to be visualized.
% Coefficients are colored using PLOTMODE.
% PLOTMODE = 'lvl' (By scale) or
% PLOTMODE = 'glb' (All scales) or
% PLOTMODE = 'abslvl' or 'lvlabs' (Absolute value and By scale) or
% PLOTMODE = 'absglb' or 'glbabs' (Absolute value and All scales)
%
% MYIMAGE(...,'plot') is equivalent to MYIMAGE(...,'absglb')
%
% You get 3-D plots (surfaces) using the same keywords listed
% above for the PLOTMODE parameter, preceded by '3D'. For example:
% MYIMAGE(...,'3Dplot') or MYIMAGE(...,'3Dlvl').
% H is the figure handle.
%----------------------------------------------------------------------
if nargin==1
plotmode='plot';
cmap=colormap(jet(240));
end
if nargin==2
cmap=colormap(jet(240));
end
NBC = 240;
if strmatch('3D',plotmode)
dim_plot = '3D';
else
dim_plot = '2D';
end
switch plotmode
case {'lvl','3Dlvl'}
lev_mode = 'row';
abs_mode = 0;
case {'glb','3Dglb'}
lev_mode = 'mat';
abs_mode = 0;
case {'abslvl','lvlabs','3Dabslvl','3Dlvlabs'}
lev_mode = 'row';
abs_mode = 1;
case {'absglb','glbabs','plot','2D','3Dabsglb','3Dglbabs','3Dplot','3D'}
lev_mode = 'mat';
abs_mode = 1;
otherwise
plotmode = 'absglb';
lev_mode = 'mat';
abs_mode = 1;
dim_plot = '2D';
end
if abs_mode , dummyCoefs = abs(dummyCoefs); end
plotPARAMS = {NBC,lev_mode,abs_mode,cmap};
switch dim_plot
case '2D'
axeAct = gca;
plotCOEFS(axeAct,dummyCoefs,plotPARAMS);
h=axeAct;
case '3D'
axeAct = gca;
surfCOEFS(axeAct,dummyCoefs,plotPARAMS);
h=axeAct;
end
%----------------------------------------------------------------------
function plotCOEFS(axeAct,coefs,plotPARAMS)
= deal(plotPARAMS{:});
coefs = wcodemat(coefs,NBC,lev_mode,abs_mode);
img = image(coefs);
set(axeAct,'YDir','normal','Box','On');
title('Matrix''s 2-D Visualization','Parent',axeAct);
xlabel('time','Parent',axeAct);
ylabel('frequency','Parent',axeAct);
colormap(cmap);
%----------------------------------------------------------------------
function surfCOEFS(axeAct,coefs,plotPARAMS)
= deal(plotPARAMS{:});
img = surf(coefs);
set(axeAct,'YDir','normal','Box','On');
title('Matrix''s 3-D Visualization','Parent',axeAct);
xlabel('time','Parent',axeAct);
ylabel('frequency','Parent',axeAct);
zlabel('amplitude','Parent',axeAct);
xl = ;
yl = ;
zl = ;
set(axeAct,'Xlim',xl,'Ylim',yl,'Zlim',zl,'view',[-30 35]);
colormap(cmap);
shading('interp')
回复 7 # zhangjian123q1 的帖子
已经贴出来了~ 你的IMF 就没有筛选好 你在筛选一下 就可以分开 用归一化IMF和原信号的相关系数 来筛选IMf
C:\Users\Administrator\Desktop 呵呵 新手 不会传图片 不是很清楚 楼主这就是筛选后的结果
{:{39}:} 我看看代码，找找原因 楼主加上这段代码试试看
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/fs;
end
figure;
plot(Cenf(1,:)*fs,bjp);
xlabel('频率 / Hz');
ylabel('幅值');
页:
[1]
2