Python/SciPy 中的峰值查找算法

我可以自己写一些东西,通过找到一阶导数的过零点之类的,但它似乎是一个足够普通的函数,可以包含在标准库中。有人知道吗?

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

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

一维峰值:

FFT output with peaks

二维峰值:

Radon transform output with circled peak

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

通常,您只关心一些强峰,所以选择它们要么是因为它们高于某个阈值,要么是因为它们是按振幅排序的有序列表中的第一个 N峰。

正如我所说,我知道如何自己写这样的东西。我只是想知道是否有一个预先存在的函数或包可以很好地工作。

更新:

翻译了一个 MATLAB 脚本和它的工作体面的一维情况下,但可以更好。

更新:

一维情况下的 创造了一个更好的版本

207409 次浏览

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

插值的样条插值和平滑效果非常好,可能有助于拟合峰值,然后找到峰值的位置。

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

以可靠的方式检测频谱中的峰值已经被研究了相当多,例如80年代所有的音乐/音频信号的正弦建模工作。在文献中寻找“正弦建模”。

如果你的信号像这个例子一样清晰,一个简单的“给我一个振幅比 N 个邻居高的东西”应该会相当不错。如果你有嘈杂的信号,一个简单而有效的方法是及时查看你的峰值,跟踪它们: 然后你检测谱线而不是谱峰。在信号的滑动窗口上计算 FFT,得到一组时间频谱(也称为光谱图)。然后观察光谱峰在时间上的演变(即在连续的窗口中)。

我正在研究一个类似的问题,我发现一些最好的参考来自化学(来自质谱数据中的峰值发现)。为了对峰值查找算法进行更好的全面回顾,请阅读 这个。这是我遇到过的最清晰的峰值寻找技术评论之一。(小波是在嘈杂数据中寻找此类峰值的最佳方法。).

看起来你们的峰值清晰可见而且没有隐藏在噪音中。在这种情况下,我建议使用平滑的 Savitiky-Golay 衍生工具来找到峰值(如果你只是区分上面的数据,你会得到一大堆假阳性。).这是一种非常有效的技术,并且非常容易实现(您确实需要一个矩阵类 w/basic 操作)。如果你只是找到第一个 S-G 导数的零交叉点,我想你会很高兴的。

有一个功能的科学命名为 scipy.signal.find_peaks_cwt,听起来是适合你的需要,但我没有经验,所以我不能推荐。.

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

对于那些不确定在 Python 中使用哪种峰值查找算法的人,这里快速概述一下替代算法: https://github.com/MonsieurV/py-findpeaks

为了让自己与 MatLab findpeaks函数等价,我发现 Marcos Duarte 的 [检测] _ 峰函数是一个不错的选择。

非常容易使用:

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

这会给你:

detect_peaks results

首先,如果没有进一步的说明,“峰值”的定义是模糊的。例如,对于接下来的系列,您会调用5-4-5一个峰值还是两个峰值?

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

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

峰值检测是极值理论文献中一个研究得很好的课题,也被称为“极值的去序”。它的典型应用包括根据对环境变量的连续读数识别危险事件,例如分析风速以探测风暴事件。

函数 scipy.signal.find_peaks,顾名思义,对此很有用。但是了解其参数 widththresholddistance尤其是 prominence对于获得较好的峰提取效果至关重要。

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

什么是 (地形学)日珥? 它是 “从山顶到达更高地形所需的最低高度”,正如我们在这里看到的:

enter image description here

我们的想法是:

日珥越高,峰值越“重要”。

测试:

enter image description here

我故意使用了一个(噪声)频率变化的正弦波,因为它显示了许多困难。我们可以看到,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()

为了同时检测正峰和负峰,峰值检测是有帮助的。

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

PeakDetection

正如在这个 呼叫的底部提到的,没有一个通用的峰值定义。因此,找到峰值的通用 算法如果不带来额外的假设(条件、参数等)就不能工作。本页提供了一些最 脱光了的建议。以上答案中列出的所有文献或多或少都是一种迂回的方式来做同样的事情,所以请随意选择。

在任何情况下,您都有责任根据您的经验和所涉及的光谱(曲线)属性(噪音、采样、带宽等) ,缩小特性需要具备的属性范围,以便将其归类为峰值