Rで周期解析(修正版) - meineko’s blog

上記のは、計算式が違ってました。
Mhhさんのアドヴァイスに従って修正しました。
#自乗して、マイナス~~で大きいの~~をプラスに変えるのか。

データ読み込む部分は同じです。

>d <- read.csv("cygchi.csv", header=F)

正しい、pergrmです。

pergrm <- function(d,start,end,div) {
    d$V2 <- d$V2 - mean(d$V2)
    p <- seq(start,end,length=div+1)
    pw <- numeric(div+1)
    for (i in 1:(div+1)) {
        tp <- p[i]
        ph <- ((d$V1/tp) %% 1)*pi*2
        ss <- sin(ph)
        cc <- cos(ph)
        s <- sum(d$V2*ss)^2 + sum(d$V2*cc)^2
        pw[i] <- s
    }
    return(data.frame(period=p,power=pw))
}