pythonnumpyconvolutiondata-fitting

Shifted result for convolution of two functions when using asymmetric arrays


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()

shifted convolution

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).


Solution

  • 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]