Encontrar el ancho máximo medio de un pico

Question

He estado tratando de calcular el ancho medio máximo (FWHM) del pico azul (ver imagen). El pico verde y el pico magenta combinados conforman el pico azul. He estado usando la siguiente ecuación para encontrar el FWHM de los picos verde y magenta: fwhm = 2*np.sqrt(2*(math.log(2)))*sd donde sd = desviación estándar. Creé los picos verde y magenta y sé la desviación estándar, que es la razón por la que puedo usar esa ecuación.

que crearon los picos verdes y magenta utilizando el siguiente código:

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) 

Si no sabía el pico azul estaba compuesto de dos picos y sólo tenía el pico azul en mis datos, ¿cómo Encuentro el FWHM?

He estado usando este código para encontrar el pico superior:

peak_top = 0.0e-1000 
for i in x_axis: 
    if i > peak_top: 
     peak_top = i 

pude dividir el peak_top por 2 para encontrar la mitad de la altura y luego tratar de encontrar los valores de y correspondientes a la media altura, pero luego me encontraría en problemas si no hay valores de x que coincidan exactamente con la mitad de la altura.

Estoy bastante seguro de que hay una solución más elegante para la que estoy tratando.

Answer 1

Puede utilizar spline para adaptarse a la [curva azul - pico/2], y luego encontrar sus raíces:

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

aquí está el resultado:

Answer 2

Si los datos tiene ruido (y siempre lo hace en el mundo real), una solución más sólida sería ajustar un Gaussian a los datos y extraer FWHM de eso:

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 

La imagen trazada se puede generar por:

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

Una mejor aproximación sería filtrar los datos ruidosos por debajo de un umbral dado antes del ajuste.

Answer 3

Aquí hay una pequeña y agradable función que utiliza el enfoque spline.

from scipy.interpolate import splrep, sproot, splev 

class MultiplePeaks(Exception): pass 
class NoPeaksFound(Exception): pass 

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. 
    """ 

    half_max = amax(y)/2.0 
    s = splrep(x, y - half_max, k=k) 
    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]) 

Answer 4

Esto funcionó para mí en IPython (rápido y sucio, se puede reducir a 3 líneas):

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) 

Algunas adiciones se pueden hacer para que la resolución más precisa, pero en el límite que hay muchas muestras a lo largo del eje X y los datos no son demasiado ruidosos, esto funciona muy bien.

Incluso cuando los datos no son gaussianos y un poco ruidosos, me funcionaron (simplemente tomo la primera y última vez que half max cruza los datos).