如何在Python中计算曲线干扰尖峰下的面积

我的需求 #

我有两个np.array：一个存储频率值（x轴），另一个存储对应的信号强度/功率谱密度（y轴）。无噪声的原始信号近似为对数曲线，形态会随数据变化。但信号中存在多个干扰尖峰，部分尖峰相互重叠，我需要：

• 计算每个尖峰下的面积；
• 若尖峰重叠，需分别计算单个尖峰的面积（理想情况要考虑重叠处的相长干扰，或分割相交区域）。

已尝试但失效的方案 #

• 提取峰值与底部宽度：当尖峰过宽或相互重叠时，该方法频繁失效；
• 拟合/滤波生成基准曲线：尝试过多种曲线拟合算法和滤波器，但始终得不到与干净信号匹配的基准线，无法通过与原曲线对比计算尖峰面积。

示例曲线 #

下图展示了data曲线、我尝试拟合的model1曲线，以及手绘的代表干净信号的黄色曲线；粉色阴影区域是需要计算面积的其中一个尖峰区域：

数据形态示例 #

频率数组（起始值为0.5）：

[ 0.5 ...  79.5  80.   80.5  81.   81.5  82.   82.5  83.   83.5  84.
 84.5  85.   85.5  86.   86.5  87.   87.5  88.   88.5  89.   89.5  90.
 90.5  91.   91.5  92.   92.5  93.   93.5  94.   94.5  95.   95.5  96.
 96.5  97.   97.5  98.   98.5  99.   99.5 100. ]

信号数组与频率数组长度一致，示例片段：

[6.83248573e-27 6.38424451e-27 4.40532611e-27 2.46641238e-27
 2.79056227e-27 1.91667602e-27 2.01585530e-27 2.81595644e-27
 ...]

解决方案思路 #

结合你的场景，这里有几个针对性的方向可以尝试：

1. 改进基准曲线生成：稳健对数拟合 #

既然原始信号是对数形态，我们可以避开尖峰区域来拟合基准线：

• 分段筛选非尖峰区域：用滑动窗口计算信号的标准差，筛选出波动较小的区间（这些是无干扰的原始信号区域），只在这些区域拟合对数模型y = a*ln(x) + b；
• 稳健拟合算法：使用带Huber损失的最小二乘拟合（比如scipy.optimize.least_squares），降低尖峰这类异常值对拟合结果的权重，生成更接近真实干净信号的基准线。

2. 重叠尖峰分离：多峰拟合或小波变换 #

如果你的尖峰符合特定模型（比如高斯、洛伦兹），可以尝试：

• 多峰组合拟合：将原始信号建模为「对数基准线 + N个尖峰函数」的组合，通过优化每个尖峰的位置、幅度、宽度参数，分离重叠的尖峰，之后对单个尖峰函数积分得到面积；
• 小波变换分解：选择合适的小波基（比如Morlet小波）对信号做多尺度分解，提取不同尺度下的尖峰成分，分离后再分别计算面积。

3. 重叠区域分割：导数或交点检测 #

对于重叠尖峰的分割：

• 导数谷底检测：计算信号的一阶导数，找到相邻尖峰之间的谷底（导数由负转正的点），以此作为两个尖峰的分界；
• 函数交点求解：如果已经拟合出单个尖峰的函数模型，直接求解两个重叠尖峰函数的交点，以交点为界分别计算各自的面积，重叠区域可根据相长干扰的物理模型调整（比如功率叠加的话，重叠区域面积需按实际情况叠加）。

代码示例：稳健拟合生成基准线 #

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

# 定义对数基准线模型
def log_baseline(params, x):
    a, b = params
    return a * np.log(x) + b

# Huber损失函数：降低异常值（尖峰）的权重
def huber_cost(params, x, y):
    residuals = y - log_baseline(params, x)
    delta = 0.8  # 调整该值控制对尖峰的容忍度
    cost = np.where(np.abs(residuals) <= delta, 
                    0.5 * residuals**2, 
                    delta * (np.abs(residuals) - 0.5 * delta))
    return cost

# 假设x_freq是频率数组，y_signal是信号数组
x_freq = np.array([0.5, 1.0, ..., 100.0])
y_signal = np.array([6.83e-27, 6.38e-27, ...])

# 初始参数猜测
initial_guess = [1.0, 0.0]

# 执行稳健拟合
fit_result = least_squares(huber_cost, initial_guess, args=(x_freq, y_signal))
baseline = log_baseline(fit_result.x, x_freq)

# 计算单个尖峰的面积（这里以第一个检测到的尖峰为例）
# 先检测尖峰区域：信号高于基准线的部分
peak_mask = y_signal > baseline
peak_x = x_freq[peak_mask]
peak_y = y_signal[peak_mask] - baseline[peak_mask]

# 用梯形法积分计算面积
single_peak_area = np.trapz(peak_y, peak_x)
print(f"单个尖峰面积：{single_peak_area}")

内容的提问来源于stack exchange，提问作者Marco Boerner