#        Math 541: Topics in Applied Mathematics: Wavelet Algorithms
#                  M.V.Wickerhauser   Sept 26, 2008
#
# Filter convolution/decimation and adjoint convolution/decimation.
#
# Assume the filter `f' is supported in the interval [0,length(f)-1]
#   and indexed by 1,2,...,length(f).
# Assume the input `u' is supported in the interval [0,length(u)-1]
#   and indexed by 1,2,...,length(u).
# In both cases, the output will be supported in [0,outlength-1] and
#   indexed by 1,2,...,outlength.

# Define some filters.
# Haar:
haarh<- c(sqrt(1/2), sqrt(1/2));
haarg<- c(sqrt(1/2), -sqrt(1/2));

# Daubechies 4:
d4h <- c(  4.82962913144534160E-01, 8.36516303737807940E-01,
	   2.24143868042013390E-01, -1.29409522551260370E-01);

d4g <- c(   1.29409522551260370E-01, 2.24143868042013390E-01,
	   -8.36516303737807940E-01, 4.82962913144534160E-01 );

# Coifman 6:
c6h <- c(  -0.072732619512526,    0.33789766245748,     0.8525720202116,   
            0.38486484686486,    -0.072732619512526,   -0.015655728135792 );

c6g <- c(  0.015655728135792,  -0.072732619512526,  -0.38486484686486, 
           0.8525720202116,    -0.33789766245748,   -0.072732619512526 );

# Discrete filter convolution/decimation
filter <- function(u, f) {
  # Compute the support interval for the output array:
  olen <- floor((length(u)+length(f))/2);   # from AWAFTTS, p.162
  # Allocate an output array of suitable length:
  out <- rep(0,olen);  
  # Pad the input array with zeroes at both ends prior to convolution:
  pad <- rep(0,length(f));
  u <- c(pad, u, pad); 
  # Compute the output values by filter convolution/decimation
  for(i in 1:olen-1) {  # index i=0,1,...,olen-1
    out[i+1] <- sum( u[2*i+1 + length(f):1] * f);  # reverses and shifts u index
  }
  return(out);
}

# Discrete filter adjoint convolution/decimation
afilter <- function(u, f) {
  # Compute the support interval for the output array:
  fhalflen <- length(f)/2
  ohalflen <- length(u)+length(f);
  olen <- 2*ohalflen;   # from AWAFTTS, p.163
  # Allocate an output array of suitable length:
  out <- rep(0,olen);  
  # Pad the input array with zeroes at both ends prior to adjoint convolution:
  pad <- rep(0,length(f)); 
  u <- c( pad, u, pad); 
  # Compute the output values by adjoint filter convolution/decimation
  half <- 1:fhalflen;   # half the f indices, starting with 1.
  fe <- f[2*half-1];    # "even"-indexed filter coefficients 
  fo <- f[2*half];      # "odd"-indexed filter coefficients
  for(j in 2*(1:ohalflen -1) ) {  # j=0,2,4,...
    startj <- j/2;
    out[j+1] <- sum( u[startj + half] * fe);  # even j case
    out[j+2] <- sum( u[startj+1 + half] * fo);  # odd j case
  }
  return(out);
}

# Plot the wavelet and scaling function for a given filter
phipsi<-function(h,g,depth) {
  phi<-afilter(1,h);  # initial scaling function approximation
  psi<-afilter(1,g);  # initial wavelet function approximation
  for(i in 1:depth) {
    phi<-afilter(phi,h);
    psi<-afilter(psi,h);
  }
  matplot(cbind(phi,psi),type="l");
  return(c(phi,psi));
}