Подтвердить что ты не робот

Надежный алгоритм определения ширины пиков

enter image description here

Я спросил как программно судить диапазоны спектра и @detly предложил использовать FWHM (полная ширина в половине максимум) для определения ширины пиков. Я обыскал вокруг и обнаружил, что FWHM можно использовать для подгонки моделей (я на самом деле это непрофессионал!), Особенно модели Guassian. В частности, 2.354 * sigma ширина для модели Guassian.

Я смотрю на гуасианскую подгонку из-за наличия плохих пиков. На этой картинке есть 4 хорошо сформированных пика. Тогда есть "двойной" пик (хотя оба они могут быть важны) и два разбросанных пика. Они докажут невозможный вызов наивному FWHM.

Можете ли вы помочь в создании гуасианского фитинга (для расчета FWHM) пика в Scip/Numpy, учитывая приблизительное местоположение его координаты x? Если гуасиан - это плохой выбор, то какая-то другая схема.

4b9b3361

ОТВЕТЫ

Ответ 1

Приспособление гауссианцев - хороший подход. И если у вас есть окие догадки к начальным значениям параметров, вы можете попробовать и угадать их все сразу. Большая проблема - шум, действительно, вы, вероятно, хотите либо подогнать каждый пик в изоляции (т.е. Учитывать только диапазон, за который данный пик за раз), либо получить кривую базовой линии на основе ваших данных и вычесть ее.

Вот некоторый код, который пытается подогнать несколько гауссиан. Я ввел некоторые довольно рыхлые параметры для того, что я считал 8 наиболее известными пиками, плюс один дополнительный очень широкий гуассиан, чтобы попытаться захватить фоновый шум. Перед обработкой я немного очистил ваше размещенное изображение, чтобы помочь получить данные из него (удалил курсор мыши и края оси и перевернул изображение).

enter image description here

код:

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()

Это настраиваемые выходные параметры Guassian:

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