28

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

我可以将除以peak_top2 以找到半高，然后尝试找到与半高相对应的 y 值，但是如果没有与半高完全匹配的 x 值，我会遇到麻烦。

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

4

8 回答 8

20

您可以使用样条曲线拟合[蓝色曲线 - 峰值/ 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()

结果如下：

在此处输入图像描述

于 2012-05-14T12:55:38.990 回答
15

这在 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 轴有很多样本并且数据不太嘈杂的限制下，这很有效。

即使数据不是高斯的并且有点嘈杂，它也对我有用（我只是第一次和最后一次半最大值穿过数据）。

于 2013-05-10T20:00:20.833 回答
12

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

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

于 2012-05-15T16:33:06.650 回答
8

这是使用样条方法的一个不错的小功能。

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])
于 2013-01-14T22:16:38.607 回答
2

您应该使用 scipy 来解决它：首先find_peaks然后peak_widths。使用rel_height (0.5)中的默认值，您正在测量峰的半峰宽度。

于 2019-11-15T14:06:37.420 回答
0

对于具有许多数据点的单调函数，如果不需要完美的准确性，我会使用：

def FWHM(X, Y):
    deltax = x[1] - x[0]
    half_max = max(Y) / 2.
    l = np.where(y > half_max, 1, 0)

    return np.sum(l) * deltax
于 2021-09-21T07:03:58.603 回答
0

我实现了一个经验解决方案，该解决方案在haggis.math.full_width_half_max. 用法非常简单：

fwhm = full_width_half_max(x, y)

该函数是稳健的：它只是使用请求的插值方案找到数据的最大值和越过“中途下降”阈值的最近点。

以下是使用其他答案中的数据的几个示例。

@HYRY的流畅数据

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

# Compute using my function
fwhm, (x1, y1), (x2, y2) = full_width_half_max(x, blue, return_points=True)

# Print comparison
print('HYRY:', r2 - r1, 'MP:', fwhm)

plt.plot(x, blue)
plt.axvspan(r1, r2, facecolor='g', alpha=0.5)
plt.plot(x1, y1, 'r.')
plt.plot(x2, y2, 'r.')

对于平滑数据，结果非常准确：

HYRY: 26.891157007233254 MP: 26.891193606203814

在此处输入图像描述

@Hooked 的嘈杂数据

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

# Compute using my function
fwhm, (x1, y1), (x2, y2) = full_width_half_max(X, Y, return_points=True)

# Print comparison
print('Hooked:', FWHM, 'MP:', fwhm)

plt.plot(X, Y)
plt.plot(X, gauss(X, p1), lw=3, alpha=.5, color='r')
plt.axvspan(fit_mu - FWHM / 2, fit_mu + FWHM / 2, facecolor='g', alpha=0.5)
plt.plot(x1, y1, 'r.')
plt.plot(x2, y2, 'r.')

对于嘈杂的数据（基线有偏差），结果并不一致。

Hooked: 1.9903193212254346 MP: 1.5039676990530118

一方面，高斯拟合对于数据不是很理想，但另一方面，选择与半最大值阈值相交的最近点的策略也可能不是最优的。

在此处输入图像描述

于 2022-01-10T09:48:07.973 回答
0

如果您更喜欢插值而不是拟合：

import numpy as np

def get_full_width(x: np.ndarray, y: np.ndarray, height: float = 0.5) -> float:
    height_half_max = np.max(y) * height
    index_max = np.argmax(y)
    x_low = np.interp(height_half_max, y[:index_max], x[:index_max])
    x_high = np.interp(height_half_max, np.flip(y[index_max:]), np.flip(x[index_max:]))

    return x_high - x_low
于 2022-02-09T16:20:13.090 回答