声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1780|回复: 6

[HHT] 求助如何用matlab语言表示时变信号,EMD端点延拓法如何使用?

[复制链接]
发表于 2013-12-4 21:01 | 显示全部楼层 |阅读模式

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

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

x
致版主:
       请不要再删了,都删了两次了,我提的问题都在论坛里找过,要不就是有相关问题没回答,要不压根就没回答。多谢了。

小弟最近研究HHT算法,现在准备用HHT 和FFT对一时变信号做比较,但时变信号不会在MATLAB中表示,
时变信号如下:
t<0.05 f(t)=sin(2*pi*50*t)+0.2*sin(2*pi*3*50*t)
t≥0.05 f(t)=sin(2*pi*50*t)+0.2*sin(2*pi*5*50*t)

是个很基本的问题,但小弟初学者,真的不会,请大家帮忙。

另外,在抑制EMD端点效应时,看到论坛里有人贴出了次端点抑制法,具体怎么操作?请懂的人说下,感激不尽!

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2013-12-6 10:38 | 显示全部楼层
x=(sin(2*pi*50*t)+sin(2*pi*150*t)).*(t>=0&t<0.05)+(sin(2*pi*50*t)+sin(2*pi*250*t)).*(t>=0.05&t<=1);
写出来了,但是结果还是不对。
 楼主| 发表于 2013-12-6 10:43 | 显示全部楼层
下面将程序列出来,希望有大神能指导下。
  1. clear all;
  2. clc;

  3. %示例程序
  4. N=1024;
  5. n=0:N-1;
  6. fs=12800;
  7. t=n/fs;
  8. x=(sin(2*pi*50*t)+sin(2*pi*150*t)).*(t>=0&t<0.5)+(sin(2*pi*50*t)+sin(2*pi*250*t)).*(t>=0.5&t<=1);
  9. data=x;
  10. imf=emd(data);                        %对输入信号进行EMD分解   
  11. [A,f,t]=hhspectrum(imf);            %对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号
  12. [E,t,Cenf]=toimage(A,f);            %将每个IMF信号合成求取Hilbert谱,E:对应的振幅值,Cenf:每个网格对应的中心频率  这里横轴为时间,纵轴为频率        
  13.                                                    %即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系:时间,中心频率,振幅)         
  14. cemd_visu(data,1:length(data),imf);   %显示每个IMF分量及残余信号--------------------------------------------
  15. disp_hhs(E);                          %希尔伯特谱----------------------------------------------------------
  16. %画出边际谱
  17. colormap(flipud(gray));               %黑白显示
  18. %N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
  19. for k=1:size(E,1)
  20.     bjp(k)=sum(E(k,:))*1/fs;
  21. end
  22. figure(3);
  23. plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
  24. xlabel('频率 / Hz');
  25. ylabel('幅值');

  26. %绘制瞬时包络图和瞬时频率图
  27. figure;
  28. subplot(221),plot(t/N,A(1,:));xlabel('时间 t/s');ylabel('幅值');title('imf1分量瞬时包络');
  29. subplot(222),plot(t/N,f(1,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf1分量瞬时频率');
  30. subplot(223),plot(t/N,A(2,:));xlabel('时间 t/s');ylabel('幅值');title('imf2分量瞬时包络');
  31. subplot(224),plot(t/N,f(2,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf2分量瞬时频率');
复制代码
下面是仿真的结果,还有我理想的结果(黑白的)。
求大神指教!

imf2.png
imf1.png
发表于 2013-12-6 11:22 | 显示全部楼层
本帖最后由 yghit08 于 2013-12-6 11:27 编辑

都提示了,生成数据段有问题,还需要怎么说。
不知道找相近的帖子回复吗?
删除后你是否考虑过自己的问题,还是一个劲的发相同内容的帖子,肯定被删除。每个新手都希望自己的帖子被大家关注,但是你的帖子里的内容只对你自己一个人有用,而且还是自己编程的问题,一方面将误导HHT用户,另一方面惯坏了你这种不思进取的发帖、学习习惯。
%%%%%%%%
clear,clc
N=1280;
fs=1280;
dt=1/fs;
t=dt:dt:N*dt;
for i=1:N
if t(i)<0.5
x(i)=sin(2*pi*50*t(i));
else
x(i)=sin(2*pi*150*t(i));
end
end
%%%%数据生成段
关于端点延拓,默认的是镜像法,运行过程中有默认,仔细看程序,单独讨论端点问题的帖子亦有,请搜索板块。
无任何新意的EMD方面的新帖都不建议开帖,跟帖回复。
 楼主| 发表于 2013-12-7 20:40 | 显示全部楼层
版主我要先感谢你,因为最后的结果还是出来了,这个程序段在你看来确实很简单,但是我真的很不明白错在哪里,整整三天都没法,这才上论坛求助。我自己也加了很多群,但是因为专业不对所以我还是没有得到我想要的结果,你给我的程序段我试过了,加进去之的结果如下: imf1.png
还是得不到截图上面的结果。

PS:我将你给我的程序段加进去后,改动了一下频率,结果又和我自己做的波形一样,我想可能是频率设置有问题,我知道不该重复发帖,也知道你审核的很严,我这里有个HHT的工具箱,是目前最全的,是我自己搜集整理的,但是因为大小问题始终无法上传,我很郁闷,因为自己明明可以帮忙的,但因为附件大小问题帮不上,请版主适当放宽我的上传大小,我一定好好为大家提供我自己在HHT编程上的问题,毕竟是初学者,我遇到的问题都是一些小错误,希望版主能理解,再次感谢!以后不会再重复发帖的。
 楼主| 发表于 2013-12-7 21:11 | 显示全部楼层
版主,你好,我又调了下,我将N和fs都调整到了2560,结果正确了,如下所示:
imf2.png

应该还是频率设置的问题,当然我数据给的也有问题,谢谢版主的指导,对我意义重大,再次感谢!以后绝不重复发帖。

点评

不是我严格,主要是你的问题主要是个人编程问题,和算法一点关系都没有,问题也不典型。且目前我个人也不希望这个板块主要在讨论HHT这个可能没有前途的算法  发表于 2013-12-7 21:13
发表于 2013-12-7 21:11 | 显示全部楼层
猫的礼拜七 发表于 2013-12-7 20:40
版主我要先感谢你,因为最后的结果还是出来了,这个程序段在你看来确实很简单,但是我真的很不明白错在哪里, ...

附件大小是论坛的设置,我没法给你调大,可以分割成好几个上传。
另:一方面你说解决了问题,另一方面又说没得到你想要的结果,我感到困惑。
建议:频率别设置太高,太高数据点数太短,还是在数据生成阶段会出问题(循环段需要改数据长度),且数据长度过长使得运行时间比较长,可能一般的电脑没法运行。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-25 09:59 , Processed in 0.077918 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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