###################################################
### Generate Wavelet Matrix in DWT              ###
###################################################
library(wavethresh)

#Note, we actually need to transpose the matrix, because
#  it is set up for the inverse transform.
t(GenW(n=4, family="DaubExPhase", filter=1))  #Haar matrix
t(GenW(n=8, family="DaubExPhase", filter=1))  #Haar matrix

t(GenW(n=8, family="DaubLeAsymm", filter=6))  #Daubechies Least Asymmetric, 6
t(GenW(n=8, family="DaubExPhase", filter=8))  #Daubechies Extremal Phase, 8



###################################################
### Generate some data for demo purposes        ###
###################################################
set.seed(14)               #for replicability

#Periodic signal
X1 <- seq(0, 4*pi, length=2^10)
Y1 <- sin(6*X1) + 1/2*cos(2*X1) + rnorm(length(X1), sd=0.2)

#Periodic signal, but not full cycle shown
X2 <- seq(0, 3.5*pi, length=2^10)
Y2 <- sin(X2) + 1/3*cos(X2/2) + 1/6*cos(3*X2) + 1/5*cos(3*X2+0.5) + rnorm(length(X2), sd=0.2)

#Changing amplitude & frequency
X3 <- seq(1, 150, length=2^12)
Y3 <- (seq(0.5, 5, length=length(X3))*sin(X3/(seq(0.5, 5, length=length(X3))))
             + seq(0.1, 1, length=length(X3))*cos(X3/(seq(0.1, 1, length=length(X3)))) 
             + rnorm(length(X3), sd=0.2) )

###################################################
### Wavelet Coefficient Structure               ###
###################################################
library(wavethresh)
Y2wd <- wd(Y2)     #DWT for Y2
accessD(Y2wd, level=3)  #level 3 coefficients
plot(Y2wd)              #plot all wavelet coefficients




###################################################
### Wavelet thresholding Motivation             ###
###################################################
Y3wd <- wd(Y3)    #DWT for Y3

par(mfrow=c(1,2))
plot(X3, Y3, type="l", main="Observed Signal", ylim=c(-6, 3), ylab="Observation", xlab="time")
plot(seq(1, 2^12-1), Y3wd$D, main="Wavelet Coefficients",
     xlab="Vector index", ylab="Coef. size", pch=20, cex=0.75)



##################################################
### Wavelet thresholding examples              ###
##################################################
Y3wdsoft <- threshold(Y3wd)                                  #soft threshold
Y3wdhard <- threshold(Y3wd, type="hard", policy="universal") #hard threshold

newY3soft <- wr(Y3wdsoft)   #get the thresholded signal, soft
newY3hard <- wr(Y3wdhard)   #get the thresholded signal, hard

plot.ts(Y3, main="Soft Thresholding vs. True Signal", col="lightgray")
lines(seq(1, 2^12), newY3soft, col="red")
lines(seq(1, 2^12), (seq(0.5, 5, length=length(X3))*sin(X3/(seq(0.5, 5, length=length(X3))))
             + seq(0.1, 1, length=length(X3))*cos(X3/(seq(0.1, 1, length=length(X3))))), col="black")

plot.ts(Y3, main="Hard Thresholding vs. True Signal", col="lightgray")
lines(seq(1, 2^12), newY3hard, col="red")
lines(seq(1, 2^12), (seq(0.5, 5, length=length(X3))*sin(X3/(seq(0.5, 5, length=length(X3))))
             + seq(0.1, 1, length=length(X3))*cos(X3/(seq(0.1, 1, length=length(X3))))), col="black")



#####
#Terre Haute Relative Humidity Data
##You need to read in the .csv file here
#THClimate <- read.csv()

setwd('C:/Users/heyman/Desktop/IRSA talk')
THClimate <- read.csv('THClimate.csv', header=TRUE)

RelHum <- THClimate$HOURLYRelativeHumidity[1:2^10]   #Grab the series we want, to make things easier

plot(seq(1,2^10), RelHum, main="Hourly Relative Humidity in Terre Haute, Dec. 2008 - Jan. 2009",
     xlab="Hours after 12:53am on 2008-12-29", ylab="Relative Humidity", type="l")

RelHumwd <- wd(RelHum)                #do the DWT
RelHumwdS <- threshold(RelHumwd)      #soft threshold relative humidity
RelHumwdH <- threshold(RelHumwd, type="hard", policy="universal") #hard threshold relative humidity

newRelHumS <- wr(RelHumwdS)  #get the soft threshold series in time domain
newRelHumH <- wr(RelHumwdH) #get the hard threshold series in the time domain


plot.ts(RelHum, main="Soft Thresholding vs. Observed Relative Humidity", col="lightgray")
lines(seq(1, 2^10), newRelHumS, col="red")


plot.ts(RelHum, main="Hard Thresholding vs. Observed Relative Humidity", col="lightgray")
lines(seq(1, 2^10), newRelHumH, col="red")




#################################################
### Ice Core Example
#################################################
##You need to read the csv file in here.  If you set the working directory 
#   in the Terre Haute example, you should be good to go here.

Vostok <- read.csv('Vostok.csv', header=TRUE)
head(Vostok)

oldestLIDIE <- rev(Vostok$LIDIE[(dim(Vostok)[1]-2^10+1):dim(Vostok)[1]]) #grab oldest 2^10 records, reverse order
oldestDepth <- rev(-Vostok$Depth[(dim(Vostok)[1]-2^10+1):dim(Vostok)[1]]) #make it depth underground, reverse order
plot(oldestDepth, oldestLIDIE, main="Deepest LIDIE in Vostok Ice Core Chronology",
     xlab="Depth Below Surface (m)",   ylab="Lock-in depth ice equivalent (m)", type="l")



##Some questions we can try to answer
# - Does wavelet thresholding seem to extract a signal in these data?
# - What wavelet coefficients hold the most "energy"?
# - Can we use the property that coefficients are de-correlated to our advantage in analysis?
# - More?