Python/SciPy的峰值查找算法

Question

我可以自己写一些代码，通过找到一阶导数的零交叉点来实现，但这似乎是一个常见的功能，应该在标准库中有提供。有人知道吗？

我具体的应用是处理一个二维数组，但通常这种功能是用来在快速傅里叶变换（FFT）中寻找峰值等。

具体来说，在这类问题中，通常会有多个强峰值，还有很多小的“峰值”是由噪声造成的，这些噪声应该被忽略。这些只是示例，不是我的实际数据：

一维峰值：

二维峰值：

这个寻找峰值的算法会找到这些峰值的位置（不仅仅是它们的数值），理想情况下还会找到真正的插值峰值，而不仅仅是最大值的索引，可能会用到二次插值之类的方法。

通常你只关心几个强峰值，所以它们要么是因为超过了某个阈值而被选中，要么是因为它们是按幅度排序的前n个峰值。

正如我所说，我知道怎么自己写这样的东西。我只是想问一下，是否有现成的函数或包是公认的效果不错。

更新：

我翻译了一个MATLAB脚本，在一维情况下效果还不错，但可以更好。

更新的更新：

sixtenbe创建了一个更好的版本，适用于一维情况。

Answer 1

在scipy这个库里，有一个叫做 scipy.signal.find_peaks_cwt 的函数，听起来很适合你的需求。不过我自己没有用过这个函数，所以不能给你推荐。

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

Answer 2

我正在研究一个类似的问题，发现一些最好的参考资料来自化学，特别是关于质谱数据中峰值的寻找。想要深入了解峰值寻找算法，可以看看这篇文章。这是我见过的关于峰值寻找技术最清晰的综述之一。（在嘈杂数据中，使用小波变换是寻找这类峰值的最佳方法。）

看起来你的峰值非常明显，并没有被噪声掩盖。在这种情况下，我建议使用平滑的Savitzky-Golay导数来寻找峰值（如果你只是对上面的数据进行微分，结果会一团糟，出现很多错误的结果）。这是一种非常有效的技术，而且实现起来也相对简单（你需要一个支持基本操作的矩阵类）。如果你只找到第一个Savitzky-Golay导数的零交叉点，我想你会满意的。

Answer 3

函数 scipy.signal.find_peaks，顾名思义，是用来找峰值的。不过，想要准确提取峰值，了解它的一些参数，比如 width（宽度）、threshold（阈值）、distance（距离），尤其是 prominence（显著性）是非常重要的。

根据我的测试和文档，显著性这个概念是“最有用的”，它能帮助我们保留真正的峰值，而去掉那些噪声峰值。

什么是 (地形)显著性？它是指“从山顶下降到任何更高地形所需的最小高度”，可以参考下面的图：

这个概念的意思是：

显著性越高，峰值就越“重要”。

测试结果：

我故意使用了一个（带噪声的）频率变化的正弦波，因为它展示了很多问题。我们可以看到，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()