## ======================================================================
## Copyright 1998--2010, Peter F. Craigmile, All Rights Reserved
## Address comments about this software to pfc@stat.osu.edu.
##
## The use of this software is permitted for academic purposes only,
## provided its use is acknowledged.  Commercial use is not allowed
## without the permission of the author.  This software is made available
## AS IS, and no warranty -- about the software, its performance, or its
## conformity to any specification -- is given or implied.
## ======================================================================


## ======================================================================
## File     : modwt.R
## Contains : Defining the maximum overlap discrete wavelet
##            transform (MODWT) and its inverse (IMODWT) in R.
## Version  : 1.0
## Updated  : pfc@stat.osu.edu, Apr 2010.
## Requires : 'dwt_modwt_cascades.R', 'dwt_filters.R' and 'filters.R'
##
## Notes    : This R code does not require extra dynamically allocated
##            C code, and so is very portable, but several times slower
##            than dedicated C code would be.
##            I could have optimized the code more, but I left the
##            filtering details intact for educational purposes.
##
## References:
##
##     D. B. Percival and A. T. Walden (2000),
##     Wavelet Methods for Time Series Analysis.
##     Cambridge, England: Cambridge University Press.
## ======================================================================




modwt <- function (x, wavelet = "LA8",
                   nlevels = floor(log2(length(x)))) {
  ## ======================================================================
  ## Purpose : Perform a MODWT decomposition of the data 'x' to
  ##           level 'nlevels' using a given 'wavelet'.
  ## ======================================================================
  
  filter <- dwt.filter(wavelet)
  
  js <- 1:nlevels
  
  N <- length(x)
  Nj <- N/2^js
  L  <- length(filter$scaling)
  Bj <- ceiling((L - 2) * (1 - 2^(-js)))
  Lj <- (2^js - 1) * (L - 1) + 1
  ## NOTE: P&W P198 contains details on coefficients affected by the boundaries.
  
  V <- x
  W <- list(nlevels)
  
  for (j in js) {
    results <- modwt.cascade.next(V, j, filter)
    V      <- results$V
    W[[j]] <- results$W
  }
  
  structure(list(x = x, W = W, V = V, nlevels = nlevels, N = N,
                 filter = filter),
            class = "modwt")
}




imodwt <- function (dx) {
  ## ======================================================================
  ## Purpose : Performs an inverse MODWT of the modwt class 'dx'.
  ## ======================================================================
  
  V <- dx$V
  for (j in dx$nlevels:1)
    V <- modwt.cascade.previous(V, dx$W[[j]], j, dx$filter)
  
  V
}



modwt.wavelets.zero <- function (dx, levels) {
  ## ======================================================================
  ## Purpose : Make the wavelet coefficients on level 'levels' of the
  ##           modwt class 'dx' zero.
  ## Missing : Code to zero out only non-boundary coefficients.
  ## ======================================================================

  dz <- dx
  for (j in levels)
    dz$W[[j]] <- rep(0,length(dz$W[[j]]))
  dz
}




## ======================================================================
## Purpose: Calculates the level 'j' detail coefficients of the MODWT
##          object 'dx'.
## Assumes: 1 <= j <= dx$nlevels.
## History: pfc@stat.osu.edu, Jul 2002.
## ======================================================================

modwt.detail <- function (dx, j)
{
  dz   <- modwt.wavelets.zero(dx, seq(dx$nlevels)[-j])
  dz$V <- dz$V * 0
  imodwt(dz)
}




## ======================================================================
## Purpose: Calculates the level 'j' smooth coefficients of the MODWT
##          object 'dx'.
## Assumes: 1 <= j <= dx$nlevels.
## History: pfc@stat.osu.edu, Jul 2002.
## ======================================================================

modwt.smooth <- function (dx, j)
{
  imodwt(modwt.wavelets.zero(dx, 1:j))
}