使用非对称数组时两函数卷积结果偏移问题求解

拟合光谱时需将特定声子线型与探测器分辨率函数卷积，使用零中心对称数组（如np.linspace(-20,20,200)）时卷积效果正常，但使用非对称数组（如np.linspace(-15,20,200)，与实际数据集形状一致）时，卷积结果出现无法修正的偏移。现有卷积实现代码如下：

import numpy as np, matplotlib.pyplot as plt

# 样品温度（开尔文）
T = 300

# 声子线型
def lineshape(x,F,xc,wL,x0):
    bose = 1/(1 - np.exp(-(x-x0)/(0.0862*T)))
    return bose * 4/np.pi * F*(x-x0)*wL / (((x-x0)**2 - (xc**2 + wL**2))**2 + 4*((x-x0)*wL)**2)

# 测量弹性线（伪沃伊特函数）
def elastic(x, A, x0):
    mu, wL, wG = 0.54811, 4.88248, 2.64723  # 探测器参数
    lorentz = 2/np.pi * wL/(4*(x-x0)**2 + wL**2)
    gauss = np.sqrt(4*np.log(2)/np.pi)/wG * np.exp(-4*np.log(2)*((x-x0)/wG)**2)
    return A * ( mu*lorentz + (1-mu)*gauss)

# 分辨率函数
def res_func(x):
    return elastic(x,1,0)

# 线型与分辨率函数的卷积
def convolve(x, F,xc,wL,x0):
    arr = lineshape(x,F,xc,wL,x0)
    kernel = res_func(x)
    npts = min(arr.size, kernel.size)
    pad = np.zeros(npts)
    a1 = np.concatenate((pad, arr, pad))
    conv = np.convolve(a1, kernel/np.sum(kernel), mode='valid')
    m = int(npts/2)
    return conv[m:npts+m]

# 测试范围
x1 = np.linspace(-20,20,200)
x2 = np.linspace(-15,20,200)

# 随机线型参数
p = (1,8,2,0)

# 绘制卷积结果
fig, ax = plt.subplots()
ax.plot(x1, convolve(x1, *p), label='conv_sym')
ax.plot(x1,lineshape(x1, *p), label='dho')
ax.plot(x2, convolve(x2, *p), label='conv_asym')
ax.legend()
plt.show()

尝试多种卷积方式和填充方法仍未解决，卷积函数定义参考lmfit的CompositeModel示例文档，需找到通用解决方案。

问题根源 #

原代码中，分辨率函数res_func(x)直接基于输入x生成以0为中心的核，但当输入x区间不对称时，核的中心在数组中的索引位置与线型的中心（x0=0）的索引位置不匹配，导致卷积后结果整体偏移。此外，手动填充和截取数组的方式容易引入索引误差，加重偏移问题。

修改方案 #

调整卷积函数，确保分辨率核的中心始终对准线型的中心x0，并使用numpy内置的卷积模式避免手动截取错误：

import numpy as np, matplotlib.pyplot as plt

T = 300

def lineshape(x,F,xc,wL,x0):
    bose = 1/(1 - np.exp(-(x-x0)/(0.0862*T)))
    return bose * 4/np.pi * F*(x-x0)*wL / (((x-x0)**2 - (xc**2 + wL**2))**2 + 4*((x-x0)*wL)**2)

def elastic(x, A, x0):
    mu, wL, wG = 0.54811, 4.88248, 2.64723
    lorentz = 2/np.pi * wL/(4*(x-x0)**2 + wL**2)
    gauss = np.sqrt(4*np.log(2)/np.pi)/wG * np.exp(-4*np.log(2)*((x-x0)/wG)**2)
    return A * ( mu*lorentz + (1-mu)*gauss)

def res_func(x):
    return elastic(x,1,0)

# 修改后的卷积函数
def convolve(x, F, xc, wL, x0):
    arr = lineshape(x, F, xc, wL, x0)
    # 生成相对于线型中心x0的位移数组，确保核中心对准x0
    dx = x - x0
    kernel = res_func(dx)
    kernel = kernel / np.sum(kernel)  # 归一化核，保证卷积后强度守恒
    # 使用mode='same'，输出与输入数组长度一致的结果，自动处理边界填充
    return np.convolve(arr, kernel, mode='same')

# 测试
x1 = np.linspace(-20,20,200)
x2 = np.linspace(-15,20,200)
p = (1,8,2,0)

fig, ax = plt.subplots()
ax.plot(x1, convolve(x1, *p), label='conv_sym')
ax.plot(x1, lineshape(x1, *p), label='dho')
ax.plot(x2, convolve(x2, *p), label='conv_asym')
ax.legend()
plt.show()

关键修改说明 #

1. 核中心对齐：通过dx = x - x0生成相对于线型中心的位移数组，用dx构建分辨率核，确保核的中心始终对应线型的中心x0，无论输入x区间是否对称。
2. 简化卷积流程：使用np.convolve(..., mode='same')，自动处理边界填充（默认补0）并输出与输入数组长度一致的结果，避免手动填充和截取带来的索引错误。
3. 核归一化：确保分辨率核的总和为1，保证卷积后线型的强度守恒，不影响后续拟合的参数准确性。

修改后，非对称区间的卷积结果将与线型中心对齐，偏移问题解决，且该方案适用于任意x区间范围。

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