测量信号的峰值检测

Question

我们使用数据采集卡从一个设备中获取读数，该设备会使其信号增加到峰值，然后回落到接近原始值。为了找到峰值，我们目前搜索数组以找到最高读数，并使用索引确定峰值的时间，该时间用于我们的计算。

如果最高值不是我们要寻找的峰值，而设备没有正常工作，则可能会看到第二个峰值，该峰值可能比初始峰值更高。我们在90秒的时间内从16个设备中每秒进行10次读取。

我的初步想法是循环读数，检查前一个和后一个点是否小于当前点，以找到峰值并构建峰值数组。也许我们应该查看当前位置两侧多个点的平均值，以允许系统中的噪声。这是继续进行的最佳方式，还是有更好的技术？

---

我们确实使用LabVIEW，我已经查看了LAVA论坛，发现有许多有趣的示例。这是我们测试软件的一部分，我们试图避免使用过多的非标准VI库，因此我希望得到有关流程/算法的反馈，而不是具体的代码。

- Swinders

Answer 1

有许多经典的峰值检测方法，其中任何一种都可能有效。您需要查看特定数据质量的限制。以下是基本描述：

1. 在数据的任意两个点之间（x(0), y(0)）和（x(n), y(n)），对于0 <= i < n, 将y(i+1)-y(i)相加，并将其称为T（“travel”），并适当设置k，将R（“rise”）设置为y(n)-y(0)+k。 T / R> 1表示峰值。如果噪声产生的大幅度波动不太可能或噪音在基础曲线形状周围对称分布，则此方法可行。对于您的应用程序，请接受早期得分高于给定阈值的峰值，或者分析每个升值值的旅行曲线以获取更有趣的属性。

2. 使用匹配滤波器来评估与标准峰形状的相似性得分（基本上，使用某些形状的归一化点积来获得余弦指标的相似性）。

3. 对标准峰值形状进行反卷积，并检查高值（尽管我通常发现简单仪器输出的噪音对2不敏感）。

4. 平滑数据并检查三个等间距点的三元组，如果x0 < x1 < x2 且 y1 > 0.5 *（y0 + y2），或者检查像这样的欧几里得距离：D（（x0，y0），（x1，y1））+ D（（x1，y1），（x2，y2））> D（（x0，y0），（x2，y2）），它依赖于三角形不等式。再次使用简单比率将为您提供得分机制。

5. 对数据进行非常简单的两个高斯混合模型拟合（例如，Numerical Recipes有一段漂亮的现成代码）。取较早的峰值。这将正确处理重叠的峰值。

6. 在数据中找到与简单的高斯，柯西，泊松或其他任何曲线最佳匹配的位置。在广泛范围内评估该曲线并从数据的副本中减去它的峰值位置之后重复此操作。采用最早的峰值，其模型参数（可能是标准差，但某些应用程序可能关心峰度或其他特征）符合某些条件。注意，在峰值从数据中减去时留下的人工效应。

我以前也做过您正在做的事情：在DNA序列数据中查找峰值，在从测量曲线估计的导数中查找峰值，并在直方图中查找峰值。

我建议您仔细进行适当的基线处理。在存在噪声的情况下，Wiener过滤或其他过滤器或简单的直方图分析通常是一种简单的基线处理方法。

最后，如果您的数据通常存在噪声，并且您从卡片上获取的数据是未参考的单端输出（甚至是参考的，但不是差分的），并且如果您将大量观测值平均到每个数据点中，请尝试对这些观测值进行排序，并且删除前后四分之一，并对剩余部分进行平均。有许多这样的异常值消除策略可以真正有用。

- Thomas Kammeyer

---

关于“+ k适当小”，k是否为所有计算的常数值？你如何选择它？适当小是什么意思？ - Loren

这只是为了防止溢出：当你将T除以R时，如果曲线是平的或几乎平的，或者第一个和最后一个y值相等或几乎相等，那么R将接近于零。适当地小就足够避免在R不接近零时扰乱比率。这主要是一个实际的、经验性的问题。它只是为了防止零除，你也可以选择用其他方式来做。 - Thomas Kammeyer

Answer 2

您可以尝试信号平均化，即对于每个点，将其值与周围的3个或更多点平均。如果噪声干扰很大，则即使这样做也可能无法解决问题。

我意识到这种方法不依赖特定编程语言，但是猜测您正在使用LabView，那么可以使用LabView提供的预打包的信号处理VI进行平滑和降噪处理。前往NI论坛获取更专业的帮助。

- Brendan

Answer 3

我愿意为这个主题提供一种算法，这是我自己开发的算法：

