In [59]:
df<-read.csv("fluxes_DE-Hai.csv")
options(repr.plot.width=8, repr.plot.height=4)
gpp<-df$GPP
N=length(gpp)
tvec<-as.POSIXct(df$time)
plot(tvec,gpp,"l")
In [60]:
fgpp<-fft(gpp)
print(head(fgpp))
print(" ")
print(tail(fgpp))
plot(Re(fgpp),log="y")

[1] 12153.39022+ 0.00000i   -95.23078+45.60007i   -32.43182- 9.70071i
[4]     7.29781+32.97094i   -68.61068-22.10979i    15.35743-33.01480i
[1] " "
[1] -10.41620- 3.26473i  15.35743+33.01480i -68.61068+22.10979i
[4]   7.29781-32.97094i -32.43182+ 9.70071i -95.23078-45.60007i

Warning message in xy.coords(x, y, xlabel, ylabel, log):
“1849 y values <= 0 omitted from logarithmic plot”
In [61]:
length(fgpp)
3652
  • half of the information is redundant
  • we can only look at the first half
  • and set the mean to zero, anyway
In [62]:
fgpp[1]<-0

fgpp_short<-fgpp[1:((N/2)+1)]

plot(1:((N/2)+1),abs(fgpp_short)^2,log="xy")

Warning message in xy.coords(x, y, xlabel, ylabel, log):
“1 y value <= 0 omitted from logarithmic plot”

But what is the x axis in this plot?

1->Inf

2->N

3->N/2

4->N/3

...

N/2+1->2

In [63]:
x<-N/(0:(N/2))
plot(x,abs(fgpp_short)^2,log="xy",xlim=c(4000,2),xlab = "Period",ylab="Power")

Warning message in xy.coords(x, y, xlabel, ylabel, log):
“1 y value <= 0 omitted from logarithmic plot”

Now try to remove annual cycle

In [64]:
inds<-which(x>362 & x<368)
inds
11
In [65]:
fgpp_filtered<-fgpp
fgpp_filtered[inds]<-0
gpp_seasonal<-fft(fgpp_filtered,inverse=T)/N
plot(tvec,Re(gpp_seasonal),"l")

What went wrong?

Ahhh - we forgot the negative frequencies

In [66]:
get_mirror_indices<-function(inds,N) {
  N-inds+2
}
inds<-c(inds,get_mirror_indices(inds,N))
inds
  1. 11
  2. 3643
In [67]:
fgpp_filtered[inds]<-0
gpp_anomalies<-fft(fgpp_filtered,inverse=T)/N
plot(tvec,Re(gpp_anomalies),"l")

Still a lot of seasonality left -> seasonal cycle has some harmonics

In [68]:
nlower<-360;nupper<-370
inds<-which((x>nlower & x<nupper) | (x>nlower/2 & x<nupper/2) | (x>nlower/3 & x<nupper/3) | (x>nlower/4 & x<nupper/4) | (x>nlower/5 & x<nupper/5) | (x>nlower/6 & x<nupper/6))
inds<-c(inds,get_mirror_indices(inds,N))
fgpp_filtered<-fgpp
fgpp_filtered[inds]<-0
gpp_anomalies<-Re(fft(fgpp_filtered,inverse=T))/N
plot(tvec,Re(gpp_anomalies),"l")

This looks good, now let's check what exactly we just filtered out

In [69]:
fgpp_seasonal <- fgpp-fgpp_filtered
gpp_seasonal  <- Re(fft(fgpp_seasonal,inverse=T))/N
plot(tvec,Re(gpp_seasonal),"l")

The sum (in real space)

In [70]:
plot(tvec,gpp_seasonal+gpp_anomalies,"l")

Alternative: EMD

In [71]:
library(EMD)
emd_decomp<-emd(gpp)
In [72]:
par(mfrow=c(1,1))
plot(tvec,emd_decomp$imf[,1],"l",lwd=0.5,ylim=c(-5,5))
for (i in 2:ncol(emd_decomp$imf)) {
    lines(tvec,emd_decomp$imf[,i],lwd=0.5)
}

The Hilbert transform

Works only on narrow frequency bands, so lets look at the seasonal cycle

