Python/SciPy的找峰算法

Question

我可以自己写一些东西，找到一阶导数的零点，但它似乎是一个常见的函数，足以包含在标准库中。有人听说过吗？

我的特殊应用程序是一个2D数组，但通常用于查找FFTs中的峰值，等等。

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

一维峰：

​

​

二维峰：

​

​

找峰算法将找到这些峰值的位置(不仅仅是它们的值)，并且理想情况下会找到真正的样本间峰值，而不仅仅是具有最大值的索引，可能会使用二次插值或其他方法。

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

就像我说的，我自己知道怎么写这样的东西。我只是想问一下，是否有一个预先存在的函数或程序包可以正常工作。

更新：

I 翻译MATLAB脚本，它很适合一维的情况，但可能更好。

更新：

1-D病例的6张创建了一个更好的版本。

Answer 1

函数scipy.signal.find_peaks (顾名思义)对此非常有用。但是，了解其参数width、threshold、distance 以及最重要的，才能得到一个很好的提取峰。

根据我的测试和文档，日珥的概念是“有用的概念”，以保持好的峰值，并丢弃噪声的峰值。

(地形)日珥是什么？它是“从山顶下降到任何较高地形所需的最低高度”，如下所示：

​

​

这样做的目的是：

日珥越高，峰越“重要”。

测试：

​

​

我故意使用(噪音)频率变化的正弦波，因为它显示出许多困难.我们可以看到，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 2

我正在研究一个类似的问题，我发现一些最好的参考资料来自化学(从大量的数据中发现峰值)。要对峰值查找算法进行彻底的检查，请阅读这。这是我所见过的最清晰的找峰技术之一。(小波是在含噪数据中找到这类峰值的最佳方法。)

看起来你的山峰是清晰的，没有隐藏在噪音中。在这种情况下，我建议使用光滑的、精明的-golay导数来找出峰值(如果您只是区分上面的数据，就会发现大量的假阳性。)这是一种非常有效的技术，很容易实现(您确实需要一个矩阵类w/ basic操作)。如果你简单地找到第一个S导数的过零点，我想你会很高兴的。

Answer 3

这里有一个叫做scipy.signal.find_peaks_cwt的函数，听起来很适合你的需要，但是我没有经验，所以我不能推荐。

cwt.html