声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 18765|回复: 62

[HHT] 部分hlbert程序解读(toimage函数)

  [复制链接]
发表于 2011-3-17 15:52 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

回复
分享到:

使用道具 举报

发表于 2011-3-17 18:17 | 显示全部楼层
回复 1 # qqzhouqianziyi 的帖子

将提示的错误贴出来
 楼主| 发表于 2011-3-18 08:51 | 显示全部楼层
回复 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);
也都运行不出来
未命名1.jpg
1.jpg
 楼主| 发表于 2011-3-18 09:37 | 显示全部楼层
回复 2 # qqvirile 的帖子

你好 你那有hhspectrum和toimage的源程序么,能不能给我一下 我想看看是不是我的原程序的问题
发表于 2011-3-18 09:50 | 显示全部楼层
toimage是把一维的变成二维的显示出来。
你下面的求边际谱的程序没有问题,你怎么不把你出错的求幅值的程序贴出来呢??
求幅值:>> [A,f,t]=hhspectrum(imf);
>> [E,t,cenf]=toimage(A,f);
>> disp_hhs(E,[],fs); 可能是你最后显示语句有问题吧,你按照我的格式改改
 楼主| 发表于 2011-3-18 15:42 | 显示全部楼层
回复 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)结果如下图:

未命名3.jpg
发表于 2011-3-19 15:42 | 显示全部楼层
回复 6 # qqzhouqianziyi 的帖子

把你最后一句的fs改为你的1000即可
 楼主| 发表于 2011-3-19 15:47 | 显示全部楼层
回复 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);
 楼主| 发表于 2011-3-19 16:00 | 显示全部楼层
回复 7 # chenlu1986 的帖子

你好 能加你QQ么 想跟你在QQ请教请教  我那程序一直弄不出来 急用
发表于 2011-3-19 22:08 | 显示全部楼层
你是要先进行EMD然后求HT谱和边际谱吗?但是你程序中的这句:y=[y1;y2;zeros(size(y1))];
并不是说把两个信号叠加,况且下面也没有分解的emd指令。我只能说按照你给的程序disp_hhs(E,[],1000);在我电脑上可以求出HT谱。但是真的我不明白你的本意,你并没有emd分解。。。。
 楼主| 发表于 2011-3-20 09:44 | 显示全部楼层
回复 10 # chenlu1986 的帖子

啊 我这个信号是用的网上的一个  y=[y1;y2;zeros(size(y1))];相当于一个IMF集  你说的那种先把信号进行EMD分解 在求HHT谱和编辑普我也做过  跟这个出来的图形一样的  只是信号不同而已  我是在是整不明白到底是哪的问题  跟MATLAB版本有关么
 楼主| 发表于 2011-3-20 10:06 | 显示全部楼层
回复 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);
谢谢
发表于 2011-3-20 19:53 | 显示全部楼层
回复 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谱和边际谱图如下:

HT谱.jpg
边际谱.jpg
 楼主| 发表于 2011-3-21 10:51 | 显示全部楼层
回复 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)请高手在帮忙看看我这问题到底出在哪  谢谢
未命名1.jpg
未命名2.jpg

新建 文本文档.rar

2.16 KB, 下载次数: 27

发表于 2011-3-21 12:09 | 显示全部楼层
本帖最后由 chenlu1986 于 2011-3-21 12:10 编辑

回复 14 # qqzhouqianziyi 的帖子

你安装libsvm工具箱了吗?或者是不是在它的工作路径下运行程序?
先检查一下,我下午找时间再测试一下告诉你。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 10:51 , Processed in 0.068595 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表