[Python] Fitting multiple Lorentzians to Brillouin Spectrum using Scipy in Python 3

I am trying to fit Brillouin Spectra (with several peaks) using scipy.optimize.curve_fit. I have have multiple spectra with several peaks and I am trying to fit them with lorentzian functions (one Lorentzian per peak). I am trying to automate the process for bulk analysis (i.e., using the peak finding algorithm of scipy to get peak positions, peak widths and peaks heights and use them as initial guesses for the fit). I am now working on one spectrum to see if the general idea works, then I will extend it to be automatic and work with all the spectra I have. So far, I have done this:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit

#import y data from linked google sheet
y_data = np.loadtxt( 'y_peo.tsv' )
#define x data
x_data = np.arange( len( y_data ) )

#find peak properties (peak position, amplitude, full width half maximum ) to use as
#initial guesses for the curve_fit function
pk, properties = find_peaks(y_data, height = 3, width = 3, prominence=0.1 ) #pk returns peaks position, properties returns
#other properties associated with the peaks
I = properties ['peak_heights'] #amplitude
fwhm = (properties['widths']) #full width half maximum

#define sum of lorentzians
def func(x, *params):
y = np.zeros_like(x)
for i in range(0, len(params), 3):
x_0 = params
I = params[i+1]
y_0 = params[i+2]
y = y + (I*y_0**2)/(y_0**2 +(x-x_0)**2)
return y

#initial guess list, to be populated with parameters found above (pk, I, fwhm)
guess = []

for i in range(len(pk)):
guess.append(pk)
guess.append(I)
guess.append(fwhm)

#convert list to array
guess=np.array(guess)

#fit
popt, pcov = curve_fit(func, x_data, y_data, p0=guess, method = 'lm', maxfev=1000000)
print(popt)
fit = func(x_data, *popt)

#plotting
fig= plt.figure(figsize=(10,5))
ax= fig.add_subplot(1,1,1)
ax.plot(pk, y_data[pk], 'o', ms=5)
ax.plot(x_data, y_data, 'ok', ms=1)
ax.plot(x_data, fit , 'r--', lw=0.5)


Where y_peo is the dependent variable (here are the values as a google sheet file: https://docs.google.com/spreadsheets/d/1UB2lqs0Jmhthed7B0U9rznqcbpEMRmRX99c-iOBHLeE/edit?usp=sharing) and pixel is the independent variable (an arbitrary array of values created in the code). This is what I get: Result of spectrum fit. Any idea why the last triplet is not fitted as expected? I checked and all the peaks are found correctly by the scipy.signal.find_peaks function - hence I think the problem is with the scipy.optimize.curve_fit, since I had to increase the number of maxfev for the fit to 'work'. Any idea on how to tackle this problem in a smarter way?

Thanks in advance,

Giuseppe.

Continue reading...