如果我有一个生成如下图的数据集,我将如何通过算法确定所示峰的 x 值(在本例中为三个):
如何在数据集中找到峰值?
一种通用的方法是平滑数据,然后通过将局部最大值滤波器与平滑滤波器进行比较来找到峰值。在R
:
argmax <- function(x, y, w=1, ...) {
require(zoo)
n <- length(y)
y.smooth <- loess(y ~ x, ...)$fitted
y.max <- rollapply(zoo(y.smooth), 2*w+1, max,
align="center")
delta <- y.max - y.smooth[-c(1:w, n+1-1:w)]
i.max <- which(delta <= 0) + w
list(x=x[i.max], i=i.max, y.hat=y.smooth)
}
它的返回值包括局部最大值 ( x
) 的参数——它回答了这个问题——以及这些局部最大值出现的 x 和 y 数组的索引 ( i
)。
有两个参数需要根据情况进行调整: w
是用于计算局部最大值的窗口的半宽。(它的值应该大大小于数据数组长度的一半。)小值会拾取微小的局部凸起,而较大的值会直接通过这些凸起。另一个——在这段代码中没有明确——是平滑器的span
参数。loess
(它通常在 0 和 1 之间;它反映了作为 x 值范围的比例的窗口宽度。)较大的值将更积极地平滑数据,使局部凹凸完全消失。
要查看此调整的效果,让我们创建一个小测试函数来绘制结果:
test <- function(w, span) {
peaks <- argmax(x, y, w=w, span=span)
plot(x, y, cex=0.75, col="Gray", main=paste("w = ", w, ",
span = ", span, sep=""))
lines(x, peaks$y.hat, lwd=2) #$
y.min <- min(y)
sapply(peaks$i, function(i) lines(c(x[i],x[i]), c(y.min,
peaks$y.hat[i]),
col="Red", lty=2))
points(x[peaks$i], peaks$y.hat[peaks$i], col="Red", pch=19,
cex=1.25)
}
以下是一些应用于一些合成的、有轻微噪声的数据的实验。
x <- 1:1000 / 100 - 5
y <- exp(abs(x)/20) * sin(2 * x + (x/5)^2) + cos(10*x) / 5 +
rnorm(length(x), sd=0.05)
par(mfrow=c(3,1))
test(2, 0.05)
test(30, 0.05)
test(2, 0.2)
宽窗口(中间图)或更激进的平滑(底部图)消除了在顶部图中检测到的局部最大值。这里的最佳组合可能是一个宽窗口并且只有温和的平滑,因为激进的平滑似乎会移动这些峰值(参见底部图中的中间点和右侧点,并将它们的位置与原始数据的明显峰值进行比较)。在这个例子中,w=50
并且span=0.05
做得很好(未显示)。
请注意,未检测到端点处的局部最大值。这些可以单独检查。(为了支持这一点,argmax
返回平滑的 y 值。)
与用于通用工作的更正式的建模相比,这种方法有几个优点:
它不采用任何先入为主的数据模型。
它可以适应数据特征。
它可以适应检测人们感兴趣的峰的种类。
正如我在评论中提到的,如果时间序列似乎是周期性拟合的,谐波回归模型提供了一种通过应用一阶和二阶导数测试来平滑函数并识别峰值的方法。Huber 指出了一种非参数检验,该检验在有多个峰值且函数不一定是周期性的情况下具有优势。但是没有免费的午餐。虽然他提到他的方法有优点,但如果参数模型合适,则可能存在缺点。这始终是使用非参数技术的另一面。尽管它避免了参数化假设,但当参数化假设合适时,参数化方法会更好。他的程序也没有充分利用数据中的时间序列结构。
我认为,虽然指出建议程序的优点是适当的,但指出潜在的缺点也很重要。我的方法和 Huber 的方法都以有效的方式找到了峰值。但是,我认为当局部最大值低于先前确定的最高峰时,他的程序需要做更多的工作。
信号处理中的经典峰值检测方法如下:
- 根据采样率和信号属性将信号过滤到某个合理的合理范围,例如对于 ECG,IIR 带通滤波器@0.5-20Hz,零相位滤波器将确保不会引入相移(和相关的时间延迟)
- 然后可以使用希尔伯特变换或小波方法来强调峰值
- 然后可以应用静态或动态阈值,其中所有高于阈值的样本都被视为峰值。在动态阈值的情况下,它通常定义为高于或低于平均值的移动平均估计值的阈值N个标准差。
另一种可行的方法是将急剧高通滤波的信号与高度平滑(低通或中值滤波)的信号进行比较,然后应用第 3 步。
希望这可以帮助。