我在拟合光谱时遇到了问题,需要我使用特定线形和探测器分辨率函数的卷积。当我使用对称数组(以零为中心的数组)时,卷积似乎运行良好。然而,当我使用不以零为中心的数组时,会出现我无法纠正的转变,这与我试图拟合的实际数据集的形状一致。卷积按以下方式完成:
import numpy as np, matplotlib.pyplot as plt
# Sample temperature in Kelvin
T = 300
# phonon lineshape
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)
# measured elastic line
def elastic(x, A, x0):
mu, wL, wG = 0.54811, 4.88248, 2.64723 # detector parameters
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) # pseudo-voigt
# resolution function
def res_func(x):
return elastic(x,1,0)
# convolution of lineshape and resolution function
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]
# testing ranges
x1 = np.linspace(-20,20,200)
x2 = np.linspace(-15,20,200)
# random lineshape parameters
p = (1,8,2,0)
# plotting the convolutions
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 部分text)中的示例。
是的,保持卷积核围绕 x=0 对称很重要。使用
mode='valid'
时,正确移动结果也很重要。我建议根据输入 x
计算内核的范围和步长,但保持对称。尝试这样的事情:
# convolution of lineshape and resolution function
def convolve(x, F,xc,wL,x0):
arr = lineshape(x,F,xc,wL,x0)
# kernel = res_func(x)
npts = len(arr)
# note that the kernel is forced to be symmetric
# and will extend slightly past the data range
kxstep = x[1]-x[0]
kxmax = max(abs(x[0]), abs(x[-1])) + kxstep*(npts//10)
nkern = 1 + 2*int(kxmax/kxstep)
kernel = res_func(np.linspace(-kxmax, kxmax, nkern))
# now pad with the length of the kernel(!) to make sure the data
# array is long enough for mode='valid'
pad = np.zeros(nkern)
a1 = np.concatenate((pad, arr, pad))
# you could also pad with value at xmin/xmax instead of 0, with:
pad = np.ones(nkern)
a1 = np.concatenate((pad*arr[0], arr, pad*arr[-1]))
conv = np.convolve(a1, kernel/np.sum(kernel), mode='valid')
# now we need a more careful shift
noff = int((len(conv) - npts)/2)
return conv[noff:npts+noff]