信号处理matlab程序示例（FFT,HHT）

燕大雪雁

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

亦可见：http://blog.sina.com.cn/s/blog_5d7954df01011jdl.html

1. clear all;
   
2. clc;
   
3. 
   
4. %给定参数
   
5. fs=1000;
   
6. N=1000;
   
7. n=0:N-1;
   
8. t=n/fs;
   
9. %构造仿真信号
   
10. 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));
    
11. %emd分解
    
12. imf=emd(x);
    
13. %绘制时域波形
    
14. figure,plot(t,x);xlabel('时间 t/s');ylabel('幅值');title('仿真信号时域波形');
    
15. %绘制FFT频谱
    
16. nfft= 2^nextpow2(length(t));%找出大于y的个数的最大的2的指数值
    
17. ff=fs*(0:nfft/2-1)/nfft;%FFT变换后对应的频率的序列
    
18. fftx=fft(x,nfft);%求FFT变换
    
19. ps=fftx.*conj(fftx)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是他本身
    
20. figure;
    
21. subplot(211),plot(ff,abs(fftx(1:nfft/2))*2/nfft);ylabel('幅值');xlabel('频率');title('FFT频谱');
    
22. subplot(212),plot(ff,ps(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
    
23. %绘制emd分解结果
    
24. plot_imf(x,t,imf);
    
25. 
    
26. %求时频谱
    
27. [A,f,t]=hhspectrum(imf(1:end-1, :));%对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号
    
28. 
    
29. %绘制瞬时包络图和瞬时频率图
    
30. figure;
    
31. subplot(221),plot(t/N,A(1,:));xlabel('时间 t/s');ylabel('幅值');title('imf1分量瞬时包络');
    
32. subplot(222),plot(t/N,f(1,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf1分量瞬时频率');
    
33. subplot(223),plot(t/N,A(2,:));xlabel('时间 t/s');ylabel('幅值');title('imf2分量瞬时包络');
    
34. subplot(224),plot(t/N,f(2,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf2分量瞬时频率');
    
35. 
    
36. %即时频图(用颜色表示第三维值的大小)和三维图（三维坐标系：时间，中心频率，振幅）
    
37. [E,t,Cenf]=toimage(A,f,t,length(t));%将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率  这里横轴为时间，纵轴为频率        
    
38. 
    
39. %绘制Hilbert-Huang谱
    
40. figure;
    
41. set(gcf,'Color','w');
    
42. imagesc(t/N,[0,0.5*fs],E);
    
43. set(gca,'YDir','normal')
    
44. colormap('jet')
    
45. colorbar;
    
46. xlabel('时间 t/s');
    
47. ylabel('频率 f/Hz');axis([0 1 0 200])
    
48. title('Hilbert-Huang Spectrum') ;
    
49. 
    
50. %画出边际谱
    
51. %N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解
    
52. for k=1:size(E,1)
    
53.     bjp(k)=sum(E(k,:))*1/fs;
    
54. end
    
55. figure;
    
56. plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率
    
57. xlabel('频率 f/Hz');
    
58. ylabel('幅值');
    
59. title('边际谱');
    
60. 
    

复制代码

恋恋依雪

沙发，学习了

905lili

原创的好

无一涯

支持一下

syxqq123

那个频谱和功率谱的量级怎么差那么多，

lx24816

好东西  要顶起啊

maybeloveu

感谢lz分享

wxy_0911

jifeiyun

顶起来

dongwenlian

dynamo05

写的很有启发性，非常感谢，

hubiyang

??? Undefined command/function 'plot_imf'.楼主运行您的程序后出现这个问题是什么原因？

tayo134

好，谢谢

简彦成

blueyese

回复 1 # 燕大雪雁 的帖子

大神 求QQ 指导了
谢谢了