周期图中众数的解释

机器算法验证 r 解释 光谱分析
2022-03-28 10:31:43

我有一个以 1000 Hz 到 3 分钟采样的数据集。所以有180000个数据点。我为这些数据绘制了周期图,并得到了一系列峰值。

  1. 强度似乎最高为 0。这是什么意思?
  2. 强度似乎相对较高,为 0.005。我该如何解释?这是否意味着0.005处的峰值表示每单位时间有0.005个周期,或者每1/0.005左右一个完整的周期,即200个单位时间?由于实际数据中200个单位/时间= 0.2s,是否意味着我的信号中最主要的频率是5 Hz?
  3. 有没有办法将周期​​图的 x 和 y 值提取为 R 中的数据框?我想提取 y 的相应值并使用ggplot.

周期图

1个回答

提前道歉,我写了一点小说。

tl;博士

1) 数据中可能存在均值或趋势。这通常以零频率表示。

2) 是的,您可以说此时有一个周期分量,周期​​为 200 个时间单位(秒)。

3)是的,你可以这样做,这取决于方法。如果您使用了 spec.pgram() 那么:

    pgram <- spec.pgram(yourDataVector)
    new.df <- data.frame(x = pgram$freq, y = pgram$spec)

但是,出于对无偏估计者的喜爱,请不要使用周期图(有很多统计上更好的方法,Rayleigh 评论了 1905-ish 中周期图的不良特性(旁注:spec.pgram 实际上并没有计算周期图,但使用 10% 余弦锥度为您提供频谱的“直接估计” - 还有更好的(就偏差和方差而言)方法))。

完整版本:

1) f = 0 或第一个频率区间中的频率表明您的数据可能具有尚未删除的趋势或均值。此外,f = 0 附近的极低频信号实际上看起来并不像在 f = 0 处。在我看来,有一个远离零的大峰值可能很有趣。仅在 f = 0 到 f = 0.002 的频带中绘制频谱可能是值得的(当然,您必须根据它的外观进行调整)。

2) 是的,您可以说您的数据中有一个周期为 200 秒的周期分量,它对您的原始时间序列贡献最大。这仅对 f = 0.0015 (ish) 以上的频率有效,因为这里肯定也有一些低功率信号。但我认为这很难卖,因为在 f = 0.005 附近有许多频率具有相似的功率,谐波 F 测试(如下)可以对此有所帮助。

3)是的,你当然可以这样做,除非我不确定你用什么来计算周期图(spec.pgram()?)。

现在,如果您允许的话,我想提倡一种理论上和经验上更好的方法。在实践中,周期图是一个严重偏差的估计并且也是不一致的(随着样本量的增加,方差不会收敛到零)。

由于您已经在使用 R,因此 multitaper 包实现了一个频谱估计器,该估计器使用一组正交数据锥度,以最佳方式将功率集中在感兴趣频率附近的频带(-W,W)中 - 换句话说,您实际上是在估计你认为你在估计什么(低偏差)。

因为您正在使用多个锥度,所以您可以对每个锥度的频谱估计进行平均,从而减少估计的方差(一致性)。

最后,该包实现了自适应加权,可在偏差和方差之间提供出色的平衡。

使用 multitaper 包非常简单,我将同时回答问题 3):

    install.packages("multitaper")
    library(multitaper)
    s <- spec.mtm(yourDataVector, deltat = yourSamplingRate, nw = 6, k = 10, Ftest = TRUE)
    s.df <- data.frame(x = s$freq, y = s$spec)
    plot(s$freq, s$mtm$Ftest, type='l')

注意: spec.mtm() 将从数据中删除一个恒定的平均值,除非您设置参数 center = "none" 。

好的,所以上面很简单。“nw”和“k”分别是时间带宽参数(设置波段(-W,W))和锥度数。通常,您会将 nw 设置在 3 到 15 之间的某个位置 - 取决于您认为数据中模式/周期分量的间距有多近。锥度的数量 k 应该始终是一个整数,其中 k <= 2*nw (在这种情况下,大多数情况下不希望使用更高阶的锥度)。

谐波 F 检验(我们在此处绘制)基本上是测试频率贡献的方差是否与白噪声过程的零假设下的预期相匹配。它遵循具有 (2, 2*k-2) 自由度的 F 分布。因此,无论频谱幅度如何,您都可以使用此方法检测重要频率。一个好的“经验法则”是将置信水平设置为 1 - 1/N(其中 N 是样本数)。或者您可以将其设置为您想要的任何内容。要获得 R 中的临界值:

qf(1-1/N, 2, 2*s$mtm$k - 2) .

最后,根据我的经验,在 log-y 轴上绘制光谱通常是一个好主意,因为高功率频率会洗掉其他所有东西。

plot(s$freq, s$spec, type='l', log = 'y')

例如在 base-R 图形中。

DJ Thomson - 1982 - “频谱估计和谐波分析”是该方法的开创性论文。