嗯，这些程序我都有，我有三个版本的程序：National Taiwan Central University, Matlab File Exchange Center，和G Rilling（2007版）的。我试过了后两个版本：Matlab File Exchange Center是分解出IMF，然后画出 时间-频率图，而幅度好像是用灰度（我不确定）表示的，没有给出画幅度-频率-时间三维图的画法，也没有给出边际谱的画法；而G rilling的程序也是一样，我在论坛上看到的画A-f-t谱图和边际谱图的的方法是：
>> x=xlsread('matlab');
>> y=x(:,2);
>> imf=emd(y);
>>[A,fa,tt]=hhspectrum(imf);
>> [E,tt1]=toimage(A,fa,tt,length(tt));
>> disp_hhs(E,tt1) %画Hilbert-Huang spectrum
>>fs=500000;%采样频率500kHz
>> for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/fs;
end
>>f=(0:(length(y)-3))/(length(y))*(fs/2);%频率
>>plot(f,bjp)； %画边际谱
其实我不理解为什么按照上面求得的bjp(k)=sum(E(k,:))*1/fs;就是边际谱？难道函数toimage的功能就是排序吗？在这里，我的数据y是9000个点，得到的E是8998*8998的矩阵，难道E的每一行就代表同一个瞬时频率下的8998个瞬时幅值吗？
还有一个问题，有解释说边际谱的纵轴并不是实际幅值，而是代表某个频率出现的概率的大小，如果是这样的话，sum(bjp)是否应该等于1呀？可是我求得的值是sum(bjp)=0.0034，完全不等于1嘛，晕掉了。。。
还有，你看我用下面的程序求瞬时幅度-时间谱和瞬时频率-时间谱可以吗？
igure;
for k=1:length(imf) %幅度-时间谱
plot(tt,A(k,:));
hold on;
end
figure;
for k=1:length(imf) %频率-时间谱
plot(tt,fa(k,:));
hold on;
end
实际上就是把所有的瞬时幅度叠加在一张图上，把所有的瞬时频率叠加在另一张图上，这涉及到我对Hilbert-Huang Spectrum的理解，我不知道对不对。
另外，我试图用disp_hhs(E,tt1)画图时，出现错误了，??? Out of memory. Type HELP MEMORY for your options.
Error in ==> disp_hhs at 67
im = 10*log10(im/M);
你知道是什么意思吗？很奇怪，好像和调用toimage时的参数length(tt)有关，如果我去掉参数length(tt)，运行[E,tt1]=toimage(A,fa,tt,);之后再运行disp_hhs(E,tt1)就可以出图，但这个时候的E变成400*8998的矩阵了。。。如果要完整的画出Hilbert-Huang Spectrum，一定要用disp_hhs吗，这画出来的图背景是蓝色的，横坐标是时间，纵坐标是normalized frequency，那幅度怎么表示呢？
唉，不懂的太多，请你不吝赐教，等你的回复。。。谢谢了！！！ |
[复制链接]
发表于 2014-6-26 13:45
显身卡
