使用Python 3中的Scipy将多个洛伦兹族拟合到布里渊谱 [英] Fitting multiple Lorentzians to Brillouin Spectrum using Scipy in Python 3

问题描述

我正在尝试使用scipy.optimize.curve_fit拟合布里渊光谱(具有多个峰).我有多个具有多个峰的光谱，并且正在尝试用洛伦兹函数(每个峰一个洛伦兹函数)拟合它们.我正在尝试自动化进行批量分析的过程(即使用scipy的峰发现算法来获取峰位置，峰宽度和峰高，并将其用作拟合的初始猜测).现在，我正在研究一个频谱，以了解总体思路是否可行，然后将其扩展为自动的，并可以使用我拥有的所有频谱.到目前为止，我已经做到了:

 将numpy导入为np导入matplotlib.pyplot作为plt从scipy.signal导入find_peaks从scipy.optimize导入curve_fit#import链接的Google工作表中的数据y_data = np.loadtxt('y_peo.tsv')#定义x数据x_data = np.arange(len(y_data))#寻找峰值属性(峰值位置，幅度，半峰全宽)用作#curve_fit函数的初始猜测pk，属性= find_peaks(y_data，高度= 3，宽度= 3，突出度= 0.1)#pk返回峰的位置，属性返回#与峰相关的其他属性I =属性['peak_heights']#振幅fwhm =(properties ['widths'])#全角一半最大值#定义洛伦兹人的总和def func(x，* params):y = np.zeros_like(x)对于范围(0，len(params)，3)中的i:x_0 =参数[i]I =参数[i + 1]y_0 =参数[i + 2]y = y +(I * y_0 ** 2)/(y_0 ** 2 +(x-x_0)** 2)返回y#initial猜测列表，将使用上面找到的参数填充(pk，I，fwhm)猜测= []对于我在范围内(len(pk)):guess.append(pk [i])guess.append(I [i])guess.append(fwhm [i])#将列表转换为数组guess = np.array(猜测)#合身popt，pcov = curve_fit(函数，x_data，y_data，p0 =猜测，方法='lm'，maxfev = 1000000)打印(弹出)适合= func(x_data，* popt)#绘图无花果= 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) 

其中 y_peo 是因变量(此处是作为Google工作表文件的值:

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]
        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[i])
    guess.append(I[i])
    guess.append(fwhm[i]) 

#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.

解决方案

With some small modifications the code of the OP runs just fine. Now it looks like:

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

def lorentzian( x, x0, a, gam ):
    return a * gam**2 / ( gam**2 + ( x - x0 )**2 )

def multi_lorentz( x, params ):
    off = params[0]
    paramsRest = params[1:]
    assert not ( len( paramsRest ) % 3 )
    return off + sum( [ lorentzian( x, *paramsRest[ i : i+3 ] ) for i in range( 0, len( paramsRest ), 3 ) ] )

def res_multi_lorentz( params, xData, yData ):
    diff = [ multi_lorentz( x, params ) - y for x, y in zip( xData, yData ) ]
    return diff

y0 = np.loadtxt( 'y_peo.tsv' )
yData = np.loadtxt( 'y_peo.tsv' )
xData = np.arange( len( yData ) )

yGround = min( yData )
yData = yData - yGround
yAmp = max( yData )
yData = yData / yAmp 

#initial properties of peaks 
pk, properties = find_peaks( yData, height = .05, width = 3 )#, prominence=0.1 )
#extract peak heights and fwhm 
I = properties [ 'peak_heights' ]
fwhm = properties[ 'widths' ]

guess = [0]
for i in range( len( pk ) ): 
    guess.append( pk[i] )
    guess.append( I[i] )
    guess.append( fwhm[i] ) 

guess=np.array( guess )

popt, pcov = leastsq( res_multi_lorentz, x0=guess, args=( xData, yData ) )
print( popt )


testData = [ multi_lorentz( x, popt ) for x in xData ]
fitData = [ yGround + yAmp * multi_lorentz( x, popt ) for x in xData ]

fig= plt.figure( figsize=( 10, 5 ) )
ax= fig.add_subplot( 2, 1, 1 )
bx= fig.add_subplot( 2, 1, 2 )
ax.plot( pk, yData[pk], 'o', ms=5 )
ax.plot( xData, yData, 'ok', ms=1 )
ax.plot( xData, testData , 'r--', lw=0.5 )
bx.plot( xData, y0, ls='', marker='o', markersize=1 )
bx.plot( xData, fitData )
plt.show()

giving

[9.39456104e-03 6.55864388e+01 5.73522507e-02 5.46727721e+00
 1.21329586e+02 2.97187567e-01 2.12738107e+00 1.76823266e+02
 1.57244131e-01 4.55424037e+00 3.51556692e+02 3.08959955e-01
 4.39114581e+00 4.02954496e+02 9.02677035e-01 2.27647259e+00
 4.53994668e+02 3.74379310e-01 4.02432791e+00 6.15694190e+02
 4.51943494e-01 4.06522919e+00 6.63307635e+02 1.05793098e+00
 2.32168786e+00 7.10644233e+02 4.33194434e-01 4.17050014e+00
 8.61276198e+02 1.79240633e-01 4.26617114e+00 9.06211127e+02
 5.03070470e-01 2.10563379e+00 9.50973864e+02 1.78487912e-01
 4.01254815e+00]

and

这篇关于使用Python 3中的Scipy将多个洛伦兹族拟合到布里渊谱的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！