1.什么是HHT？
HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。
2.EMD分解的步骤。
EMD分解的流程图如下：
3.实例演示。
给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t)+5*sin(2*pi*35t)
(1)为了对比，先用fft对求上述信号的幅频和相频曲线。
1. 2. 3. 4. 5. 6. 7. 8.
function fftfenxi clear;clc; N=2048;
?t默认计算的信号是从0开始的
t=linspace(1,2,N);deta=t(2)-t(1);1/deta x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
% N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi; %
x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t); 9. y = x; 10. m=0:N-1;
11. f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的 12. %下面计算的Y就是x(t)的傅里叶变换数值
13. %Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移
后[-2,2]之间的频谱值 14. Y=fft(y);
15. z=sqrt(Y.*conj(Y));
16. plot(f(1:100),z(1:100)); 17. title('幅频曲线') 18. xiangwei=angle(Y);
19. figure(2)
20. plot(f,xiangwei) 21. title('相频曲线') 22. figure(3) 23. plot(t,y,'r')
24. %axis([-2,2,0,1.2]) 25. title('原始信号')
复制代码
（2）用Hilbert变换直接求该信号的瞬时频率
1. clear;clc;clf;
2. %假设待分析的函数是z=t^3 3. N=2048;
4. ?t默认计算的信号是从0开始的
5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta; 6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t); 7. z=x;
8. hx=hilbert(z);
9. xr=real(hx);xi=imag(hx); 10. %计算瞬时振幅
11. sz=sqrt(xr.^2+xi.^2); 12. %计算瞬时相位 13. sx=angle(hx); 14. %计算瞬时频率 15. dt=diff(t); 16. dx=diff(sx); 17. sp=dx./dt;
18. plot(t(1:N-1),sp) 19. title('瞬时频率') 20.
复制代码
小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。Hilbert变换是求取瞬时频率的方法，但如果只用Hilbert变换求出来的瞬时频率也不准确。（出现负频，实际上负频没有意义！）
（3）用HHT求取信号的时频谱与边际谱
1. function HHT 2. clear;clc;clf; 3. N=2048;
4. ?t默认计算的信号是从0开始的
5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta; 6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t); 7. z=x;
8. c=emd(z); 9.
10. %计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性 11. [m,n]=size(c); 12. for i=1:m;
13. a=corrcoef(c(i,:),z); 14. xg(i)=a(1,2); 15. end 16. xg; 17.
18. for i=1:m-1
19. %--------------------------------------------------------------------
20. %计算各IMF的方差贡献率
21. %定义：方差为平方的均值减去均值的平方 22. %均值的平方
23. %imfp2=mean(c(i,:),2).^2 24. %平方的均值
25. %imf2p=mean(c(i,:).^2,2) 26. %各个IMF的方差
27. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; 28. end;
29. mmse=sum(mse); 30. for i=1:m-1
31. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; 32. %方差百分比，也就是方差贡献率 33. mseb(i)=mse(i)/mmse*100; 34. %显示各个IMF的方差和贡献率 35. end;
36. %画出每个IMF分量及最后一个剩余分量residual的图形 37. figure(1) 38. for i=1:m-1
39. disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]); 40. end;
41. subplot(m+1,1,1) 42. plot(t,z)
43. set(gca,'fontname','times New Roman') 44. set(gca,'fontsize',14.0)
45. ylabel(['signal','Amplitude']) 46.
47. for i=1:m-1
48. subplot(m+1,1,i+1); 49. set(gcf,'color','w') 50. plot(t,c(i,:),'k')
51. set(gca,'fontname','times New Roman') 52. set(gca,'fontsize',14.0) 53. ylabel(['imf',int2str(i)]) 54. end
55. subplot(m+1,1,m+1); 56. set(gcf,'color','w') 57. plot(t,c(m,:),'k')
58. set(gca,'fontname','times New Roman') 59. set(gca,'fontsize',14.0) 60. ylabel(['r',int2str(m-1)]) 61.
62. %画出每个IMF分量及剩余分量residual的幅频曲线 63. figure(2)
64. subplot(m+1,1,1) 65. set(gcf,'color','w') 66. [f,z]=fftfenxi(t,z); 67. plot(f,z,'k')
68. set(gca,'fontname','times New Roman') 69. set(gca,'fontsize',14.0)
70. ylabel(['initial signal',int2str(m-1),'Amplitude']) 71.
72. for i=1:m-1
73. subplot(m+1,1,i+1); 74. set(gcf,'color','w')
75. [f,z]=fftfenxi(t,c(i,:)); 76. plot(f,z,'k')
77. set(gca,'fontname','times New Roman') 78. set(gca,'fontsize',14.0)
79. ylabel(['imf',int2str(i),'Amplitude']) 80. end
81. subplot(m+1,1,m+1); 82. set(gcf,'color','w')
83. [f,z]=fftfenxi(t,c(m,:)); 84. plot(f,z,'k')
85. set(gca,'fontname','times New Roman') 86. set(gca,'fontsize',14.0)
87. ylabel(['r',int2str(m-1),'Amplitude']) 88.
89. hx=hilbert(z);
90. xr=real(hx);xi=imag(hx); 91. %计算瞬时振幅
92. sz=sqrt(xr.^2+xi.^2); 93. %计算瞬时相位 94. sx=angle(hx); 95. %计算瞬时频率 96. dt=diff(t); 97. dx=diff(sx); 98. sp=dx./dt; 99. figure(6) 100. plot(t(1:N-1),sp) 101. title('瞬时频率') 102. 103. %计算HHT时频谱和边际谱 104. [A,fa,tt]=hhspectrum(c); 105. [E,tt1]=toimage(A,fa,tt,length(tt)); 106. figure(3) 107. disp_hhs(E,tt1) %二维图显示HHT时频谱，E是求得的HHT谱 108. pause 109. figure(4) 110. for i=1:size(c,1) 111. faa=fa(i,:); 112. [FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图 113. surf(FA,TT1,E) 114. title('HHT时频谱三维显示') 115. hold on 116. end 117. hold off 118. E=flipud(E); 119. for k=1:size(E,1) 120. bjp(k)=sum(E(k,:))*1/fs; 121. end 122. f=(1:N-2)/N*(fs/2); 123. figure(5) 124. plot(f,bjp); 125. xlabel('频率 / Hz'); 126. ylabel('信号幅值'); 127. title('信号边际谱')%要求边际谱必须先对信号进行EMD分解 128. 129. function [A,f,tt] = hhspectrum(x,t,l,aff) 130.