python - 错误的 Voigt 输出/具有不对称 x 输入的卷积

Question

目前，我正在对我的数据拟合高斯或洛伦兹，但两者拟合得不够好，我想切换到 Voigt 拟合，即它们的卷积。

我从 https://www.originlab.com/doc/Origin-Help/Voigt-FitFunc 检索了 Voigt 函数并将其写入Python。

import numpy as np
import matplotlib.pyplot as plt

def lorentz_voigt(x, A, xc, wL):
    return (2*A / np.pi) * (wL / ((4*(x.astype(float) - xc)**2) + wL**2))

def gauss_voigt(x, wG):
    return np.sqrt((4*np.log(2)) / np.pi) * ((np.exp(-(((4*np.log(2)) / (wG**2))*(x.astype(float))**2))) / (wG))

def Voigt(x, xc, A, wG, wL):
    return np.convolve(lorentz_voigt(x, A, xc, wL), gauss_voigt(x, wG), 'same')

symx = np.linspace(-100, 100, 1001)
asymx = np.linspace(0, 100, 1001)
symy = Voigt(symx, 50, 1, 5, 5)
asymy = Voigt(asymx, 50, 1, 5, 5)
plt.clf()
plt.plot(symx, symy)
plt.plot(asymx, asymy)

如下图所示，具有不对称 x 轴输入的 Voigt 函数(橙色)无法再现正确的 Voigt 轮廓(如蓝色所示)。

我的数据的波数为 600-4000 cm^-1，我想知道是否需要在 -4000 到 600 cm^{-1<的数据中添加零} 因为这是卷积的纯粹数学限制，或者我的代码/解决方案中是否有错误？

Answer 1

卷积对 x 轴一无所知。最好让它返回完整的卷积，然后将其缩小到原始的 x 范围，例如:

import numpy as np
import matplotlib.pyplot as plt

def lorentz_voigt(x, A, xc, wL):
    return (2*A / np.pi) * (wL / ((4*(x.astype(float) - xc)**2) + wL**2))

def gauss_voigt(x, xc, wG):
    return np.sqrt((4*np.log(2)) / np.pi) * ((np.exp(-(((4*np.log(2)) / (wG**2))*((x-xc).astype(float))**2))) / (wG))

def Voigt(x, xc, A, wG, wL):
    return np.convolve(lorentz_voigt(x, A, xc, wL), gauss_voigt(x, xc, wG), 'full')

symx = np.linspace(-100, 100, 1001)
asymx = np.linspace(0, 200, 1001)
symy = Voigt(symx, 50, 1, 5, 5)[::2]
asymy = Voigt(asymx, 50, 1, 5, 5)[::2]

plt.clf()
plt.plot(symx, symy,'r')
plt.plot(asymx, asymy,'b--')

另请注意，您对 gauss_voigt 的定义错过了中心坐标 xc。

正如我在上面的评论中所说，还有另一种不基于卷积的计算 voigt 函数的方法。卷积在计算时间方面是昂贵的，当用作拟合模型时可能会变得烦人。以下示例代码不需要卷积。

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import wofz

def voigt(x, amp, pos, fwhm, shape):
    tmp = 1/wofz(np.zeros((len(x))) +1j*np.sqrt(np.log(2.0))*shape).real
    return tmp*amp*wofz(2*np.sqrt(np.log(2.0))*(x-pos)/fwhm+1j*np.sqrt(np.log(2.0))*shape).real

x = np.linspace(0, 100, 1001)
y0 = voigt(x, 1, 50, 5, 0)
y1 = voigt(x, 1, 50, 5, 1)
plt.plot(x,y0,'k',label='shape = 0')
plt.plot(x,y1,'b',label='shape = 1')
plt.legend(loc=0)
plt.show()