Statistics and Data Visualizations in Climate Science

with R and Python

 

A Cambridge University Press Book by

SSP Shen and GR North

 

 

Version 1.0 released in July 2023 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
https://shen.sdsu.edu
Email:

 

 

 

 

Chapter 8: Spectral Analysis and Filtering of Time Series

R Plot Fig. 8.1: Waves of Different Amplitudes, Periods and Phases

#setwd("~/climstats")
#setEPS()
#postscript("fig0801.eps", height=6, width=10)
par(mar = c(4.5, 4.5, 2.5, 0.5))
par(mfrow=c(2,2))
t = seq(0, 2, len = 1000)
y1 = sin(2*pi*t) #A=1, tau=1, phi=0
y2 = sin(2*pi*t/0.5) #A=1, tau=0.5, phi=0
y3 = 0.5*sin(2*pi*t  - pi/2) #A=0.5, tau=1, phi=0
y4 = 0.5*sin(2*pi*t/0.5 + pi/6) #A=0.5, tau=1, phi=pi/6
plot(t, y1, 
     type = 'l', lwd = 2,
     xlab = " ", ylab = "T(t)",
     cex.lab =1.5, cex.axis = 1.5,
     cex.main = 1.5,
     main=expression(paste(
       'Function T(t): A=1, ', tau, '=1, ', phi,'=0')))
plot(t, y2, 
     type = 'l', lwd = 2,
     xlab = " ", ylab = "T(t)",
     cex.lab =1.5, cex.axis = 1.5,
     cex.main = 1.5, col = 'red',
     main=expression(paste(
       'Function T(t): A=1, ', tau, '=0.5, ', phi,'=0')))
plot(t, y3, 
     type = 'l', lwd = 2,
     xlab = "Time t", ylab = "T(t)",
     cex.lab =1.5, cex.axis = 1.5,
     cex.main = 1.5, col = 'blue',
     main=expression(paste(
       'Function T(t): A=0.5, ', tau, '=1, ', phi,'=', - pi/2)))
plot(t, y4, 
     type = 'l', lwd = 2,
     xlab = "Time t", ylab = "T(t)",
     cex.lab =1.5, cex.axis = 1.5,
     cex.main = 1.5, col = 'purple',
     main=expression(paste(
       'Function T(t): A=0.5, ', tau, '=0.5, ', phi,'=', pi/6)))

#dev.off()

 

R Plot Fig. 8.2: Wave Superposition

#setEPS()
#postscript("fig0802.eps", height=6, width=10)
par(mar = c(4.5, 4.8, 2.5, 0.5))
plot(t, y1 + y2 + 2*y3 + y4, 
     type = 'l', lwd = 4,
     xlab = "Time t", ylab = "T(t)",
     cex.lab =1.5, cex.axis = 1.5,
     cex.main = 1.5, col = 'blue',
     main='Superposition of several harmonics')

#dev.off()

 

R Code for the Discrete Sine Transform (DST)

M = 5
time = 1:M - 1/2
freq = 1:M - 1/2
time_freq = outer(time, freq)
#Construct a sine transform matrix
Phi = sqrt(2/M)*sin((pi/M)*time_freq)
#Verify Phi is an orthogonal matrix
round((t(Phi)%*%Phi), digits = 2)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1
#The output is an identity matrix

ts = time^2  #Given time series data
ts
## [1]  0.25  2.25  6.25 12.25 20.25
#[1]  0.25  2.25  6.25 12.25 20.25
ts_dst = t(Phi)%*%ts #DST transform
Recon = Phi%*%ts_dst #inverse DST 
t(Recon) #Verify Recon = ts
##      [,1] [,2] [,3]  [,4]  [,5]
## [1,] 0.25 2.25 6.25 12.25 20.25
#[1,] 0.25 2.25 6.25 12.25 20.25

 

R Plot Fig. 8.3: Polar Expression of a Complex Number

#setwd("~/climstats")
#setEPS()
#postscript("fig0803.eps", height=5.5, width=5.5)
par(mar = c(0,0,0,0))
r=10
bb=1.5*r
t=seq(0, 2*pi, length=200)
x=r*cos(t)
y=r*sin(t)
plot(x,y,type="l", lwd=3, asp=1, 
     xlim=c(-bb + 3, bb),ylim=c(-bb + 3, bb),
     xaxt="n",yaxt="n",ann=FALSE,bty="n")
aa=1.4*r
x1=r*cos(pi/3)
y1=r*sin(pi/3)
arrows(c(-aa + 2, 0), c(0, -aa + 2), c(aa,0), c(0, aa), 
       length=0.1, code=2, lwd=2)
arrows(0,0,0.98*x1,0.98*y1, 
       length=0.15,code=2, lwd=1, angle=15)
t2=seq(0,pi/3,length=10)
x2=0.22*r*cos(t2)
y2=0.22*r*sin(t2)
lines(x2,y2, type="l")
points(x1,y1,pch=19,cex=1)
segments(x1,0, x1,y1)
text(1.1*r,0.1*r, "1", cex=1.3)
text(-0.1*r, 1.1*r, "i", cex=1.5)
text(1.2*x1,1.1*y1, "(x,y)", cex=1.3)
text(0.3*r,-0.1*r, expression(paste(cos,theta)), 
     cex=1.3)
text(1.35*x1,0.5*y1, expression(paste(sin,theta)), 
     cex=1.3)
text(1.35*r,-0.1*r, "Real Axis", cex=1.3)
text(0.5*r, 1.35*r, "Imaginary Axis", cex=1.3)
text(-0.1*r, -0.1*r, "0", cex=1.5)
text(0.3*r*cos(pi/6),0.3*r*sin(pi/6), 
     expression(paste(theta)), cex=1.3)
text(1.9*x1,1.3*y1, 
expression(paste(e^{i*theta},"=", "cos",theta,"+ i sin",theta)), 
     cex=1.5)