LMD算法Python程序

陈波

自己写的一个简单的LMD算法程序，需要 Python 3.0 或以上版本。
注意：未作计算优化，速度较慢。

欢迎大家提意见。

1. 
   
2. """
   
3. This module implements Local Mean Decomposition signal processing algorithm
   
4. as described in [1].
   
5. 
   
6. 
   
7. References:
   
8. 
   
9. [1] Jonathan S. Smith. The local mean decomposition and its application to EEG
   
10. perception data. Journal of the Royal Society Interface, 2005, 2(5):443-454
    
11. """
    
12. 
    
13. #
    
14. # 算法控制参数
    
15. #
    
16. 
    
17. # 是否把信号的端点看作伪极值点
    
18. INCLUDE_ENDPOINTS = True
    
19. 
    
20. # 滑动平均算法: 最多迭代的次数
    
21. MAX_SMOOTH_ITERATION = 12
    
22. 
    
23. # 分离局部包络信号时最多迭代的次数
    
24. MAX_ENVELOPE_ITERATION = 200
    
25. 
    
26. ENVELOPE_EPSILON = 0.01
    
27. CONVERGENCE_EPSILON = 0.01
    
28. 
    
29. # 最多生成的积函数个数
    
30. MAX_NUM_PF = 8
    
31. 
    
32. # for debugging
    
33. DBG_FIG = None
    
34. 
    
35. def is_monotonous(signal):
    
36.     """判断信号是否是(非严格的)单调序列"""
    
37. 
    
38.     def is_monotonous_increase(signal):
    
39.         y0 = signal[0]
    
40.         for y1 in signal:
    
41.             if y1 < y0:
    
42.                 return False
    
43.             y0 = y1
    
44.         return True
    
45. 
    
46.     def is_monotonous_decrease(signal):
    
47.         y0 = signal[0]
    
48.         for y1 in signal:
    
49.             if y1 > y0:
    
50.                 return False
    
51.             y0 = y1
    
52.         return True
    
53. 
    
54.     if len(signal) <= 0:
    
55.         return True
    
56.     else:
    
57.         return is_monotonous_increase(signal) \
    
58.             or is_monotonous_decrease(signal)
    
59. 
    
60. def find_extrema(signal):
    
61.     """找出信号的所有局部极值点位置."""
    
62.    
    
63.     n = len(signal)
    
64. 
    
65.     # 局部极值的位置列表
    
66.     extrema = [ ]
    
67. 
    
68.     lookfor = None
    
69.     for x in range(1, n):
    
70.         y1 = signal[x-1]; y2 = signal[x]
    
71.         if lookfor == "min":
    
72.             if y2 > y1:
    
73.                 extrema.append(x-1)
    
74.                 lookfor = "max"
    
75.         elif lookfor == "max":
    
76.             if y2 < y1:
    
77.                 extrema.append(x-1)
    
78.                 lookfor = "min"
    
79.         else:
    
80.             if y2 < y1:
    
81.                 lookfor = "min"
    
82.             elif y2 > y1:
    
83.                 lookfor = "max"
    
84. 
    
85.     # 优化(微调)极值点的位置:
    
86.     # 当连续多个采样值相同时，取最中间位置作为极值点位置
    
87. 
    
88.     def micro_adjust(x):
    
89.         assert(0 <= x < n)
    
90.         y = signal[x]; lo = hi = x
    
91.         while (lo-1 >= 0) and (signal[lo-1] == y): lo -= 1
    
92.         while (hi+1 < n) and (signal[hi+1] == y): hi += 1
    
93.         if INCLUDE_ENDPOINTS:
    
94.             if lo == 0: return 0
    
95.             if hi == n-1: return n-1
    
96.         return (lo + hi) // 2
    
97. 
    
98.     extrema = [micro_adjust(x) for x in extrema]
    
99. 
    
100.     if extrema and INCLUDE_ENDPOINTS:
     
101.         if extrema[0] != 0:
     
102.             extrema.insert(0, 0)
     
103.         if extrema[-1] != n-1:
     
104.             extrema.append(n-1)
     
105. 
     
106.     return extrema
     
107. 
     
108. def moving_average_smooth(signal, window):
     
109.     """使用移动加权平均算法平滑方波信号."""
     
110. 
     
111.     n = len(signal)
     
112. 
     
113.     # at least one nearby sample is needed for average
     
114.     if window < 3:
     
115.         window = 3
     
116. 
     
117.     # adjust length of sliding window to an odd number for symmetry
     
118.     if (window % 2) == 0:
     
119.         window += 1
     
120.         
     
121.     half = window // 2
     
122. 
     
123.     weight = list(range(1, half+2)) + list(range(half, 0, -1))
     
124.     assert(len(weight) == window)
     
125. 
     
126.     def sliding(signal, window):
     
127.         for x in range(n):
     
128.             x1 = x - half; x2 = x1 + window
     
