经验模态分解（EMD）

Pseudo-lover

今天看了些EMD信号分解方面的东西，matlab官网上有个Hilbert-Huang Transform的代码，代码效率极高啊，人家3句语句就解决了一个大问题，很牛啊！还有一个GRilling的EMD工具箱，好多文件，功能应该相当强大。

    这里研究了研究matlab官网的代码，加了些注释、功能演示，效果如下

    原始信号由3个正弦信号加噪声组成，如下

下面为做EMD分解的结果

第三次分解信号的瞬时频率如下

第四次分解信号的Hilbert分析

具体代码如下

1. <font color="#000000">test.m文件
   
2. clc
   
3. clear all
   
4. close all
   
5. % [x, Fs] = wavread('Hum.wav');
   
6. % Ts = 1/Fs;
   
7. % x = x(1:6000);
   
8. Ts = 0.001;
   
9. Fs = 1/Ts;
   
10. t=0:Ts:1;
    
11. x = sin(2*pi*10*t) + sin(2*pi*50*t) + sin(2*pi*100*t) + 0.1*randn(1, length(t));
    
12. imf = emd(x);
    
13. plot_hht(x,imf,1/Fs);
    
14. k = 4;
    
15. y = imf{k};
    
16. N = length(y);
    
17. t = 0:Ts:Ts*(N-1);
    
18. [yenvelope, yfreq, yh, yangle] = HilbertAnalysis(y, 1/Fs);
    
19. yModulate = y./yenvelope;
    
20. [YMf, f] = FFTAnalysis(yModulate, Ts);
    
21. Yf = FFTAnalysis(y, Ts);
    
22. figure
    
23. subplot(321)
    
24. plot(t, y)
    
25. title(sprintf('IMF%d', k))
    
26. xlabel('Time/s')
    
27. ylabel(sprintf('IMF%d', k));
    
28. subplot(322)
    
29. plot(f, Yf)
    
30. title(sprintf('IMF%d的频谱', k))
    
31. xlabel('f/Hz')
    
32. ylabel('|IMF(f)|');
    
33. subplot(323)
    
34. plot(t, yenvelope)
    
35. title(sprintf('IMF%d的包络', k))
    
36. xlabel('Time/s')
    
37. ylabel('envelope');
    
38. subplot(324)
    
39. plot(t(1:end-1), yfreq)
    
40. title(sprintf('IMF%d的瞬时频率', k))
    
41. xlabel('Time/s')
    
42. ylabel('Frequency/Hz');
    
43. subplot(325)
    
44. plot(t, yModulate)
    
45. title(sprintf('IMF%d的调制信号', k))
    
46. xlabel('Time/s')
    
47. ylabel('modulation');
    
48. subplot(326)
    
49. plot(f, YMf)
    
50. title(sprintf('IMF%d调制信号的频谱', k))
    
51. xlabel('f/Hz')
    
52. ylabel('|YMf(f)|');
    
53. 
    
54. findpeaks.m文件
    
55. function n = findpeaks(x)
    
56. % Find peaks. 找极大值点，返回对应极大值点的坐标
    
57. n    = find(diff(diff(x) > 0) < 0); % 相当于找二阶导小于0的点
    
58. u    = find(x(n+1) > x(n));
    
59. n(u) = n(u)+1;                      % 加1才真正对应极大值点
    
60. % 图形解释上述过程
    
61. % figure
    
62. % subplot(611)
    
63. % x = x(1:100);
    
64. % plot(x, '-o')
    
65. % grid on
    
66. %
    
67. % subplot(612)
    
68. % plot(1.5:length(x), diff(x) > 0, '-o')
    
69. % grid on
    
70. % axis([1,length(x),-0.5,1.5])
    
71. %
    
72. % subplot(613)
    
73. % plot(2:length(x)-1, diff(diff(x) > 0), '-o')
    
74. % grid on
    
75. % axis([1,length(x),-1.5,1.5])
    
76. %
    
77. % subplot(614)
    
78. % plot(2:length(x)-1, diff(diff(x) > 0)<0, '-o')
    
79. % grid on
    
80. % axis([1,length(x),-1.5,1.5])
    
81. %
    
82. % n    = find(diff(diff(x) > 0) < 0);
    
83. % subplot(615)
    
84. % plot(n, ones(size(n)), 'o')
    
85. % grid on
    
86. % axis([1,length(x),0,2])
    
87. %
    
88. % u    = find(x(n+1) > x(n));
    
89. % n(u) = n(u)+1;
    
90. % subplot(616)
    
91. % plot(n, ones(size(n)), 'o')
    
92. % grid on
    
93. % axis([1,length(x),0,2])
    
94. 
    
95. plot_hht.m文件
    
96. function plot_hht(x,imf,Ts)
    
97. % Plot the HHT.
    
98. % :: Syntax
    
99. %    The array x is the input signal and Ts is the sampling period.
    
100. %    Example on use: [x,Fs] = wavread('Hum.wav');
     
101. %                    plot_hht(x(1:6000),1/Fs);
     
102. % Func : emd
     
103. % imf = emd(x);
     
104. for k = 1:length(imf)
     
105.     b(k) = sum(imf{k}.*imf{k});
     
106.     th   = unwrap(angle(hilbert(imf{k})));  % 相位
     
107.     d{k} = diff(th)/Ts/(2*pi);          % 瞬时频率
     
