# Wavelets: There are several packages in R for wavelet analysis:
# wavelets, waveslim, wavethresh.
# I show here the use of some simple functions (wavelet transform and
# its inverse, plotting the wavelet coefficients and the
# wavelet decomposition with the "least asymmetric Daubechies L=8 wavelets"


x <- ts(scan(url("http://stat.ethz.ch/Teaching/Datasets/ecg.dat")))
x <- ts(scan("ecg.dat"))
# Electrocardiagram measurements during the normal sinus rhythm of a patient
# who occasionally experiences arythmia. Data are in units of millivolts
# and the sampling interval is 1/180 seconds. Aus Percival&Walden (2000).
par(mfrow=c(2,1))
plot(x)
spec.pgram(x)
spec.pgram(x,spans=c(7,7))
spec.pgram(x,spans=c(15,15))
f.haar(x,6)


#filter coefficients of the "least asymmetric Daubechies L=8 wavelets"
g <- c(-0.0757657147893407,-0.0296355276459541,0.4976186676324578,
       0.8037387518052163,0.2978577956055422,-0.0992195435769354,
       -0.0126039672622612,0.0322231006040713)
h <- g[8:1]*(-1)^(0:7)
# Checking the orthogonality relations
sum(g)
sum(g[1:6]*g[3:8])
sum(g[1:4]*g[5:8])
sum(g[1:2]*g[7:8])
sum(h)
sum(h[1:6]*h[3:8])
sum(h[1:4]*h[5:8])
sum(h[1:2]*h[7:8])
sum(h*g)
sum(h[1:6]*g[3:8])
sum(g[1:6]*h[3:8])
sum(h[1:4]*g[5:8])
sum(g[1:4]*h[5:8])
sum(h[1:2]*g[7:8])
sum(g[1:2]*h[7:8])

#Transfer functions 
f.transfer <- function(coef,freq)
{
  ## Purpose: Computing the transfer function \sum coeff(u)*exp(-i*freq*u) 
  ## -------------------------------------------------------------------------
  ## Arguments: coef=filter coefficients, freq=frequency
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  2 Dec 96, 11:11
  p <- length(coef)-1
  z <- (cos(freq*2*pi)+sin(freq*2*pi)*1i)
  a  <- outer(z,0:p,FUN="^")
  a%*%coef
}
freq <- seq(0,0.5,0.01)
par(mfrow=c(2,1))
htransf <- f.transfer(h,freq)
gtransf <- f.transfer(g,freq) 
plot(freq,Mod(htransf),type="l",xlab="Frequenz",ylab="Amplitude")
lines(freq,Mod(gtransf),col=2)
plot(freq,Arg(htransf),type="l",xlab="Frequenz",ylab="Phase")
lines(freq,Arg(gtransf), col=2)
# Filtering with coefficients g dampens the high frequencies 
# and amplifies the low frequencies (low-pass filter), and
# vice versa for h.
# Both have a linear phase shift. 

# wavelet transformation 
f.wavelet <- function(x,j0,h,g)
{
  ## Purpose: Computing the wavelet transform
  ## -------------------------------------------------------------------------
  ## Arguments: x=data (if length is not a power of 2, values at the end are
  ##            removed)
  ##            j0=coarsest level to be computed
  ##            h,g=coefficients for the wavelet transform
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 29 Apr 96, modified 13 Dec 2000
  jstar <- floor(log(length(x),2))
  n <- 2^jstar
  x <- x[1:n]
  m <- n
  m2 <- n
  for (i in (1:j0)){
    m <- m/2
    a <- filter(x[1:m2],g,circular=T)[2*(1:m)]
    d <- filter(x[1:m2],h,circular=T)[2*(1:m)]
    x[1:m2] <- c(a,d)
    m2 <- m}
  x
}

w <- f.wavelet(x,6,h,g)

#Plotting of the wavelet coefficients
f.waveplot <- function(x,w,j0)
{
  ## Purpose: plotting of wavelet coefficients and original series
  ## -------------------------------------------------------------------------
  ## Arguments: x = series, w = wavelet coefficients, j0=coarsest scale
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 14 Dec 2000, 10:15
  par(mfrow=c(j0+2,1),mar=c(2,6,2,2))
  n <- length(w)
  l <- n/2^j0
  plot(w[1:l],type="h", xlab="",ylab="appr")
  abline(h=0)
  m <- l
  for (i in (j0:1)){
    plot(w[(m+1):(m+l)],type="h",xlab="t",ylab="det")
    abline(h=0)
    m <- m+l
    l <- 2*l}
  plot(x,type="l",xlab="",ylab="x")
}

