Python/SciPy的峰值查找算法

Question

我可以通过查找一阶导数或其他东西的零交叉来自己编写一些东西,但它似乎是一个通用的函数,可以包含在标准库中.谁知道一个？

我的特定应用是2D阵列,但通常它将用于在FFT等中查找峰值.

具体而言,在这些类型的问题中,存在多个强峰,然后是许多较小的"峰值",这些"峰值"仅由应该忽略的噪声引起.这只是例子; 不是我的实际数据:

一维峰值:

二维峰值:

峰值寻找算法将找到这些峰值的位置(不仅仅是它们的值),理想情况下会找到真正的样本间峰值,而不仅仅是具有最大值的索引,可能使用二次插值等.

通常,您只关心一些强峰,因此它们要么被选中,要么是因为它们高于某个阈值,要么是因为它们是有序列表的前n个峰值,按振幅排列.

正如我所说,我知道如何自己写这样的东西.我只是问是否有一个已知的功能或包已知可以正常工作.

更新:

我翻译了一个MATLAB脚本,它适用于1-D案例,但可能更好.

更新更新:

sixtenbe 为1-D案例创造了更好的版本.

Answer 1

scipy.signal.find_peaks顾名思义,该功能对此非常有用.但是,要理解以及它的参数是非常重要的width,threshold,distance 和高于一切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 2

我正在研究一个类似的问题,我发现一些最好的参考资料来自化学(来自质谱数据中的峰值).有关峰值查找算法的详细审查,请阅读此内容.这是我遇到的最佳峰值发现技术最清晰的评论之一.(小波最适合在噪声数据中找到此类峰值.).

看起来你的峰值是清晰的,并没有隐藏在噪音中.在这种情况下,我建议使用平滑的savtizky-golay衍生物来找到峰值(如果你只是区分上面的数据,那么你就会有一堆误报.).这是一种非常有效的技术,并且很容易实现(你需要一个矩阵类w /基本操作).如果您只是找到第一个SG衍生物的过零点,我想您会很高兴.

Answer 3

在scipy中有一个功能scipy.signal.find_peaks_cwt,听起来像是适合您的需要,但我没有经验,所以我不能推荐..

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

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

哪个会给你:

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

Answer 6

已经对可靠地检测频谱中的峰值进行了相当多的研究,例如80年代对音乐/音频信号的正弦建模的所有工作.在文献中寻找"正弦模型".

如果您的信号与示例一样干净,那么简单的"给我一个振幅高于N个邻居的东西"应该可以很好地工作.如果您有嘈杂的信号,一种简单但有效的方法是查看您的峰值,跟踪它们:然后检测谱线而不是谱峰.IOW,您可以在信号的滑动窗口上计算FFT,以获得一组光谱(也称为频谱图).然后,您可以看到光谱峰值的演变(即在连续的窗口中).