matplotlib 求峰的半高全宽

u4vypkhs  于 2023-08-06  发布在  其他
关注(0)|答案(9)|浏览(173)

我一直试图计算出蓝色峰值的半高宽(FWHM)(见图)。绿色峰值和洋红峰值组合构成蓝色峰值。我一直在使用以下公式来计算绿色和洋红峰值的半高宽:其中sd =标准偏差。我创建了绿色和洋红色的峰值,我知道标准差,这就是为什么我可以使用这个方程。
我使用以下代码创建了绿色和洋红色峰值:

  1. def make_norm_dist(self, x, mean, sd):
  2. import numpy as np
  3. norm = []
  4. for i in range(x.size):
  5. norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))]
  6. return np.array(norm)

字符串
如果我不知道蓝峰是由两个峰组成的,并且我的数据中只有蓝峰,那么我如何找到FWHM?


的数据
我一直在用这段代码来寻找峰顶:

  1. peak_top = 0.0e-1000
  2. for i in x_axis:
  3. if i > peak_top:
  4. peak_top = i


我可以将peak_top除以2以找到半高,然后尝试找到对应于半高的y值,但是如果没有x值完全匹配半高,我就会遇到麻烦。
我很确定有一个更优雅的解决方案,我正在尝试的一个。

omtl5h9j

omtl5h9j1#

您可以使用样条曲线拟合[blue curve - peak/2],然后找到它的根:

  1. import numpy as np
  2. from scipy.interpolate import UnivariateSpline
  3. def make_norm_dist(x, mean, sd):
  4. return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2))
  5. x = np.linspace(10, 110, 1000)
  6. green = make_norm_dist(x, 50, 10)
  7. pink = make_norm_dist(x, 60, 10)
  8. blue = green + pink
  9. # create a spline of x and blue-np.max(blue)/2
  10. spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0)
  11. r1, r2 = spline.roots() # find the roots
  12. import pylab as pl
  13. pl.plot(x, blue)
  14. pl.axvspan(r1, r2, facecolor='g', alpha=0.5)
  15. pl.show()

字符串
结果如下:
x1c 0d1x的数据

展开查看全部
nhaq1z21

nhaq1z212#

这在iPython中对我有效(快速而肮脏,可以减少到3行):

  1. def FWHM(X,Y):
  2. half_max = max(Y) / 2.
  3. #find when function crosses line half_max (when sign of diff flips)
  4. #take the 'derivative' of signum(half_max - Y[])
  5. d = sign(half_max - array(Y[0:-1])) - sign(half_max - array(Y[1:]))
  6. #plot(X[0:len(d)],d) #if you are interested
  7. #find the left and right most indexes
  8. left_idx = find(d > 0)[0]
  9. right_idx = find(d < 0)[-1]
  10. return X[right_idx] - X[left_idx] #return the difference (full width)

字符串
可以进行一些添加以使分辨率更准确,但是在X轴沿着有许多样本并且数据没有太多噪声的限制下,这很有效。
即使当数据不是高斯的并且有点嘈杂时,它也对我有效(我只是取第一次和最后一次半最大值穿过数据)。

ltskdhd1

ltskdhd13#

如果您的数据有噪声(在真实的世界中总是如此),则更可靠的解决方案是将高斯拟合到数据并从中提取FWHM:

  1. import numpy as np
  2. import scipy.optimize as opt
  3. def gauss(x, p): # p[0]==mean, p[1]==stdev
  4. return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))
  5. # Create some sample data
  6. known_param = np.array([2.0, .7])
  7. xmin,xmax = -1.0, 5.0
  8. N = 1000
  9. X = np.linspace(xmin,xmax,N)
  10. Y = gauss(X, known_param)
  11. # Add some noise
  12. Y += .10*np.random.random(N)
  13. # Renormalize to a proper PDF
  14. Y /= ((xmax-xmin)/N)*Y.sum()
  15. # Fit a guassian
  16. p0 = [0,1] # Inital guess is a normal distribution
  17. errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
  18. p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y))
  19. fit_mu, fit_stdev = p1
  20. FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
  21. print "FWHM", FWHM

字符串


的数据
可通过以下方式生成打印图像:

  1. from pylab import *
  2. plot(X,Y)
  3. plot(X, gauss(X,p1),lw=3,alpha=.5, color='r')
  4. axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
  5. show()


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

展开查看全部
lc8prwob

lc8prwob4#

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

  1. from scipy.interpolate import splrep, sproot, splev
  2. class MultiplePeaks(Exception): pass
  3. class NoPeaksFound(Exception): pass
  4. def fwhm(x, y, k=10):
  5. """
  6. Determine full-with-half-maximum of a peaked set of points, x and y.
  7. Assumes that there is only one peak present in the datasset. The function
  8. uses a spline interpolation of order k.
  9. """
  10. half_max = amax(y)/2.0
  11. s = splrep(x, y - half_max, k=k)
  12. roots = sproot(s)
  13. if len(roots) > 2:
  14. raise MultiplePeaks("The dataset appears to have multiple peaks, and "
  15. "thus the FWHM can't be determined.")
  16. elif len(roots) < 2:
  17. raise NoPeaksFound("No proper peaks were found in the data set; likely "
  18. "the dataset is flat (e.g. all zeros).")
  19. else:
  20. return abs(roots[1] - roots[0])

