关于破凰的边界谱程序

serendipity_w

关于这段程序，有点儿不明白啊，f=(0:N-3)/N*(fs/2);
这里fs/2是什么意思啊？
还有做出来的图怎么显示频率是比450还大
可是这两个正弦信号的频率分别是30和45啊
HHT.fig (3.92 KB, 下载次数: 40)

2007-7-30 16:40 上传

点击文件名下载附件

%破凰版边界谱分析
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)
    b(k)=sum(E(k,:))*1/fs;  
end
f=(0:N-3)/N*(fs/2);
plot(f,b);
xlabel('频率 / Hz');
ylabel('幅值');

[ 本帖最后由 serendipity_w 于 2007-7-30 16:40 编辑 ]

破凰

新版的不用E=flipud(E);
fs/(2*N)是时频图上y方向的分辨率

serendipity_w

分辨率为什么要处以2N呢？还有是Y方向的分辨率吗？不是X？
还有画出来的谱图上的频率对应真是频率吗？

serendipity_w

哦，明白了，看来就是flipud的问题，谢谢破凰的解答
但是有一点我还是不大明白为什么要除以2N而不是N呢？

zhlong

f=(0:N-3)/N*(fs/2); fs为采样频率，已经对应真实频率了。

频率轴的频率范围是0~fs/2，总共应该有N点的，但是由于diff的原因只能得到N-2个点。这个和FFT的频率刻度设置是一样的。

serendipity_w

明白了，谢谢zhlong的回答！

江湖夜雨

哪位有rilling完整的HHT画三维图，时频谱，边际谱的程序？信号版我没找到完整的，而且用nspab等函数画出来的毕竟不是rilling的正品，哪位热心版友能贴出来让我们参考一下啊

zhlong

http://zhlong.ys168.com

江湖夜雨

我看见了，但是我想要的是使用rilling的程序画图的主程序，不是程序包，这个程序包我有

zhlong

http://forum.vibunion.com/forum/ ... 6amp%3Btypeid%3D187

江湖夜雨

好像没有边际谱，还有三维谱的，我看边际谱帖子的效果也不是特别好，尤其是对实际数据啊

zhlong

实测信号EMD分解效果通常都不是太好，所以边际谱效果可想而知。
至于三维谱，只不过看图的角度不一样，你可以用别的函数如surf、mesh等来画。

zhangnan3509

咱们版也应该把这个问题解决一下了！:handshake