我正在分析一些脑电图波形并基于它们开发分类器。在我拥有的数据中，我通常会在非常低的频率处找到一个强峰值（我对此峰值不感兴趣），然后是更高频率的另一个峰值（称为f0），我有一种预感，在我想用我的分类器预测的两个类之间的f0 。一个示例功率谱（将原始谱分成 200 个区间后）如下所示：
谁能告诉一种可靠地检测和估计峰值出现在此类光谱中的频率的好技术（即，我想估计蓝线指示的频率）？我正在使用 R 进行数据分析，并且需要 R 中的解决方案。
在 R 中查找信号的功率谱中的峰值
机器算法验证
r
光谱分析
2022-03-21 17:15:30
1个回答
您可以将其视为分析定位以图形方式观察到的局部最大值的问题。下面的代码通过一维优化算法定位这些最大值，该算法应用于覆盖整个 x 轴的重叠间隔。在每个间隔，目标函数由样本周期图的样条插值定义。
我认为在包或帖子中的某个地方有一个更完整的实现这个想法，但我找不到它。
例如，我为“AirPassengers”数据的日志获取样本周期图：
sp <- spectrum(log(AirPassengers), span = c(3, 5))
接下来，我定义目标函数、检查是否存在局部最大值的区间以及梯度和 Hessian 的一些阈值，这些阈值将用于决定是否接受作为局部最大值的解。
# objective function, spline interpolation of the sample spectrum
f <- function(x, q, d) spline(q, d, xout = x)$y
x <- sp$freq
y <- log(sp$spec)
nb <- 10 # choose number of intervals
iv <- embed(seq(floor(min(x)), ceiling(max(x)), len = nb), 2)[,c(2,1)]
# make overlapping intervals to avoid problems if the peak is close to
# the ends of the intervals (two modes could be found in each interval)
iv[-1,1] <- iv[-nrow(iv),2] - 2
# The function "f" is maximized at each of these intervals
iv
# [,1] [,2]
# [1,] 0.0000000 0.6666667
# [2,] -1.3333333 1.3333333
# [3,] -0.6666667 2.0000000
# [4,] 0.0000000 2.6666667
# [5,] 0.6666667 3.3333333
# [6,] 1.3333333 4.0000000
# [7,] 2.0000000 4.6666667
# [8,] 2.6666667 5.3333333
# [9,] 3.3333333 6.0000000
# choose thresholds for the gradient and Hessian to accept
# the solution is a local maximum
gr.thr <- 0.001
hes.thr <- 0.03
现在为每个间隔运行优化算法。梯度和 Hessian 是通过函数grad和hessian来自 package的数值计算的numDeriv。
require("numDeriv")
vals <- matrix(nrow = nrow(iv), ncol = 3)
grd <- hes <- rep(NA, nrow(vals))
for (j in seq(1, nrow(iv)))
{
opt <- optimize(f = f, maximum = TRUE, interval = iv[j,], q = x, d = y)
vals[j,1] <- opt$max
vals[j,3] <- exp(opt$obj)
grd[j] <- grad(func = f, x = vals[j,1], q = x, d = y)
hes[j] <- hessian(func = f, x = vals[j,1], q = x, d = y)
if (abs(grd[j]) < gr.thr && abs(hes[j]) > hes.thr)
vals[j,2] <- 1
}
# it is convenient to round the located peaks in order to avoid
# several local maxima that essentially the same point
vals[,1] <- round(vals[,1], 2)
if (anyNA(vals[,2])) {
peaks <- unique(vals[-which(is.na(vals[,2])),1])
} else
peaks <- unique(vals[,1])
在这种情况下，定位的峰如下：
peaks
# [1] 0.06 1.00 2.01 2.99 4.00 4.98
从图形上我们可以看到，在这个例子中，它们与周期图中观察到的峰值很好地匹配：
plot(sp$freq, sp$spec, log = "y", type = "l")
abline(v = peaks, lty = 2)
其它你可能感兴趣的问题