Supplement to time series decomposition
Dave Eslinger, Marcus Beck
November 10, 2016
This document demonstrates some generic approaches to time series decomposition using Fourier transformations (fft function in the base stats package. The first examples uses generic data and the second uses SWMP data.
Generic data
#Domain setup
T=2^8.
dt <- 1/(2**4) #s
n <- T/dt
#F F/2
#freq = .01
df <- 1/T
t <- seq(0,T,by=dt) #also try ts function
#CREATE OUR FREQUENCY ARRAY
f <- 1:length(t)/T
period <- 1/f
#CREATE OUR TIME SERIES DATA
#y <- 5*sin(2*pi*1/2*t ) + 2*sin(2*pi*5/16*t) + 6*sin(2*pi*3/32*t) + 1*sin(2*pi*1/256*t)
y <- 2.5 *sin(2*pi*1*t) +
+ 1*sin(2*pi*1/2*t ) +
+ 2*sin(2*pi*5/16*t) +
+ 3*sin(2*pi*1/128*t)
#FOURIER TRANSFORM WORK
Y <- fft(y)
mag <- sqrt(Re(Y)^2+Im(Y)^2)*2/n
phase <- atan(Im(Y)/Re(Y))
Yr <- Re(Y)
Yi <- Im(Y)
mag_sorted <- sort(mag[1:(length(mag)/2)], decreasing = TRUE, index.return = TRUE)
head(mag_sorted$x,10)## [1] 3.00040640 2.48358408 1.99986860 0.99829432 0.16662343 0.14515756
## [7] 0.08112848 0.07429189 0.05384002 0.04970293#PLOTTING
plot(t[1:length(f)/2],y[1:length(f)/2],'l')
plot(f[1:length(f)/8],mag[1:length(f)/8],'l')
#Square wave, magnitude -1 to 1, based on frequency of y
#sq = (2.0 * (sin(2*pi*(1/8)*t) > 0)) - 1
sq <- 1.0 + sin(2*pi*(1/8)*t)
sq <- ifelse(sq > 1, 2.0,sq)
sq <- ifelse(sq < 1, 0.,sq)
# if(foo > 0) {
# sq = 1
# } elseif (sq < 0) {
# sq = -1
# }
#FOURIER TRANSFORM WORK
Ysq <- fft(sq)
magsq <- sqrt(Re(Ysq)^2+Im(Ysq)^2)*2/n
#PLOTTING
plot(t[1:length(f)/4],sq[1:length(f)/4],'l')
plot(f[1:length(f)/8],magsq[1:length(f)/8],'l')
#plot(t[1:length(f)/32],magsq[1:length(f)/32],'l')
# Time series that is not exact mulitple of sampling:
#CREATE OUR TIME SERIES DATA
y1 <- 3*sin(2*pi*1/4.1*t ) +5* sin(2*pi*1/16.1*t)
#FOURIER TRANSFORM WORK
Y1 <- fft(y1)
mag1 <- sqrt(Re(Y1)^2+Im(Y1)^2)*2/n
phase1 <- atan(Im(Y1)/Re(Y1))
Yr <- Re(Y1)
Yi <- Im(Y1)
mag_sorted1 <- sort(mag1[1:(length(mag1)/2)], decreasing = TRUE, index.return = TRUE)
head(mag_sorted1$x,10)## [1] 4.9303855 2.0897591 1.7235769 0.6589699 0.6030792 0.5175817 0.4346006
## [8] 0.3942738 0.3629837 0.2829360#PLOTTING
plot(t[1:length(f)/4],y1[1:length(f)/4],'l')
plot(f[1:length(f)/128],mag1[1:length(f)/128],'l')
SWMP data
# Time series test
#
## load the SWMPr package, you must have it installed
library(SWMPr)
## get data for apacpwq, all years
## import and assign to 'dat'
dat <- apacpwq
## view first six rows
head(dat)## datetimestamp temp f_temp spcond f_spcond sal f_sal do_pct
## 1 2012-01-01 00:00:00 16.9 <0> 45.97 <0> 29.9 <0> 89.3
## 2 2012-01-01 00:15:00 17.0 <0> 45.99 <0> 29.9 <0> 88.2
## 3 2012-01-01 00:30:00 17.1 <0> 44.51 <0> 28.8 <0> 88.8
## 4 2012-01-01 00:45:00 17.2 <0> 43.19 <0> 27.9 <0> 89.0
## 5 2012-01-01 01:00:00 17.2 <0> 35.30 <0> 22.3 <0> 87.7
## 6 2012-01-01 01:15:00 17.3 <0> 36.30 <0> 23.0 <0> 87.6
## f_do_pct do_mgl f_do_mgl depth f_depth cdepth f_cdepth level f_level
## 1 <0> 7.2 <0> 1.24 <0> 1.15 <3> NA <-1>
## 2 <0> 7.1 <0> 1.22 <0> 1.13 <3> NA <-1>
## 3 <0> 7.2 <0> 1.21 <0> 1.12 <3> NA <-1>
## 4 <0> 7.3 <0> 1.19 <0> 1.10 <3> NA <-1>
## 5 <0> 7.4 <0> 1.18 <0> 1.09 <3> NA <-1>
## 6 <0> 7.3 <0> 1.16 <0> 1.07 <3> NA <-1>
## clevel f_clevel ph f_ph turb f_turb chlfluor f_chlfluor
## 1 NA NA 7.9 <0> 2 <0> NA <-1>
## 2 NA NA 8.0 <0> 4 <0> NA <-1>
## 3 NA NA 8.0 <0> 4 <0> NA <-1>
## 4 NA NA 7.9 <0> 4 <0> NA <-1>
## 5 NA NA 7.9 <0> 1 <0> NA <-1>
## 6 NA NA 7.9 <0> 1 <0> NA <-1>dat <- qaqc(dat)
plot(temp ~ datetimestamp, dat,type='l')
plot(sal ~ datetimestamp, dat,type='l')
plot(do_mgl ~ datetimestamp, dat,type='l')
plot(do_pct ~ datetimestamp, dat,type='l')
plot(ph ~ datetimestamp, dat,type='l')
# Define the needed paramters to create a time series:
# Time step, i.e.sampling interval, IN DAYS:
delta_T = 15. / (60 * 24) # 15 minute data, converted into fraction of a day: 15 min/(minutes/hour * hours/ day)
# Start "time" is not really a time, but the start of the time series:
T_Start = 0.0
# End time is really the number of days worth of data:
T_End = (length(dat$temp)-1.)/(4. * 24.) # = number of intervals / (4 obs/hour * 24 hours/day) gives # of days
# As an aside: Data series that are powers of two make this more efficient. Let's not worry about that at present.
# Create an array of times for each data point:
Time <- seq(T_Start,T_End,by=delta_T)
# CREATE OUR FREQUENCY ARRAY
f <- seq(from=0.0, to=length(Time)-1.0,by=1.0)/(T_End-T_Start)
#Lets look at Water Temperature:
ts <- na.approx(dat$temp) # Fill missing (NA) values with interpolated points
bar <- matrix(c(Time,ts),ncol=2)
foobar <- subset(bar, bar[,2] < 35)
plot(Time,ts,'l',col='black',xlab="Time (Days)",ylab="Temperature") # BECAUSE YOU ALWAYS LOOK AT THE DATA FIRST!!!
ts_avg = mean(ts)
abline(h = ts_avg, col='blue') # Add a line for the mean
# Hmm, clearly a trend, so lets find and plot the linear portion:
ts_regr <- lm(ts ~ Time)
ts_regr##
## Call:
## lm(formula = ts ~ Time)
##
## Coefficients:
## (Intercept) Time
## 2.233e+01 9.813e-04abline(ts_regr,'l',col='red')
# Now, calculate that trend and remove it:
T_trend = predict.lm(ts_regr)
plot(Time,T_trend,xlim=c(T_Start,T_End)) # Plot this to check that it is correct
abline(ts_regr, 'l', col='red') #compare to slope & intercept
# And it looks correct, so remove it:
ts_detrend = ts - T_trend
plot(Time,ts_detrend,'l',xlim=c(T_Start,T_End))
# Now the data are STATIONARY, i.e., the MEAN and varibility are constant throu time.
# So lets find what frequencies are dominating the signal!
#FOURIER TRANSFORM WORK
Y <- fft(ts_detrend)
mag <- sqrt(Re(Y)^2+Im(Y)^2)*2/length(ts_detrend)
phase <- atan(Im(Y)/Re(Y))
Yr <- Re(Y)
Yi <- Im(Y)
#Plot the results
plot(f[1:length(f)/2000],mag[1:length(f)/2000],'l',
xlab="Freq., day^-1",ylab= "Power")
abline(h=1,col='red')
abline(h=2,col='red')
abline(h=5,col='red')
mag_sorted <- sort(mag[1:(length(mag)/2)], decreasing = TRUE, index.return = TRUE)
head(mag_sorted$x,10)## 3 2 13 5 7 4 12
## 7.5111968 1.0862413 0.8591905 0.8402766 0.8248885 0.7912128 0.6465453
## 6 16 15
## 0.4773758 0.4558474 0.4497400head(mag,10)## 1 2 3 4 5
## 5.609761e-12 1.086241e+00 7.511197e+00 7.912128e-01 8.402766e-01
## 6 7 8 9 10
## 4.773758e-01 8.248885e-01 2.259119e-01 4.103275e-01 2.271429e-01topMag <- mag_sorted$ix[1:10]
f[topMag]## [1] 0.002736017 0.001368009 0.016416103 0.005472034 0.008208051
## [6] 0.004104026 0.015048094 0.006840043 0.020520128 0.019152120f[4]## [1] 0.004104026f[2]## [1] 0.001368009f[18]## [1] 0.02325615print(mag_sorted$x[0:10],topMag,f[topMag])## 3 2 13 5 7 4 12 6 16 15
## 7.511 1.086 0.859 0.840 0.825 0.791 0.647 0.477 0.456 0.450