Python/SciPy的峰值查找算法

23 浏览
0 Comments

Python/SciPy的峰值查找算法

我可以通过找到一阶导数的零点或其他方法自己编写代码，但这似乎是一个常见的函数，应该包含在标准库中。大家有知道的吗？

我的特定应用是一个2D数组，但通常它会用于在FFT等中找到峰值。

具体来说，在这些问题中，有多个强峰，然后有很多小的“峰”，只是由噪声引起的，应该被忽略。这些只是示例；不是我的实际数据：

一维峰：

\"FFT

二维峰：

\"Radon

峰值查找算法会发现这些峰的位置（而不仅仅是它们的值），理想情况下会找到真实的样本峰值，而不仅仅是使用最大值的索引，可能是使用二次插值等方法。

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

正如我所说的，我知道如何自己编写这样的代码。我只是想问一下是否有已存在的函数或包已知道能很好地工作。

更新：

我翻译了一个MATLAB脚本，对于1-D情况还行，但还有改进空间。

更新更新：

sixtenbe为一维情况创建了一个更好的版本。

admin 更改状态以发布 2023年5月24日
0

2 答案

0 Comments

我正在研究一个类似的问题，我发现一些最好的参考资料来自于化学方面（从质谱数据中寻找峰）。如果想要了解峰值查找算法的彻底回顾，请阅读这份文献。这是我遇到过的最清晰的峰值查找技术的最佳回顾之一。（小波在嘈杂的数据中找到这种峰值是最好的方法。）

看起来您的峰已经很明确地定义了并没有隐藏在噪声中。在这种情况下，我建议使用平滑的Savitzky-Golay导数来查找峰值（如果只是对数据进行微分，你就会有一堆错误的正例）。这是一种非常有效的技术，而且相当容易实现（你需要一个具有基本操作的矩阵类）。如果你仅仅找到第一个S-G导数的零交点，我认为你会很满意。

0
0 Comments

如其名，函数scipy.signal.find_peaks能够很好地完成这个任务。但了解其参数width、threshold、distance以及最重要的prominence非常重要，才能得到良好的峰值提取。

根据我的测试和文档，prominence的概念是过滤掉嘈杂峰值并保留正确峰值的“有用概念”。

(地形) prominence是什么？它是“从山顶下降到任何更高的地形所需的最小高度”，如下图所示：

enter image description here

prominence越高，峰值越“重要”。

测试：

enter image description here

我特意使用了一个（嘈杂的）变频正弦波，因为它显示出了许多困难。我们可以看到，在高频部分，如果将最小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()

0