Hilbert边际谱程序（matlab）

form

原帖由 破凰 于 2007-4-27 21:51 发表

                               
登录/注册后可看大图

form求边际谱时所用程序是没有问题的，用得是矩形积分公式。
他所的结果不正确的原因是：输入的应是调用了toimage后的结果，而不是调用了hhspectrum后的结果。
下面给一段程序，大家可以去试下。边际谱的分析 ...

谢谢破凰，总算有人给以肯定的答案了，我有些疑问：
1）按照上面程序运行，出来的谱线反而在右侧，去掉E=flipud(E);才对应了那副图，是否因为E的频率是反的，但size（E，1）算得是大小，因此这句for k=1:size(E,1)，变量k是从小到大顺过来了。所以不用翻转了。
2）IMF合集是边际谱正确与否的关键。生成的仿真信号，自然可以直接把各个成分当作IMF集直接作边际谱，结果是很好的，但经过emd分解后难免会有虚假的IMF分量，此时对IMF合集再作边际谱肯定是有出入的；如果是实际数据，经emd分解后有许多IMF分量，如何界定对那部分IMF作边际谱？要去掉虚假成分？而且边际谱的幅值有些问题。

liming982

请问xiaoshi509，能量谱你调试了吗。有程序吗？分享下啊

junzifei

（flipud是一个使矩阵上下翻转的函数。在Grilling 提供的程序toimage中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用flipud将矩阵翻转后，更便于我们阅读时频图。对于边际谱来说，如果不对E翻转，边际谱图中的频率将是从从右往左递增的。）

好像不用flipud才对，
我这边是那样的

版主用的是新的程序吗？

[ 本帖最后由 zhangnan3509 于 2007-6-11 11:13 编辑 ]

zhangnan3509

你指的是什么新程序？

zhangnan3509

看来我的帖子是白发了！呵呵！你都认为了fft作出的结果是正确的，别的做法就一定要和它的结果一样，那就没办法了！边际谱和它不一样就不可以吗？这两个的理论都不一样，结果能一样吗？

zhangnan3509

在nspab里面有解释，上面说的很清楚，而且在小议里面有做法，用法也很明确。

zhangnan3509

您说的这个问题form已经解释过了，不用反转，但是我个人意见还是要反转的！看一下省略E=flipud(E)的效果图;

[ 本帖最后由 zhangnan3509 于 2007-6-11 11:24 编辑 ]

破凰

新版的确实不用翻转了！

zhangnan3509

你说的是3月份的吧

破凰

是啊

zhlong

原帖由 破凰 于 2007-6-12 12:52 发表

                               
登录/注册后可看大图

是啊

2007年3月的确实不用再翻转了。还有一个问题请教，你给的下面程序中求边际谱的频率坐标系列时是不是应该用fs/2而不是fs？

form求边际谱时所用程序是没有问题的，用的是矩形积分公式。
他所的结果不正确的原因是：输入的应是调用了toimage后的结果，而不是调用了hhspectrum后的结果。
下面给一段程序，大家可以去试下。边际谱的分析结果是完全正确的。
clear;
fs=1000;    %fs为采样频率;
N=1000;     %采样点数
t=1/fs:1/fs:1;
y1=2*sin(60*pi*t);
y2=5*sin(90*pi*t);
y=[y1;y2;zeros(size(y1))]; %IMF集
%%%%%%%%%%%%%求边际谱
[A,fa,tt]=hhspectrum(y);
[E,tt1]=toimage(A,fa,tt,length(tt));
E=flipud(E);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;  
end
f=(1:N-2)/N*fs;
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');

破凰

是的，求频率坐标系列时应该用fs/2而不是fs。
f=(0:N-3)/N*fs/2;
当不是对信号进行整周期采样时，效果就没那么好了，尤其是幅值很不准！

zhlong

原帖由 破凰 于 2007-6-12 16:59 发表

                               
登录/注册后可看大图

是的，求频率坐标系列时应该用fs/2而不是fs。
f=(0:N-3)/N*fs/2;
当不是对信号进行整周期采样时，效果就没那么好了，尤其是幅值很不准！

谢谢你的回复！
上面帖子的3是2的笔误吧？
还有这个幅值不稳定问题我也发现了。另外我在你发程序的帖子中将fs改为fs/2，要不还会有人问这个得出的频率是理论频率的两倍的问题。

[ 本帖最后由 zhlong 于 2007-8-10 15:44 编辑 ]

破凰

是N-3，用hhspectrum求瞬时频率和瞬时幅值后，就少了两个数据，所以应是f=(0:N-3)/N*fs/2

wy558558558

本人有个问题用ＨＨＴ是否比小波分析好很多呢．他们的优缺点是什么？谢谢．请大虾们指点．:lol ．因为本文以前是用小波的，如果ＨＨＴ好的话就换个算法．