In [73]:
fgpp_seasonal_hilbert<-fgpp_seasonal
fgpp_seasonal_hilbert[1:(N/2)]<-fgpp_seasonal_hilbert[1:(N/2)]*complex(imaginary=-1)
fgpp_seasonal_hilbert[(N/2+1):N]<-fgpp_seasonal_hilbert[(N/2+1):N]*complex(imaginary=1)
gpp_seasonal_hilbert<-Re(fft(fgpp_seasonal_hilbert,inverse=T))/N
plot(tvec,Re(gpp_seasonal),"l",ylim=c(-6,6),xlab = "",ylab = "Signal + Hilbert Transform")
lines(tvec,Re(gpp_seasonal_hilbert),col=2)

We can extract the instantaneous amplitude

In [74]:
amp <- sqrt(abs(gpp_seasonal)^2+abs(gpp_seasonal_hilbert)^2)
plot(tvec,Re(gpp_seasonal),"l",ylim=c(-6,6))
lines(tvec,Re(gpp_seasonal_hilbert),col=2)
lines(tvec,amp,col=3)
In [75]:
plot(tvec,Re(gpp_seasonal),"l",ylim=c(-6,6))
lines(tvec,amp,col=3)
lines(tvec,atan2(Re(gpp_seasonal_hilbert),Re(gpp_seasonal)),col=4)

AR models

In [76]:
plot(gpp_anomalies[1:(N-1)],gpp_anomalies[2:N])
In [77]:
myfit <- lm(gpp_anomalies[1:(N-1)] ~ gpp_anomalies[2:N])
myfit

Call:
lm(formula = gpp_anomalies[1:(N - 1)] ~ gpp_anomalies[2:N])

Coefficients:
       (Intercept)  gpp_anomalies[2:N]  
         2.271e-05           3.422e-01
In [89]:
plot(tvec[2:N],residuals(myfit),"l")
In [27]:
#Create fast AR(1) function
library(rustinr)
rust('
// #[rustr_export]
pub fn generate_AR(phi: f64, noise: Vec<f64>) -> Vec<f64> {
    let mut lastval = noise[0];
    let mut v_out = vec!(lastval);
    for i in 1..noise.len() {
      let newval = lastval * phi + noise[i];
      v_out.push(newval);
      lastval = newval
    }
    return v_out
}
')
In [84]:
options(repr.plot.width=8, repr.plot.height=8)
y1<-generate_AR(0,rnorm(N))
y2<-generate_AR(0.5,rnorm(N))
y3<-generate_AR(0.9,rnorm(N))
y4<-generate_AR(1.0,rnorm(N))

par(mfrow=c(2,2))
plot(tvec,y1,"l",main="Phi=0",lwd=0.2)
plot(tvec,y2,"l",main="Phi=0.5",lwd=0.2)
plot(tvec,y3,"l",main="Phi=0.9",lwd=0.2)
plot(tvec,y4,"l",main="Phi=1.0",lwd=0.2)
In [85]:
print(cor(y3,y4))
print(cor(y1,y2))

[1] 0.06669395
[1] 0.02724317

Spectrum (acf) of an AR Process

In [86]:
par(mfrow=c(2,2))
plot(x,abs(fft(y1)[1:(N/2+1)])^2,xlim = c(N,2),log="x",main = "Spectrum Phi=0")
plot(x,abs(fft(y2)[1:(N/2+1)])^2,xlim = c(N,2),log="x",main = "Spectrum Phi=0.5")
plot(x,abs(fft(y3)[1:(N/2+1)])^2,xlim = c(N,2),log="x",main = "Spectrum Phi=0.9")
plot(x,abs(fft(y4)[1:(N/2+1)])^2,xlim = c(N,2),log="xy",main = "Spectrum Phi=1.0")
In [87]:
par(mfrow=c(2,2))
acf(y1)
acf(y2)
acf(y3)
acf(y4)
In [39]:
options(repr.plot.width=8, repr.plot.height=4)
ysurr<-generate_AR(myfit$coefficients[2],rnorm(N,sd = sd(residuals(myfit))))
plot(tvec,ysurr,"l",lwd=0.2)
lines(tvec,gpp_anomalies,col=4,lwd=0.2)