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的数据中添加零，因为这是卷积的纯数学限制，或者是否有错误我的代码的/解决方案？

标签： pythonpython-3.xmathcurve-fittingconvolution

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 的定义错过了中心 corrdinate 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()