Tìm chiều rộng tối đa một nửa đầy đủ của một đỉnh

Question

Tôi đã cố gắng để tìm ra đầy đủ chiều rộng tối đa một nửa (FWHM) của đỉnh màu xanh (xem ảnh). Đỉnh xanh và đỉnh đỏ tươi kết hợp tạo nên đỉnh màu xanh dương. Tôi đã sử dụng phương trình sau để tìm FWHM của các đỉnh xanh và đỏ tươi: fwhm = 2*np.sqrt(2*(math.log(2)))*sd trong đó sd = độ lệch chuẩn. Tôi tạo ra các đỉnh xanh và đỏ tươi và tôi biết độ lệch chuẩn là lý do tại sao tôi có thể sử dụng phương trình đó.

Tôi tạo ra các đỉnh núi xanh và đỏ tươi bằng cách sử dụng đoạn mã sau:

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) 

Nếu tôi không biết đỉnh màu xanh được tạo thành từ hai đỉnh và tôi chỉ có đỉnh màu xanh trong dữ liệu của tôi, làm thế nào sẽ Tôi tìm thấy FWHM?

Tôi đã được sử dụng mã này để tìm đỉnh đầu:

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

tôi có thể chia peak_top bởi 2 để tìm chiều cao nửa và sau đó cố gắng và tìm giá trị y tương ứng với chiều cao một nửa, nhưng sau đó tôi sẽ gặp rắc rối nếu không có giá trị x chính xác phù hợp với chiều cao một nửa.

Tôi khá chắc chắn có một giải pháp thanh lịch hơn với người tôi đang cố gắng.

Answer 1

Bạn có thể sử dụng spline để phù hợp với [xanh đường cong - đỉnh cao/2], và sau đó tìm thấy nó là rễ:

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

Dưới đây là kết quả:

Answer 2

Nếu dữ liệu của bạn có tiếng ồn (và nó luôn luôn trong thế giới thực), một giải pháp mạnh mẽ hơn là để phù hợp với Gaussian với dữ liệu và giải nén FWHM từ đó:

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 

Hình ảnh vẽ có thể được tạo ra bởi:

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

Một xấp xỉ thậm chí tốt hơn sẽ lọc ra các dữ liệu ồn ào dưới một ngưỡng cho trước phù hợp.

Answer 3

Dưới đây là một chức năng nhỏ đẹp bằng cách sử dụng phương pháp 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

này đã làm việc cho tôi trong ipython (nhanh chóng và dơ bẩn, có thể được giảm xuống còn 3 lines):

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) 

Một số bổ sung có thể được thực hiện để làm cho độ phân giải chính xác hơn, nhưng trong giới hạn mà có nhiều mẫu dọc theo trục X và dữ liệu không quá ồn ào, điều này hoạt động rất tốt.

Ngay cả khi dữ liệu không phải là Gaussian và một chút ồn ào, nó làm việc cho tôi (tôi chỉ dành thời gian nửa đầu và kéo dài tối đa vượt qua các dữ liệu).