Perform Fourier-Analyse zu einer Zeitreihe in R

Question

Ich möchte fourier ausführen zu einer Zeitreihe verwandeln R. mit möchte ich mag:

1. Erhalten Sie die Summe aus dem 5. bis 18. Harmonischen
2. Grundstück Jede Welle
3. und als csv-Datei ausgegeben.

Hier ist der Link zu den Daten: Link to Data

Hier ist mein ursprünglicher Code.

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) 

Gibt es eine Möglichkeit, diesen Code zu vereinfachen? Wenn ich zu den 18. Harmonischen kommen muss, wird die Gleichung sehr lang. Ich weiß auch noch nicht, wie ich die Harmonischen hier hinzufügen soll.

Vielen Dank im Voraus,

Answer 1

eine viel einfachere Lösung ist die schnelle Fourier-Transformation (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) 
} 

Dann wir res anrufen müssen verwenden, um die Zeitreihen als x vorbei, die Zahl der Harmonische wie n und das Upsampling (so dass wir Zeitpunkte neben den ursprünglichen auftragen) als up.

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

---

Um die Summe der 5. bis 18. Harmonischen zu erhalten, es ist einfach ein Unterschied zwischen Serie

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() 

---

Hinzufügen die Argumente add und col ermöglichen es uns, mehrere Wellen als auch mit spezifischen Farben

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() 

---

Gibt es eine Möglichkeit, so extrahiert die Daten jeder Reihe dann speichern als csv plotten Datei. In diesem Beispiel sollte ich 18 CSV-Dateien für die 18 Wellen haben.

I edited den Code ein 0 Harmonischen (im Grunde ein Mittelwert) zu ermöglichen, so dass nun extrahieren Sie die einzelnen Wellen als:

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 
} 

Dann können Sie write.table verwenden, um eine CSV-Datei zu schreiben.