Hilbert边际谱程序（matlab）

drifterwang

实际的效果根本不可能每个imf分量里只有一个频率分量，并且也不可能没有干扰

dacyzhang

E=flipud(E);
是否应该去掉这个才能得到38楼的图？？

drifterwang

2007版以后的matlab要去掉

csutszy

谢谢，新手期待深入学习

土木年华

回复 simon21 的帖子
谢谢分享

   

pandacute

下面是站内有人提供的hht边际谱的程序：
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=(0:N-3)/N*(fs/2);
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');

请问：为什么f=(0:N-3)/N*(fs/2);不是f=(0:N-1)/N*(fs/2);呢，请高手赐教，谢谢。

shinhwaforever

感谢楼主的这个帖子还有楼上破凰给出的相关资料,让我知道了不少东西.

linqus

pandacute 发表于 2010-9-4 09:58

                               
登录/注册后可看大图

下面是站内有人提供的hht边际谱的程序：
clear;
fs=1000;    %fs为采样频率;

其实你将前面的N＝1000去掉，
在求边际谱前用N=size(E,1)
就要用f=(0:N-1)/N*(fs/2);了

pandacute

回复 linqus 的帖子

你好，我还不是很明白，为什么另N=1000后，要写成f=(0:N-3)/N*(fs/2);这样呢？还有bjp(k)=sum(E(k,:))*1/fs;改成bjp(k)=sum(E(k,1:end-1))*1/fs;或bjp(k)=sum(E(k,2:end))*1/fs;会不会更妥当一些，因为这是在每一个频率点处求离散情况下的面积啊！谢谢指教

WJzhengdong

谁懂啊，我也被这个问题困扰。给解释一下吧 。为什么另N=1000后，要写成f=(0:N-3)/N*(fs/2);这样呢

wangyanfang678

回复145楼
  整周期采样问题，如果采样点数不是整周期采样，会引起能量泄露问题，而绘图时要求f,bjp同样长度，所以为保证同样长度，根据bjp调整f的长度。可以看workplace中的数据长度来调整f的长度。不知道是不是这样？

hivdg

回复143楼

其实你将前面的N＝1000去掉，
在求边际谱前用N=size(E,1)
就要用f=(0:N-1)/N*(fs/2);了

f=(0:N-1)/N*(fs/2)应该改为f=(0:N-1)/（N+2)*(fs/2)才对吧

a1runner

学习了，懂了很多

bearkingfang

程序直接复制入matlab，建立m文件运行后，为什么输出的边际谱的频率出现在450Hz附近？按照程序内设定的三个imf分量，不是应该是30和45么？有所不解，是不是我对边际谱本身理解不够啊？求高手解释E:\ss.fig

mimijun

我的EMD分解是利用Rilling的程序进行的 根据上面给出的我做了实验可为什么边际谱会这样呢？
clc;close all;clear all
L=1000;
Fs=500;T=1/Fs;
n=(0:L-1)/Fs;
x1=sin(2*pi*84*n);
x2=0.4*sin(2*pi*12*n);
x3=0.9*sin(2*pi*26*n);
sig=x1+x2+x3;
imf=emd(sig,'MAXMODES',6);
%%%%%%%%%%%%%求边际谱
[A,fa,tt]=hhspectrum(imf);
[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=(0:L-3)/L*(Fs/2);
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值')