测量信号的峰值检测

StackOverflow https://stackoverflow.com/questions/3260

  •  08-06-2019
  •  | 
  •  

题

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

如果最高值是我们正在寻找的峰值，那么这种方法效果很好，但如果设备工作不正常，我们可以看到第二个峰值，该峰值可能高于初始峰值。我们在 90 秒的时间内每秒从 16 个设备获取 10 个读数。

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

我们确实使用 LabVIEW，我已经检查过 熔岩论坛 还有很多有趣的例子。这是我们测试软件的一部分，我们试图避免使用太多非标准 VI 库，因此我希望获得有关所涉及的流程/算法的反馈，而不是具体的代码。

有帮助吗？

解决方案

您可以尝试信号平均，即对于每个点，将其值与周围 3 个或更多点进行平均。如果噪声信号很大，那么即使这样也可能无济于事。

我意识到这与语言无关，但猜测您正在使用 LabView，LabView 附带有许多预打包的信号处理 VI，您可以使用它们来进行平滑和降噪。这 NI 论坛 是在此类事情上获得更专业帮助的好地方。

其他提示

有很多很多经典的峰值检测方法，其中任何一种都可能有效。您必须特别了解是什么限制了数据的质量。以下是基本描述：

  1. 在数据中的任意两点之间， (x(0), y(0)) 和 (x(n), y(n)), ， 加起来 y(i + 1) - y(i) 为了 0 <= i < n 并称之为 T （“旅行”）并设置 R （“升起到 y(n) - y(0) + k 对于适当小的 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. 将一个非常简单的 2 高斯混合模型拟合到您的数据（例如，Numerical Recipes 有一个很好的现成代码块）。以较早的高峰为例。这将正确处理重叠峰。

  6. 在数据中找到与简单高斯曲线、柯西曲线、泊松曲线或什么曲线的最佳匹配。在广泛的范围内评估该曲线，并在注意到其峰值位置后从数据副本中减去它。重复。取模型参数（可能是标准差，但某些应用程序可能关心峰度或其他特征）满足某些标准的最早峰值。注意从数据中减去峰值时留下的伪影。最佳匹配可能由上面 #2 中建议的匹配评分类型来确定。

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

我鼓励您仔细关注正确的基线。维纳滤波或其他滤波或简单的直方图分析通常是在存在噪声的情况下确定基线的简单方法。

最后，如果您的数据通常有噪声，并且您从卡上获取数据作为未引用的单端输出（甚至引用，只是不是差分），并且如果您将大量观察结果平均到每个数据点中，请尝试对这些数据进行排序观察并丢弃第一个和最后一个四分位数并对剩下的进行平均。有许多此类非常有用的异常值消除策略。

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

有一组非常最新的实现 TSpectrum* 类的 根 （核/粒子物理分析工具）。该代码适用于一维到三维数据。

ROOT 源代码是可用的，因此您可以根据需要获取此实现。

来自 频谱 类文档：

此类中使用的算法已发布在以下参考文献中：

[1] M.Morhac 等人：多维重合伽马射线光谱的背景消除方法。物理研究中的核仪器和方法A 401（1997）113-132。

[2] M.Morhac 等人：有效的一维黄金反卷积及其在伽玛射线光谱分解中的应用。物理研究中的核仪器和方法A 401（1997）385-408。

[3] M.Morhac 等人：多维重合伽马射线光谱中峰的鉴定。研究物理学中的核仪器和方法A 443（2000），108-125。

对于那些没有 NIM 在线订阅的人来说，这些论文来自课程文档的链接。

所做工作的简短版本是，将直方图展平以消除噪声，然后在展平的直方图中通过强力检测局部最大值。

我想为这个线程贡献一个算法 我已经发展了自己:

它基于以下原理 分散: ：如果新数据点与某个移动平均值相差给定的 x 个标准差，则算法会发出信号（也称为 z 分数）。该算法非常稳健，因为它构造了一个 分离 移动均值和偏差，使得信号不会破坏阈值。因此，无论先前信号的数量如何，未来信号的识别精度都大致相同。该算法需要 3 个输入： lag = the lag of the moving window, threshold = the z-score at which the algorithm signals 和 influence = the influence (between 0 and 1) of new signals on the mean and standard deviation. 。例如，一个 lag of 5 将使用最后 5 个观测值来平滑数据。A threshold 如果数据点与移动平均值相差 3.5 个标准差，则 3.5 将发出信号。和 influence 0.5 给出信号 一半 正常数据点的影响。同样，一个 influence 0 完全忽略信号以重新计算新阈值：因此，影响力为 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

演示

Demonstration of robust thresholding algorithm

> 原答案

这个方法基本上来自David Marr的书《Vision》

高斯模糊您的信号与预期的峰值宽度。这消除了噪声尖峰，并且您的相位数据不会受到损坏。

然后边缘检测（LOG即可）

那么你的边缘就是特征的边缘（比如峰）。在边缘之间查找峰值，按大小对峰值进行排序，然后就完成了。

我已经使用了它的变体并且效果非常好。

我想你想要 互相关 您的信号与预期的示例信号。但是，自从我学习信号处理以来已经很久了，即使那时我也没有太在意。

我对仪器不太了解，所以这可能完全不切实际，但话又说回来，这可能是一个有用的不同方向。如果您知道读数如何会失败，并且在出现此类失败的情况下，峰值之间有一定的间隔，为什么不在每个间隔进行梯度下降。如果下降带你回到之前搜索过的区域，你可以放弃它。根据采样表面的形状，这也可能帮助您比搜索更快地找到峰值。

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

你可以应用一些 标准差 按照您的逻辑并注意超过 x% 的峰值。

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top