在频谱中找到峰值位置 numpy

Question

我有一个 TOF 谱，我想使用 python (numpy) 实现一个算法，它可以找到谱的所有最大值并返回相应的 x 值。
我在网上查了一下，发现下面报道的算法。

这里的假设是，在最大值附近，之前的值与最大值之间的差值大于数字 DELTA。问题是我的频谱是由均匀分布的点组成的，甚至接近最大值，因此 DELTA 永远不会超过并且函数 peakdet 返回一个空数组。

您知道如何解决这个问题吗？我非常感谢 cmets 能够更好地理解代码，因为我是 python 的新手。

谢谢！

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('Input vectors v and x must have same length')
   if not isscalar(delta):
      sys.exit('Input argument delta must be a scalar')
   if delta <= 0:
      sys.exit('Input argument delta must be positive')

   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)

下面显示了频谱的一部分。实际上，我的峰比这里显示的要多。

Answer 1

我认为这可以作为一个起点。我不是信号处理专家，但我在生成的信号Y 上尝试了这个，它看起来很像你的信号，但噪音更大：

from scipy.signal import convolve
import numpy as np
from matplotlib import pyplot as plt
#Obtaining derivative
kernel = [1, 0, -1]
dY = convolve(Y, kernel, 'valid') 

#Checking for sign-flipping
S = np.sign(dY)
ddS = convolve(S, kernel, 'valid')

#These candidates are basically all negative slope positions
#Add one since using 'valid' shrinks the arrays
candidates = np.where(dY < 0)[0] + (len(kernel) - 1)

#Here they are filtered on actually being the final such position in a run of
#negative slopes
peaks = sorted(set(candidates).intersection(np.where(ddS == 2)[0] + 1))

plt.plot(Y)

#If you need a simple filter on peak size you could use:
alpha = -0.0025
peaks = np.array(peaks)[Y[peaks] < alpha]

plt.scatter(peaks, Y[peaks], marker='x', color='g', s=40)

样本结果： 对于嘈杂的，我用alpha过滤了峰值：

如果alpha 需要更多复杂性，您可以尝试从使用例如发现的峰值动态设置 alpha。假设它们是混合高斯（我最喜欢的是 Otsu 阈值，存在于 cv 和 skimage）或某种聚类（k-means 可以工作）。

作为参考，这是我用来生成信号的：

Y = np.zeros(1000)

def peaker(Y, alpha=0.01, df=2, loc=-0.005, size=-.0015, threshold=0.001, decay=0.5):  
    peaking = False
    for i, v in enumerate(Y):
        if not peaking:
            peaking = np.random.random() < alpha
            if peaking:
                Y[i] = loc + size * np.random.chisquare(df=2)
                continue
        elif Y[i - 1] < threshold:
            peaking = False

        if i > 0:
            Y[i] = Y[i - 1] * decay

peaker(Y)

编辑：支持降级基线

我通过这样做模拟了一条倾斜的基线：

Z = np.log2(np.arange(Y.size) + 100) * 0.001
Y = Y + Z[::-1] - Z[-1]

然后用固定的 alpha 进行检测（注意我在 alpha 上更改了符号）：

from scipy.signal import medfilt

alpha = 0.0025
Ybase = medfilt(Y, 51) # 51 should be large in comparison to your peak X-axis lengths and an odd number.
peaks = np.array(peaks)[Ybase[peaks] - Y[peaks] > alpha] 

导致以下结果（基线绘制为黑色虚线）：

编辑 2：简化和注释

正如@skymandr 评论的那样，我简化了代码，为convolves 使用一个内核。这也消除了调整收缩的幻数，因此任何大小的内核都应该这样做。

选择"valid" 作为convolve 的选项。它可能与"same" 一样有效，但我选择了"valid"，因此我不必考虑边缘条件以及算法是否可以检测到那里的伪峰。

Answer 2

从 SciPy 1.1 版开始，您还可以使用find_peaks：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

np.random.seed(0)

Y = np.zeros(1000)

# insert @deinonychusaur's peaker function here

peaker(Y)

# make data noisy
Y = Y + 10e-4 * np.random.randn(len(Y))
# find_peaks gets the maxima, so we multiply our signal by -1
Y *= -1 
# get the actual peaks
peaks, _ = find_peaks(Y, height=0.002)
# multiply back for plotting purposes
Y *= -1
plt.plot(Y)
plt.plot(peaks, Y[peaks], "x")
plt.show()

这将绘制（注意我们使用height=0.002，它只会找到高于 0.002 的峰值）：

除了height，我们还可以设置两个峰之间的最小距离。如果您使用distance=100，则绘图如下所示：

你可以使用

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

在上面的代码中。

Answer 3

在查看了答案和建议后，我决定提供一个我经常使用的解决方案，因为它简单且易于调整。 它使用滑动窗口并计算局部峰值出现的次数作为最大值，因为窗口沿 x 轴移动。正如@DrV 建议的那样，不存在“局部最大值”的通用定义，这意味着一些调整参数是不可避免的。此功能使用“窗口大小”和“频率”来微调结果。窗口大小以自变量 (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])))

加上一个sn-p用来产生下图：

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 数据，我使用了上面提供的 @deinonychusaur 函数，并从 @Cleb 的回答中添加了一些噪音。

我希望有些人会觉得这很有用，但它的效率主要取决于实际的峰形和距离。

Answer 4

找到最小值或最大值并不是那么简单，因为“局部最大值”没有通用的定义。

您的代码似乎在寻找一个混合最小值，然后如果信号低于最大值减去某个增量值之后的最大值下降，则将其作为最大值接受。之后，它开始寻找具有类似标准的最小值。您的数据是缓慢下降还是缓慢上升并不重要，因为当达到最大值时会记录最大值，并在水平低于滞后阈值时附加到最大值列表中。

这是一种找到局部最小值和最大值的可能方法，但它有几个缺点。其中之一是该方法不是对称的，即如果相同的数据反向运行，结果不一定相同。

很遗憾，我无法提供更多帮助，因为正确的方法实际上取决于您正在查看的数据、其形状和噪音。如果您有一些样品，那么我们也许可以提出一些建议。