声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 19013|回复: 90

[综合] 信号处理matlab程序示例(FFT,HHT)

  [复制链接]
发表于 2012-4-26 19:56 | 显示全部楼层 |阅读模式

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

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

x
在研究和学习利用matlab进行信号分析的过程中,自己写的一个程序,贴出来大家一块学习一下。

亦可见:http://blog.sina.com.cn/s/blog_5d7954df01011jdl.html
  1. clear all;
  2. clc;

  3. %给定参数
  4. fs=1000;
  5. N=1000;
  6. n=0:N-1;
  7. t=n/fs;
  8. %构造仿真信号
  9. x=2*sin(2*pi*120*t)+(1+sin(2*pi*5*t)).*sin(2*pi*50*t+0.5*cos(2*pi*5*t));
  10. %emd分解
  11. imf=emd(x);
  12. %绘制时域波形
  13. figure,plot(t,x);xlabel('时间 t/s');ylabel('幅值');title('仿真信号时域波形');
  14. %绘制FFT频谱
  15. nfft= 2^nextpow2(length(t));%找出大于y的个数的最大的2的指数值
  16. ff=fs*(0:nfft/2-1)/nfft;%FFT变换后对应的频率的序列
  17. fftx=fft(x,nfft);%求FFT变换
  18. ps=fftx.*conj(fftx)/nfft;%conj()函数是求y函数的共轭复数,实数的共轭复数是他本身
  19. figure;
  20. subplot(211),plot(ff,abs(fftx(1:nfft/2))*2/nfft);ylabel('幅值');xlabel('频率');title('FFT频谱');
  21. subplot(212),plot(ff,ps(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
  22. %绘制emd分解结果
  23. plot_imf(x,t,imf);

  24. %求时频谱
  25. [A,f,t]=hhspectrum(imf(1:end-1, :));%对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号

  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分量瞬时频率');

  32. %即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系:时间,中心频率,振幅)
  33. [E,t,Cenf]=toimage(A,f,t,length(t));%将每个IMF信号合成求取Hilbert谱,E:对应的振幅值,Cenf:每个网格对应的中心频率  这里横轴为时间,纵轴为频率        

  34. %绘制Hilbert-Huang谱
  35. figure;
  36. set(gcf,'Color','w');
  37. imagesc(t/N,[0,0.5*fs],E);
  38. set(gca,'YDir','normal')
  39. colormap('jet')
  40. colorbar;
  41. xlabel('时间 t/s');
  42. ylabel('频率 f/Hz');axis([0 1 0 200])
  43. title('Hilbert-Huang Spectrum') ;

  44. %画出边际谱
  45. %N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
  46. for k=1:size(E,1)
  47.     bjp(k)=sum(E(k,:))*1/fs;
  48. end
  49. figure;
  50. plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
  51. xlabel('频率 f/Hz');
  52. ylabel('幅值');
  53. title('边际谱');

复制代码


1.jpg
2.jpg
3.jpg
4.jpg
5.jpg
6.jpg

点评

赞成: 5.0
赞成: 5
  发表于 2014-6-30 23:48

评分

3

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2012-4-27 09:00 | 显示全部楼层
沙发,学习了
发表于 2012-4-27 09:05 | 显示全部楼层
原创的好
发表于 2012-4-27 23:14 | 显示全部楼层
支持一下
发表于 2012-4-28 14:45 | 显示全部楼层
那个频谱和功率谱的量级怎么差那么多,
发表于 2012-5-2 20:49 | 显示全部楼层
好东西  要顶起啊
发表于 2012-5-3 08:49 | 显示全部楼层
感谢lz分享
发表于 2012-5-3 10:40 | 显示全部楼层
发表于 2012-5-3 14:57 | 显示全部楼层
顶起来
发表于 2012-5-3 18:07 | 显示全部楼层
发表于 2012-5-8 03:25 | 显示全部楼层
写的很有启发性,非常感谢,
发表于 2012-5-9 14:38 | 显示全部楼层
??? Undefined command/function 'plot_imf'.楼主运行您的程序后出现这个问题是什么原因?
发表于 2012-5-9 16:26 | 显示全部楼层
好,谢谢
发表于 2012-5-9 19:04 | 显示全部楼层
发表于 2012-5-9 21:14 | 显示全部楼层
回复 1 # 燕大雪雁 的帖子

大神 求QQ 指导了
谢谢了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-24 21:30 , Processed in 0.083030 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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