Python Workshop I: Fitting a single symmetric peak

Time to play with real data! In this workshop, we will fit a peak with a single symmetric function. The script for the entire analysis is given below. Copy the whole thing into an editor and save it as ‘name_of_file.py’. I have used the Raman spectrum of silicon centered at the dominant peak (520 cm-1). This data is available here. Please download it and put it into the same directory as the script. The code is pretty self-explanatory but here is an outline of the analysis routine:

1. Import data

2. Background removal. In this section, the regions outside the peak are isolated and fit to a linear background. To do this, we first describe the regions of the required data. The variables ‘ind_bg_low’ and ‘ind_bg_high’ are actually boolean arrays created by the specified limits. So, ‘ind_bg_low’ will be ‘True’ for all x greater than the minimum value of x and less than 450. These two variables are then used as ‘indices’ to splice the full x and y datasets.

3. Fitting. As with any optimization routine, one needs to define a function that can be fitted to the data. In this example, we have chosen a Lorentzian peak with three parameters that can be fitted: the half-width at half-maximum (hwhm), peak center and the peak intensity. These parameters are assigned to a list of parameters called ‘p’. The items of the list, e.g. item 1(index=0), correspond to a given parameter. We have defined the list to have the parameters in this order: hwhm, peak center, intensity. In addition to defining a function, we also need to define a minimisation criterion. The leastsq function minimises the difference between the y values of the data and the fit. The function ‘residuals’  is defined such that it returns this difference and hence, it is the function that is given to the leastsq function call. The first part of the results outputted by the leastsq function contain the parameters of best fit. These parameters are then put into the variable ‘best_parameters’. The fit is evaluated by substituting ‘best_parameters’ into the Lorentzian function.

4. Plotting data and fit

This particular routine is an introduction to simple scenarios of each step. In the next workshop, we will look at more sophisticated background correction, fitting and plotting options.

#########################################################################
####################### IMPORTING REQUIRED MODULES ######################

import numpy
import pylab
from scipy.optimize import leastsq # Levenberg-Marquadt Algorithm #

#########################################################################

x = a[:,0]
y = a[:,1]

#########################################################################
########################### DEFINING FUNCTIONS ##########################

def lorentzian(x,p):
numerator =  (p[0]**2 )
denominator = ( x - (p[1]) )**2 + p[0]**2
y = p[2]*(numerator/denominator)
return y

def residuals(p,y,x):
err = y - lorentzian(x,p)
return err

#########################################################################
######################## BACKGROUND SUBTRACTION #########################

# defining the 'background' part of the spectrum #
ind_bg_low = (x > min(x)) & (x < 450.0)
ind_bg_high = (x > 590.0) & (x < max(x))

x_bg = numpy.concatenate((x[ind_bg_low],x[ind_bg_high]))
y_bg = numpy.concatenate((y[ind_bg_low],y[ind_bg_high]))
#pylab.plot(x_bg,y_bg)

# fitting the background to a line #
m, c = numpy.polyfit(x_bg, y_bg, 1)

# removing fitted background #
background = m*x + c
y_bg_corr = y - background
#pylab.plot(x,y_bg_corr)

#########################################################################
############################# FITTING DATA ## ###########################

# initial values #
p = [5.0,520.0,12e3]  # [hwhm, peak center, intensity] #

# optimization #
pbest = leastsq(residuals,p,args=(y_bg_corr,x),full_output=1)
best_parameters = pbest[0]

# fit to data #
fit = lorentzian(x,best_parameters)

#########################################################################
############################## PLOTTING #################################

pylab.plot(x,y_bg_corr,'wo')
pylab.plot(x,fit,'r-',lw=2)
pylab.xlabel(r'$\omega$ (cm$^{-1}$)', fontsize=18)
pylab.ylabel('Intensity (a.u.)', fontsize=18)

pylab.show()

Post contributed by Shrividya Ravi, VUW

5 thoughts on “9. Python Workshop I: Fitting a single symmetric peak”

1. Sebastian says:

Thanks for your explanation. But what if I don’t have a numerical set of data but a microscope image with a laser groove in it?
First I did otsu-segmentation and output the intensity profile of the groove with .

profile = im0.sum(axis=0) #
plt.figure()
plt.plot(profile)

I then get an image with a peak. How to determine its FWHM?

Thanks, Sebastian

• Mithil says:

Hi Sebastian,

Please post here if you could find how to deteremine FWHM.

@James, very useful documentation. Thanks.

-Mithil

2. Matthew says:

Thanks for this great example!

Just a (maybe not so minor) heads-up: The numerator of the Lorentzian in the example above is squared, but according to both MathWorld and Wikipedia it shouldn’t be.

3. Richard Dai says: