Python/SciPy 的寻峰算法

Question

我可以通过查找一阶导数的过零或其他东西自己编写一些东西，但这似乎是一个足够通用的函数，可以包含在标准库中。有人知道吗？

我的特定应用是二维数组，但通常它会用于在 FFT 等中查找峰值。

具体来说，在这类问题中，有多个强峰值，然后是许多较小的“峰值”，这些“峰值”只是由应该忽略的噪声引起的。这些只是例子；不是我的实际数据：

一维峰：

二维峰：

寻峰算法会找到这些峰值的位置（不仅仅是它们的值），并且理想情况下会找到真正的样本间峰值，而不仅仅是具有最大值的索引，可能使用quadratic interpolation 或其他东西。

通常您只关心几个强峰，因此选择它们要么是因为它们高于某个阈值，要么是因为它们是有序列表的前 n 个峰，按幅度排序。

正如我所说，我知道如何自己写这样的东西。我只是想问一下，是否有一个预先存在的已知运行良好的函数或包。

更新：

我translated a MATLAB script，它适用于一维情况，但可能会更好。

更新更新：

sixtenbe created a better version 用于一维情况。

Answer 1

我不认为您正在寻找的东西是由 SciPy 提供的。在这种情况下，我会自己编写代码。

scipy.interpolate 的样条插值和平滑非常好，可能有助于拟合峰值，然后找到它们的最大值。

Answer 2

有一些标准的统计函数和方法可以找到数据的异常值，这可能是您在第一种情况下所需要的。使用导数可以解决您的第二个问题。但是，我不确定是否有一种方法可以同时解决连续函数和采样数据。

Answer 3

已经对以可靠方式检测频谱中的峰值进行了大量研究，例如在 80 年代对音乐/音频信号进行正弦建模的所有工作。在文献中查找“Sinusoidal Modeling”。

如果您的信号与示例一样干净，那么简单的“给我一些幅度高于 N 个邻居的信号”应该可以很好地工作。如果您有噪声信号，一个简单但有效的方法是及时查看峰值并跟踪它们：然后您检测谱线而不是谱峰。 IOW，您在信号的滑动窗口上计算 FFT，以便及时获得一组频谱（也称为频谱图）。然后，您可以查看光谱峰值在时间上的演变（即在连续窗口中）。

Answer 4

我正在研究一个类似的问题，我发现一些最好的参考来自化学（来自质谱数据中的峰）。要全面了解峰值查找算法，请阅读this。这是我遇到过的关于峰值发现技术的最清晰的评论之一。 （小波最适合在噪声数据中找到此类峰值。）。

您的峰值看起来很清晰，没有隐藏在噪音中。在这种情况下，我建议使用平滑的 savtizky-golay 导数来找到峰值（如果你只是区分上面的数据，你会有一堆误报。）。这是一种非常有效的技术，并且很容易实现（您确实需要一个带有基本操作的矩阵类）。如果你只是找到第一个 S-G 导数的过零，我想你会很高兴的。

Answer 5

scipy 中有一个名为scipy.signal.find_peaks_cwt 的函数，听起来很适合您的需求，但是我没有这方面的经验，所以我不能推荐..

http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks_cwt.html

Answer 6

对于那些不确定在 Python 中使用哪种寻峰算法的人，这里是替代方案的快速概述：https://github.com/MonsieurV/py-findpeaks

希望自己与 MatLab findpeaks 函数等效，我发现 Marcos Duarte 的 detect_peaks function 是一个不错的选择。

非常容易使用：

import numpy as np
from vector import vector, plot_peaks
from libs import detect_peaks
print('Detect peaks with minimum height and distance filters.')
indexes = detect_peaks.detect_peaks(vector, mph=7, mpd=2)
print('Peaks are: %s' % (indexes))

这会给你：

Answer 7

首先，如果没有进一步的说明，“峰值”的定义是模糊的。例如，对于以下系列，您将 5-4-5 称为一峰还是二峰？

1-2-1-2-1-1-5-4-5-1-1-5-1

在这种情况下，您至少需要两个阈值：1) 一个高阈值，只有高于该阈值才能将极值记录为峰值； 2) 一个低阈值，这样由低于它的小值分隔的极值将成为两个峰值。

峰值检测是极值理论文献中深入研究的主题，也称为“极值去簇”。其典型应用包括根据环境变量的连续读数识别危险事件，例如分析风速以检测风暴事件。

Answer 8

函数scipy.signal.find_peaks，顾名思义，对此很有用。但重要的是要很好地理解它的参数width、threshold、distance以及最重要的是prominence，以获得良好的峰值提取。

根据我的测试和文档，prominence 的概念是“有用的概念”，可以保留好的峰值，并丢弃嘈杂的峰值。

(topographic) prominence 是什么？这是“从山顶下降到任何更高地形所需的最低高度”，如下所示：

想法是：

突出度越高，峰越“重要”。

测试：

我故意使用（嘈杂的）频率变化正弦曲线，因为它显示出很多困难。我们可以看到width 参数在这里不是很有用，因为如果您将最小width 设置得太高，那么它将无法在高频部分跟踪非常接近的峰值。如果您将width 设置得太低，信号左侧会出现许多不需要的峰值。 distance 也有同样的问题。 threshold 只与直接邻居比较，这里没有用。 prominence 是提供最佳解决方案的那个。请注意，您可以组合其中的许多参数！

代码：

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

x = np.sin(2*np.pi*(2**np.linspace(2,10,1000))*np.arange(1000)/48000) + np.random.normal(0, 1, 1000) * 0.15
peaks, _ = find_peaks(x, distance=20)
peaks2, _ = find_peaks(x, prominence=1)      # BEST!
peaks3, _ = find_peaks(x, width=20)
peaks4, _ = find_peaks(x, threshold=0.4)     # Required vertical distance to its direct neighbouring samples, pretty useless
plt.subplot(2, 2, 1)
plt.plot(peaks, x[peaks], "xr"); plt.plot(x); plt.legend(['distance'])
plt.subplot(2, 2, 2)
plt.plot(peaks2, x[peaks2], "ob"); plt.plot(x); plt.legend(['prominence'])
plt.subplot(2, 2, 3)
plt.plot(peaks3, x[peaks3], "vg"); plt.plot(x); plt.legend(['width'])
plt.subplot(2, 2, 4)
plt.plot(peaks4, x[peaks4], "xk"); plt.plot(x); plt.legend(['threshold'])
plt.show()

Answer 9

要检测正峰和负峰，PeakDetect 很有帮助。

from peakdetect import peakdetect

peaks = peakdetect(data, lookahead=20) 
# Lookahead is the distance to look ahead from a peak to determine if it is the actual peak. 
# Change lookahead as necessary 
higherPeaks = np.array(peaks[0])
lowerPeaks = np.array(peaks[1])
plt.plot(data)
plt.plot(higherPeaks[:,0], higherPeaks[:,1], 'ro')
plt.plot(lowerPeaks[:,0], lowerPeaks[:,1], 'ko')