找到一个峰的半峰全宽

Question

我一直在试图找出蓝峰的全峰半峰值（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 

我可以除以2 peak_top找到一半的高度，然后试图找到对应的y值半高，但如果没有完全匹配半高的x值，我会遇到麻烦。

我很确定有一个更优雅的解决方案，我正在尝试。

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

下面是结果：

Answer 2

如果您的数据有噪声（并且它总是在现实世界中），更可靠的解决方案是将高斯拟合为数据并从中提取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 

所绘制的图像可以通过产生：

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

一个更好的近似拟合之前将过滤出低于给定的阈值的噪声数据。

Answer 3

这是一个很好的使用样条方法的小函数。

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

这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轴有很多样本，数据不会太嘈杂，这很好。

即使数据不是高斯和有点嘈杂，它为我工作（我只是采取第一次和最后一次半最大穿越数据）。