scipy 信号 find_peaks_cwt 没有准确找到峰值？

Posted 2023-03-12

【中文标题】scipy 信号 find_peaks_cwt 没有准确找到峰值？

【英文标题】：scipy signal find_peaks_cwt not finding the peaks accurately?

【发布时间】：2014-10-23 15:23:54

【问题描述】：

我有一个一维信号，我试图在其中找到峰值。我希望能完美地找到它们。

我目前正在做：

import scipy.signal as signal
peaks = signal.find_peaks_cwt(data, np.arange(100,200))

下图是带有红点的图表，显示了find_peaks_cwt() 发现的峰值位置。

如您所见，计算的峰值不够准确。真正重要的是右边的三个。

我的问题：如何使这个更准确？

更新：数据在这里：http://pastebin.com/KSBTRUmW

对于某些背景，我想做的是在图像中找到手指之间的空间。绘制的是手周围轮廓的 x 坐标。青色斑点 = 峰。如果有更可靠/更强大的方法，请发表评论。

【问题讨论】：

如果我理解正确，峰位精度受峰宽参数中最小数字的限制。所以不是np.arange(100,200)，也许np.array([10, 50, 100, 200])会给出更好的结果。

所以，我刚刚尝试了这个，它似乎确实更好，但它也发现了很多局部颠簸，而不仅仅是大的颠簸，这使得输出充满误报。

我也遇到了类似的问题，想问一下如何正确确定峰宽参数？我想测试您的数据但我无法打开链接（无法检索到请求的 url）

再试一次链接-我刚刚检查过它确实有效。

【参考方案1】：

已解决，解决方案：

先过滤数据：

  window = signal.general_gaussian(51, p=0.5, sig=20)
  filtered = signal.fftconvolve(window, data)
  filtered = (np.average(data) / np.average(filtered)) * filtered
  filtered = np.roll(filtered, -25)

然后按照rapelpy的回答使用angrelextrema。

结果：

【讨论】：

能否添加查找峰的代码？我认为rapelpy的答案已经过编辑

cjm2671，请详细解释您的答案，包括与拉普拉斯核的卷积、缩放和滚动。

Below我根据这个答案发布了一个工作函数。

【参考方案2】：

使用这个函数有一个更简单的解决方案： https://gist.github.com/endolith/250860 这是http://billauer.co.il/peakdet.html的改编版

我刚刚尝试了您提供的数据，结果如下。无需预过滤...

享受:-)

【讨论】：

基于 peakdet 函数，我希望看到噪声峰值出现。你怎么没有在最左边、最宽的峰中检测到两个峰？

@aerijman：噪声量由delta 值（也称为峰/谷阈值）控制。为了达到这个结果，delta 可能是 ~75。要查看您提到的嘈杂峰，它需要明显更小（可能约为 5）。

【参考方案3】：

获取原始数据后进行编辑。

argelmax 和 arglextrma 退出比赛。

曲线非常嘈杂，因此您必须使用较小的峰宽（如 pv. 所述）和噪音。

我找到的最好的看起来不是很好。

import numpy as np
import scipy.signal as signal

peakidx = signal.find_peaks_cwt(y_array, np.arange(10,15), noise_perc=0.1)
print peakidx

[10, 100, 132, 187, 287, 351, 523, 597, 800, 1157, 1451, 1673, 1742, 1836]

【讨论】：

我试过这个，不幸的是我的结果比你的要差imgur.com/E4QbRXY。我可以先对数据进行平滑处理，尽管我觉得这可能会扭曲结果。

也许这会起作用：signal.argrelmax(y_array, order=5) 其中 order 是每边用于比较的点数。

