声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 55413|回复: 189

[HHT] 完整的EMD分解全过程,有Hilbert谱和边际谱

  [复制链接]
发表于 2010-4-14 11:06 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
%示例程序
N=1000;
n=1:N;
fs=1000;
t=n/fs;
fx=10;
fy=50;
x=cos(2*pi*fx*t);
y=10*cos(2*pi*fy*t);
z=x+y;
data=z;
imf=emd(data);                        %对输入信号进行EMD分解   
[A,f,t]=hhspectrum(imf);            %对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号
[E,t,Cenf]=toimage(A,f);            %将每个IMF信号合成求取Hilbert谱,E:对应的振幅值,Cenf:每个网格对应的中心频率  这里横轴为时间,纵轴为频率        
                                                   %即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系:时间,中心频率,振幅)         
cemd_visu(data,1:length(data),imf);   %显示每个IMF分量及残余信号--------------------------------------------
disp_hhs(E);                          %希尔伯特谱----------------------------------------------------------
%画出边际谱
%N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure(3);
plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
xlabel('频率 / Hz');
ylabel('幅值');

分解后的IMF

分解后的IMF

希尔伯特谱

希尔伯特谱

边际谱

边际谱

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2010-4-15 18:50 | 显示全部楼层

回复 地板 dlingling 的帖子

我上面就是源程序,你直接新建一个空白文件,然后复制进去就可以了
发表于 2010-4-19 14:31 | 显示全部楼层
按照楼主的程序得到楼主一样的结果。不过在Hilbert谱中得到的幅值(10HZ\50HZ)却与信号原始的幅值对不上,10HZ的差不多是10而50hz差不多是1。边际谱很好地对应上了。如何解释。
发表于 2010-4-19 16:35 | 显示全部楼层
我安装的工具箱用来运行你这个程序会出现这样的错误啊

??? Undefined function or method 'instfreq' for input arguments of type 'double'.

Error in ==> hhspectrum at 79
  f(i,:)=instfreq(an(i,:)',tt,l)';
发表于 2010-4-19 17:12 | 显示全部楼层
??? Error using ==> instfreq at 66
T can not be equal to 1 neither to the last element of X

Error in ==> hhspectrum at 45
  f(i,:)=instfreq(an(i,:)',tt,l)';

出现如此错误信息!!
请解答
 楼主| 发表于 2010-4-19 20:16 | 显示全部楼层

回复 12楼 pro123 的帖子

请贴出你的源程序
 楼主| 发表于 2010-4-19 20:20 | 显示全部楼层

回复 10楼 franciscowu 的帖子

你好,Hilbert谱主要是观察频率随时间的变化规律,以及能量的相对变化,而且这个HILBERT谱的幅值是对颜色进行归一化后的值,颜色的相对变化来看出它的相对高低变化的,所以你说的值并不是它的原始幅值,不知道这样和你解释你能懂吗?
发表于 2010-4-19 21:54 | 显示全部楼层

回复 楼主 fancy78066994 的帖子

楼主,这个程序在这些试验中,还是有不错的结果的,但是在实际的应用中,确总是有各种个样的问题。特别想问下,如果去求每一个imf分量的频谱呢?比如是经过分解后,并通过[A,f,t]=hhspectrum(imf);   求得各个分量的频谱。如果想看看第一分量的频谱图,是否用      plot(f(1,:));就可以看出?
 楼主| 发表于 2010-4-19 22:33 | 显示全部楼层

回复 15楼 aprilcat 的帖子

呵呵,你这个问题问得相当好,我在我专业上就用了这个思想的,用HHT变换来分解某些信号,而只需要分析单个IMF所表示的那个频段的信号,就需要分析单个分量的频谱,有两种途径,1)将单个IMF可以进行傅里叶变换得到其频谱,2)利用边际谱来分析它的频谱;plot(f(1,:))只是这个IMF的瞬时频率,即时频图,并不是频谱图,往往我们用HHT时,未分析IMF的阶数与原始信号的关系,结合你自己的专业,来找出它们之间的对应关系,是解决某些问题很好的一个途径,有些时候,小波变换和HHT有很多相似的地方~!你大可以看看小波的很多文章~!!

评分

1

查看全部评分

发表于 2010-4-20 09:44 | 显示全部楼层

回复 14楼 fancy78066994 的帖子

楼主的意思我能够明白。就是坐标单位的问题。另,如果我要求幅值是非对数显示的幅值,频率要求是非归一化得频率,这个怎么处理,请求楼主解答
 楼主| 发表于 2010-4-20 09:51 | 显示全部楼层

回复 17楼 franciscowu 的帖子

关于显示的问题,我没有怎么深入看过,你确实有兴趣的话,可以去查看toimage函数,这只是实现的问题!
发表于 2010-4-20 11:35 | 显示全部楼层

回复 16楼 fancy78066994 的帖子

根据你的意思,我对单一的imf做了傅立叶变换,但是得到的确实一堆很乱的线段,这一般是什么问题呢?比如我要求第四层的频谱,总共是2500个点f4=imf(4,:);f4=fft(f4,2500);
2.能否详细的说下imf阶数与原始信号的关系呢?运用什么的标准呢?
发表于 2010-4-21 22:05 | 显示全部楼层
有关论文谈到运用相关系数筛选的办法,楼上可以看看。
 楼主| 发表于 2010-4-22 09:07 | 显示全部楼层

回复 19楼 aprilcat 的帖子

这要看你是分析什么问题了,你分析的信号的主要频带是在哪?信号杂乱只能说明你对要分析的问题的目标不明确,做EMD分解时往往会出现很多虚假分量,这时你可利用相关系数法去剔除这些你所不需要的分量,总之,具体的问题具体分析,你要多看文献,多思考~!
发表于 2010-4-23 10:54 | 显示全部楼层
那个EMD.m文件去那里下载啊!
回复 支持 2 反对 0

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-28 23:14 , Processed in 0.156822 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表