它基于分散原理：如果新数据点与某个移动平均值相距x个标准差，则该算法会发出信号（也称为z-score）。该算法非常强大，因为它构建了一个独立的移动平均和偏差，使得信号不会破坏阈值。因此，未来的信号将以大致相同的精度标识，而不受先前信号量的影响。该算法需要3个输入：lag = 移动窗口的滞后期、threshold = 算法信号的z-score和influence = 新信号对均值和标准差的影响（0到1之间）。例如，lag为5将使用最后5次观测结果平滑数据。当数据点与移动平均值相差3.5个标准差时，threshold将会发出信号。当influence为0.5时，信号的影响力只有普通数据点的一半。同样地，influence为0将完全忽略信号以重新计算阈值：因此，influence为0是最强大的选择。

它的工作原理如下：

伪代码

# Let y be a vector of timeseries data of at least length lag+2
# Let mean() be a function that calculates the mean
# Let std() be a function that calculates the standard deviaton
# Let absolute() be the absolute value function

# Settings (the ones below are examples: choose what is best for your data)
set lag to 5;          # lag 5 for the smoothing functions
set threshold to 3.5;  # 3.5 standard deviations for signal
set influence to 0.5;  # between 0 and 1, where 1 is normal influence, 0.5 is half

# Initialise variables
set signals to vector 0,...,0 of length of y;   # Initialise signal results
set filteredY to y(1,...,lag)                   # Initialise filtered series
set avgFilter to null;                          # Initialise average filter
set stdFilter to null;                          # Initialise std. filter
set avgFilter(lag) to mean(y(1,...,lag));       # Initialise first value
set stdFilter(lag) to std(y(1,...,lag));        # Initialise first value

for i=lag+1,...,t do
  if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then
    if y(i) > avgFilter(i-1)
      set signals(i) to +1;                     # Positive signal
    else
      set signals(i) to -1;                     # Negative signal
    end
    # Adjust the filters
    set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  else
    set signals(i) to 0;                        # No signal
    # Adjust the filters
    set filteredY(i) to y(i);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  end
end

---

演示

> 欲了解更多信息，请参见原始答案

- Jean-Paul

---

如果我不仅想检测峰值，而且还想分离形成峰值的单个波，我该如何使用这个算法？我不仅想要峰值，还想知道检测到的尖峰的宽度/持续时间和上升时间/斜率。 - Brandon Brown

@BrandonBrown 这应该不难吧？只需修改原始算法，使每次出现信号时保存所需的详细信息。例如，实现一个变量来“标记”峰值开始的时间，跟踪您需要的值，并在峰值停止时保存它们。 - Jean-Paul

Answer 4

这个问题已经被详细研究过了。

在ROOT（一个核/粒子物理分析工具）的TSpectrum*类中有一系列最新的实现。该代码适用于一到三维数据。

ROOT源代码是可用的，如果您想要，可以获取此实现。

从TSpectrum类文档中可以看到：

该类中使用的算法已经在以下参考文献中发表：

[1] M.Morhac et al.: 多维符合伽马射线谱的背景消除方法。《物理研究A：核仪器和物理研究方法》401 (1997) 113-132。

[2] M.Morhac et al.: 高效一维和二维Gold反褶积及其在伽马射线谱分解中的应用。《物理研究A：核仪器和物理研究方法》401 (1997) 385-408。

[3] M.Morhac et al.: 多维符合伽马射线谱峰位的鉴定。《物理研究A：核仪器和物理研究方法》443(2000), 108-125。

对于没有NIM在线订阅的人，这些论文链接可以在课程文档中找到。

---

简单来说，所做的是将直方图压平以消除噪声，然后在压平的直方图中通过蛮力法检测局部极大值。

- dmckee

Answer 5

这种方法基本上来自David Marr的书“视觉”。

使用期望峰宽对信号进行高斯模糊处理。这样可以消除噪声尖峰，而您的相位数据不会受到损害。

然后进行边缘检测（LOG即可）。

然后您的边缘就是特征（如峰）的边缘。在边缘之间寻找峰值，按大小排序，完成。

我已经使用了这种方法的变体，效果非常好。

- Tim Williscroft

Answer 6

我认为你想用期望的模板信号对你的信号进行交叉相关。但是，自从我学习信号处理以来已经过了很长时间，那时我也没有特别关注。请参考交叉相关。

- Anthony Cramp

Answer 7

你可以在逻辑中应用一些标准差，并注意超过x%的峰值。

- GateKiller

Answer 8

我对仪器仪表了解不多，所以这可能完全不切实际，但另一方面它可能是一个有用的不同方向。如果您知道读数可能失败的方式，并且在这种故障情况下给定峰值之间存在一定的时间间隔，为什么不在每个间隔处进行梯度下降。如果下降将您带回到以前搜索过的区域，则可以放弃它。根据采样表面的形状，这也可能比搜索更快地找到峰值。

- John with waffle

Answer 9

所需峰值和不需要的第二个峰值之间是否存在定性差异？如果两个峰值都是“尖锐”的——即时间持续短暂——当在频域中查看信号（通过执行FFT）时，您将在大多数频段上获得能量。但是，如果“好”的峰值可靠地在不存在于“坏”峰值中的频率上存在能量，或者反之亦然，则可以通过这种方式自动区分它们。

- Adam Hollidge