部分hlbert程序解读（toimage函数）

qqzhouqianziyi

我在做hilbert谱和边际谱   但每次到下面两句的时候就运行不出来  有哪位高手帮忙解读一下 急需  谢谢
[E, tt1] = toimage(A,fa,tt,length(tt));
disp_hhs(E, tt1);
toimage函数究竟是求什么的  还有下面那句总是不显示图形

qqvirile

回复 1 # qqzhouqianziyi 的帖子

将提示的错误贴出来

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=[y1;y2;zeros(size(y1))];
[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);
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);
也都运行不出来

qqzhouqianziyi

回复 2 # qqvirile 的帖子

你好 你那有hhspectrum和toimage的源程序么，能不能给我一下 我想看看是不是我的原程序的问题

chenlu1986

toimage是把一维的变成二维的显示出来。
你下面的求边际谱的程序没有问题，你怎么不把你出错的求幅值的程序贴出来呢？？
求幅值：>> [A,f,t]=hhspectrum(imf);
>> [E,t,cenf]=toimage(A,f);
>> disp_hhs(E,[],fs); 可能是你最后显示语句有问题吧，你按照我的格式改改

qqzhouqianziyi

回复 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=[y1;y2;zeros(size(y1))];
[A,f,t]=hhspectrum(y);
[E,t,cenf]=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)结果如下图：

chenlu1986

回复 6 # qqzhouqianziyi 的帖子

把你最后一句的fs改为你的1000即可

qqzhouqianziyi

回复 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);

qqzhouqianziyi

回复 7 # chenlu1986 的帖子

你好 能加你QQ么 想跟你在QQ请教请教  我那程序一直弄不出来 急用

chenlu1986

你是要先进行EMD然后求HT谱和边际谱吗？但是你程序中的这句：y=[y1;y2;zeros(size(y1))];
并不是说把两个信号叠加，况且下面也没有分解的emd指令。我只能说按照你给的程序disp_hhs(E,[],1000);在我电脑上可以求出HT谱。但是真的我不明白你的本意，你并没有emd分解。。。。

qqzhouqianziyi

回复 10 # chenlu1986 的帖子

啊 我这个信号是用的网上的一个  y=[y1;y2;zeros(size(y1))];相当于一个IMF集  你说的那种先把信号进行EMD分解 在求HHT谱和编辑普我也做过  跟这个出来的图形一样的  只是信号不同而已  我是在是整不明白到底是哪的问题  跟MATLAB版本有关么

qqzhouqianziyi

回复 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);
谢谢

chenlu1986

回复 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')
>> [A,f,t]=hhspectrum(imf);
>> [E,t,cenf]=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谱和边际谱图如下：

qqzhouqianziyi

回复 13 # chenlu1986 的帖子

你好 你那个程序我试了一下 还是不对  找了好几次原因  还是没找出来  前面的还好 就是后面那几句：
[A,f,t]=hhspectrum(imf);
[E,t,cenf]=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

回复 14 # qqzhouqianziyi 的帖子

你安装libsvm工具箱了吗？或者是不是在它的工作路径下运行程序？
先检查一下，我下午找时间再测试一下告诉你。