检测噪声时间序列中的周期最大值（峰值）（在 R 中？）

Question

这个问题是关于确定数字序列中最大值的数量和位置的算法。因此，这个问题有统计的味道，但它更倾向于编程，因为我对具体的统计属性不感兴趣，并且解决方案需要在R中。使用统计来回答这个问题是可以的，但不是要求。

我想提取时间序列数据（即数字的有序序列）中的最大周期。此类数据的一个示例是太阳耀斑时间序列（~11 年周期，9 至 14 年之间）。循环不会以完美的间隔重复，并且峰值并不总是相同的高度。

我发现最近有一篇论文描述了这种算法，该论文实际上使用了太阳耀斑作为示例（图 5，Scholkmann 等人，2012 年，算法）。我希望这个算法或同样有效的算法可以作为 R 包提供。

链接到 Scholkmann 论文“基于自动多尺度的峰值检测” http://www.mdpi.com/1999-4893/5/4/588

我已经尝试过“pastecs”包中的“转折点”功能，但它似乎太敏感（即检测到太多峰值）。我想首先尝试平滑时间序列，但我不确定这是否是最好的方法（我不是专家）。

感谢您的指点。

Answer 1

这是一个涉及 R 中的包的解决方案。我添加了自己的小函数，以方便在接近最大值wmtsa时搜索最大值。wmtsa::wavCWTPeaks

PeakCycle <- function(Data=as.vector(sunspots), SearchFrac=0.02){
    # using package "wmtsa"
    #the SearchFrac parameter just controls how much to look to either side 
    #of wavCWTPeaks()'s estimated maxima for a bigger value
    #see dRange
    Wave <- wavCWT(Data)
    WaveTree <- wavCWTTree(Wave)
    WavePeaks <- wavCWTPeaks(WaveTree, snr.min=5)
    WavePeaks_Times <- attr(WavePeaks, which="peaks")[,"iendtime"]

    NewPeakTimes <- c()
    dRange <- round(SearchFrac*length(Data))
    for(i in 1:length(WavePeaks_Times)){
        NewRange <- max(c(WavePeaks_Times[i]-dRange, 1)):min(c(WavePeaks_Times[i]+dRange, length(Data)))
        NewPeakTimes[i] <- which.max(Data[NewRange])+NewRange[1]-1
    }

    return(matrix(c(NewPeakTimes, Data[NewPeakTimes]), ncol=2, dimnames=list(NULL, c("PeakIndices", "Peaks"))))
}

dev.new(width=6, height=4)
par(mar=c(4,4,0.5,0.5))
plot(seq_along(as.vector(sunspots)), as.vector(sunspots), type="l")
Sunspot_Ext <- PeakCycle()
points(Sunspot_Ext, col="blue", pch=20)