HHT分析实例；请高人指导

zhangshun5233

对一简单的信号做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);
[A,fa,tt]=hhspectrum(imf);
[E,tt1,Cenf]=toimage(A,fa);
figure;
myimage(E,'2D');   % 画出二维希尔伯特谱
figure;
myimage(E,'3D');   %三维


                                                           图1 原始信号

                                                          图2 经验模态分解

                                                      图3 二维HHT频谱

xiyangcc98

应该可以分开吧，IMF余量能量不大

zhangshun5233

可以分出来，imf1的频率为35Hz，imf2的频率为10Hz；但是在第三个图中，红色线条的纵坐标应该有很明显的两个：一个是35，一个是10。怎么做呢？

zhangnan3509

频率太接近了，换一个差距大一点的会好一些

zhangshun5233

回复 4 # zhangnan3509 的帖子

假如是s=5*sin(2*pi*50*t)+5*sin(2*pi*200*t);的话，出来的怎么一个是对应20，一个对应80呢？

summerxt404

回复 1 # zhangshun5233 的帖子

建议楼主把myimage的代码发出来看看，这样才能更好的解决问题，呵呵！~

zhangjian123q1

需要myimage的函数，楼主发我一个，zhangjianjohnny@gmail.com

zhangshun5233

“myimage”，我曾经在这个论坛里面找的，顺便谢谢那个楼主哦，忘记是谁了哈~

1. function h=myimage(dummyCoefs,plotmode,cmap)
   
2. % MYIMAGE 3-D visualization of coefficients.
   
3. % DUMMYCOEFS is the coefficients matrix to be visualized.
   
4. % Coefficients are colored using PLOTMODE.
   
5. % PLOTMODE = 'lvl' (By scale) or
   
6. % PLOTMODE = 'glb' (All scales) or
   
7. % PLOTMODE = 'abslvl' or 'lvlabs' (Absolute value and By scale) or
   
8. % PLOTMODE = 'absglb' or 'glbabs' (Absolute value and All scales)
   
9. %
   
10. % MYIMAGE(...,'plot') is equivalent to MYIMAGE(...,'absglb')
    
11. %
    
12. % You get 3-D plots (surfaces) using the same keywords listed
    
13. % above for the PLOTMODE parameter, preceded by '3D'. For example:
    
14. % MYIMAGE(...,'3Dplot') or MYIMAGE(...,'3Dlvl').
    
15. % H is the figure handle.
    
16. %----------------------------------------------------------------------
    
17. if nargin==1
    
18. plotmode='plot';
    
19. cmap=colormap(jet(240));
    
20. end
    
21. if nargin==2
    
22. cmap=colormap(jet(240));
    
23. end
    
24. NBC = 240;
    
25. if strmatch('3D',plotmode)
    
26. dim_plot = '3D';
    
27. else
    
28. dim_plot = '2D';
    
29. end
    
30. switch plotmode
    
31. case {'lvl','3Dlvl'}
    
32. lev_mode = 'row';
    
33. abs_mode = 0;
    
34. case {'glb','3Dglb'}
    
35. lev_mode = 'mat';
    
36. abs_mode = 0;
    
37. case {'abslvl','lvlabs','3Dabslvl','3Dlvlabs'}
    
38. lev_mode = 'row';
    
39. abs_mode = 1;
    
40. case {'absglb','glbabs','plot','2D','3Dabsglb','3Dglbabs','3Dplot','3D'}
    
41. lev_mode = 'mat';
    
42. abs_mode = 1;
    
43. otherwise
    
44. plotmode = 'absglb';
    
45. lev_mode = 'mat';
    
46. abs_mode = 1;
    
47. dim_plot = '2D';
    
48. end
    
49. if abs_mode , dummyCoefs = abs(dummyCoefs); end
    
50. plotPARAMS = {NBC,lev_mode,abs_mode,cmap};
    
51. switch dim_plot
    
52. case '2D'
    
53. axeAct = gca;
    
54. plotCOEFS(axeAct,dummyCoefs,plotPARAMS);
    
55. h=axeAct;
    
56. case '3D'
    
57. axeAct = gca;
    
58. surfCOEFS(axeAct,dummyCoefs,plotPARAMS);
    
59. h=axeAct;
    
60. end
    
61. %----------------------------------------------------------------------
    
62. function plotCOEFS(axeAct,coefs,plotPARAMS)
    
63. [NBC,lev_mode,abs_mode,cmap] = deal(plotPARAMS{:});
    
64. coefs = wcodemat(coefs,NBC,lev_mode,abs_mode);
    
65. img = image(coefs);
    
66. set(axeAct,'YDir','normal','Box','On');
    
67. title('Matrix''s 2-D Visualization','Parent',axeAct);
    
68. xlabel('time','Parent',axeAct);
    
69. ylabel('frequency','Parent',axeAct);
    
70. colormap(cmap);
    
71. %----------------------------------------------------------------------
    
72. function surfCOEFS(axeAct,coefs,plotPARAMS)
    
73. [NBC,lev_mode,abs_mode,cmap] = deal(plotPARAMS{:});
    
74. img = surf(coefs);
    
75. set(axeAct,'YDir','normal','Box','On');
    
76. title('Matrix''s 3-D Visualization','Parent',axeAct);
    
77. xlabel('time','Parent',axeAct);
    
78. ylabel('frequency','Parent',axeAct);
    
79. zlabel('amplitude','Parent',axeAct);
    
80. xl = [1 size(coefs,2)];
    
81. yl = [1 size(coefs,1)];
    
82. zl = [min(min(coefs)) max(max(coefs))];
    
83. set(axeAct,'Xlim',xl,'Ylim',yl,'Zlim',zl,'view',[-30 35]);
    
84. colormap(cmap);
    
85. shading('interp')

复制代码

zhangshun5233

回复 7 # zhangjian123q1 的帖子

已经贴出来了~

ghliu1986

你的IMF 就没有筛选好 你在筛选一下 就可以分开 用归一化IMF和原信号的相关系数 来筛选IMf

ghliu1986

C:\Users\Administrator\Desktop

ghliu1986

呵呵 新手 不会传图片 不是很清楚 楼主这就是筛选后的结果

天涯青驹

jiaguangfei

我看看代码，找找原因

safin0524

楼主加上这段代码试试看

for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/fs;
end
figure;
plot(Cenf(1,:)*fs,bjp);
xlabel('频率 / Hz');
ylabel('幅值');