用于检测峰宽的鲁棒算法

Jes*_*ose 15 python scipy

在此输入图像描述

我问如何以编程方式判断频段@detly使用FWHM(半高全宽),以确定峰的宽度建议.我四处搜索,发现FWHM可用于拟合模型(我实际上是一个门外汉!),特别是高斯模型.具体来说,2.354 * sigma是高斯模型的宽度.

由于存在不良峰值,我正在寻找高斯适合.这张照片中有4个结构良好的峰.然后是"双峰"(虽然两者都很重要)和两个展开峰.他们将证明对天真的FWHM来说是一个不可能的挑战.

你能帮助生成Scip/Numpy峰值的guassian拟合(用于FWHM计算),因为它是x坐标的近似位置吗?如果Guassian是一个糟糕的选择,那么其他一些计划.

fra*_*xel 21

拟合高斯是一种很好的方法.如果你对初始参数值有一些好的猜测,你可以尝试一次猜测它们.一个很大的问题是噪声,实际上你可能想要孤立地拟合每个峰值(即,只考虑给定峰值一次的范围),或者在数据中获得基线噪声曲线并减去它.

以下是一些适合多个高斯人的代码.我已经输入了一些非常宽松的参数,我认为它是8个最突出的峰值,另外还有一个非常广泛的Guassian试图捕捉背景噪音.在处理之前,我稍微清理了您发布的图像,以帮助从中获取数据(移除鼠标光标和轴边缘,并反转图像).

在此输入图像描述

码:

import Image
from scipy import *
from scipy.optimize import leastsq
import numpy as np

im = Image.open("proc.png")
im = im.convert('L')
h, w = im.size
#Extract data from the processed image:
im = np.asarray(im)
y_vals, junk  = np.mgrid[w:0:-1, h:0:-1]
y_vals[im < 255] = 0
y_vals = np.amax(y_vals,axis=0)

def gaussian(x, A, x0, sig):
    return A*exp(-(x-x0)**2/(2.0*sig**2))

def fit(p,x):
    return np.sum([gaussian(x, p[i*3],p[i*3+1],p[i*3+2]) 
                   for i in xrange(len(p)/3)],axis=0)

err = lambda p, x, y: fit(p,x)-y

#params are our intitial guesses for fitting gaussians, 
#(Amplitude, x value, sigma):
params = [[50,40,5],[50,110,5],[100,160,5],[100,220,5],
          [50,250,5],[100,260,5],[100,320,5], [100,400,5],   
          [30,300,150]]  # this last one is our noise estimate
params = np.asarray(params).flatten()

x  = xrange(len(y_vals))
results, value = leastsq(err, params, args=(x, y_vals))

for res in results.reshape(-1,3):
    print "amplitude, position, sigma", res

import pylab
pylab.subplot(211, title='original data')
pylab.plot(y_vals)
pylab.subplot(212, title='guassians fit')
y = fit(results, x)
pylab.plot(x, y)
pylab.savefig('fig2.png')
pylab.show()
Run Code Online (Sandbox Code Playgroud)

这些是拟合输出高斯参数:

#Output:
amplitude, position, sigma [ 23.47418352  41.27086097   5.91012897]
amplitude, position, sigma [  16.26370263  106.39664543    3.67827824]
amplitude, position, sigma [  59.74500239  163.11210316    2.89866545]
amplitude, position, sigma [  57.91752456  220.24444939    2.91145375]
amplitude, position, sigma [  39.78579841  251.25955921    2.74646334]
amplitude, position, sigma [  86.50061756  260.62004831    3.08367483]
amplitude, position, sigma [  79.26849867  319.64343319    3.57224402]
amplitude, position, sigma [ 127.97009966  399.27833126    3.14331212]
amplitude, position, sigma [  20.21224956  379.41066063  195.47122312]
Run Code Online (Sandbox Code Playgroud)

  • @BandGap - 首先我使用Gimp编辑图像以移除鼠标光标和轴边缘,然后将其反转并对图像进行阈值处理以使其成为二进制.然后在代码中提取数据的#Extract数据中的4行提取数据:通过创建行值网格,如果图像中相应的位置为零则将所有这些值设置为零,则需要列中的最大行索引,导致第一个绘图"原始数据"对应于已发布的图像. (5认同)