Python/SciPy的寻峰算法

Question

我可以通过寻找一阶导数的零交叉点等方法自己编写算法，但这似乎是一个常见的函数应该包含在标准库中。 有人知道这个吗？

我的特定应用程序是一个二维数组，但通常它将用于查找FFT等中的峰值。

具体来说，在这些问题的情况下，存在多个强峰，然后有许多较小的“峰”，这些峰仅由应忽略的噪声引起。 这些只是示例； 不是我的实际数据：

1维峰值：

2维峰值：

峰值查找算法将找到这些峰值的位置（不仅仅是其价值），理想情况下将找到真正的样本间隔峰值，而不仅仅是具有最大值的索引，可能使用二次插值等方法。

通常，您只关心几个强峰，因此它们将被选择，因为它们在某个阈值以上，或者因为它们是一个有序列表的前n个峰，按幅度排序。

正如我所说，我知道如何自己编写这样的代码。 我只是想问是否有已知效果良好的预先存在的函数或软件包。

更新：

我翻译了MATLAB脚本，对于1-D情况可以使用，但可以更好。

更新更新：

sixtenbe创建了更好的版本 用于1-D情况。

- endolith

---

@endolith，你有翻译成Python的这个MATLAB文件的原始文件吗？谢谢！ - Spacey

@Mohammad：http://billauer.co.il/peakdet.html https://gist.github.com/250860#file_peakdet.m - endolith

这个怎么样：http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks_cwt.html - dashesy

这个回答在这里值得注意。 - Gulzar

Answer 1

函数 scipy.signal.find_peaks 的作用是查找峰值，但是理解其参数 width、threshold、distance 和最重要的 prominence 对于获得良好的峰值提取非常重要。

根据我的测试和文档，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()

- Basj

---

这就是我想要的。但是你是否知道有没有在二维数组中找到突出点的实现？ - Jason

@Jason 我刚刚看到了在2D数组中检测峰值，值得一读！ - Basj

@jpmorr，您能否发布一个可运行的示例，展示wlen + prominence？这将非常有用！ - Basj

文档（或其他几个解释）没有很清楚地说明“prominence”这个参数的含义，这真的很烦人。它似乎是最重要的参数。感谢您的解释。 - SuperCodeBrah

@thentangler，突出性与下降的高度有关，因此通过观察信号的振幅（这里从-1到1），您可以使用类似于“1”的正确数量级。 - Basj

Answer 2

我遇到了类似的问题，并且发现一些最好的参考资料来自化学领域（针对质谱数据中的峰值寻找）。如果想要了解峰值寻找算法的全面评估，请参考这篇文章。这是我所见过的最好、最清晰的峰值寻找技术综述之一（小波变换在嘈杂数据中寻找此类峰值最为有效）。

看起来你的峰值非常明显，不会被噪声掩盖。在这种情况下，我建议使用平滑的Savitzky-Golay导数来寻找峰值（如果仅仅对数据进行微分，会得到许多错误的结果）。这是一种非常有效的技术，并且很容易实现（你需要一个具有基本操作的矩阵类）。如果你只是寻找第一阶段S-G导数的零交叉点，我相信你会很满意。

- Paul

---

我正在寻找一个通用的解决方案，而不是只适用于特定图像的解决方案。我将一个MATLAB脚本改编成Python，并且它运行得相当不错。 - endolith

为什么导数的零点比仅测试三个点中间的一个是另外两个中较大还是较小更好？我已经应用了SG变换，似乎是额外的成本。 - kirill_igum

Answer 3

scipy中有一个名为scipy.signal.find_peaks_cwt的函数，听起来很适合你的需求，但我没有使用过，所以无法推荐。

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

- Hanan Shteingart

Answer 4

如果对于在Python中使用哪种峰值查找算法还不确定，这里提供一份快速概述备选方案：https://github.com/MonsieurV/py-findpeaks

我自己想找一个和MatLab findpeaks 函数相等的东西，后来发现Marcos Duarte的detect_peaks函数挺不错的。

很容易使用：

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))

这将为您提供：

- Yoan Tournade

---

自从这篇文章写出来以后，scipy 添加了 find_peaks 函数。 - onewhaleid

Answer 5

为了检测正负峰值，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')

- Timoth Dev A

---

谢谢。这是最简单和最有效的解决方案。它轻松地找到了我的数据中的峰值和谷值。 - StuckInPhDNoMore

Answer 6

可靠地检测频谱中的峰值已经得到了相当多的研究，例如80年代对于音乐/音频信号的正弦建模的所有工作。在文献中寻找“正弦建模”。如果您的信号像示例一样干净，那么一个简单的“给我一个振幅比N个邻居高的东西”应该能够相当好地工作。如果您有嘈杂的信号，则一种简单但有效的方法是查看时间上的峰值并跟踪它们：然后您会检测到频谱线而不是频谱峰值。换句话说，您在信号的滑动窗口上计算FFT，以获取一组随时间变化的频谱（也称为谱图）。然后，您查看频谱峰值随时间的演变（即在连续的窗口中）。

- David Cournapeau

---

查看时间峰值？检测光谱线？我不确定这是什么意思。这对方波有效吗？ - endolith

哦，你说的是使用STFT而不是FFT。这个问题并不特别涉及FFT；那只是一个例子。它是关于在任何一般的1D或2D数组中找到峰值的问题。 - endolith

Answer 7

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

- nullpointer

Answer 8

我认为 SciPy 并没有提供你所需要的内容。在这种情况下，我会自己编写代码。

scipy.interpolate 中的样条插值和平滑非常好，可能对拟合峰值并找到其最大值的位置非常有帮助。

- Eric O. Lebigot

---

抱歉，我认为这应该是一条评论，而不是答案。它只是建议自己写代码，并提出了一些可能有用的功能（顺便提一句，Paul答案中提到的更相关）。 - Ami Tavory

Answer 9

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

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

在这种情况下，您需要至少两个阈值：1）高阈值仅在其上方才能将极值记录为峰值；2）低阈值，以便在其下面由小值分隔的极值将成为两个峰值。

峰值检测是极值理论文献中研究得比较充分的主题，也被称为“极值去聚类”。其典型应用包括基于环境变量的连续读数识别危险事件，例如分析风速以检测风暴事件。

- Ian Liu

Answer 10

如页面底部所述，峰值没有普遍的定义。因此，一个能够找到峰值的通用算法不能在不引入其他假设（条件、参数等）的情况下工作。本页面提供了一些最简化的建议。上面回答中列出的所有文献都是以更或多或少迂回的方式完成相同的任务，因此请随意选择。

无论如何，根据您的经验和所涉及光谱（曲线）的特性（噪声、采样、带宽等），您有责任缩小特征需要具备的属性范围，以便将其归类为峰值。

- Dusan Kojic