对简单仿真信号进行EMD分解和后续Hilbert谱分析时的一些问题
本人这个学期的毕业设计与HHT在轨道检测中的应用有关,刚接触这个方面的东西感觉还蛮有兴趣的,本人不久前将于德介老师的《机械故障诊断的Hibert-Huang变换方法》大致看了下,对EMD、IMF和HHT等概念有了一些基本的了解,这几天天天在泡论坛,想通过Matlab将上书中的一些简单仿真信号的处理自己实现一遍,例如书中P29的仿真信号:x(t)=2*sin(30*pi*t)+4*sin(20*pi*t).*sin(2*pi*1/10*t)+sin(2*pi*5*t),t∈
我的matlab程序如下(很多参照了论坛中达人提供的程序):
T=0:0.01:1;
x=2*sin(30*pi*T)+4*sin(20*pi*T).*sin(2*pi*1/10*T)+sin(2*pi*5*T);
plot(T,x)%绘出仿真波形图,如图1
OPTIONS.DISPLAY=2;
OPTIONS.T=0:0.01:1;
imf=emd(x,OPTIONS);
emd_visu(x,0:0.01:1,imf)%绘出各个IMF波形,如图2
=hhspectrum(imf);
=toimage(A,f);
disp_hhs(im,[],100);
colormap(flipud(gray))%绘出Hilbert谱,如图3
从仿真信号的函数可以看到理论上对信号EMD分解出来应该是一个15hz的正弦信号、一个5hz的正弦信号和一个10hz的调幅信号,可是在图2中一共得到了4个IMF和一个residue,而且前3个IMF的效果很差,比书上给的EMD分解图效果差很多,第2个IMF几乎没有调幅的特征了,我不知道是我做的EMD分解过程有误还是这就是EMD的端点效应问题?要是是端点效应的问题话,不知道用Matlab该如何实现抑制端点效应?
还有我用disp_hhs(im,[],100)语句倒是把hilbert谱的纵坐标转换成了实际频率(采样频率100hz),或者可以单独把hilbert谱的横坐标转换成实际时间,但不知道怎么同时把两个坐标都转换成实际的,我想通过disp_hhs(im,tt,[],100)(tt=0:0.01:1)同时将两坐标转换过来,但会有错误,提示说输入参数数量太多了。。。搞了半天不知道该怎么解决这个问题。。。
[ 本帖最后由 Snikas 于 2008-4-12 19:50 编辑 ] 我也很奇怪这个问题,不过你可以修改disp_hhs.m这个函数的代码,将其中
if fs == 0
imagesc(t,,im,);
ylabel('normalized frequency')
else
imagesc(t,,im,);改为imagesc(t/fs,,im,);
ylabel('frequency')
end
然后通过disp_hhs(im,[],fs);调用,则可以同时在图中表现对应的实际时间和频率了
回复 4楼 的帖子
确实可行,谢谢哈:victory: 原程序确实有点问题哈,自己说明里都说了有四个输入参量,可是实际又不行。。。回复 楼主 的帖子
disp_hhs(im,tt,[],100)中的tt好像不是0:0.01:1,因为程序中对时间进行了截取,可以这样处理(t-1)/fs 为什么我在调试程序时总是出现下面的错误呢???? Undefined command/function 'instfreq'.
Error in ==> hhspectrum at 79
f(i,:)=instfreq(an(i,:)',tt,l)';
Error in ==> Untitled5 at 8
=hhspectrum(imf);
是没有instfreq函数吗?
回复 5楼 的帖子
没有instfreq这个函数 将error(nargchk(1,3,nargin));改为error(nargchk(1,4,nargin));调用时用disp_hhs(im,t/fs,[],fs),也可以解决问题的! 很感谢楼上的几位,让我学到了好多东西。知道了,调用的函数能修改。谢谢! 这个与设置的迭代次数及采样点数多少有关。
我分解的是下面的图。
对简单仿真信号进行EMD分解和后续Hilbert谱分析时的一些问题
请问楼上参数是怎么设置的? 参看相关文章,MSSP杂志上 ,有专门的文章讲述,参数的选择问题,楼上可以看一下 回复 1 # Snikas 的帖子请问楼主如何才能让每一个imf图上的横纵坐标的值显示出来呢,谢谢哈,急!!! 请问楼主如何才能让每一个imf图上的横纵坐标的值显示出来呢,谢谢哈,急!!!{:{39}:} 谢谢,有收获 回复 7 # jinnian 的帖子
确实好用,谢谢!
页:
[1]
2