algorithm -测量信号峰值检测

Translate

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

如果最大值是我们正在寻找的峰值,则此方法效果很好,但是如果设备无法正常工作,我们可以看到第二个峰值,该峰值可能高于初始峰值。在90秒钟的时间内,我们从16台设备中每秒读取10个读数。

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


我们确实使用LabVIEW,并且已经检查了LAVA论坛还有很多有趣的例子。这是我们测试软件的一部分,我们试图避免使用过多的非标准VI库,因此我希望能就所涉及的过程/算法而不是特定的代码获得反馈。

This question and all comments follow the "Attribution Required."

所有的回答

Translate

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

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

来源
Translate

有很多经典的峰值检测方法,其中任何一种都可以使用。您必须查看特别是什么会限制数据质量。以下是基本说明:

  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高斯混合模型拟合到您的数据(例如,“数值配方”有一个很好的现成代码块)。采取较早的高峰。这将正确处理重叠的峰。

  6. 在数据中找到最匹配的简单高斯曲线,柯西曲线,泊松曲线或您拥有的曲线。在宽范围内评估该曲线,并在注意到其峰值位置后从数据副本中减去它。重复。取其模型参数(可能为标准偏差,但某些应用程序可能关心峰度或其他特征)满足某些标准的最早峰。当从数据中减去峰时,请注意是否留有伪影。最佳匹配可能取决于上面#2中建议的匹配得分类型。

我已经完成了您之前的工作:在DNA序列数据中查找峰,在根据测量曲线估算的导数中查找峰,并在直方图中查找峰。

我鼓励您仔细参加适当的基准调查。在存在噪声的情况下,维纳滤波或其他滤波或简单的直方图分析通常是基线的简单方法。

最后,如果您的数据通常比较嘈杂,并且您从卡中获取的数据是未引用的单端输出(或者甚至是被引用的,只是没有差分),并且如果您对每个数据点平均了很多观察结果,请尝试对这些数据进行排序观察并丢弃第一个和最后一个四分位数,然后对剩余的值求平均。有许多这样的异常消除策略可能非常有用。

来源
Translate

已经对该问题进行了详细的研究。

在其中有一组最新的实现TSpectrum*(核/粒子物理分析工具)。该代码适用于一维到三维数据。

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

来自频谱类文档:

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

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

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

[3] M.Morhac等人:多维重合伽马射线光谱中峰的识别。 《研究物理中的核仪器与方法》,A 443(2000),108-125。

这些文章是从班级文档中链接的,这些课程没有您的NIM在线订阅。


这样做的简短形式是直方图被展平以消除噪声,然后在展平的直方图中用蛮力检测局部最大值。

来源
Translate

我想为此算法贡献一个我发展了自己:

它基于以下原则分散:如果新数据点是相对于某个移动平均值的给定x倍标准偏差,则算法会发出信号(也称为分数)。该算法非常健壮,因为它可以构造一个分离移动均值和偏差,以使信号不会破坏阈值。因此,无论先前信号的数量如何,都以大致相同的精度识别未来信号。该算法有3个输入:lag = the lag of the moving window, threshold = the z-score at which the algorithm signalsinfluence = the influence (between 0 and 1) of new signals on the mean and standard deviation。例如,一个lag中的5个将使用最后5个观测值来平滑数据。一种threshold如果数据点与移动平均值相差3.5个标准偏差,则信号3.5将发出信号。和influence0.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

>原始答案

来源
Translate

此方法基本上来自David Marr的著作“ Vision”

高斯模糊会以预期的峰宽模糊信号。这消除了噪声尖峰,并且相位数据没有损坏。

然后边缘检测(LOG将执行)

然后,您的边缘就是特征(如峰)的边缘。在峰的边缘之间查找峰,按大小对峰进行排序,就完成了。

我对此使用了各种变体,它们工作得很好。

来源
Translate

我想你想互相关您的信号具有预期的示例性信号。但是,自从我学习信号处理以来已经有很长的时间了,即使那时我也没有引起太大的注意。

来源
Translate

我对仪器不是很了解,因此这可能是完全不切实际的,但同样,它可能是一个有用的不同方向。如果您知道读数可能失败的原因,并且给出此类失败的峰之间有一定的间隔,那么为什么不每个间隔都进行梯度下降呢?如果下降使您回到之前搜索过的区域,则可以放弃。根据采样表面的形状,这也可能会帮助您比搜索更快地找到峰。

来源
Translate

所需峰与不需要的第二峰之间是否存在质的差异?如果两个峰值都是“尖锐的”(即持续时间较短),则在频域中查看信号(通过执行FFT)时,您将在大多数频带上获得能量。但是,如果“好”峰确实有能量存在于“坏”峰中不存在的频率上,反之亦然,则可以通过这种方式自动区分它们。

来源
Translate

你可以申请一些标准偏差根据您的逻辑,并注意超过x%的峰值。

来源