如何使用Python Scipy对带积分限的离散数据进行数值积分

你的原有实现逻辑是可行的，以下是几种更优的实现方式：

1. 高效索引定位优化（推荐日常使用） #

原来的np.where需要遍历全数组匹配条件，对于welch输出的严格单调升序频率数组，用np.searchsorted二分查找的效率高得多，代码也更简洁，逻辑和你原有实现完全一致：

def psd_integrate(f, p, f0, fN):
    # 增加边界判断避免索引错误
    if f0 >= f[-1] or fN <= f[0]:
        return 0.
    # 二分查找定位上下限对应索引
    k0 = np.searchsorted(f, f0, side='left')
    kN = np.searchsorted(f, fN, side='right')
    return simpson(p[k0:kN], x=f[k0:kN])

该方案在频率数组长度较大时性能优势非常明显，且输出结果和你原有函数完全一致。

2. 带边界插值的高精度实现 #

如果积分上下限刚好不在现有频率采样点上，且对精度要求较高，可以先对边界点做线性插值补全，减少边界截断误差：

from scipy.interpolate import interp1d

def psd_integrate_high_precision(f, p, f0, fN):
    if f0 >= f[-1] or fN <= f[0]:
        return 0.
    k0 = np.searchsorted(f, f0, side='left')
    kN = np.searchsorted(f, fN, side='right')
    # 提取插值所需的相邻点
    f_seg = f[k0-1:kN+1] if k0>0 else f[:kN+1]
    p_seg = p[k0-1:kN+1] if k0>0 else p[:kN+1]
    interp_func = interp1d(f_seg, p_seg, kind='linear')
    # 拼接插值后的边界点与原有效区间
    f_integrate = np.concatenate([[f0], f[k0:kN], [fN]])
    p_integrate = np.concatenate([[interp_func(f0)], p[k0:kN], [interp_func(fN)]])
    return simpson(p_integrate, x=f_integrate)

该方案适合频率分辨率较低、积分区间较小的场景，能有效降低边界截断带来的误差。

3. 均匀频率下的简化实现 #

由于welch返回的频率数组是均匀间隔的，你也可以直接传入采样间隔dx参数省略频率数组的传递，进一步提升计算速度：

def psd_integrate_uniform(f, p, f0, fN):
    if f0 >= f[-1] or fN <= f[0]:
        return 0.
    df = f[1] - f[0]
    k0 = np.searchsorted(f, f0, side='left')
    kN = np.searchsorted(f, fN, side='right')
    return simpson(p[k0:kN], dx=df)

测试示例 #

用你提供的MWE测试0.5~6Hz区间的功率，结果会接近理论值1（两个振幅为1的正弦波平均功率和为0.5+0.5=1）：

f,p = welch(y,fs=freq,nperseg=2**13)
print(psd_integrate(f, p, 0.5, 6))
print(psd_integrate_high_precision(f, p, 0.5, 6))

你提到的使用seaborn修改绘图样式的操作不影响上述计算逻辑。

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