f.waveplot(x,w,6)

# Inverse wavelet transform
f.waveletinv <- function(w,j0,h,g)
{
  ## Purpose: computing the inverse wavelet transform
  ## -------------------------------------------------------------------------
  ## Arguments: w=sequence, containing the wavelet coefficients
  ##            j0=coarsest level used in the transform
  ##            h,g=filter coefficients of the wavelet transforms
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 30 Apr 96, modified Dec 14, 2000
  jstar <- floor(log(length(w),2))
  n <- 2^jstar
  w <- w[1:n]
  m <- 2^(jstar-j0)
  m2 <- m
  for (i in (1:j0)){
    m2 <- 2*m
    vp <- rep(0,m2)
    vp[2*(1:m)] <- w[1:m]
    vp <- vp[m2:1]
    wp <- rep(0,m2)
    wp[2*(1:m)] <- w[(m+1):m2]
    wp <- wp[m2:1]
    w[m2:1] <- filter(vp,g,circular=T)+filter(wp,h,circular=T)
    m <- m2}
  w
}

# Testing that one is the inverse of the other
w <- f.wavelet(x,6,h,g)
max(abs(x-f.waveletinv(w,6,h,g)))

# Plotting the mother wavelet by setting all wavelet coefficients except one
# to zero. 
y <- f.waveletinv(c(rep(0,11),1,rep(0,2036)),8,h,g)
par(mfrow=c(1,1))
plot(y,type="l")

# Decomposition of a time series into components with different scales.
f.wavemra <- function(x,j0,h,g)
{
  ## Purpose: Computing and plotting of a multiresolution analysis 
  ## -------------------------------------------------------------------------
  ## Arguments: x=series of length 2^jstar
  ##            j0=scale of the coarsest approximation
  ##            h,g=filter coefficients of the wavelet
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 14 Dec 2000, 09:46
  jstar <- floor(log(length(x),2))
  n <- 2^jstar
  x <- x[1:n]
  w <- f.wavelet(x,j0,h,g)
  l <- n
  l2 <- n
  d <- matrix(0,nrow=j0+1,ncol=n)
  for (i in (1:j0)){
    l <- l/2
    v <- c(rep(0,l),w[(l+1):l2],rep(0,n-l2))
    d[i,] <- f.waveletinv(v,i,h,g)
    l2 <- l}
  v <- c(w[1:l],rep(0,n-l))
  d[j0+1,] <- f.waveletinv(v,j0,h,g)
  par(mfrow=c(j0+2,1),mar=c(2,6,2,2))
  plot(d[j0+1,],type="l",xlab="",ylab="appr")
  for (i in (j0:1)) {
    plot(d[i,],type="l",xlab="",ylab="det")}
  plot(x[1:n],type="l",xlab="",ylab="x")     
}

f.wavemra(x,6,h,g)

#same thing with Haar wavelets
x11()
f.haar(x,6)

f.haar <- function(x,j0)
{
  ## Purpose: Computing and plotting the Haar approximation of scale j0
  ##          and the details for finer scales
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 13 Dec 2000, 16:37
  n <- (length(x)%/%2^j0)*2^j0
  a <- matrix(0,nrow=j0+1,ncol=n)
  a[1,] <- x[1:n]
  y <- a[1,]
  l <- n
  b <- 1
    for (i in (1:j0)){
    l <- l/2
    b <- 2*b
    y <- 0.5*(y[2*(1:l)-1]+y[2*(1:l)])
    a[i+1,] <- rep(y,rep(b,l))}
  par(mfrow=c(j0+2,1),mar=c(2,6,2,2))
  plot(a[j0+1,],type="l",xlab="",ylab="appr")
  for (i in (j0:1)) {
    plot(a[i,]-a[i+1,],type="l",xlab="",ylab="det")}
  plot(a[1,],type="l",xlab="",ylab="x")
}