129.             w1 = 0; w2 = window
     
130.             if x1 < 0: w1 = -x1; x1 = 0
     
131.             if x2 > n: w2 = n-x2; x2 = n
     
132.             yield signal[x1:x2], weight[w1:w2]
     
133. 
     
134.     def weighted(signal, weight):
     
135.         assert(len(signal) == len(weight))
     
136.         return sum(y*w for y,w in zip(signal, weight)) / sum(weight)
     
137. 
     
138.     def is_smooth(signal):
     
139.         for x in range(1, n):
     
140.             if signal[x] == signal[x-1]:
     
141.                 return False
     
142.         return True
     
143. 
     
144.     smoothed = signal
     
145.     for i in range(MAX_SMOOTH_ITERATION):
     
146.         smoothed = [weighted(s,w) for s,w in sliding(smoothed, window)]
     
147.         if is_smooth(smoothed): break
     
148. 
     
149.     return smoothed
     
150. 
     
151. def local_mean_and_envelope(signal, extrema):
     
152.     """根据极值点位置计算局部均值函数和局部包络函数."""
     
153.     n = len(signal); k = len(extrema); assert(1 < k <= n)
     
154.     # construct square signal
     
155.     mean = [ ]; enve = [ ]
     
156.     prev_mean = (signal[extrema[0]] + signal[extrema[1]]) / 2
     
157.     prev_enve = abs(signal[extrema[0]] - signal[extrema[1]]) / 2
     
158.     e = 1
     
159.     for x in range(n):
     
160.         if (x == extrema[e]) and (e + 1 < k):
     
161.             next_mean = (signal[extrema[e]] + signal[extrema[e+1]]) / 2
     
162.             mean.append((prev_mean + next_mean) / 2)
     
163.             prev_mean = next_mean
     
164.             next_enve = abs(signal[extrema[e]] - signal[extrema[e+1]]) / 2
     
165.             enve.append((prev_enve + next_enve) / 2)
     
166.             prev_enve = next_enve
     
167.             e += 1
     
168.         else:
     
169.             mean.append(prev_mean)
     
170.             enve.append(prev_enve)
     
171.     # smooth square signal
     
172.     window = max(extrema[i]-extrema[i-1] for i in range(1, k)) // 3
     
173.     return mean, moving_average_smooth(mean, window), \
     
174.            enve, moving_average_smooth(enve, window)
     
175. 
     
176. def extract_product_function(signal):
     
177.     s = signal; n = len(signal)
     
178.     envelopes = [ ] # 每次迭代得到的包络信号
     
179. 
     
180.     def component():
     
181.         c = [ ]
     
182.         for i in range(n):
     
183.             y = s[i]
     
184.             for e in envelopes:
     
185.                 y = y * e[i]
     
186.             c.append(y)
     
187.         return c
     
188. 
     
189.     for i in range(MAX_ENVELOPE_ITERATION):
     
190.         extrema = find_extrema(s)
     
191.         if len(extrema) <= 3: break
     
192.         m0, m, a0, a = local_mean_and_envelope(s, extrema)
     
193.         for y in a:
     
194.             if y <= 0: raise ValueError("invalid envelope signal")
     
195.         # 对原信号进行调制
     
196.         h = [y1-y2 for y1,y2 in zip(s, m)]
     
197.         t = [y1/y2 for y1,y2 in zip(h, a)]
     
198.         if DBG_FIG:
     
199.             DBG_FIG.plot(i, component(), s, m0, m, a0, a, h, t)
     
200.         # 得到纯调频信号时终止迭代
     
201.         err = sum(abs(1-y) for y in a) / n
     
202.         if err <= ENVELOPE_EPSILON: break
     
203.         # 调制信号收敛时终止迭代
     
204.         err = sum(abs(y1-y2) for y1,y2 in zip(s, t)) / n
     
205.         if err <= CONVERGENCE_EPSILON: break
     
206.         envelopes.append(a); s = t
     
207. 
     
208.     return component()
     
209. 
     
210. def lmd(signal):
     
211.     pf = [ ]
     
212.     # 每次迭代得到一个PF分量, 直到残余函数接近为一个单调函数为止
     
213.     residue = signal[:]
     
214.     while (len(pf) < MAX_NUM_PF) and \
     
215.           (not is_monotonous(residue)) and \
     
216.           (len(find_extrema(residue)) >= 5):
     
217.         component = extract_product_function(residue)
     
218.         residue = [y1-y2 for y1,y2 in zip(residue, component)]
     
219.         pf.append((component, residue))
     
220.     return pf
     

复制代码

sorry

谢谢您的分享，先收藏了

紫云轩8023

LZ，请问您有LMD算法的MATLAB代码么？谢谢！！！

犟牛

紫云轩8023 发表于 2015-10-23 13:44
LZ，请问您有LMD算法的MATLAB代码么？谢谢！！！

matlab版有人发过,可以搜索一下