python - 我可以拟合质心和峰高随机分布的多个高斯谱吗?
问题描述
新手在这里,但我已经尝试在发布之前进行尽职调查。为任何无意的失礼道歉。
我正在以电压与时间序列的形式从示波器获取数据。时间箱的宽度为 0.8 纳秒。我运行多个“数据捕获”周期。单次捕获将具有固定数量的样本,以及 5 到 15 个高斯峰,确切的峰数未知。高斯峰具有相对受限的 FWHM(在 2 到 3 纳秒之间)、变化的峰高和随机到达时间(即质心位置不是周期性的)。
我一直在使用 Python 对这些数据进行高斯拟合,并且使用 scipy.optimise 库和 astropy 库取得了一些成功。下面包含使用 scipy.optimise 的代码。我可以拟合多个高斯,但我的代码中的一个关键步骤是提供峰数的“猜测”,并为每个峰估计质心位置、峰高和峰宽。有没有办法概括这段代码而不必提供“猜测”?如果我放松“猜测”中的条件,那么拟合就会失去质量。我知道峰将是宽度受限的高斯峰,但希望概括代码以适应任何给定数据捕获中的峰质心和峰高。
import ctypes
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#Get data from file
with open('test3.txt') as f:
w, h = [float(x) for x in next(f).split()]
print(w, h)
array = [[float(x) for x in line.split()] for line in f]
#Separate
x, z = zip(*array)
#Change sign since fitting routine seems to
#prefer positive numbers
y=[ -p for p in z]
def func(x, *params):
y = np.zeros_like(x)
for i in range(0, len(params), 3):
ctr = params[i]
amp = params[i+1]
wid = params[i+2]
y = y + amp * np.exp( -((x - ctr)/wid)**2)
return y
#Guess the peak positions, heights, and widths
guess = [16, 5, 2, 75, 5, 2, 105, 5, 2, 139, 5, 2, 225, 5, 2, 315, 5, 2, 330, 5, 2]
#Fit and print parameters to screen
popt, pcov = curve_fit(func, x, y, p0=guess)
print(popt)
fit = func(x, *popt)
#Plot stuff
plt.plot(x, y)
plt.plot(x, fit , 'r-')
plt.show()
结果如下所示: 数据和拟合图
数据文件在这里:https ://spaces.hightail.com/receive/5MY7Vc7r9R
这类似于如何在 Python 中将多个高斯曲线拟合到质谱数据?并将多个高斯拟合到 python 中的数据,但这两个依赖于拟合周期性数据集或具有已知峰值位置、宽度和高度的数据集。他们让我走到这一步很有用,但我现在被困住了。有什么想法或建议我可以跟进吗?
谢谢,阿迪
解决方案
我的想法是我们将曲线的值与其平均值进行比较。
该multiplier
变量表示该值必须大于平均值多少倍才能让我们了解这是峰值之一。超过该值的峰值的第一个点被认为是逼近该峰值平均值的起点。
我还将列表替换为 x 和 y 的数组。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#Get data from file
with open('test3.txt') as f:
w, h = [float(x) for x in next(f).split()]
array = [[float(x) for x in line.split()] for line in f]
#Separate
x, z = zip(*array)
x = np.array(x)
y = np.array([ -p for p in z])
#Change sign since fitting routine seems to
#prefer positive numbers
def func(x, *params):
y = np.zeros_like(x)
for i in range(0, len(params), 3):
ctr = params[i]
amp = params[i+1]
wid = params[i+2]
y = y + amp * np.exp( -((x - ctr)/wid)**2)
return y
#Guess the peak positions, heights, and widths
# guess = [16, 5, 2, 75, 5, 2, 105, 5, 2, 139, 5, 2, 225, 5, 2, 315, 5, 2, 330, 5, 2]
def getPeaks(x, y, multiplier):
x_peaks = []
isPeak = False
for i, j in zip(x, y):
if j > y.mean() * multiplier:
if not isPeak:
isPeak = True
x_peaks.append(i)
else:
isPeak = False
return x_peaks
multiplier = 3
x_peaks = getPeaks(x, y, multiplier)
guess = []
for i in range(len(x_peaks)):
guess.append(x_peaks[i])
guess.append(5)
guess.append(2)
#Fit and print parameters to screen
popt, pcov = curve_fit(func, x, y, p0=guess)
print(popt)
fit = func(x, *popt)
#Plot stuff
plt.plot(x, y)
plt.plot(x, fit , 'r--')
# plt.plot(popt[::3], np.ones_like(popt[::3]) * multiplier, 'ko')
plt.show()
推荐阅读
- java - 使用 @SpringAware MapLoader 正确使用 Spring Boot Hazelcast 自动配置
- javascript - discord.js if date.now() 函数
- javascript - 如何获取 async 的值并传递给函数?
- string - 是否可以批量操作文件名?
- typescript - 如何从您的 Ionic 应用程序中访问任何自定义端口
- angular - 从服务调用组件 - 角度
- html - 如何在降价中禁止列表项之间的换行
- oracle - 在 WINDOWS 10 机器上安装 ORACLE 11 G
- java - 获取 CtClass 超类名称或 CtClass 超类实例
- python - 保存具有当前日期和时间的 csv