【问题标题】:Counting Peaks in a Time Series计算时间序列中的峰值
【发布时间】:2019-07-31 14:00:29
【问题描述】:

我正在计算一个 numpy 数组中的波峰和波谷的数量。

我有一个像这样的 numpy 数组:

stack = np.array([0,0,5,4,1,1,1,5,1,1,5,1,1,1,5,1,1,5,1,1,5,1,1,5,1,1,5,1,1])

绘制,这个数据看起来像这样:

我正在寻找这个时间序列中的峰值数量:

这是我的代码,它适用于像这样在时间序列表示中有明显波峰和波谷的示例。我的代码返回找到峰值的数组的索引。

#example
import numpy as np
from scipy.signal import argrelextrema

stack = 
np.array([0,0,5,4,1,1,1,5,1,1,5,1,1,1,5,1,1,5,1,1,5,1,1,5,1,1,5,1,1])

# for local maxima
y = argrelextrema(stack, np.greater)

print(y)

结果:

(array([ 2,  7, 10, 14, 17, 20, 23, 26]),)

已经找到了8个清晰的峰,可以正确计数。

我的解决方案似乎不适用于不太清晰和更混乱的数据。

下面的一个数组不好用,找不到我需要的峰:

array([ 0.        ,  5.70371806,  5.21210157,  3.71144767,  3.9020162 ,
    3.87735984,  3.89030171,  6.00879918,  4.91964227,  4.37756275,
    4.03048542,  4.26943028,  4.02080471,  7.54749062,  3.9150576 ,
    4.08933851,  4.01794766,  4.13217794,  4.15081972,  8.11213474,
    4.6561735 ,  4.54128693,  3.63831552,  4.3415324 ,  4.15944019,
    8.55171441,  4.86579459,  4.13221943,  4.487663  ,  3.95297979,
    4.35334706,  9.91524674,  4.44738182,  4.32562141,  4.420753  ,
    3.54525697,  4.07070637,  9.21055852,  4.87767969,  4.04429321,
    4.50863677,  3.38154581,  3.73663523,  3.83690315,  6.95321174,
    5.11325128,  4.50351938,  4.38070175,  3.20891173,  3.51142661,
    7.80429569,  3.98677631,  3.89820773,  4.15614576,  3.47369797,
    3.73355768,  8.85240649,  6.0876192 ,  3.57292324,  4.43599135,
    3.77887259,  3.62302175,  7.03985076,  4.91916556,  4.22246518,
    3.48080777,  3.26199699,  2.89680969,  3.19251448])

绘制,这个数据看起来像:

同样的代码返回:

(array([ 1,  4,  7, 11, 13, 15, 19, 23, 25, 28, 31, 34, 37, 40, 44, 50, 53,
   56, 59, 62]),)

此输出错误地将数据点计为峰值。

理想输出

理想的输出应该返回清晰峰的数量,在这种情况下为 11,它们位于索引处:

[1,7,13,19,25,31,37,44,50,56,62]

我相信我的问题是由于 argrelextrema 函数的聚合性质而出现的。

【问题讨论】:

  • 这取决于您如何定义“明确定义的峰值”。从表面上看,这似乎很简单,但您希望在程序中使用什么样的规则?
  • 这可能与Peak-finding algorithm for Python/SciPy 重复。请注意,此解决方案确实提供了 prominence 输入,可让您指定要计算的峰值的突出程度。
  • 看起来您需要定义一些阈值并搜索大于它的值。

标签: python arrays scipy


【解决方案1】:

看起来argrelextrema 可以帮助您完成大部分工作。它具有您想要的所有峰,但也有一些额外的峰。您需要提出一个适合您情况的标准并过滤掉您不想要的峰。

例如,如果您不希望峰值小于 5,您可以这样做:

In [17]: result = argrelextrema(a, np.greater)                                                           

In [18]: result                                                                                          
Out[18]: 
(array([ 1,  4,  7, 11, 13, 15, 19, 23, 25, 28, 31, 34, 37, 40, 44, 50, 53,
        56, 59, 62]),)

In [19]: result[0][a[result[0]] > 5]                                                                     
Out[19]: array([ 1,  7, 13, 19, 25, 31, 37, 44, 50, 56, 62])

【讨论】:

    【解决方案2】:

    您可以使用一些阈值来找到峰值:

    prev = stack[0] or 0.001
    threshold = 0.5
    peaks = []
    
    for num, i in enumerate(stack[1:], 1):
        if (i - prev) / prev > threshold:
            peaks.append(num)
        prev = i or 0.001
    
    print(peaks)
    # [1, 7, 13, 19, 25, 31, 37, 44, 50, 56, 62]
    

    【讨论】:

    • 这个答案似乎运作良好。谢谢! 0.001 代表什么?您是如何选择 0.5 的阈值的?
    • @Murray 不客气。 i or 0.001 解决了除以零的问题。如果 i 为 0,则此表达式的计算结果为 0.001(您也可以使用 0.1,具体取决于您的数据)。我刚刚取了第一个值 0.5,它起作用了。
    【解决方案3】:

    您应该在 scipy.signal 模块中尝试 find_peaks

    数据示例

    import numpy as np
    from scipy.signal import find_peaks
    import matplotlib.pyplot as plt
    
    a = np.array([ 0.        ,  5.70371806,  5.21210157,  3.71144767, 3.9020162 ,  3.87735984,  3.89030171,  6.00879918,  4.91964227,  4.37756275,
    4.03048542,  4.26943028,  4.02080471,  7.54749062,  3.9150576 ,
    4.08933851,  4.01794766,  4.13217794,  4.15081972,  8.11213474,
    4.6561735 ,  4.54128693,  3.63831552,  4.3415324 ,  4.15944019,
    8.55171441,  4.86579459,  4.13221943,  4.487663  ,  3.95297979,
    4.35334706,  9.91524674,  4.44738182,  4.32562141,  4.420753  ,
    3.54525697,  4.07070637,  9.21055852,  4.87767969,  4.04429321,
    4.50863677,  3.38154581,  3.73663523,  3.83690315,  6.95321174,
    5.11325128,  4.50351938,  4.38070175,  3.20891173,  3.51142661,
    7.80429569,  3.98677631,  3.89820773,  4.15614576,  3.47369797,
    3.73355768,  8.85240649,  6.0876192 ,  3.57292324,  4.43599135,
    3.77887259,  3.62302175,  7.03985076,  4.91916556,  4.22246518,
    3.48080777,  3.26199699,  2.89680969,  3.19251448])
    
    # Here you should fine tune parameters to get what you want
    peaks = find_peaks(a, prominence=1.5)
    
    # Plotting
    plt.plot(a)
    plt.title("Finding Peaks")
    
    [plt.axvline(p, c='C3', linewidth=0.3) for p in peaks[0]]
    
    plt.show()
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2012-08-28
      • 1970-01-01
      • 1970-01-01
      • 2020-06-16
      • 2020-10-08
      • 1970-01-01
      • 1970-01-01
      • 2018-06-08
      相关资源
      最近更新 更多