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

Question

我正在尝试使用 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")

Answer 1

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

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