声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2459|回复: 3

[其他] LMD算法Python程序

[复制链接]
发表于 2015-10-13 13:52 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 陈波 于 2015-10-13 13:52 编辑

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

欢迎大家提意见。

  1. """
  2. This module implements Local Mean Decomposition signal processing algorithm
  3. as described in [1].


  4. References:

  5. [1] Jonathan S. Smith. The local mean decomposition and its application to EEG
  6. perception data. Journal of the Royal Society Interface, 2005, 2(5):443-454
  7. """

  8. #
  9. # 算法控制参数
  10. #

  11. # 是否把信号的端点看作伪极值点
  12. INCLUDE_ENDPOINTS = True

  13. # 滑动平均算法: 最多迭代的次数
  14. MAX_SMOOTH_ITERATION = 12

  15. # 分离局部包络信号时最多迭代的次数
  16. MAX_ENVELOPE_ITERATION = 200

  17. ENVELOPE_EPSILON = 0.01
  18. CONVERGENCE_EPSILON = 0.01

  19. # 最多生成的积函数个数
  20. MAX_NUM_PF = 8

  21. # for debugging
  22. DBG_FIG = None

  23. def is_monotonous(signal):
  24.     """判断信号是否是(非严格的)单调序列"""

  25.     def is_monotonous_increase(signal):
  26.         y0 = signal[0]
  27.         for y1 in signal:
  28.             if y1 < y0:
  29.                 return False
  30.             y0 = y1
  31.         return True

  32.     def is_monotonous_decrease(signal):
  33.         y0 = signal[0]
  34.         for y1 in signal:
  35.             if y1 > y0:
  36.                 return False
  37.             y0 = y1
  38.         return True

  39.     if len(signal) <= 0:
  40.         return True
  41.     else:
  42.         return is_monotonous_increase(signal) \
  43.             or is_monotonous_decrease(signal)

  44. def find_extrema(signal):
  45.     """找出信号的所有局部极值点位置."""
  46.    
  47.     n = len(signal)

  48.     # 局部极值的位置列表
  49.     extrema = [ ]

  50.     lookfor = None
  51.     for x in range(1, n):
  52.         y1 = signal[x-1]; y2 = signal[x]
  53.         if lookfor == "min":
  54.             if y2 > y1:
  55.                 extrema.append(x-1)
  56.                 lookfor = "max"
  57.         elif lookfor == "max":
  58.             if y2 < y1:
  59.                 extrema.append(x-1)
  60.                 lookfor = "min"
  61.         else:
  62.             if y2 < y1:
  63.                 lookfor = "min"
  64.             elif y2 > y1:
  65.                 lookfor = "max"

  66.     # 优化(微调)极值点的位置:
  67.     # 当连续多个采样值相同时,取最中间位置作为极值点位置

  68.     def micro_adjust(x):
  69.         assert(0 <= x < n)
  70.         y = signal[x]; lo = hi = x
  71.         while (lo-1 >= 0) and (signal[lo-1] == y): lo -= 1
  72.         while (hi+1 < n) and (signal[hi+1] == y): hi += 1
  73.         if INCLUDE_ENDPOINTS:
  74.             if lo == 0: return 0
  75.             if hi == n-1: return n-1
  76.         return (lo + hi) // 2

  77.     extrema = [micro_adjust(x) for x in extrema]

  78.     if extrema and INCLUDE_ENDPOINTS:
  79.         if extrema[0] != 0:
  80.             extrema.insert(0, 0)
  81.         if extrema[-1] != n-1:
  82.             extrema.append(n-1)

  83.     return extrema

  84. def moving_average_smooth(signal, window):
  85.     """使用移动加权平均算法平滑方波信号."""

  86.     n = len(signal)

  87.     # at least one nearby sample is needed for average
  88.     if window < 3:
  89.         window = 3

  90.     # adjust length of sliding window to an odd number for symmetry
  91.     if (window % 2) == 0:
  92.         window += 1
  93.         
  94.     half = window // 2

  95.     weight = list(range(1, half+2)) + list(range(half, 0, -1))
  96.     assert(len(weight) == window)

  97.     def sliding(signal, window):
  98.         for x in range(n):
  99.             x1 = x - half; x2 = x1 + window
  100.             w1 = 0; w2 = window
  101.             if x1 < 0: w1 = -x1; x1 = 0
  102.             if x2 > n: w2 = n-x2; x2 = n
  103.             yield signal[x1:x2], weight[w1:w2]

  104.     def weighted(signal, weight):
  105.         assert(len(signal) == len(weight))
  106.         return sum(y*w for y,w in zip(signal, weight)) / sum(weight)

  107.     def is_smooth(signal):
  108.         for x in range(1, n):
  109.             if signal[x] == signal[x-1]:
  110.                 return False
  111.         return True

  112.     smoothed = signal
  113.     for i in range(MAX_SMOOTH_ITERATION):
  114.         smoothed = [weighted(s,w) for s,w in sliding(smoothed, window)]
  115.         if is_smooth(smoothed): break

  116.     return smoothed

  117. def local_mean_and_envelope(signal, extrema):
  118.     """根据极值点位置计算局部均值函数和局部包络函数."""
  119.     n = len(signal); k = len(extrema); assert(1 < k <= n)
  120.     # construct square signal
  121.     mean = [ ]; enve = [ ]
  122.     prev_mean = (signal[extrema[0]] + signal[extrema[1]]) / 2
  123.     prev_enve = abs(signal[extrema[0]] - signal[extrema[1]]) / 2
  124.     e = 1
  125.     for x in range(n):
  126.         if (x == extrema[e]) and (e + 1 < k):
  127.             next_mean = (signal[extrema[e]] + signal[extrema[e+1]]) / 2
  128.             mean.append((prev_mean + next_mean) / 2)
  129.             prev_mean = next_mean
  130.             next_enve = abs(signal[extrema[e]] - signal[extrema[e+1]]) / 2
  131.             enve.append((prev_enve + next_enve) / 2)
  132.             prev_enve = next_enve
  133.             e += 1
  134.         else:
  135.             mean.append(prev_mean)
  136.             enve.append(prev_enve)
  137.     # smooth square signal
  138.     window = max(extrema[i]-extrema[i-1] for i in range(1, k)) // 3
  139.     return mean, moving_average_smooth(mean, window), \
  140.            enve, moving_average_smooth(enve, window)

  141. def extract_product_function(signal):
  142.     s = signal; n = len(signal)
  143.     envelopes = [ ] # 每次迭代得到的包络信号

  144.     def component():
  145.         c = [ ]
  146.         for i in range(n):
  147.             y = s[i]
  148.             for e in envelopes:
  149.                 y = y * e[i]
  150.             c.append(y)
  151.         return c

  152.     for i in range(MAX_ENVELOPE_ITERATION):
  153.         extrema = find_extrema(s)
  154.         if len(extrema) <= 3: break
  155.         m0, m, a0, a = local_mean_and_envelope(s, extrema)
  156.         for y in a:
  157.             if y <= 0: raise ValueError("invalid envelope signal")
  158.         # 对原信号进行调制
  159.         h = [y1-y2 for y1,y2 in zip(s, m)]
  160.         t = [y1/y2 for y1,y2 in zip(h, a)]
  161.         if DBG_FIG:
  162.             DBG_FIG.plot(i, component(), s, m0, m, a0, a, h, t)
  163.         # 得到纯调频信号时终止迭代
  164.         err = sum(abs(1-y) for y in a) / n
  165.         if err <= ENVELOPE_EPSILON: break
  166.         # 调制信号收敛时终止迭代
  167.         err = sum(abs(y1-y2) for y1,y2 in zip(s, t)) / n
  168.         if err <= CONVERGENCE_EPSILON: break
  169.         envelopes.append(a); s = t

  170.     return component()

  171. def lmd(signal):
  172.     pf = [ ]
  173.     # 每次迭代得到一个PF分量, 直到残余函数接近为一个单调函数为止
  174.     residue = signal[:]
  175.     while (len(pf) < MAX_NUM_PF) and \
  176.           (not is_monotonous(residue)) and \
  177.           (len(find_extrema(residue)) >= 5):
  178.         component = extract_product_function(residue)
  179.         residue = [y1-y2 for y1,y2 in zip(residue, component)]
  180.         pf.append((component, residue))
  181.     return pf
复制代码



评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2015-10-21 14:59 | 显示全部楼层
谢谢您的分享,先收藏了
发表于 2015-10-23 13:44 | 显示全部楼层
LZ,请问您有LMD算法的MATLAB代码么?谢谢!!!

点评

matlab版有人发过,可以搜索一下  详情 回复 发表于 2015-10-23 14:56
发表于 2015-10-23 14:56 | 显示全部楼层
紫云轩8023 发表于 2015-10-23 13:44
LZ,请问您有LMD算法的MATLAB代码么?谢谢!!!

matlab版有人发过,可以搜索一下
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-12-23 21:12 , Processed in 0.064056 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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