-2

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

Zen
  • 83
  • 2
  • 9
  • 1
    You are not putting a Voigt but a pseudo-Voigt and I am not sure if it is a good one, see [here](https://en.wikipedia.org/wiki/Voigt_profile#Pseudo-Voigt_approximation). Most importantly, you are not fitting anything. And finally, you should consider to make the peak-position a variable of your function. BTW, why is the gaussian shifted by `0.1` but the Lorentzian by `0.5`. This is somewhat incorrect. – mikuszefski Mar 20 '18 at 09:02
  • As Voigt you might want to use the [Faddeeva function](https://en.wikipedia.org/wiki/Faddeeva_function), which can be expressed in terms of `erf`, which is part of `scipy.special`. Might need some data scaling to work properly, though. – mikuszefski Mar 20 '18 at 09:08
  • I have edited my code again. Could you pls check it and tell me whats wrong there? – Zen Mar 20 '18 at 09:22
  • Please edit again. First of all, most of the `import`s are not required for your example, check [How to create a Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve). The the the line where `amp` etc is set has the wrong indent. the two lines with `plt` afterwards should be removed. – mikuszefski Mar 20 '18 at 09:41
  • Tip, do not mix the order of variables and parameters. Better to have `V( x, sigma, gamma, amp, a, b )` – mikuszefski Mar 20 '18 at 09:44
  • I have edited again. Could you please tell me the whats wrong with the program part – Zen Mar 20 '18 at 09:50
  • Finally, what is actually your problem? You are still not fitting anything. – mikuszefski Mar 20 '18 at 09:50
  • I need to fit the peak with the function that I have defined. But its not fitting. I do not understand why. – Zen Mar 20 '18 at 09:51
  • I saw that you are a physicist. This is the laser spectra measured by Fabry Perot Interferometer. – Zen Mar 20 '18 at 09:57
  • First, your code does not run but produces error messages. Second, what do you mean by fitting? You are just putting some values by hand, right? There is not fit algorithm you call....and thanks for the background info. – mikuszefski Mar 20 '18 at 09:57
  • I have attached a link to my text file. It does not produce error message. It is producing a graph but it does not produce a fit that I have defined. – Zen Mar 20 '18 at 10:02
  • See the edit. I have attached a plot that it produces – Zen Mar 20 '18 at 10:03

1 Answers1

3

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

mikuszefski
  • 3,943
  • 1
  • 25
  • 38
  • Thank you. This is exactly what I want. – Zen Mar 20 '18 at 11:53
  • BTW, what are those numbers 2.69269, 2.42843, 4.47163, ... that you have used in f to define the pseudo_voigt? – Zen Mar 20 '18 at 16:28
  • @Zen THe numbers are given in the Wiki article on the Voigt profile I linked above. As Voigt is the convolution, you may approximate by a sum, but the way this sum depends on gamma and sigma is not trivial. The numbers are the according expansion terms. And I guess a lot of people figured to what order this expansion makes sense in typical applications like spectral line analysis. – mikuszefski Mar 21 '18 at 07:32
  • Is it also possible to get the chisq value of the fits in python? I want to present the chisq value of the fits in the legend. – Zen Mar 22 '18 at 15:05
  • @Zen you should get the square deviations by calling `res( sol, x, y )**2` or the total `chi2` by `sum( res( sol, x, y )**2 )` – mikuszefski Mar 22 '18 at 15:30
  • @Zen Note, you also need this to calculate the errors as described [here](https://stackoverflow.com/a/14857441/803359) – mikuszefski Mar 22 '18 at 15:43