使用 Python 计算图中的波数答案

【问题标题】：Count Waves in Plot Using Python使用 Python 计算图中的波数
【发布时间】：2020-11-24 04:08:06
【问题描述】：

我刚刚生成了一个如下所示的图；我试图计算这张图表中有多少个全波。如注释所述，应该有 7 个全波。我试图平滑图表，但最终无处可去。我还在 scipy 中研究了 find_peaks，但它似乎不适用于这种情况，因为一个波可能有多个峰值。

编辑：

这个问题来自现实中的一个问题：搞清楚图中货架的数量（见下图）：

阈值处理后，我得到以下数据文件：

原文：file-1

行平均值：file-2

file-2 的结果是开头所示的绘图。

【问题讨论】：

你如何定义“全波”？高于然后低于某个阈值？
@StackTracer 这是棘手的部分。我在图片中标记了 7 个全波。
好像你可以将值设置为1，低于平均值，0高于平均值，然后使用find_peaks
@lovechillcool 是的，在尝试分析真实数据时，这往往是主要问题。
这样的峰很难平滑，但也许你可以使用 Savitzky–Golay 过滤器（来自 scipy.signal import savgol_filter）。它易于使用。

标签： python scipy matlab-figure smoothing

【解决方案1】：

我认为一般来说这是一个很难回答的问题，因为有些数据集可以解决方案，而有些数据集则不能。你的数据是倒置的，所以第一步是将y翻转到-y，所以最小值将被解释为最大值（并且可能取绝对值以避免处理负数。）

第一个选项是使用scipy.signal.find_peaks。在这里了解您的数据，您可以利用一些参数：根据我的经验，高度、距离和突出度是最有用的。 There is a nice explanation关于find_peaks的参数。在大多数情况下，这将正确识别峰值，但需要时间来适当地设置参数。

与scipy.signal.find_peaks_cwt 类似的解决方案，在这里（在大多数情况下）您需要调整宽度参数，即（来自文档）：

用于计算 CWT 矩阵的一维宽度数组。一般来说，此范围应涵盖目标峰的预期宽度。

但同样，这需要对您的数据有一些先验知识。

因为你有周期性数据，也许你可以利用FFT找到特征频率来调整find_peaks和find_peaks_cwt里面的参数。由于您没有提供数据集，因此我只有合成数据要处理。请注意，我返回len(peaks) - 1，因为通常在边界上会有一个额外的时间段被计算在内。

import numpy as np
from scipy.signal import find_peaks, find_peaks_cwt
import matplotlib.pyplot as plt

# some generic data
x = np.linspace(0, 1000, 10000)
y = 250 + 100 * np.sin(0.08 * x) - np.random.normal(30, 20, 10000)

def count_waves_1(x, y):
    peaks, props = find_peaks(y, prominence=120, height= np.max(y) / 10, distance=200)
    # here you can make use of props to filter the peaks by different properties,
    # for example extract only the n largest prominence peak:
    #
    # ind = np.argpartition(props["prominences"], -n_largest)[-n_largest:]
    # peaks = peaks[ind]

    plt.plot(x, y)
    plt.plot(x[peaks], y[peaks], 'ro')
    return len(peaks) - 1

first_solution = count_waves_1(x, y)


def count_waves_2(x, y):
    peaks = find_peaks_cwt(y, widths=np.arange(100, 200))
    plt.plot(x, y)
    plt.plot(x[peaks], y[peaks], 'ro')
    return len(peaks) - 1

second_solution = count_waves_2(x, y)

print(first_solution, second_solution)

【讨论】：