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

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了从Scipy Signal对象获取峰的宽度和面积相关的知识,希望对你有一定的参考价值。

如何使用cwt get peak方法从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)

这只是返回一个数组?

答案

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

用于计算CWT矩阵的1-D宽度数组。通常,该范围应涵盖目标峰的预期宽度。

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)

这里有阵列areaspeak_beginspeak_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()

picture

以上是关于从Scipy Signal对象获取峰的宽度和面积的主要内容,如果未能解决你的问题,请参考以下文章

修改 scipy.signal.welch 方法以在平均之前拒绝某些光谱

如何在 scipy.signal.filter_design.ellip 中使用 rp、rs 和 Wn 参数?

使用 scipy.signal 在 Python 中进行卷积和反卷积

scipy.signal.wiener 总是显示警告

scipy.signal信号处理的库(笔记06)

了解 scipy.signal.spectrogram() 的输出