我正在尝试使用 scipy.optimize.curve_fit 拟合布里渊光谱(有几个峰)。我有多个具有多个峰的光谱,我正在尝试将它们与洛伦兹函数(每个峰一个洛伦兹)拟合。我正在尝试自动化批量分析的过程(即,使用 scipy 的峰值查找算法来获取峰值位置、峰值宽度和峰值高度,并将它们用作拟合的初始猜测)。我现在正在研究一个光谱,看看总体思路是否有效,然后我会将其扩展为自动并使用我拥有的所有光谱。到目前为止,我已经这样做了:
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)
其中y_peo是因变量(这里是谷歌表格文件中的值:https ://docs.google.com/spreadsheets/d/1UB2lqs0Jmhthed7B0U9rznqcbpEMRmRX99c-iOBHLeE/edit?usp=sharing )和像素是自变量(任意在代码中创建的值数组)。这就是我得到的:频谱拟合的结果。知道为什么最后一个三胞胎没有按预期安装吗?我检查并通过 scipy.signal.find_peaks 函数正确找到了所有峰值 - 因此我认为问题出在 scipy.optimize.curve_fit 上,因为我必须增加 maxfev 的数量才能“工作”。关于如何以更聪明的方式解决这个问题的任何想法?
慕田峪4524236
相关分类