找到一個峯的半峯全寬

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軸有很多樣本，數據不會太嘈雜，這很好。

即使數據不是高斯和有點嘈雜，它爲我工作（我只是採取第一次和最後一次半最大穿越數據）。