在python中用voigt函数拟合数据

Question

我正在尝试用 Voigt 函数拟合我的数据。我使用了下面给出的代码。但是合适的范围不合适..而且我不知道如何设置合适的范围。谁能帮帮我？

import numpy as np
import matplotlib.pyplot as plt
from scipy import asarray as exp
from numpy import genfromtxt

data= genfromtxt ('calibration.txt')
x=data[:,0]
y=data[:,1]
plt.xlim(0,1)
plt.ylim(0,1.25)
plt.xlabel("Voltage [V]")
plt.ylabel("Intensity")


def V(amp,x, sigma, gamma,a,b):
"""
Return the Voigt line shape at x with Lorentzian component HWHM gamma
and Gaussian component sigma, a&b as the center.

"""

    return amp*np.exp(-(x-a)**2/(2*(sigma)**2))+gamma/np.pi/((x-b)**2+(gamma)**2)
amp,sigma, gamma,a,b =0.9, 0.1,0.04, 0.5,0.5
plt.plot(x,y,'b.',x, V(amp, x, sigma, gamma,a,b))
plt.show()

这是我的数据的链接 https://www.dropbox.com/s/vm9ta6samnlc0s2/calibration.txt?dl=0 感谢您的任何帮助。 PS：该程序产生下面给出的情节： https://www.dropbox.com/s/3rbuq4v7gcc92m7/figure_1.png?dl=0

Answer 1

我不太确定你在做什么或试图做什么，但这就是我会做的事情（假设所有峰的 sigma 和 gamma 都是相同的。如果这在一个有意义的法布里-珀罗）

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def cauchy(x, x0, g):
    return 1. / ( np.pi * g * ( 1 + ( ( x - x0 )/ g )**2 ) )

def gauss( x, x0, s):
    return 1./ np.sqrt(2 * np.pi * s**2 ) * np.exp( - (x-x0)**2 / ( 2 * s**2 ) )

def pseudo_voigt( x, x0, s, g, a ):
    fg = 2 * s * np.sqrt( 2 * np.log(2) )
    fl = 2 * g
    f = ( fg**5 +  2.69269 * fg**4 * fl + 2.42843 * fg**3 * fl**2 + 4.47163 * fg**2 * fl**3 + 0.07842 * fg * fl**4+ fl**5)**(1./5.)
    eta = 1.36603 * ( fl / f ) - 0.47719 * ( fl / f )**2 + 0.11116 * ( f / fl )**3
    return a * ( eta * cauchy( x, x0, f) + ( 1 - eta ) * gauss( x, x0, f ) )

def all_peaks(x, mus, amps,  s, g ):
    out = 0
    for m, a in zip( mus, amps ):
        out += pseudo_voigt( x, m, s, g, a )
    return out

def res( params, xData, yData):
    mus = params[0:5]
    amp = params[5:10]

    sig = params[-3]
    gam = params[-2]
    off = params[-1]
    yth = np.fromiter( ( abs( off ) + all_peaks( x , mus, amp, sig, gam) for x in xData ), np.float )
    diff = yth - yData
    return diff

sigma, gamma = 0.007, 0.02
offset = 0.01
muList = [ 0.5, 2.6, 4.8, 6.8,  8.9 ]
ampList = [ .135 ] * 5

data = np.loadtxt( 'calibration.txt' )
x = data[ :,0 ]
y = data[ :,1 ]

sol, err = leastsq( res, muList + ampList + [sigma , gamma, offset ], args=(x, y) )
print sol
plt.xlabel( "Voltage [V]" )
plt.ylabel( "Intensity" )

plt.plot( x,y,ls='', marker='o' )
plt.plot( x, sol[-1] + all_peaks( x, sol[0:5],sol[5:10], sol[-3], sol[-2]) )
plt.show()

给了

[
    4.97681822e-01 2.63788309e+00 4.74796088e+00 6.83620027e+00 8.90127524e+00 
    1.28754082e-01 1.35709531e-01 1.34679136e-01 1.35460544e-01 1.39491029e-01 
    5.61700040e-03 1.93814469e-02 9.99057213e-03
]

还有下图