字符串

展开查看全部
kxxlusnw

kxxlusnw5#

你应该使用scipy来解决它:先find_peaks,然后peak_widths。默认值为 rel_height(0.5),您测量的是峰的半峰宽。

ryhaxcpt

ryhaxcpt6#

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

  1. import numpy as np
  2. def get_full_width(x: np.ndarray, y: np.ndarray, height: float = 0.5) -> float:
  3. height_half_max = np.max(y) * height
  4. index_max = np.argmax(y)
  5. x_low = np.interp(height_half_max, y[:index_max+1], x[:index_max+1])
  6. x_high = np.interp(height_half_max, np.flip(y[index_max:]), np.flip(x[index_max:]))
  7. return x_high - x_low

字符串

nhjlsmyf

nhjlsmyf7#

对于具有许多数据点的单调函数,如果不需要完美的精度,我会用途:

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

字符串

6ovsh4lw

6ovsh4lw8#

我在haggis.math.full_width_half_max中实现了一个经验解决方案,它可以很好地处理噪声和非高斯数据。用法非常简单:

  1. fwhm = full_width_half_max(x, y)

字符串
该函数具有鲁棒性:它使用所请求的内插方案简单地找到数据的最大值和与“向下一半”阈值交叉的最近点。
这里有几个例子使用其他答案的数据。

@HYRY的平滑数据

  1. def make_norm_dist(x, mean, sd):
  2. return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2))
  3. x = np.linspace(10, 110, 1000)
  4. green = make_norm_dist(x, 50, 10)
  5. pink = make_norm_dist(x, 60, 10)
  6. blue = green + pink
  7. # create a spline of x and blue-np.max(blue)/2
  8. spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0)
  9. r1, r2 = spline.roots() # find the roots
  10. # Compute using my function
  11. fwhm, (x1, y1), (x2, y2) = full_width_half_max(x, blue, return_points=True)
  12. # Print comparison
  13. print('HYRY:', r2 - r1, 'MP:', fwhm)
  14. plt.plot(x, blue)
  15. plt.axvspan(r1, r2, facecolor='g', alpha=0.5)
  16. plt.plot(x1, y1, 'r.')
  17. plt.plot(x2, y2, 'r.')


对于平滑数据,结果非常精确:

  1. HYRY: 26.891157007233254 MP: 26.891193606203814


x1c 0d1x的数据

@Hooked's Noisy Data

  1. def gauss(x, p): # p[0]==mean, p[1]==stdev
  2. return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))
  3. # Create some sample data
  4. known_param = np.array([2.0, .7])
  5. xmin,xmax = -1.0, 5.0
  6. N = 1000
  7. X = np.linspace(xmin,xmax,N)
  8. Y = gauss(X, known_param)
  9. # Add some noise
  10. Y += .10*np.random.random(N)
  11. # Renormalize to a proper PDF
  12. Y /= ((xmax-xmin)/N)*Y.sum()
  13. # Fit a guassian
  14. p0 = [0,1] # Inital guess is a normal distribution
  15. errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
  16. p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y))
  17. fit_mu, fit_stdev = p1
  18. FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
  19. # Compute using my function
  20. fwhm, (x1, y1), (x2, y2) = full_width_half_max(X, Y, return_points=True)
  21. # Print comparison
  22. print('Hooked:', FWHM, 'MP:', fwhm)
  23. plt.plot(X, Y)
  24. plt.plot(X, gauss(X, p1), lw=3, alpha=.5, color='r')
  25. plt.axvspan(fit_mu - FWHM / 2, fit_mu + FWHM / 2, facecolor='g', alpha=0.5)
  26. plt.plot(x1, y1, 'r.')
  27. plt.plot(x2, y2, 'r.')


对于噪声数据(具有有偏基线),结果并不一致。

  1. Hooked: 1.9903193212254346 MP: 1.5039676990530118


一方面,高斯拟合对于数据不是非常优的,但另一方面,拾取与半最大阈值相交的最近点的策略可能也不是最优的。


展开查看全部
pbossiut

pbossiut9#

一个有效的python实现是values是一个列表:

  1. def calculate_FWHM(values):
  2. # Find the maximum value and its index
  3. max_value = max(values)
  4. max_index = values.index(max_value)
  5. # Find the value that is half the maximum
  6. half_max = max_value / 2
  7. # Find the indices where the values are closest to half the maximum on both sides of the peak
  8. left_index = next((i for i in range(max_index, -1, -1) if values[i] <= half_max), 0)
  9. right_index = next((i for i in range(max_index, len(values)) if values[i] <= half_max), len(values) - 1)
  10. # return the FWHM
  11. return right_index - left_index

字符串

展开查看全部

相关问题