很遗憾没有；它产生的分数更少，但似乎为坏分数牺牲了好分数:(。还有其他想法吗？

你能把你的数据上传到某个地方吗？

好的，我已经用多种方法解决了！我使用了 angrelextrema，但首先将信号与拉普拉斯窗口进行卷积以消除噪声；然后它是完美的。它几乎是像素完美的！谢谢！

【参考方案4】：

根据@cjm2671 的回答，这是一个在噪声信号中查找相对最大值和最小值的工作示例：

import numpy as np
import matplotlib.pyplot as plt

from scipy.ndimage.filters import gaussian_filter1d
from scipy import signal

data =np.array([5.14,5.22,5.16,4.82,4.46,4.36,4.4,4.35,4.13,3.83,3.59,3.51,3.46,3.27,3.08,3.03,2.95,2.96,2.98,3.02,3.09,3.14,3.06,2.84,2.68,2.72,2.92,3.23,3.44,3.5,3.28,3.34,3.73,3.97,4.26,4.48,4.5,5.06,6.02,6.68,7.09,7.58,8.6,9.85,10.7,11.3,11.3,11.6,12.3,12.6,12.8,12.8,12.5,12.4,12.2,12.2,12.3,11.9,11.2,10.6,10.3,10.3,10.,9.53,8.97,8.55,8.49,8.41,8.09,7.71,7.34,7.26,7.42,7.47,7.37,7.17,7.05,7.02,7.09,7.23,7.18,7.16,7.47,7.92,8.55,8.68,8.31,8.52,9.11,9.59,9.83,9.73,10.2,11.1,11.6,11.7,11.7,12.,12.6,13.1,13.3,13.2,13.,12.6,12.3,12.2,12.3,12.,11.6,11.1,10.9,10.9,10.7,10.3,9.83,9.64,9.63,9.37,8.88,8.39,8.14,8.12,7.92,7.48,7.06,6.87,6.87,6.63,6.17,5.71,5.45,5.45,5.34,5.05,4.78,4.57,4.47,4.37,4.16,3.95,3.88,3.83,3.69,3.64,3.57,3.5,3.51,3.33,3.14,3.09,3.06,3.12,3.11,2.94,2.83,2.76,2.74,2.77,2.75,2.73,2.72,2.59,2.47,2.53,2.54,2.63,2.76,2.78,2.75,2.69,2.54,2.42,2.58,2.79,2.83,2.78,2.71,2.77,2.88,2.97,2.97,2.9,2.92,3.16,3.29,3.28,3.49,3.97,4.32,4.49,4.82,5.08,5.48,6.03,6.52,6.72,7.16,8.18,9.52,10.9,12.1,12.6,12.9,13.3,13.3,13.6,13.9,13.9,13.6,13.3,13.2,13.2,12.8,12.,11.4,11.,10.9,10.4,9.54,8.83,8.57,8.61,8.24,7.54,6.82,6.46,6.43,6.26,5.78,5.29,5.,5.08,5.14,5.,4.84,4.56,4.38,4.52,4.84,5.33,5.52,5.56,5.82,6.54,7.27,7.74,7.64,8.14,8.96,9.7,10.2,10.2,10.5,11.3,12.,12.4,12.5,12.3,12.,11.8,11.8,11.9,11.6,11.,10.3,10.,9.98,9.6,8.87,8.16,7.76,7.74,7.54,7.03,6.54,6.25,6.26,6.09,5.66,5.31,5.08,5.19,5.4,5.38,5.38,5.22,4.95,4.9,5.02,5.28,5.44,5.93,6.77,7.63,8.48,8.89,8.97,9.49,10.3,10.8,11.,11.1,11.,11.,10.9,11.1,11.1,11.,10.7,10.5,10.4,10.3,10.4,10.3,10.2,10.1,10.2,10.4,10.4,10.5,10.7,10.8,11.,11.2,11.2,11.2,11.3,11.4,11.4,11.3,11.2,11.2,11.,10.7,10.4,10.3,10.3,10.2,9.9,9.62,9.47,9.46,9.35,9.12,8.82,8.48,8.41,8.61,8.83,8.77,8.48,8.26,8.39,8.84,9.2,9.31,9.18,9.11,9.49,9.99,10.3,10.5,10.4,10.2,10.,9.91,10.,9.88,9.47,9.,8.78,8.84,8.8,8.55,8.17,8.02,8.03,7.78,7.3,6.8,6.54,6.53,6.35,5.94,5.54,5.33,5.32,5.14,4.76,4.43,4.28,4.3,4.26,4.11,4.,3.89,3.81,3.68,3.48,3.35,3.36,3.47,3.57,3.55,3.43,3.29,3.19,3.2,3.17,3.21,3.33,3.37,3.33,3.37,3.38,3.26,3.34,3.62,3.86,3.92,3.83,3.69,4.2,4.78,5.03,5.13,5.07,5.4,6.,6.42,6.5,6.45,6.48,6.55,6.66,6.79,7.06,7.33,7.53,7.9,8.17,8.29,8.6,9.05,9.35,9.51,9.69,9.88,10.2,10.6,10.8,10.6,10.7,10.9,11.2,11.3,11.3,11.4,11.5,11.6,11.8,11.7,11.3,11.1,10.9,11.,11.2,11.1,10.6,10.3,10.1,10.2,10.,9.6,9.03,8.73,8.73,8.7,8.53,8.26,8.06,8.03,8.03,7.97,7.94,7.77,7.64,7.85,8.29,8.65,8.68,8.61,9.08,9.66,9.86,9.9,9.71,10.,10.9,11.4,11.6,11.8,11.8,11.9,11.9,12.,12.,11.7,11.3,10.9,10.8,10.7,10.4,9.79,9.18,8.89,8.87,8.55,7.92,7.29,6.99,6.98,6.73,6.18,5.65,5.35,5.35,5.22,4.89,4.53,4.28,4.2,4.05,3.83,3.67,3.61,3.61,3.48,3.27,3.05,2.9,2.93,2.99,2.99,2.98,2.94,2.88,2.89,2.92,2.86,2.97,3.,3.02,3.03,3.11,3.07,3.46,3.96,4.09,4.25,4.3,4.67,5.7,6.33,6.68,6.9,7.09,7.66,8.25,8.75,8.87,8.97,9.78,10.9,11.6,11.8,11.8,11.9,12.3,12.6,12.8,12.9,12.7,12.4,12.1,12.,12.,11.9,11.5,11.1,10.9,10.9,10.7,10.5,10.1,9.91,9.84,9.63,9.28,9.,8.86,8.95,8.87,8.61,8.29,7.99,7.95,7.96,7.92,7.87,7.77,7.78,7.9,7.73,7.51,7.43,7.6,8.07,8.62,9.06,9.24,9.13,9.14,9.46,9.76,9.8,9.78,9.73,9.82,10.2,10.6,10.8,10.8,10.9,11.,10.9,11.,11.,10.9,10.9,11.,10.9,10.8,10.5,10.2,10.2,10.2,9.94,9.51,9.08,8.88,8.88,8.62,8.13,7.64,7.37,7.37,7.23,6.91,6.6,6.41,6.42,6.29,5.94,5.57,5.43,5.46,5.4,5.17,4.95,4.84,4.87,4.9,4.69,4.4,4.24,4.26,4.35,4.34,4.19,3.96,3.97,4.42,5.03,5.34,5.15,4.73,4.86,5.35,5.88,6.35,6.52,6.81,7.26,7.62,7.66,8.01,8.91,10.,10.9,11.3,11.1,10.9,10.9,10.8,10.9,11.,10.7,10.2,9.68,9.43,9.42,9.17,8.66,8.13,7.83,7.81,7.62,7.21,6.77,6.48,6.44,6.31,6.06,5.72,5.47,5.45,5.42,5.31,5.23,5.22,5.3,5.32,5.16,4.96,4.82,4.73,4.9,4.95,4.91,4.92,5.41,6.04,6.34,6.8,7.08,7.26,7.95,8.57,8.78,8.95,9.06,9.14,9.2,9.33,9.53,9.65,9.69,9.53,9.18,9.02,9.,8.82,8.42,8.05,7.85,7.84,7.79,7.58,7.28,7.09,7.07,6.94,6.68,6.35,6.09,6.2,6.27,6.24,6.16,5.91,5.86,6.02,6.19,6.45,6.92,7.35,7.82,8.4,8.87,9.,9.09,9.61,9.99,10.4,10.8,10.7,10.7,11.1,11.4,11.5,11.5,11.3,11.3,11.4,11.7,11.8,11.5,11.,10.5,10.4,10.3,9.94,9.23,8.52,8.16,8.15,7.86,7.23,6.59,6.26,6.25,6.04,5.55,5.06,4.81,4.78,4.62,4.28,3.98,3.84,3.92,3.93,3.68,3.46,3.31,3.16,3.11,3.18,3.19,3.14,3.28,3.3,3.16,3.19,3.04,3.07,3.59,3.83,3.82,3.95,4.06,4.71,5.39,5.89,6.06,6.08,6.45,6.97,7.57,8.1,8.25,8.55,8.92,9.09,9.2,9.32,9.36,9.45,9.65,9.73,9.7,9.82,9.94,9.92,9.97,9.93,9.78,9.63,9.48,9.49,9.48,9.2,8.81,8.34,8.,8.06,7.98,7.63,7.47,7.37,7.24,7.2,7.05,6.93,6.83,6.59,6.44,6.42,6.33,6.18,6.37,6.29,6.1,6.34,6.57,6.54,6.77,7.21,7.58,7.86,8.11,8.57,9.07,9.45,9.67,9.68,9.87,10.2,10.4,10.4,10.4,10.4,10.4,10.5,10.6,10.7,10.4,9.98,9.58,9.45,9.51,9.44,9.09,8.68,8.46,8.36,8.17,7.88,7.55,7.34,7.3,7.17,6.97,6.88,6.69,6.69,6.77,6.77,6.81,6.67,6.5,6.57,6.99,7.4,7.59,7.8,8.45,9.47,10.4,10.8,10.9,10.9,11.,11.4,11.8,12.,11.9,11.4,10.9,10.8,10.8,10.5,9.76,8.99,8.59,8.58,8.43,8.05,7.61,7.26,7.16,6.99,6.58,6.15,5.98,5.93,5.71,5.48,5.22,5.06,5.08,4.95,4.78,4.62,4.45,4.48,4.65,4.66,4.69])

dataFiltered = gaussian_filter1d(data, sigma=5)
tMax = signal.argrelmax(dataFiltered)[0]
tMin = signal.argrelmin(dataFiltered)[0]

plt.plot(data, label = 'raw')
plt.plot(dataFiltered, label = 'filtered')
plt.plot(tMax, dataFiltered[tMax], 'o', mfc= 'none', label = 'max')
plt.plot(tMin, dataFiltered[tMin], 'o', mfc= 'none', label = 'min')
plt.legend()
plt.savefig('fig.png', dpi = 300)

高斯滤波器已经使用高斯窗口实现了卷积。我们只需要将窗口的标准差作为参数。

在这种情况下，这种方法比使用signal.find_peaks_cwt 效果更好。