求峰值的全宽度半最大值

Question

我一直在试图计算出蓝峰的半高宽(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) 

如果我不知道蓝峰是由两个峰组成的，而我的数据中只有蓝峰，我如何找到半高宽呢？

​

​

我一直在使用下面的代码来寻找峰顶：

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

我可以将peak_top除以2来找到半高，然后尝试找到对应于半高的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

这对我在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[0:len(d)],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轴有许多样本且数据不太嘈杂的限制下，这非常有效。

即使数据不是高斯的，而且有一点噪声，它对我也是有效的(我只需要第一次和最后一次半最大值通过数据)。

Answer 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

​

​

打印的图像可以通过以下方式生成：

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

一种更好的近似方法是在拟合之前过滤掉低于给定阈值的噪声数据。