测量信号的峰值检测

Question

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

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

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

---

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

Answer 1

有许多经典的峰值检测方法,其中任何一种都可行.您必须看到特别限制数据质量的内容.以下是基本描述:

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序列数据中找到峰值,在测量曲线中估算衍生物的峰值,并在直方图中找到峰值.

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

最后,如果您的数据通常是嘈杂的,并且您从卡中获取数据作为未引用的单端输出(或甚至引用,只是非差分),并且如果您将大量观察结果平均到每个数据点,请尝试对这些数据进行排序观察并扔掉第一个和最后一个四分位数并平均剩下的四分位数.有许多这样的异常消除策略可以真正有用.

Answer 2

您可以尝试信号平均,即对于每个点,使用周围的3个或更多点对该值求平均值.如果噪音很大,那么即使这样也许没有用.

我意识到这是语言不可知的,但猜测你正在使用LabView,LabView附带了许多预先打包的信号处理VI,你可以使用它来进行平滑和降噪.在NI论坛是一个伟大的地方得到这样的事情更专业的帮助.

Answer 3

已经对该问题进行了一些详细研究.

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

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

从TSpectrum类文档:

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

[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在线订阅的人,这些文件与课程文档相关联.

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

Answer 4

我想为这个线程贡献一个我自己开发的算法:

它基于分散原理:如果新的数据点是距离某些移动平均值的给定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.例如,a lag的5将使用最后5个观察值来平滑数据.threshold如果数据点与移动平均值相差3.5个标准差,则3.5的A 将发出信号.influence0.5的一个信号给出正常数据点具有一半影响的信号.同样,influence0为完全忽略信号以重新计算新阈值: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

演示

>原始答案