私が希望のフーリエ変換Rを使用して、時系列に変換を実行したいR

Question

に時系列にフーリエ解析を実行するには：

1. は18日、高調波
2. プロットに第5回の合計を取得します。各波
3. とし、csvファイルとして出力します。

ここではデータへのリンクです： Link to Data

ここに私の最初のコードです。

dat <- read.csv("Baguio.csv",header=FALSE) 
y  <- dat$V1 
ssp <-spectrum(y) 
t  <- 1:73 
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)] 
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)) 
rg <- diff(range(y)) 

#blue dashed line 
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg)) 
lines(fitted(reslm)~t,col=4,lty=2) 

#green line 2nd harmonics 
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t)) 
lines(fitted(reslm2)~t,col=3) 

このコードを簡素化する方法はありますか？ 18次高調波にする必要がある場合、方程式は非常に長くなります。また、ここで高調波を追加する方法はまだ分かりません。事前に

多くのおかげで、

Answer 1

はるかに簡単に解決策は、高速フーリエ変換を使用することです（fft）

dat <- read.csv("Baguio.csv", header=FALSE) 
y  <- dat$V1 
t  <- 1:73 
rg <- diff(range(y)) 

nff = function(x = NULL, n = NULL, up = 10L, plot = TRUE, add = FALSE, main = NULL, ...){ 
    #The direct transformation 
    #The first frequency is DC, the rest are duplicated 
    dff = fft(x) 
    #The time 
    t = seq(from = 1, to = length(x)) 
    #Upsampled time 
    nt = seq(from = 1, to = length(x)+1-1/up, by = 1/up) 
    #New spectrum 
    ndff = array(data = 0, dim = c(length(nt), 1L)) 
    ndff[1] = dff[1] #Always, it's the DC component 
    if(n != 0){ 
    ndff[2:(n+1)] = dff[2:(n+1)] #The positive frequencies always come first 
    #The negative ones are trickier 
    ndff[length(ndff):(length(ndff) - n + 1)] = dff[length(x):(length(x) - n + 1)] 
    } 
    #The inverses 
    indff = fft(ndff/73, inverse = TRUE) 
    idff = fft(dff/73, inverse = TRUE) 
    if(plot){ 
    if(!add){ 
     plot(x = t, y = x, pch = 16L, xlab = "Time", ylab = "Measurement", 
     main = ifelse(is.null(main), paste(n, "harmonics"), main)) 
     lines(y = Mod(idff), x = t, col = adjustcolor(1L, alpha = 0.5)) 
    } 
    lines(y = Mod(indff), x = nt, ...) 
    } 
    ret = data.frame(time = nt, y = Mod(indff)) 
    return(ret) 
} 

はその後、我々はx、数として時系列を渡し、resを呼び出す必要があります高調波をnとし、アップサンプリング（元々の時点の時点をプロットする）をupとします。

png("res_18.png") 
res = nff(x = y, n = 18L, up = 100L, col = 2L) 
dev.off() 

---

それは単にシリーズ

追加

sum5to18 = nff(x = y, n = 18L, up = 10L, plot = FALSE) 
sum5to18$y = sum5to18$y - nff(x = y, n = 4L, up = 10L, plot = FALSE)$y 
png("sum5to18.png") 
plot(sum5to18, pch = 16L, xlab = "Time", ylab = "Measurement", main = "5th to 18th harmonics sum", type = "l", col = 2) 
dev.off() 

---

の違いだ18高調波に第5回の合計を取得するには引数addとcolので、その後、CSVとして保存し、各シリーズのデータ​​を抽出方法はあり、私たちは特定の色

colors = rainbow(36L, alpha = 0.3) 
nff(x = y, n = 36L, up = 100L, col = colors[1]) 
png("all_waves.png") 
for(i in 1:18){ 
    ad = ifelse(i == 1, FALSE, TRUE) 
    nff(x = y, n = i, up = 100L, col = colors[i], add = ad, main = "All waves up to 18th harmonic") 
} 
dev.off() 

---

で、同様に複数の波をプロットすることができますファイル。この例では、18波のための18のcsvファイルが必要です。

私は0高調波（基本的には平均値）を許可するようにコードを編集したので、今はと別の波を抽出します。そして、あなたはcsvファイルを書き込むためにwrite.tableを使用することができます

sep = array(data = NA_real_, dim = c(7300L, 2 + 18), dimnames = list(NULL, c("t", paste0("H", 0:18)))) 
sep[,1:2] = as.matrix(nff(x = y, n = 0, up = 100L, plot = FALSE)) 

for(i in 1:18L){ 
    sep[,i+2] = nff(x = y, n = i, up = 100L, plot = FALSE)$y - nff(x = y, n = i-1, up = 100L, plot = FALSE)$y 
} 

。