I have been having problems fitting a spectrum that requires me to use a convolution of a certain lineshape and a detector resolution function. When I use a symmetric array (array centered at zero) the convolution seems to work well. However a shift that I just can't correct appears when I use an array not centered at zero, which coincides with the shape for the actual dataset I'm trying to fit. The convolution is done in the following way:
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()
Does anyone know how to solve this problem? I have tried different types of convolution and various paddings, but a general solution has eluded me. I'm fitting with lmfit and the definition of the convolution function was taken from the example in the CompositeModel section text).
Yes, it is important to keep the convolution kernel symmetric around x=0. When using mode='valid' it is also important to shift the result correctly. I would suggest calculating the range and step size for the kernel based on the input x, but keeping it symmetric. Try something like:
# 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]