部分hlbert程序解读（toimage函数）
我在做hilbert谱和边际谱 但每次到下面两句的时候就运行不出来有哪位高手帮忙解读一下 急需谢谢= toimage(A,fa,tt,length(tt));
disp_hhs(E, tt1);
toimage函数究竟是求什么的还有下面那句总是不显示图形
回复 1 # qqzhouqianziyi 的帖子
将提示的错误贴出来
回复 2 # qqvirile 的帖子
比如 我运行的下面的程序：
fs=1000;
N=1000;
t=1/fs:1/fs:1;
y1=2*sin(60*pi*t);
y2=5*sin(90*pi*t);
y=;
=hhspectrum(y);
=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);
figure
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('信号幅值');
title('信号边际谱')%要求边际谱必须先对信号进行EMD分解
运行后E：998x998 double 打开E矩阵显示 Cannot display variables with more than 524288 elements.
正确的图形显示为（1）：
得到下面的图形（未命名1）：
而且每次运行disp_hhs(E, tt1);
也都运行不出来 回复 2 # qqvirile 的帖子
你好 你那有hhspectrum和toimage的源程序么，能不能给我一下 我想看看是不是我的原程序的问题 toimage是把一维的变成二维的显示出来。
你下面的求边际谱的程序没有问题，你怎么不把你出错的求幅值的程序贴出来呢？？
求幅值：>> =hhspectrum(imf);
>> =toimage(A,f);
>> disp_hhs(E,[],fs); 可能是你最后显示语句有问题吧，你按照我的格式改改 回复 5 # chenlu1986 的帖子
你好 你说的意思是这样的么：
fs=1000;
N=1000;
t=1/fs:1/fs:1;
y1=2*sin(60*pi*t);
y2=5*sin(90*pi*t);
y=;
=hhspectrum(y);
=toimage(A,f);
disp_hhs(E,[],fs);我运行了一下，结果显示：
??? Error using ==> disp_hhs
inf doit etre < 0
Error in ==> Untitled at 9
disp_hhs(E,[],fs);
我把最后一句改成disp_hhs(E,t)结果如下图：
回复 6 # qqzhouqianziyi 的帖子
把你最后一句的fs改为你的1000即可 回复 7 # chenlu1986 的帖子
我刚才试了一下 还是不行啊
??? Error using ==> disp_hhs
inf doit etre < 0
Error in ==> bianjipu at 9
disp_hhs(E,1000);
??? Error using ==> disp_hhs
inf doit etre < 0
Error in ==> bianjipu at 9
disp_hhs(E,[],fs);
??? Error using ==> disp_hhs
inf doit etre < 0
Error in ==> bianjipu at 9
disp_hhs(E,[],1000);
回复 7 # chenlu1986 的帖子
你好 能加你QQ么 想跟你在QQ请教请教我那程序一直弄不出来 急用 你是要先进行EMD然后求HT谱和边际谱吗？但是你程序中的这句：y=;
并不是说把两个信号叠加，况且下面也没有分解的emd指令。我只能说按照你给的程序disp_hhs(E,[],1000);在我电脑上可以求出HT谱。但是真的我不明白你的本意，你并没有emd分解。。。。 回复 10 # chenlu1986 的帖子
啊 我这个信号是用的网上的一个y=;相当于一个IMF集你说的那种先把信号进行EMD分解 在求HHT谱和编辑普我也做过跟这个出来的图形一样的只是信号不同而已我是在是整不明白到底是哪的问题跟MATLAB版本有关么 回复 10 # chenlu1986 的帖子
你好 你能帮我做一下HHT时频谱（二维时频和三维时频）和边际谱么
信号为：N=512;
%fft默认计算的信号是从0开始的
t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;
x=2*sin(2*pi*20*t)+6*sin(2*pi*80*t);
谢谢 回复 12 # qqzhouqianziyi 的帖子
这样吧，我把以前做的两个信号叠加然后EMD分解后求HT谱和边际谱的程序附上吧。你可以参照。
clear;
fs=1000;
tspan=2;
t=1/fs:1/fs:tspan;
N=length(t);
x=sin(2*pi*20*t);
y=sin(2*pi*60*t+140);
z=x+y;
plot(t,z);
imf=emd(z);
cemd_visu(z,1:length(z),imf);
%%%imf1的Hilbert变换
xn1=hilbert(imf(1,:));
xr1=real(xn1);
xi1=imag(xn1);
A1=sqrt(xr1.^2+xi1.^2);
figure,subplot(2,1,1);plot(t,A1);
xlabel('时间（t）');ylabel('瞬时振幅');title('imf1')
%%%imf2的Hilbert变换
xn2=hilbert(imf(2,:));
xr2=real(xn2);
xi2=imag(xn2);
A2=sqrt(xr2.^2+xi2.^2);
subplot(2,1,2);plot(t,A2);
xlabel('时间（t）');ylabel('瞬时振幅');title('imf2')
%%%imf1的瞬时相位
P1=atan2(xi1,xr1);
figure,subplot(2,1,1);plot(t,P1);
xlabel('时间（t）'); ylabel('瞬时相位');title('imf1')
%%%imf2的瞬时相位
P2=atan2(xi2,xr2);
subplot(2,1,2);plot(t,P2);
xlabel('时间（t）'); ylabel('瞬时相位');title('imf2')
%%%%imf1瞬时频率
xh1=unwrap(P1);
fs=1000;
xhd1=fs*diff(xh1)/(2*pi);
figure,subplot(2,1,1);plot(t(1:1999),xhd1);
xlabel('时间（t）');ylabel('瞬时频率');title('imf1')
%%%%imf2瞬时频率
xh2=unwrap(P2);
fs=1000;
xhd2=fs*diff(xh2)/(2*pi);
subplot(2,1,2);plot(t(1:1999),xhd2);
xlabel('时间（t）');ylabel('瞬时频率');title('imf2')
>> =hhspectrum(imf);
>> =toimage(A,f);
>> disp_hhs(E,[],1000);
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/tspan*1/fs;
end
H=size(E,1);
f=(0:H-1)/H*(fs/2);
figure ();
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');title('边际谱')
出来的HT谱和边际谱图如下：
回复 13 # chenlu1986 的帖子
你好 你那个程序我试了一下 还是不对找了好几次原因还是没找出来前面的还好 就是后面那几句：
=hhspectrum(imf);
=toimage(A,f);
disp_hhs(E,[],1000);
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/tspan*1/fs;
end
H=size(E,1);
f=(0:H-1)/H*(fs/2);
figure ();
plot(f,bjp);
运行出的图形如下（未命名1、未命名2）：
运行用的源程序见新建文本文档
另外我用的MATLAB7.8.0（R2009a）请高手在帮忙看看我这问题到底出在哪谢谢
本帖最后由 chenlu1986 于 2011-3-21 12:10 编辑
回复 14 # qqzhouqianziyi 的帖子
你安装libsvm工具箱了吗？或者是不是在它的工作路径下运行程序？
先检查一下，我下午找时间再测试一下告诉你。