在一个频谱中找到峰值的位置，使用numpy。

19 浏览
0 Comments

在一个频谱中找到峰值的位置，使用numpy。

我有一个TOF谱图，我想使用Python（numpy）实现一个算法，找出谱图中的所有极大值，并返回对应的x值。
\n我在网上查找了下面报告的算法。
\n这里的假设是在最大值附近，前一个值和最大值之间的差异大于一个数DELTA。问题是，我的谱图由等间距的点组成，即使在最大值附近也是如此，所以DELTA从来没有超过过，函数peakdet返回一个空数组。\n您有任何想法如何解决这个问题吗？我真的很感谢您对代码的评论，因为我在Python方面还是新手。\n谢谢！ \n

import sys
from numpy import NaN, Inf, arange, isscalar, asarray, array
def peakdet(v, delta, x = None): 
   maxtab = []
   mintab = []
   if x is None:
      x = arange(len(v))
      v = asarray(v)
   if len(v) != len(x):
      sys.exit('输入向量v和x的长度必须相等')
   if not isscalar(delta):
      sys.exit('输入参数delta必须是标量')
   if delta <= 0:
      sys.exit('输入参数delta必须是正数')
   mn, mx = Inf, -Inf
   mnpos, mxpos = NaN, NaN
   lookformax = True
   for i in arange(len(v)):
      this = v[i]
      if this > mx:
         mx = this
         mxpos = x[i]
      if this < mn:
         mn = this
         mnpos = x[i]
      if lookformax:
         if this < mx-delta:
            maxtab.append((mxpos, mx))
            mn = this
            mnpos = x[i]
            lookformax = False
      else:
         if this > mn+delta:
            mintab.append((mnpos, mn))
            mx = this
            mxpos = x[i]
            lookformax = True
return array(maxtab), array(mintab)

\n下面显示了部分谱图。实际上，我有比这里显示的更多的峰值。\n\"enter

0

3 答案

0 Comments

在查看了答案和建议之后，我决定提供一个我经常使用的解决方案，因为它简单直观且更易于调整。它使用滑动窗口，并计算当窗口沿x轴移动时，局部峰值作为最大值出现的次数。正如建议的那样，“局部最大值”没有普遍的定义，这意味着某些调整参数是不可避免的。该函数使用“窗口大小”和“频率”来微调结果。窗口大小以自变量（x）的数据点数来衡量，频率表示峰值检测的灵敏度（也表示为数据点数；频率值越低，产生的峰值越多，反之亦然）。主要函数如下所示：

def peak_finder(x0, y0, window_size, peak_threshold):
    # extend x, y using window size
    y = numpy.concatenate([y0, numpy.repeat(y0[-1], window_size)])
    x = numpy.concatenate([x0, numpy.arange(x0[-1], x0[-1]+window_size)])
    local_max = numpy.zeros(len(x0))
    for ii in range(len(x0)):
        local_max[ii] = x[y[ii:(ii + window_size)].argmax() + ii]
    u, c = numpy.unique(local_max, return_counts=True)
    i_return = numpy.where(c>=peak_threshold)[0]
    return(list(zip(u[i_return], c[i_return])))

同时，以下代码片段用于生成下面显示的图像：

import numpy
from matplotlib import pyplot
def plot_case(axx, w_f):
    p = peak_finder(numpy.arange(0, len(Y)), -Y, w_f[0], w_f[1])
    r = .9*min(Y)/10
    axx.plot(Y)
    for ip in p:
        axx.text(ip[0], r + Y[int(ip[0])], int(ip[0]),
                 rotation=90, horizontalalignment='center')
    yL = pyplot.gca().get_ylim()
    axx.set_ylim([1.15*min(Y), yL[1]])
    axx.set_xlim([-50, 1100])
    axx.set_title(f'window: {w_f[0]}, count: {w_f[1]}', loc='left', fontsize=10)
    return(None)
window_frequency = {1:(15, 15), 2:(100, 100), 3:(100, 5)}
f, ax = pyplot.subplots(1, 3, sharey='row', figsize=(9, 4),
                        gridspec_kw = {'hspace':0, 'wspace':0, 'left':.08,
                                       'right':.99, 'top':.93, 'bottom':.06})
for k, v in window_frequency.items():
    plot_case(ax[k-1], v)
pyplot.show()

三个案例展示了不同的参数值，分别为（从左到右）：(1)峰值过多，(2)峰值过少，(3)峰值数量适中。为了生成Y数据，我使用了上面给出的函数，并从答案中添加了一些噪声。希望对一些人有用，但它的效率主要取决于实际的峰值形状和距离。

0
0 Comments

问题的出现原因是：需要在一个频谱中找到峰值的位置。

解决方法是使用SciPy库中的find_peaks函数。首先，导入所需的库和模块。然后，生成一个长度为1000的零数组Y。通过插入峰值函数到Y中，使其具有峰值。接下来，给Y添加噪声，使其变得有噪声。为了找到峰值，我们需要将Y乘以-1，并使用find_peaks函数来获取实际的峰值位置。最后，为了绘制图形，我们需要将Y乘以-1。通过调用plot函数，我们可以绘制Y和峰值的位置。

在find_peaks函数中，可以设置height参数来限制只找到高于某个阈值的峰值。通过设置distance参数，可以限制峰值之间的最小距离。

因此，我们可以使用以下代码来找到高于0.002阈值且峰值之间距离大于100的峰值位置：

peaks, _ = find_peaks(Y, height=0.002, distance=100)

0
0 Comments

原因：这篇文章主要是为了解决在给定的频谱中找到峰值位置的问题。作者使用了一种基于信号处理的方法来解决这个问题。

解决方法：作者首先使用了scipy库中的convolve函数来计算频谱的导数。然后，通过检查导数的正负号来确定信号的斜率变化情况。接下来，作者找到了所有负斜率的位置，并筛选出最终的峰值位置。最后，作者使用matplotlib库将找到的峰值位置在频谱上进行了可视化。

文章中还提到了一些优化的方法，如根据峰值的大小设置一个阈值来筛选峰值，以及根据峰值的分布情况动态调整阈值。作者还介绍了如何处理基线变化的情况，以及如何简化代码以提高效率。

通过使用信号处理的方法，作者成功地解决了在给定频谱中找到峰值位置的问题。他提供了一种简单而有效的方法，并给出了一些优化的建议。这个方法可以在实际应用中用来分析具有峰值的频谱数据。

0