Search code examples
pythondata-fitting

Fitting the data with a voigt function in python

I am trying to fit my data with a Voigt function. I used the code given below. But the fit is not there with a right range.. and I do not know how to set the range to fit. Can anyone help me pls?

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

and here is the link to my data https://www.dropbox.com/s/vm9ta6samnlc0s2/calibration.txt?dl=0 Thank you for any help. PS: The program produces the plot given below: https://www.dropbox.com/s/3rbuq4v7gcc92m7/figure_1.png?dl=0

Solution

  • I am not really sure what you are doing or trying to do, but this is how I would do it ( assuming that the sigma and gamma are the same for all peaks. Did not think to much if this makes sense in a Fabry-Perot)

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

    which gives

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

    and the following graph

    myfit