Поиск полной ширины максимума пика
Я пытался выяснить половину максимума ширины (FWHM) синего пика (см. изображение). Зеленый пик и пик пурпурного цвета составляют синий пик. Я использовал следующее уравнение, чтобы найти FWHM пиков зеленого и пурпурного цветов: fwhm = 2*np.sqrt(2*(math.log(2)))*sd
где sd = стандартное отклонение. Я создал зеленый и пурпурный пики, и я знаю стандартное отклонение, поэтому я могу использовать это уравнение.
Я создал зеленый и пурпурный пики, используя следующий код:
def make_norm_dist(self, x, mean, sd):
import numpy as np
norm = []
for i in range(x.size):
norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))]
return np.array(norm)
Если бы я не знал, что синий пик состоял из двух пиков, и у меня был только синий пик в моих данных, как бы я нашел FWHM?
![]()
Я использую этот код, чтобы найти вершину пика:
peak_top = 0.0e-1000
for i in x_axis:
if i > peak_top:
peak_top = i
Я мог бы разделить peak_top
на 2, чтобы найти половину высоты, а затем попытаться найти y-значения, соответствующие половине высоты, но тогда я столкнулся бы с проблемой, если нет значений x, точно совпадающих с половинной высотой,
Я уверен, что есть более элегантное решение, которое я пытаюсь.
Ответы
Ответ 1
Вы можете использовать сплайн для соответствия [синей кривой - пик /2], а затем найти его корнями:
import numpy as np
from scipy.interpolate import UnivariateSpline
def make_norm_dist(x, mean, sd):
return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2))
x = np.linspace(10, 110, 1000)
green = make_norm_dist(x, 50, 10)
pink = make_norm_dist(x, 60, 10)
blue = green + pink
# create a spline of x and blue-np.max(blue)/2
spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0)
r1, r2 = spline.roots() # find the roots
import pylab as pl
pl.plot(x, blue)
pl.axvspan(r1, r2, facecolor='g', alpha=0.5)
pl.show()
Вот результат:
![enter image description here]()
Ответ 2
Это сработало для меня в iPython (быстрый и грязный, можно уменьшить до 3 строк):
def FWHM(X,Y):
half_max = max(Y) / 2.
#find when function crosses line half_max (when sign of diff flips)
#take the 'derivative' of signum(half_max - Y[])
d = sign(half_max - array(Y[0:-1])) - sign(half_max - array(Y[1:]))
#plot(X,d) #if you are interested
#find the left and right most indexes
left_idx = find(d > 0)[0]
right_idx = find(d < 0)[-1]
return X[right_idx] - X[left_idx] #return the difference (full width)
Некоторые дополнения могут быть сделаны, чтобы сделать разрешение более точным, но в пределе, что на оси X имеется много выборок, а данные не слишком шумные, это отлично работает.
Даже когда данные не являются гауссовыми и немного шумными, это сработало для меня (я просто беру первый и последний раз половину макс, пересекая данные).
Ответ 3
Если ваши данные имеют шум (и это всегда происходит в реальном мире), более надежным решением было бы привязать гауссов к данным и извлечь из них FWHM:
import numpy as np
import scipy.optimize as opt
def gauss(x, p): # p[0]==mean, p[1]==stdev
return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))
# Create some sample data
known_param = np.array([2.0, .7])
xmin,xmax = -1.0, 5.0
N = 1000
X = np.linspace(xmin,xmax,N)
Y = gauss(X, known_param)
# Add some noise
Y += .10*np.random.random(N)
# Renormalize to a proper PDF
Y /= ((xmax-xmin)/N)*Y.sum()
# Fit a guassian
p0 = [0,1] # Inital guess is a normal distribution
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y))
fit_mu, fit_stdev = p1
FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
print "FWHM", FWHM
![enter image description here]()
Настроенное изображение может быть сгенерировано с помощью:
from pylab import *
plot(X,Y)
plot(X, gauss(X,p1),lw=3,alpha=.5, color='r')
axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
show()
Еще лучшее приближение отфильтровывает шумные данные ниже заданного порога до пригонки.
Ответ 4
Вот небольшая функция, использующая сплайн-подход.
from scipy.interpolate import splrep, sproot, splev
def fwhm(x, y, k=10):
"""
Determine full-with-half-maximum of a peaked set of points, x and y.
Assumes that there is only one peak present in the datasset. The function
uses a spline interpolation of order k.
"""
class MultiplePeaks(Exception): pass
class NoPeaksFound(Exception): pass
half_max = amax(y)/2.0
s = splrep(x, y - half_max)
roots = sproot(s)
if len(roots) > 2:
raise MultiplePeaks("The dataset appears to have multiple peaks, and "
"thus the FWHM can't be determined.")
elif len(roots) < 2:
raise NoPeaksFound("No proper peaks were found in the data set; likely "
"the dataset is flat (e.g. all zeros).")
else:
return abs(roots[1] - roots[0])