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)