108. end
     
109. [u,v] = sort(-b);
     
110. b     = 1-b/max(b);                     % 后面绘图的亮度控制
     
111. % Hilbert瞬时频率图
     
112. N = length(x);
     
113. c = linspace(0,(N-2)*Ts,N-1);           % 0:Ts:Ts*(N-2)
     
114. for k = v(1:2)                          % 显示能量最大的两个IMF的瞬时频率
     
115.     figure
     
116.     plot(c,d{k});
     
117.     xlim([0 c(end)]);
     
118.     ylim([0 1/2/Ts]);
     
119.     xlabel('Time/s')
     
120.     ylabel('Frequency/Hz');
     
121.     title(sprintf('IMF%d', k))
     
122. end
     
123. % 显示各IMF
     
124. M = length(imf);
     
125. N = length(x);
     
126. c = linspace(0,(N-1)*Ts,N);             % 0:Ts:Ts*(N-1)
     
127. for k1 = 0:4:M-1
     
128.     figure
     
129.     for k2 = 1:min(4,M-k1)
     
130.         subplot(4,2,2*k2-1)
     
131.         plot(c,imf{k1+k2})
     
132.         set(gca,'FontSize',8,'XLim',[0 c(end)]);
     
133.         title(sprintf('第%d个IMF', k1+k2))
     
134.         xlabel('Time/s')
     
135.         ylabel(sprintf('IMF%d', k1+k2));
     
136.       
     
137.         subplot(4,2,2*k2)
     
138.         [yf, f] = FFTAnalysis(imf{k1+k2}, Ts);
     
139.         plot(f, yf)
     
140.         title(sprintf('第%d个IMF的频谱', k1+k2))
     
141.         xlabel('f/Hz')
     
142.         ylabel('|IMF(f)|');
     
143.     end
     
144. end
     
145. figure
     
146. subplot(211)
     
147. plot(c,x)
     
148. set(gca,'FontSize',8,'XLim',[0 c(end)]);
     
149. title('原始信号')
     
150. xlabel('Time/s')
     
151. ylabel('Origin');
     
152. subplot(212)
     
153. [Yf, f] = FFTAnalysis(x, Ts);
     
154. plot(f, Yf)
     
155. title('原始信号的频谱')
     
156. xlabel('f/Hz')
     
157. ylabel('|Y(f)|');
     
158. emd.m文件
     
159. function imf = emd(x)
     
160. % Empiricial Mode Decomposition (Hilbert-Huang Transform)
     
161. % EMD分解或HHT变换
     
162. % 返回值为cell类型，依次为一次IMF、二次IMF、...、最后残差
     
163. x   = transpose(x(:));
     
164. imf = [];
     
165. while ~ismonotonic(x)
     
166.     x1 = x;
     
167.     sd = Inf;
     
168.     while (sd > 0.1) || ~isimf(x1)
     
169.         s1 = getspline(x1);         % 极大值点样条曲线
     
170.         s2 = -getspline(-x1);       % 极小值点样条曲线
     
171.         x2 = x1-(s1+s2)/2;
     
172.       
     
173.         sd = sum((x1-x2).^2)/sum(x1.^2);
     
174.         x1 = x2;
     
175.     end
     
176.    
     
177.     imf{end+1} = x1;
     
178.     x          = x-x1;
     
179. end
     
180. imf{end+1} = x;
     
181. % 是否单调
     
182. function u = ismonotonic(x)
     
183. u1 = length(findpeaks(x))*length(findpeaks(-x));
     
184. if u1 > 0
     
185.     u = 0;
     
186. else
     
187.     u = 1;
     
188. end
     
189. % 是否IMF分量
     
190. function u = isimf(x)
     
191. N  = length(x);
     
192. u1 = sum(x(1:N-1).*x(2:N) < 0);                     % 过零点的个数
     
193. u2 = length(findpeaks(x))+length(findpeaks(-x));    % 极值点的个数
     
194. if abs(u1-u2) > 1
     
195.     u = 0;
     
196. else
     
197.     u = 1;
     
198. end
     
199. % 据极大值点构造样条曲线
     
200. function s = getspline(x)
     
201. N = length(x);
     
202. p = findpeaks(x);
     
203. s = spline([0 p N+1],[0 x(p) 0],1:N);
     
204. FFTAnalysis.m文件
     
205. % 频谱分析
     
206. function [Y, f] = FFTAnalysis(y, Ts)
     
207. Fs = 1/Ts;
     
208. L = length(y);
     
209. NFFT = 2^nextpow2(L);
     
210. y = y - mean(y);
     
211. Y = fft(y, NFFT)/L;
     
212. Y = 2*abs(Y(1:NFFT/2+1));
     
213. f = Fs/2*linspace(0, 1, NFFT/2+1);
     
214. end
     
215. HilbertAnalysis.m文件
     
216. % Hilbert分析
     
217. function [yenvelope, yf, yh, yangle] = HilbertAnalysis(y, Ts)
     
218. yh = hilbert(y);
     
219. yenvelope = abs(yh);                % 包络
     
220. yangle = unwrap(angle(yh));         % 相位
     
221. yf = diff(yangle)/2/pi/Ts;          % 瞬时频率
     
222. end</font>

复制代码

转自：http://blog.sina.com.cn/s/blog_6163bdeb0102e2b3.html