使用 trapz 进行峰值积分,但不是从基线和可视化开始

问题描述 投票:0回答:1

我正在尝试使用 numpy.trapz 或 scipy.trapz 积分信号峰值。我手动选择峰值最小值和最大值的范围。我面临的问题是我认为这些函数将 AUC 积分到基线。有时我的信号开始或结束高于基线。我愿意使用基线校正方法,但我需要在不进行任何校正的情况下可视化积分区域,以确保它们是正确的。

图片来源:https://terpconnect.umd.edu/~toh/spectrum/Integration.html

生成数据的示例

import numpy
import peakutils
from peakutils.plot import plot as pplot
from matplotlib import pyplot


centers = (30.5, 72.3)
x = numpy.linspace(0, 120, 121)
y = (peakutils.gaussian(x, 5, centers[0], 3) +
    peakutils.gaussian(x, 7, centers[1], 10) +
    numpy.random.rand(x.size))
pyplot.figure(figsize=(10,6))
pyplot.plot(x, y)
pyplot.title("Data with noise")

python numpy scipy signal-processing peak-detection
1个回答
0
投票

要获取峰值下的面积,您可以使用

scipy.interpolate.make_interp_spline
插入曲线,该曲线具有
integrate
方法。这将是
x=a
x=b
之间曲线下的整个面积,因此您需要减去通过连接
x=a
x=b
处的点和 x 轴创建的梯形面积。在我下面分享的代码中,我假设
f(a) >= 0
f(b) >= 0

注意:我添加了另一个峰值并降低了噪声水平,以使端点更加突出/明显。

import numpy as np
import peakutils
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.interpolate import make_interp_spline

plt.close("all")

rng = np.random.default_rng(42)

centers = (30.5, 72.3, 112.1)
x = np.arange(0, 120)
y = (peakutils.gaussian(x, 5, centers[0], 3)
     + peakutils.gaussian(x, 7, centers[1], 10)
     + peakutils.gaussian(x, 6, centers[2], 5)
     + rng.random(x.size)*0.2)
     
peak_indices, _ = find_peaks(-y, prominence=4)

fig, ax = plt.subplots()
ax.plot(x, y, label="Data")
ax.plot(x[peak_indices], y[peak_indices], "x", color="r", label="End points")
# ax.title("Data with noise")


x1, x2 = x[peak_indices]

interp = make_interp_spline(x, y, k=1)
m = (interp(x2) - interp(x1))/(x2 - x1)
x = np.linspace(x1, x2, 200)
ax.fill_between(x, interp(x), m*(x - x1) + interp(x1), alpha=0.5, 
                label="Desired area")
ax.fill_between(x, m*(x - x1) + interp(x1), alpha=0.5, 
                label="Trapezoid area")

ax.legend()

# Integral under the peak minus the area of the trapezoid connecting the two
# ends of the peak.
# NOTE: EDGE CASE OF PEAK ENDS GOING BELOW THE X-AXIS IS NOT BEING CONSIDERED
integral = interp.integrate(x1, x2) - 0.5*(x2-x1)*(interp(x1) + interp(x2))
print(integral)  # 163.8743118056535

© www.soinside.com 2019 - 2024. All rights reserved.