找到光谱中的峰值位置numpy

Question

我有一个TOF频谱,我想用python(numpy)实现一个算法,它找到频谱的所有最大值并返回相应的x值.
我在网上查了一下,发现下面报告的算法.

这里的假设是,在最大值附近,之前的值与最大值之间的差值大于数字DELTA.问题是我的光谱由均匀分布的点组成,即使在最大值附近,也不会超过DELTA,函数peakdet返回一个空数组.

你知道如何克服这个问题吗？我非常感谢评论,以便更好地理解代码,因为我是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:简化和评论

我简化了代码,使用一个内核convolve作为@skymandr评论的s.这也消除了调整收缩率的神奇数字,以便任何大小的内核都应该这样做.

选择"valid"作为选项convolve.它可能也有效"same",但我选择"valid"这样我就不必考虑边缘条件,如果算法可以检测出那里的spurios峰值.

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)

在上面的代码中。