从 Scipy Signal 对象获取峰的宽度和面积

Question

如何使用 cwt get peaks 方法从 Scipy Signal 函数中获取具有位置、峰面积、峰宽等属性的峰对象：

def CWT(trace):
x = []
y = []
for i in range(len(trace)):
    x.append(trace[i].Position)
    y.append(trace[i].Intensity)
x = np.asarray(x)
y = np.asarray(y)
return signal.find_peaks_cwt(x,y)

这只是返回一个数组？

Answer 1

首先，您似乎错误地使用了 find_peaks_cwt。它的两个位置参数不是数据点的 x 和 y 坐标。第一个参数是 y 值。根本不采用 x 值，它们被假定为 0,1,2,.... 第二个参数是您感兴趣的峰宽列表；

用于计算 CWT 矩阵的一维宽度数组。通常，此范围应涵盖预期的目标峰宽度。

width 参数没有理由与数据数组的大小相同。在下面的示例中，数据有 500 个值，但我使用的宽度是 30...99。

其次，这个方法只找到峰的位置（你得到的数组有峰的索引）。没有分析它们的宽度和面积。你要么不得不去别处寻找（博客文章Peak Detection in the Python World 列出了一些替代方案，尽管它们都没有返回你想要的数据），或者想出你自己的方法来估计这些东西。

我的尝试如下。它执行以下操作：

1. 通过峰之间的中点截断信号
2. 对于每个部分，使用其中值的中值作为基线
3. 声明峰值包含大于 0.5*（峰值 + 基线）的所有值，即中值和最大值之间的中间值。
4. 查找峰值的起点和终点。 （宽度就是这些的区别）
5. 将峰面积声明为在步骤 4 中找到的区间内 (y - 基线) 的总和。

完整示例：

t = np.linspace(0, 4.2, 500)
y = np.sin(t**2) + np.random.normal(0, 0.03, size=t.shape)  # simulated noisy signal
peaks = find_peaks_cwt(y, np.arange(30, 100, 10))

cuts = (peaks[1:] + peaks[:-1])//2    # where to cut the signal 
cuts = np.insert(cuts, [0, cuts.size], [0, t.size])
peak_begins = np.zeros_like(peaks)
peak_ends = np.zeros_like(peaks)
areas = np.zeros(peaks.shape)
for i in range(peaks.size):
  peak_value = y[peaks[i]]
  y_cut = y[cuts[i]:cuts[i+1]]           # piece of signal with 1 peak
  baseline = np.median(y_cut)
  large = np.where(y_cut > 0.5*(peak_value + baseline))[0]
  peak_begins[i] = large.min() + cuts[i]
  peak_ends[i] = large.max() + cuts[i]
  areas[i] = np.sum(y[peak_begins[i]:peak_ends[i]] - baseline)

数组areas、peak_begins 和peak_ends 在这里很有趣。宽度为[84 47 36]，表示峰变细（回想一下，这些是以索引为单位的，宽度是峰中数据点的数量）。我使用这些数据将峰着色为红色：

widths = peak_ends - peak_begins
print(widths, areas)
plt.plot(t, y)
for i in range(peaks.size):
  plt.plot(t[peak_begins[i]:peak_ends[i]], y[peak_begins[i]:peak_ends[i]], 'r')
plt.show()