#18.12.2000 Wavelets. Die Analyse mit Haar wavelets sollten Sie aufgrund
#der letzten Stunde verstehen. Zum allgemeinen Fall werden am Montag
#noch Unterlagen abgegeben. Vergleichen Sie aber wenn moeglich bereits
#das Ergebnis von f.haar(x,6) mit dem Ergebnis von f.wavemra(x,6,h,g).
#Fuer das leztere brauchen Sie die Koeffizienten h und g sowie die
#Funktionen f.wavelet und f.waveletinv

library(ts)

#Haar wavelets
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")
}

x <- ts(scan.url("http://stat.ethz.ch/Teaching/Datasets/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)


#Filterkoeffizienten des "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)
#Nachpruefen der Orthogonalitaetsrelationen
sum(h*g)
sum(h[1:6]*h[3:8])
sum(h[1:6]*g[3:8])
sum(g[1:6]*h[3:8])
sum(g[1:6]*g[3:8])
sum(h[1:4]*h[5:8])
sum(h[1:4]*g[5:8])
sum(g[1:4]*h[5:8])
sum(g[1:4]*g[5:8])
sum(h[1:2]*h[7:8])
sum(h[1:2]*g[7:8])
sum(g[1:2]*h[7:8])
sum(g[1:2]*g[7:8])
#Transferfunktionen
f.transfer <- function(coef,freq)
{
  ## Purpose: Berechnung von \sum coeff(u)*exp(-i*freq*u) 
  ##          (sogenannte Transferfunktion)
  ## -------------------------------------------------------------------------
  ## Arguments: coef=Koeffizienten, freq=Frequenz
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  2 Dec 96, 11:11
  p <- length(coef)-1
  z <- (cos(freq)+sin(freq)*1i)
  a  <- outer(z,0:p,FUN="^")
  a%*%coef
}
freq <- seq(0,pi,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=3)
plot(freq,Arg(htransf),type="l",xlab="Frequenz",ylab="Phase")
lines(freq,Arg(gtransf), col=3)
#Man sieht: Filtern mit Koeffizienten g daempft die hohen Frequenzen und
#verstaerkt die niedrigen (low-pass filter), bei h ist es umgekehrt.
#Beide haben eine lineare Phasenverschiebung

#Wavelet Transformation allgemein
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)

#Plotten der wavelet Koeffizienten
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 Transformation
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
}

#Test, dass die eine Transformation die Inverse der andern ist.
w <- f.wavelet(x,6,h,g)
max(abs(x-f.waveletinv(w,6,h,g)))
#Plot des skalierten Wellenpakets indem man alle ausser einem wavelet
#Koeffizienten gleich null setzt
y <- f.waveletinv(c(rep(0,11),1,rep(0,2036)),8,h,g)
par(mfrow=c(1,1))
plot(y,type="l")

#Multiresolution-Analyse (Zerlegung in Approximation und Details)
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)