.packageName <- "waveslim"
wavelet.filter <- function(wf.name, filter.seq = "L", n = 512)
{
  cascade <- function(f, x, j)
    {
      L <- length(f)
      N <- length(x)
      M <- (L - 1) * 2^j
      M1 <- M - L + 2
      M2 <- 2 * M - L + 2
      if(N > M1)
        stop("x is too long\n")
      else x <- c(x, rep(0, M1 - N))
      xj <- c(rep(0, M), x, rep(0, M))
      yj <- rep(0, M2)
      for(i in 1:L)
        yj <- yj + f[L - i + 1] * xj[1:M2 + (i - 1) * 2^j]
      yj
    }
  if(is.character(wf.name))
    wf <- wave.filter(wf.name)
  else
    wf <- wf.name
  J <- nchar(filter.seq)
  key <- rev(substring(filter.seq, 1:J, 1:J))
  f <- 1
  fl <- wf$lpf
  fh <- wf$hpf
  for(k in 1:J) {
    if(key[k] == "H")
      f <- cascade(fh, f, k - 1)
    else if(key[k] == "L")
      f <- cascade(fl, f, k - 1)
    else stop("Invalid filter.seq\n")
  }
  f
}

squared.gain <- function(wf.name, filter.seq = "L", n = 512)
{
  cascade <- function(f, x, j)
    {
      L <- length(f)
      N <- length(x)
      M <- (L - 1) * 2^j
      M1 <- M - L + 2
      M2 <- 2 * M - L + 2
      if(N > M1)
        stop("x is too long\n")
      else x <- c(x, rep(0, M1 - N))
      xj <- c(rep(0, M), x, rep(0, M))
      yj <- rep(0, M2)
      for(i in 1:L)
        yj <- yj + f[L - i + 1] * xj[1:M2 + (i - 1) * 2^j]
      yj
    }
  if(is.character(wf.name))
    wf <- wave.filter(wf.name)
  else 
    wf <- wf.name
  J <- nchar(filter.seq)
  key <- rev(substring(filter.seq, 1:J, 1:J))
  f <- 1
  fl <- wf$lpf
  fh <- wf$hpf
  for(k in 1:J) {
    if(key[k] == "H")
      f <- cascade(fh, f, k - 1)
    else if(key[k] == "L")
      f <- cascade(fl, f, k - 1)
    else stop("Invalid filter.seq\n")
  }
  Mod(fft(c(f, rep(0, n - length(f))))[1:(n/2 + 1)])^2
}
## Wavelet Variance

wave.variance <- function(x, type="eta3", p=0.025)
{
  ci.gaussian <- function(x, y, p) {
    find.first <- function(v) {
      na.length <- sum(is.na(v))
      v[na.length + 1]
    }
    x.acf <- lapply(x, FUN = my.acf)
    Aj <- unlist(lapply(x.acf, FUN = function(v) sum(v*v, na.rm = TRUE))) -
      unlist(lapply(x.acf, FUN = find.first))^2 / 2
    wv.var <- 2 * Aj / unlist(lapply(x, FUN = function(v) sum(!is.na(v))))
    return(data.frame(wavevar = y, lower = y - qnorm(1-p) * sqrt(wv.var),
                      upper = y + qnorm(1-p) * sqrt(wv.var)))
  }
  
  ci.eta1 <- function(x, y, p) {
    x4 <- lapply(x, FUN = function(v) sum(v^4, na.rm = TRUE))
    eta1 <- x4.ss * unlist(lapply(x, FUN = function(v) sum(!is.na(v))))
    return(0)
  }
  
  ci.eta2 <- function(x, y, p) {
    return(0)
  }
  
  ci.eta3 <- function(x, y, p) {
    x.length <- unlist(lapply(x, FUN=function(v)sum(!is.na(v))))
    eta3 <- pmax(x.length / 2^(1:length(x)), 1)
    return(data.frame(wavevar = y, lower = eta3 * y / qchisq(1-p, eta3),
                      upper = eta3 * y / qchisq(p, eta3)))
  }

  ci.nongaussian <- function(x, y, p) {
    K <- 5
    J <- length(x)
    x.length <- unlist(lapply(x, FUN=function(v)sum(!is.na(v))))
    x.ss <- unlist(lapply(x, FUN=function(v)v[!is.na(v)]^2))
    mt.var <- numeric(J)
    for(j in 1:J) { 
      x.dpss <- dpss.taper(x.length[j], K, 4)
      V <- apply(x.dpss, 2, sum)
      J <- apply(x.dpss * x.ss[[j]], 2, sum)
      mt.var[j] <- sum((J - y[j] * V)^2) / K / x.length[j]
    }
    return(data.frame(wavevar = y, lower = y - qnorm(1-p) * sqrt(mt.var),
                      upper = y + qnorm(1-p) * sqrt(mt.var)))
  }

  x.ss <- unlist(lapply(x, FUN = function(v) sum(v*v, na.rm=TRUE)))
  x.length <- unlist(lapply(x, FUN = function(v) sum(!is.na(v))))
  y <- x.ss / x.length

  switch(type,
    gaussian = ci.gaussian(x, y, p),
    eta1 = ci.eta1(x, y, p),
    eta2 = ci.eta2(x, y, p),
    eta3 = ci.eta3(x, y, p),
    nongaussian = ci.nongaussian(x, y, p),
    stop("Invalid selection of \"type\" for the confidence interval"))
}

##plot.var <- function(x, y=NA, ylim=range(x, y, na.rm=TRUE))
##{
##  n <- dim(x)[1]
##  plot(2^(0:(n-1)), x[,1], axes=FALSE, type="n", log="xy", ylim=ylim)
##  axis(1, at=2^(0:(n-1)))
##  axis(2)
##  polyci(x[,1], x[,-1], -1)
##  if(any(!is.na(y))) { polyci(y[,1], y[,-1], 1, color=5) }
##  abline(h=0, lty=2)
##}

## Wavelet Covariance

wave.covariance <- function(x, y)
{
  my.acf.na <- function(v) {
    v <- v[!is.na(v)]
    my.acf(v)
  }
  my.ccf.na <- function(u, v) {
    u <- u[!is.na(u)]
    v <- v[!is.na(v)]
    n <- length(u)
    u <- c(u, rep(0, n))
    v <- c(v, rep(0, n))
    n <- length(u)
    x <- Re(fft(fft(u) * Conj(fft(v)), inv = TRUE)) / 2 / n^2
    x[c((n %/% 2):n, 1:(n %/% 2 - 1))]
  }
  compute.sum.xy.ccvs <- function(x, y) {
    l <- length(x)
    xy <- numeric(l)
    for(i in 1:l)
      xy[i] <- sum(my.ccf.na(x[[i]], y[[i]])^2)
    xy
  }
  compute.xy.acvs <- function(x, y) {
    l <- length(x)
    xy <- vector("list", l)
    for(i in 1:l) {
      z <- x[[i]] * y[[i]]
      xy[[i]] <- c(rev(z), z[-1])
    }
    xy
  }
  per <- function (z) {
    n <- length(z)
    (Mod(fft(z))^2/n)[1:(n%/%2 + 1)]
  }
  per2 <- function(x, y) {
    n <- length(x)
    fft.x <- fft(x)
    fft.y <- fft(y)
    ((Conj(fft.x) * fft.y)/n)[1:(n %/% 2 + 1)]
  }
  
  l <- length(x)
  xy <- vector("list", l)
  for(i in 1:l)
    xy[[i]] <- as.vector(x[[i]] * y[[i]])
  z.ss <- unlist(lapply(xy, sum, na.rm=TRUE))
  x.na <- lapply(x, is.na)
  for(i in 1:l)
    x.na[[i]] <- !x.na[[i]]
  z.length <- unlist(lapply(x.na, sum))

  zz <- z.ss / z.length
  names(zz) <- names(x)

  x.acvs <- lapply(x, my.acf.na)
  y.acvs <- lapply(y, my.acf.na)
  sum.xy.acvs <- unlist(lapply(compute.xy.acvs(x.acvs, y.acvs), sum))
  sum.squared.xy.ccvs <- compute.sum.xy.ccvs(x, y)
  var.gamma <- (sum.xy.acvs + sum.squared.xy.ccvs) / 2 / z.length

  out <- data.frame(wavecov = zz, lower = zz - qnorm(.975) * sqrt(var.gamma),
                    upper = zz + qnorm(.975) * sqrt(var.gamma))
  return(as.matrix(out))
}

##polyci <- function(x, xci, sp, color=2)
##{
##  n <- length(x)
##  y <- 2^(0:(n-1)+sp*.05)
##  delta <- y - 2^(0:(n-1))
##  for(i in 1:n){
##     polygon(c(y[i] + .6*delta[i], y[i] + .6*delta[i], y[i] - .6*delta[i],
##        y[i] - .6*delta[i]), c(xci[i,], xci[i,2:1]), border=FALSE,
##        col=color, lty=1)
##  }
##  points(y, x, pch="-")
##}

##plot.cov <- function(x, ylim=range(x,0))
##{
##  n <- dim(x)[1]
##  plot(2^(0:(n-1)), x[,1], axes=FALSE, type="n", log="x", ylim=ylim)
##  axis(1, at=2^(0:(n-1)))
##  axis(2)
##  polyci(x[,1], x[,-1], 1)
##  abline(h=0, lty=2)
##}

## Wavelet Correlation

wave.correlation <- function(x, y, N, p = .975)
{
  sum.of.squares <- function(x) { sum(x^2, na.rm=TRUE) / sum(!is.na(x)) }
  sum.of.not.squares <- function(x) { sum(x, na.rm=TRUE) / sum(!is.na(x)) }

  l <- length(x)
  xy <- vector("list", l); xy.abs <- vector("list", l)
  for(i in 1:l) {
    xy[[i]] <- as.vector(x[[i]] * y[[i]])
    xy.abs[[i]] <- as.vector(abs(x[[i]] * y[[i]]))
  }
  xy.cov <- unlist(lapply(xy, sum.of.not.squares))
  
  x.var <- unlist(lapply(x, sum.of.squares))
  y.var <- unlist(lapply(y, sum.of.squares))
  
  xy.cor <- xy.cov / sqrt(x.var * y.var)
  n <- trunc(N/2^(1:l))
  out <- data.frame(wavecor=xy.cor,
                    lower=tanh(atanh(xy.cor)-qnorm(p)/sqrt(n-3)),
                    upper=tanh(atanh(xy.cor)+qnorm(p)/sqrt(n-3)))
  return(as.matrix(out))
}

##plot.cor <- function(x, ylim=c(-1,1), cex=NULL)
##{
##  n <- dim(x)[1]
##  plot(2^(0:(n-1)), x[,1], axes=FALSE, type="n", log="x", ylim=ylim, cex=cex)
##  axis(1, at=2^(0:(n-1)), cex=cex)
##  axis(2, cex=cex)
##  polyci(x[,1], x[,-1], 1)
##  abline(h=0, lty=2)
##}

## Plotting functions for wavelet variance and covariance

## Wavelet cross-covariance 

spin.covariance <- function(x, y, lag.max = NA)
{
  xx <- zz <- x[!is.na(x)]
  yy <- y[!is.na(y)]
  n.length <- length(xx)
  xx.length <- min(length(xx)-1, lag.max, na.rm=TRUE)

  lag1 <- numeric(xx.length + 1)
  lag2 <- numeric(xx.length + 1)
  for(i in 1:(xx.length+1)) {
      lag1[i] <- sum(xx * yy, na.rm=TRUE) / n.length
      lag2[i] <- sum(zz * yy, na.rm=TRUE) / n.length
      xx <- c(xx[2:n.length], NA)
      zz <- c(NA, zz[1:(n.length-1)])
    }
  c(rev(lag2[-1]), lag1)
}

spin.correlation <- function(x, y, lag.max = NA)
{
  xx <- zz <- x[!is.na(x)]
  yy <- y[!is.na(y)]
  n.length <- length(xx)
  xx.length <- min(length(xx)-1, lag.max, na.rm=TRUE)
  xx.var <- mean(xx^2)
  yy.var <- mean(yy^2)

  lag1 <- numeric(xx.length + 1)
  lag2 <- numeric(xx.length + 1)
  for(i in 1:(xx.length+1)) {
      lag1[i] <- sum(xx * yy, na.rm=TRUE) / sqrt(xx.var * yy.var) / n.length
      lag2[i] <- sum(zz * yy, na.rm=TRUE) / sqrt(xx.var * yy.var) / n.length
      xx <- c(xx[2:n.length], NA)
      zz <- c(NA, zz[1:(n.length-1)])
    }
  c(rev(lag2[-1]), lag1)
}

##edof <- function(x) {
##  x <- x[!is.na(x)]
##  n <- length(x)
##  x.acf <- my.acf(x)
##  n * x.acf[1]^2 /
##    sum((1 - abs(seq(-n+1,n-1))/n) * c(rev(x.acf[-1]), x.acf)^2)
##}
manual.thresh <- function(wc, max.level=4, value, hard=TRUE)
{
  wc.fine <- wc[["d1"]]
  factor <- median(abs(wc.fine)) / .6745

  wc.shrink <- wc

  if(hard) {
    # Hard thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * value
      wc.shrink[[i]] <- wci * (abs(wci) > unithresh)
    }
  }
  else {
    # Soft thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * value
      wc.shrink[[i]] <- sign(wci) * (abs(wci) - unithresh) * 
        (abs(wci) > unithresh)
    }
  }
  wc.shrink
}

universal.thresh <- function(wc, max.level=4, hard=TRUE)
{
  n <- length(idwt(wc))

  wc.fine <- wc[["d1"]]
  factor <- median(abs(wc.fine)) / .6745

  wc.shrink <- wc

  if(hard) {
    # Hard thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * sqrt(2 * log(n))
      wc.shrink[[i]] <- wci * (abs(wci) > unithresh)
    }
  }
  else {
    # Soft thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * sqrt(2 * log(n))
      wc.shrink[[i]] <- sign(wci) * (abs(wci) - unithresh) *
        (abs(wci) > unithresh)
    }
  }
  wc.shrink
}

universal.thresh.modwt <- function(wc, max.level=4, hard=TRUE)
{
  n <- length(wc[[1]])
  wc.fine <- wc[["d1"]]
  factor <- sqrt(2) * median(abs(wc.fine)) / .6745
  wc.shrink <- wc
  j <- 1
  if(hard) {
    ## Hard thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * sqrt(2 * log(n)) / 2^(j/2)
      wc.shrink[[i]] <- wci * (abs(wci) > unithresh)
      j <- j+1
    }
  }
  else {
    ## Soft thresholding
    for(i in names(wc)[1:max.level]) {
      wci <- wc[[i]]
      unithresh <- factor * sqrt(2 * log(n)) / 2^(j/2)
      wc.shrink[[i]] <- sign(wci) * (abs(wci) - unithresh) *
        (abs(wci) > unithresh)
      j <- j+1
    }
  }
  wc.shrink
}

sure.thresh <- function(wc, max.level=4, hard=TRUE)
{
  wc.shrink <- wc

  sure <- function(t, x) {
    ax <- sort(abs(x))    
    num <- match(FALSE, ax <= t, nomatch = length(ax) + 1) - 1
    length(ax) - 2 * num + sum(pmin(ax, t)^2)
  }

  for(i in names(wc)[1:max.level]) {
    wci <- wc[[i]]
    ni <- length(wci)
    factor <- median(abs(wci)) / .6745
    xi <- wci / factor
    sxi <- sort(abs(xi))^2
    s <- cumsum(sxi) + ((ni - 1):0) * sxi
    risk <- (ni - (2 * (1:ni)) + s) / ni
    surethresh <- sqrt(sxi[order(risk)[1]])
    if(hard) {
      ## Hard thresholding
      wc.shrink[[i]] <- wci * (abs(xi) > surethresh)
    }
    else {
      ## Soft thresholding
      wc.shrink[[i]] <- sign(wci) * (abs(wci) - factor*surethresh) *
        (abs(xi) > surethresh)
    }
  }
  return(wc.shrink)
}

hybrid.thresh <- function(wc, max.level = 4, verbose = FALSE, seed = 0)
{
  shrinkit <- function(coeffs, thresh)
    sign(coeffs) * pmax(abs(coeffs) - thresh, 0)

  sure <- function(t, x) {
    ax <- sort(abs(x))    
    num <- match(FALSE, ax <= t, nomatch = length(ax) + 1) - 1
    length(ax) - 2 * num + sum(pmin(ax, t)^2)
  }

  wc.shrink <- wc
  n <- length(unlist(wc))
  nlev <- log(n + 1, 2) - 1
  i <- 1
  iloc <- 1
  while(i <= max.level) {
    ## Extract current level coefficients from all wavelet coefficients
    raw <- wc[[names(wc)[i]]]
    d <- length(raw)
    ## Test:  if the variance is small enough, just use threshold sqrt(2logd)
    if((sum(raw^2) - d)/d <= sqrt(i^3/2^i)) {
      if(verbose)
        cat(paste("At level ", i, " the threshhold is sqrt(2log(d)): ", 
                  sqrt(2 * log(d)), "\n", sep = ""))
      wc.shrink[[names(wc)[i]]] <- shrinkit(wc[[names(wc)[i]]], sqrt(2*log(d)))
    }
      else {
        ## Generate random subset
        if(length(seed) != 1)
          .Random.seed <- seed
        Iset <- sort(sample(d, d/2))
        rawI <- raw[Iset] / (median(abs(raw[Iset])) / .6745)
        rawIp <- raw[ - Iset] / (median(abs(raw[ - Iset])) / .6745)
        ggI <- sort(abs(rawI))
        ggIp <- sort(abs(rawIp))        
        ## Calculate SURE for all possible thresholds
        surevecI <- sapply(c(ggI[ggI < sqrt(2 * log(d))], 0,
                             sqrt(2 * log(d))), sure, ggI)
        surevecIp <- sapply(c(ggIp[ggI < sqrt(2 * log(d))], 0, 
                              sqrt(2 * log(d))), sure, ggIp)    
        ## Threshold that minimizes risk
        llI <- length(surevecI)
        llIp <- length(surevecIp)       
        ## The minimum occurs either at sqrt(2logd), 
        if(min(surevecI) == surevecI[llI])
          threshI <- sqrt(2 * log(d))
        else if(min(surevecI) == surevecI[llI - 1])
          threshI <- 0
        else threshI <- ggI[match(min(surevecI), surevecI)]     
        ## or at 0, 
        if(min(surevecIp) == surevecIp[llIp])
          threshIp <- sqrt(2 * log(d))
        else if(min(surevecIp) == surevecI[llIp - 1])
          threshIp <- 0
        else
          threshIp <- ggIp[match(min(surevecIp), surevecIp)]
        ## or at 0, 
        if(verbose) {
          cat(paste("At level ", i, ", threshold1 is ", threshI, "\n",
                    sep = ""))
          cat(paste("At level ", i, ", threshold2 is ", threshIp,
                    "\n", sep = ""))
        }
        ## Perform shrinking
        newI <- shrinkit(rawI, threshIp)
        newIp <- shrinkit(rawIp, threshI)
        new <- rep(0, d)
        new[Iset] <- newI
        new[ - Iset] <- newIp
        wc.shrink[[names(wc)[i]]] <- new
      }
    ## Otherwise, go through all this stuff
    iloc <- iloc + 2^i
    i <- i + 1
  }
  wc.shrink
}

da.thresh <- function(wc, alpha=.05, max.level=4, verbose=FALSE,
                      return.thresh=FALSE) {
  onebyone2 <- function(dat, alpha) {
    kolsmi.chi2 <- function(dat) {
      n <- length(dat)
      return(max(abs(cumsum(dat)-(1:n)*sum(dat)/n))/sqrt(2*n))
    }
    crit <- c(seq(0.28,1.49,by=.01), seq(1.50,2.48,by=.02))
    alph <- c(.999999,.999996,.999991,.999979,.999954,.999909,.999829,
              .999697,.999489,.999174,.998715,.998071,.997192,.996028,
              .994524,.992623,.990270,.987410,.983995,.979978,.975318,
              .969983,.963945,.957186,.949694,.941466,.932503,.922817,
              .912423,.901344,.889605,.877240,.864282,.850771,.836775,
              .822247,.807323,.792013,.776363,.760418,.744220,.727811,
              .711235,.694529,.677735,.660887,.644019,.627167,.610360,
              .593628,.576998,.560495,.544143,.527959,.511970,.496192,
              .480634,.465318,.450256,.435454,.420930,.406684,.392730,
              .379072,.365714,.352662,.339918,.327484,.315364,.303556,
              .292060,.280874,.270000,.259434,.249174,.239220,.229566,
              .220206,.211140,.202364,.193872,.185658,.177718,.170050,
              .162644,.155498,.148606,.141962,.135558,.129388,.123452,
              .117742,.112250,.106970,.101896,.097028,.092352,.087868,
              .083568,.079444,.075495,.071712,.068092,.064630,.061318,
              .058152,.055128,.052244,.049488,.046858,.044350,.041960,
              .039682,.037514,.035448,.033484,.031618,.029842,.028154,
              .026552,.025030,.023588,.022218,.019690,.017422,.015390,
              .013574,.011952,.010508,.009223,.008083,.007072,.006177,
              .005388,.004691,.004078,.003540,.003068,.002654,.002293,
              .001977,.001703,.001464,.001256,.001076,.000921,.000787,
              .000671,.000572,.000484,.000412,.000350,.000295,.000250,
              .000210,.000178,.000148,.000126,.000104,.000088,.000074,
              .000060,.000051,.000042,.000035,.000030,.000024,.000020,
              .000016,.000013,.000011,.000009)
    if(alpha < min(alph) || alpha > max(alph))
      stop("alpha =",alpha,"is out of range")
    ind <- match(TRUE, alpha > alph)
    critval <- crit[ind-1]+(alph[ind-1]-alpha)*(crit[ind]-crit[ind-1]) / 
      (alph[ind-1]-alph[ind])
    i <- length(dat)
    cc <- kolsmi.chi2(dat)
    while(cc[length(cc)] > critval && i > 1) {
      i <- i-1
      cc <- c(cc,kolsmi.chi2(dat[sort(order(dat)[1:i])]))
    }
    return(cc)
  }
  getthrda2 <- function(dat, alpha) {
    a <- onebyone2(dat, alpha)
    if(length(a) == length(dat))
      if(1 - pchisq(min(dat),1) < alpha)
        ggg <- 0
      else
        ggg <- sqrt(min(dat))
    else
      ggg <- sqrt(max(dat[sort(order(dat)[1:(length(dat)-length(a)+1)])]))
    return(ggg)
  }
  shrinkit <- function(coeffs, thresh)
    sign(coeffs) * pmax(abs(coeffs) - thresh, 0)

  if(alpha <= .000009 || alpha >= .999999)
    stop("alpha out of range")
  ans <- wc
  n <- length(unlist(wc))
  nlev <- log(n+1, 2)-1
  i <- 1
  iloc <- 1
  while(i <= max.level) {
    gg <- wc[[names(wc)[i]]]
    thresh <- getthrda2(gg^2,alpha)
    if(verbose)
      cat(paste("At level ",i,", the threshold is ",thresh, "\n",sep=""))
    if(return.thresh)
      if(i == nlev)
        rt <- thresh
      else
        rt <- c(thresh, rt)
    else
      ans[[names(wc)[i]]] <- shrinkit(wc[[names(wc)[i]]], thresh)
    iloc <- iloc + 2^i
    i <- i+1
  }
  if(return.thresh)
    return(rt)
  else
    return(ans)
}
dualfilt1 <- function() {

  af1 <- c(0.03516384000000, 0,
           0, 0,
           -0.08832942000000, -0.11430184000000,
           0.23389032000000, 0,
           0.76027237000000, 0.58751830000000,
           0.58751830000000, -0.76027237000000,
           0, 0.23389032000000,
           -0.11430184000000, 0.08832942000000,
           0, 0,
           0, -0.03516384000000)
  af1 <- matrix(af1, ncol=2, byrow=TRUE)
  af2 <- c(0, -0.03516384000000,
           0, 0,
           -0.11430184000000, 0.08832942000000,
           0, 0.23389032000000,
           0.58751830000000, -0.76027237000000,
           0.76027237000000, 0.58751830000000,
           0.23389032000000, 0,
           -0.08832942000000, -0.11430184000000,
           0, 0,
           0.03516384000000, 0)
  af2 <- matrix(af2, ncol=2, byrow=TRUE)
  sf1 <- af1[nrow(af1):1, ]
  sf2 <- af2[nrow(af2):1, ]
  list(af = list(af1, af2), sf = list(sf1, sf2))
}

FSfarras <- function() {

  af1 <- c(0, 0,
           -0.08838834764832,  -0.01122679215254,
           0.08838834764832,   0.01122679215254,
           0.69587998903400,   0.08838834764832,
           0.69587998903400,   0.08838834764832,
           0.08838834764832,  -0.69587998903400,
           -0.08838834764832,   0.69587998903400,
           0.01122679215254,  -0.08838834764832,
           0.01122679215254,  -0.08838834764832,
           0, 0)
  af1 <- matrix(af1, ncol=2, byrow=TRUE)
  sf1 <- af1[nrow(af1):1, ]
  af2 <- c(0.01122679215254, 0,
           0.01122679215254, 0,
           -0.08838834764832, -0.08838834764832,
           0.08838834764832, -0.08838834764832,
           0.69587998903400, 0.69587998903400,
           0.69587998903400, -0.69587998903400,
           0.08838834764832, 0.08838834764832,
           -0.08838834764832, 0.08838834764832,
           0, 0.01122679215254,
           0, -0.01122679215254)
  af2 <- matrix(af2, ncol=2, byrow=TRUE)
  sf2 <- af2[nrow(af2):1, ]
  list(af = list(af1, af2), sf = list(sf1, sf2))
}

farras <- function() {

  af <- c(0, -0.01122679215254, 0, 0.01122679215254,
          -0.08838834764832, 0.08838834764832,
          0.08838834764832, 0.08838834764832,
          0.69587998903400, -0.69587998903400,
          0.69587998903400, 0.69587998903400,
          0.08838834764832, -0.08838834764832,
          -0.08838834764832, -0.08838834764832,
          0.01122679215254, 0, 0.01122679215254, 0)
  af <- matrix(af, nrow=10, byrow=TRUE)
  sf <- af[nrow(af):1, ]
  list(af = af, sf = sf)
}

cshift <- function(x, m) {

  N <- length(x)
  n <- 0:(N-1)
  n <- (n-m) %% N
  y <- x[n+1]
  y
}

afb <- function(x, af) {

  N <- length(x)
  L <- nrow(af)/2
  x <- cshift(x,-L)
  
  ## lowpass filter
  lo <- convolve(x, af[,1], conj=FALSE, type="open")
  lo <- cshift(lo,-(2*L-1))
  lo <- lo[seq(1, length(lo), by=2)]
  lo[1:L] <- lo[N/2+(1:L)] + lo[1:L]
  lo <- lo[1:(N/2)]
  
  ## highpass filter
  hi <- convolve(x, af[,2], conj=FALSE, type="open")
  hi <- cshift(hi,-(2*L-1))
  hi <- hi[seq(1, length(hi), by=2)]
  hi[1:L] <- hi[N/2+(1:L)] + hi[1:L]
  hi <- hi[1:(N/2)]

  list(lo = lo, hi = hi)
}

dualtree <- function(x, J, Faf, af) {

  ## normalization
  x <- x/sqrt(2)
  w <- vector("list", J+1)

  ## Tree 1
  w[[1]] <- vector("list", 2)
  temp <- afb(x, Faf[[1]])
  x1 <- temp$lo
  w[[1]][[1]] <- temp$hi
  if(J > 1) {
    for(j in 2:J) {
      w[[j]] <- vector("list", 2)
      temp <- afb(x1, af[[1]])
      x1 <- temp$lo
      w[[j]][[1]] <- temp$hi
    }
  }
  w[[J+1]] <- vector("list", 2)
  w[[J+1]][[1]] <- x1

  ## Tree 2
  temp <- afb(x, Faf[[2]])
  x2 <- temp$lo
  w[[1]][[2]] <- temp$hi
  if(J > 1) {
    for(j in 2:J) {
      temp <- afb(x2, af[[2]])
      x2 <- temp$lo
      w[[j]][[2]] <- temp$hi
    }
  }
  w[[J+1]][[2]] <- x2
  w
}

sfb <- function(lo, hi, sf) {

  N <- 2*length(lo)
  L <- nrow(sf)
  ## lo <- upfirdn(lo, sf[,1], 2, 1)
  lo <- c(matrix(c(rep(0, N/2), lo), nrow=2, byrow=TRUE))
  lo <- convolve(lo, sf[,1], conj=FALSE, type="open")
  lo <- cshift(lo, -L)
  ## hi <- upfirdn(hi, sf[,2], 2, 1)
  hi <- c(matrix(c(rep(0, N/2), hi), nrow=2, byrow=TRUE))
  hi <- convolve(hi, sf[,2], conj=FALSE, type="open")
  hi <- cshift(hi, -L)
  
  y <- lo + hi
  y[1:(L-2)] <- y[1:(L-2)] + y[N+1:(L-2)]
  y <- y[1:N]
  ## y = cshift(y, 1-L/2);
  y <- cshift(y, 1-L/2)
  y
}

idualtree <- function(w, J, Fsf, sf) {

  ## Tree 1
  y1 <- w[[J+1]][[1]]
  if(J > 1) {
    for(j in J:2) {
      y1 <- sfb(y1, w[[j]][[1]], sf[[1]])
    }
  }
  y1 <- sfb(y1, w[[1]][[1]], Fsf[[1]])
  
  ## Tree 2
  y2 <- w[[J+1]][[2]]
  if(J > 1) {
    for(j in J:2) {
      y2 <- sfb(y2, w[[j]][[2]], sf[[2]])
    }
  }
  y2 <- sfb(y2, w[[1]][[2]], Fsf[[2]])

  ## normalization
  y <- (y1 + y2)/sqrt(2)
  y
}

dualtree2D <- function(x, J, Faf, af) {

  ## normalization
  x <- x/sqrt(2)
  w <- vector("list", J+1)
  
  ## Tree 1
  w[[1]] <- vector("list", 2)
  temp <- afb2D(x, Faf[[1]])              # stage 1
  x1 <- temp$lo
  w[[1]][[1]] <- temp$hi
  if(J > 1) {
    for(j in 2:J) {
      w[[j]] <- vector("list", 2)
      temp <- afb2D(x1, af[[1]])          # remaining stages
      x1 <- temp$lo
      w[[j]][[1]] <- temp$hi
    }
  }
  w[[J+1]] <- vector("list", 2)
  w[[J+1]][[1]] <- x1                     # lowpass subband

  ## Tree 2
  temp <- afb2D(x, Faf[[2]])              # stage 1
  x2 <- temp$lo
  w[[1]][[2]] <- temp$hi
  if(J > 1) {
    for(j in 2:J) {
      temp <- afb2D(x2, af[[2]])          # remaining stages
      x2 <- temp$lo
      w[[j]][[2]] <- temp$hi
    }
  }
  w[[J+1]][[2]] <- x2                     # lowpass subband

  ## sum and difference
  for(j in 1:J) {
    for(m in 1:3) {
      A <- w[[j]][[1]][[m]]
      B <- w[[j]][[2]][[m]]
      w[[j]][[1]][[m]] <- (A+B)/sqrt(2)
      w[[j]][[2]][[m]] <- (A-B)/sqrt(2)
    }
  }
  w
}

afb2D <- function(x, af1, af2=NULL) {

  if(is.null(af2))
    af2 <- af1

  ## filter along columns
  temp <- afb2D.A(x, af1, 1)
  L <- temp$lo
  H <- temp$hi
  
  ## filter along rows
  hi <- vector("list", 3)
  temp <- afb2D.A(L, af2, 2)
  lo <- temp$lo
  hi[[1]] <- temp$hi
  temp <- afb2D.A(H, af2, 2)
  hi[[2]] <- temp$lo
  hi[[3]] <- temp$hi
  list(lo = lo, hi = hi)
}

afb2D.A <- function(x, af, d) {

  lpf <- af[, 1]     # lowpass filter
  hpf <- af[, 2]     # highpass filter

  if(d == 2)
    x <- t(x)

  ## x <- matrix(1:32, 32, 64)
  N <- nrow(x)
  L <- nrow(af)/2
  x <- cshift2D(x,-L)
  ## image(x, col=rainbow(16))

  ## lo <- upfirdn(x, lpf, 1, 2)
  lo <- convolve2D(x, lpf, conj=FALSE, type="open")
  lo <- cshift2D(lo,-(2*L-1))
  lo <- lo[seq(1, nrow(lo), by=2), ]
  lo[1:L, ] <- lo[1:L, ] + lo[(1:L)+N/2, ]
  lo <- lo[1:(N/2), ]

  ## hi <- upfirdn(x, hpf, 1, 2)
  hi <- convolve2D(x, hpf, conj=FALSE, type="open")
  hi <- cshift2D(hi,-(2*L-1))
  hi <- hi[seq(1, nrow(hi), by=2), ]
  hi[1:L, ] <- hi[1:L, ] + hi[(1:L)+N/2, ]
  hi <- hi[1:(N/2), ]

  if(d == 2) {
    lo <- t(lo)
    hi <- t(hi)
  }
  list(lo = lo, hi = hi)
}

cshift2D <- function(x, m) {

  N <- nrow(x)
  n <- 0:(N-1)
  n <- (n-m) %% N
  y <- x[n+1, ]
  y
}

convolve2D <- function(x, y, conj = TRUE, type = c("circular", "open")) {
  ## Generalize convolve to handle vector arrays by calling mvfft()
  type <- match.arg(type)
  n <- nrow(x)
  ny <- length(y)
  Real <- is.numeric(x) && is.numeric(y)
  if (type == "circular") {
    if (ny != n) 
      stop("length mismatch in convolution")
  }
  else {
    n1 <- ny - 1
    x <- rbind(matrix(0, n1, ncol(x)), x)
    n <- length(y <- c(y, rep.int(0, n - 1)))
  }
  x <- mvfft(mvfft(x) * (if (conj) Conj(fft(y)) else fft(y)), inv = TRUE)
  (if (Real) Re(x) else x)/n
}

idualtree2D <- function(w, J, Fsf, sf) {

  ## sum and difference
  for(k in 1:J) {
    for(m in 1:3) {
      A <- w[[k]][[1]][[m]]
      B <- w[[k]][[2]][[m]]
      w[[k]][[1]][[m]] <- (A+B)/sqrt(2)
      w[[k]][[2]][[m]] <- (A-B)/sqrt(2)
    }
  }

  ## Tree 1
  y1 <- w[[J+1]][[1]]
  if(J > 1) {
    for(j in J:2) 
      y1 <- sfb2D(y1, w[[j]][[1]], sf[[1]])
  }
  y1 <- sfb2D(y1, w[[1]][[1]], Fsf[[1]])

  ## Tree 2
  y2 <- w[[J+1]][[2]]
  if(J > 1) {
    for(j in J:2)
      y2 <- sfb2D(y2, w[[j]][[2]], sf[[2]])
    y2 <- sfb2D(y2, w[[1]][[2]], Fsf[[2]])
  }

  ## normalization
  y <- (y1 + y2)/sqrt(2)
  y
}

sfb2D <- function(lo, hi, sf1, sf2=NULL) {

  if(is.null(sf2))
    sf2 <- sf1

  ## filter along rows
  lo <- sfb2D.A(lo, hi[[1]], sf2, 2)
  hi <- sfb2D.A(hi[[2]], hi[[3]], sf2, 2)

  ## filter along columns
  y <- sfb2D.A(lo, hi, sf1, 1)
  y
}

sfb2D.A <- function(lo, hi, sf, d) {

  lpf <- sf[, 1]     # lowpass filter
  hpf <- sf[, 2]     # highpass filter

  if(d == 2) {
   lo <- t(lo)
   hi <- t(hi)
 }

  N <- 2*nrow(lo)
  M <- ncol(lo)
  L <- nrow(sf)
  ## y = upfirdn(lo, lpf, 2, 1) + upfirdn(hi, hpf, 2, 1);
  lo <- c(matrix(c(rep(0, length(lo)), c(lo)), nrow=2, byrow=TRUE))
  lo <- matrix(lo, N, M)
  lo <- convolve2D(lo, lpf, conj=FALSE, type="open")
  lo <- cshift2D(lo, -L)
  hi <- c(matrix(c(rep(0, length(hi)), c(hi)), nrow=2, byrow=TRUE))
  hi <- matrix(hi, N, M)
  hi <- convolve2D(hi, hpf, conj=FALSE, type="open")
  hi <- cshift2D(hi, -L)
  y <- lo + hi

  y[1:(L-2), ] <- y[1:(L-2), ] + y[N+1:(L-2), ]
  y <- y[1:N, ]
  y <- cshift2D(y, 1-L/2)

  if(d == 2)
    y <- t(y)
  y
}

pm <- function(a, b) {

  u <- (a + b)/sqrt(2)
  v <- (a - b)/sqrt(2)
  list(u = u, v = v)
}

dwpt <- function(x, wf="la8", n.levels=4, boundary="periodic") {
  N <- length(x)
  J <- n.levels
  if(log(N, 2) != trunc(log(N, 2)))
    stop("Sample size is not a power of 2")
  if(2^J > N)
    stop("wavelet transform exceeds sample size in dwt")

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  y <- vector("list", sum(2^(1:J)))
  crystals1 <- rep(1:J, 2^(1:J))
  crystals2 <- unlist(apply(as.matrix(2^(1:J) - 1), 1, seq, from=0))
  names(y) <- paste("w", crystals1, ".", crystals2, sep="")

  for(j in 1:J) {
    jj <- min((1:length(crystals1))[crystals1 == j])
    for(n in 0:(2^j/2-1)) {
      if(j > 1)
        x <- y[[(1:length(crystals1))[crystals1 == j-1][n+1]]]
      W <- V <- numeric(N/2^j)
      if(n %% 2 == 0) {
        z <- .C("dwt", as.double(x), as.integer(N/2^(j-1)), L, h, g, 
	  W=as.double(W), V=as.double(V), PACKAGE="waveslim")
        y[[jj + 2*n + 1]] <- z$W
        y[[jj + 2*n]] <- z$V
      }
      else {
        z <- .C("dwt", as.double(x), as.integer(N/2^(j-1)), L, h, g,
                W=as.double(W), V=as.double(V), PACKAGE="waveslim")
        y[[jj + 2*n]] <- z$W
        y[[jj + 2*n + 1 ]] <- z$V
      }
    }
  }
  attr(y, "wavelet") <- wf
  return(y)
}

idwpt <- function(y, y.basis)
{
  J <- trunc(log(length(y), 2))

  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  for(j in J:1) {
    a <- min((1:length(rep(1:J, 2^(1:J))))[rep(1:J, 2^(1:J)) == j])
    b <- max((1:length(rep(1:J, 2^(1:J))))[rep(1:J, 2^(1:J)) == j])
    n <- a
    while(n <= b) {
      if(y.basis[n]) {
        m <- length(y[[n]])
        XX <- numeric(2 * m)
        if(floor((n-a)/2) %% 2 == 0)
          X <- .C("idwt", as.double(y[[n+1]]), as.double(y[[n]]),
                  as.integer(m), L, h, g, out=as.double(XX),
                  PACKAGE="waveslim")$out
        else
          X <- .C("idwt", as.double(y[[n]]), as.double(y[[n+1]]), 
                  as.integer(m), L, h, g, out=as.double(XX),
                  PACKAGE="waveslim")$out
        if(j != 1) {
          y[[a-(b-a+1)/2 + (n-a)/2]] <- X
          y.basis[[a-(b-a+1)/2 + (n-a)/2]] <- 1
        }
        n <- n + 2
      }
      else { n <- n + 1 }
    }
  }
  return(X)
}

##plot.dwpt <- function(x, n.levels, pgrid=TRUE)
##{
##  J <- n.levels
##  scales <- rep(1:J, 2^(1:J))
##  y <- matrix(NA, 2*length(x[[1]]), J)
##  for(j in 1:J) {
##    a <- min((1:length(scales))[scales == j])
##    b <- max((1:length(scales))[scales == j])
##    y[, j] <- unlist(x[a:b])
##    x.length <- length(y[, j])
##  }
##  plot(ts(y), ylim=c(-.45,.45))
##  if(pgrid) {
##    lines(x.length * c(0,1), c(0,0), lty=2)
##    for(j in 1:J) {
##      lines(x.length * c(0,1), c(-j,-j), lty=2)
##      for(n in 0:2^j) lines(x.length * c(n/2^j, n/2^j), c(-j,-(j-1)), lty=2)
##    }
##  }
##  title(ylab="Level")
##}

basis <- function(x, basis.names)
{
  m <- length(x)
  n <- length(basis.names)
  y <- numeric(m)
  for(i in 1:n) { y <- y + as.integer(names(x) == basis.names[i]) }
  return(y)
}

ortho.basis <- function(xtree) {
  J <- trunc(log(length(xtree), 2))
  X <- vector("list", J)
  X[[1]] <- xtree[rep(1:J, 2^(1:J)) == 1]
  for(i in 2:J) {
    for(j in i:J) {
      if(i == 2) X[[j]] <- xtree[rep(1:J, 2^(1:J)) == j]
        X[[j]] <- X[[j]] + 2 * c(apply(matrix(xtree[rep(1:J, 2^(1:J)) == i-1]),
                                       1, rep, 2^(j-i+1)))
    }
  }
  X[[J]][X[[J]] == 0] <- 1
  ifelse(unlist(X) == 1, 1, 0)
}

##plot.basis <- function(xtree)
##{
##  J <- trunc(log(length(xtree), base=2))
##  j <- rep(1:J, 2^(1:J))
##  n <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0))
##  basis <- ifelse(xtree, paste("w", j, ".", n, sep=""), NA)
##  pgrid.plot(basis[basis != "NA"])
##  invisible()
##}

phase.shift.packet <- function(z, wf, inv=FALSE)
{
  ## Center of energy
  coe <- function(g)
    sum(0:(length(g)-1) * g^2) / sum(g^2)

  J <- length(x) - 1
  g <- wave.filter(wf)$lpf
  h <- wave.filter(wf)$hpf

  if(!inv) {
    for(j in 1:J) {
      ph <- round(2^(j-1) * (coe(g) + coe(h)) - coe(g), 0)
      Nj <- length(x[[j]])
      x[[j]] <- c(x[[j]][(ph+1):Nj], x[[j]][1:ph])
    }
    ph <- round((2^J-1) * coe(g), 0)
    J <- J + 1
    x[[J]] <- c(x[[J]][(ph+1):Nj], x[[J]][1:ph])
  } else {
    for(j in 1:J) {
      ph <- round(2^(j-1) * (coe(g) + coe(h)) - coe(g), 0)
      Nj <- length(x[[j]])
      x[[j]] <- c(x[[j]][(Nj-ph+1):Nj], x[[j]][1:(Nj-ph)])
    }
    ph <- round((2^J-1) * coe(g), 0)
    J <- J + 1
    x[[J]] <- c(x[[j]][(Nj-ph+1):Nj], x[[j]][1:(Nj-ph)])
  }
  return(x)
}

modwpt <- function(x, wf="la8", n.levels=4, boundary="periodic")
{
  N <- length(x); storage.mode(N) <- "integer"
  J <- n.levels
  if(2^J > N) stop("wavelet transform exceeds sample size in modwt")

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  ht <- dict$hpf/sqrt(2)
  storage.mode(ht) <- "double"
  gt <- dict$lpf/sqrt(2)
  storage.mode(gt) <- "double"

  y <- vector("list", sum(2^(1:J)))
  yn <- length(y)
  crystals1 <- rep(1:J, 2^(1:J))
  crystals2 <- unlist(apply(as.matrix(2^(1:J) - 1), 1, seq, from=0))
  names(y) <- paste("w", crystals1, ".", crystals2, sep="")

  W <- V <- numeric(N)
  storage.mode(W) <- storage.mode(V) <- "double"
  for(j in 1:J) {
    index <- 0
    jj <- min((1:yn)[crystals1 == j])
    for(n in 0:(2^j / 2 - 1)) {
      index <- index + 1
      if(j > 1)
        x <- y[[(1:yn)[crystals1 == j-1][index]]]
      if(n %% 2 == 0) {
        z <- .C("modwt", as.double(x), N, as.integer(j), L, ht, gt, 
                W = W, V = V, PACKAGE = "waveslim")[7:8]
        y[[jj + 2*n + 1]] <- z$W
        y[[jj + 2*n]] <- z$V
      }
      else {
        z <- .C("modwt", as.double(x), N, as.integer(j), L, ht, gt, 
                W = W, V = V, PACKAGE = "waveslim")[7:8]
        y[[jj + 2*n]] <- z$W
        y[[jj + 2*n + 1 ]] <- z$V
      }
    }
  }
  attr(y, "wavelet") <- wf
  return(y)
}

dwpt.brick.wall <- function(x, wf, n.levels, method="modwpt")
{
  N <- length(x[[1]])
  m <- wave.filter(wf)$length
  J <- n.levels
  crystals1 <- rep(1:J, 2^(1:J))
  crystals2 <- unlist(apply(as.matrix(2^(1:J) - 1), 1, seq, from=0))

  if(method=="dwpt") {
    ## for DWPT
    for(j in 1:J) {
      jj <- min((1:length(crystals1))[crystals1 == j])
      L <- switch(j,
                  (m-2)/2,
                  ((m-2)/2 + floor(m/4)),
                  ((m-2)/2 + floor((m/2 + floor(m/4))/2)))
      if(is.null(L)) L <- (m-2)
      for(n in 0:(2^j-1))
        x[[jj+n]][1:L] <- NA
    }
  }
  else {
    ## for MODWPT
    for(j in 1:J) {
      jj <- min((1:length(crystals1))[crystals1 == j])
      L <- min((2^j - 1) * (m - 1), N)
      for(n in 0:(2^j-1))
        x[[jj+n]][1:L] <- NA
    }
  }
  return(x)
}

css.test <- function(y) 
{
  K <- length(y)
  test <- numeric(K)

  for(k in 1:K) {
    x <- y[[k]]
    x <- x[!is.na(x)]
    n <- length(x)
    plus <- 1:n/(n - 1) - cumsum(x^2)/sum(x^2)
    minus <- cumsum(x^2)/sum(x^2) - 0:(n - 1)/(n - 1)
    D <- max(abs(plus), abs(minus))
    if(D < 1.224/(sqrt(n) + 0.12 + 0.11/sqrt(n))) test[k] <- 1
  }
  return(test)
}

entropy.test <- function(y)
{
  K <- length(y)
  test <- numeric(K)

  for(k in 1:K) {
    x <- y[[k]]
    test[k] <- sum(x^2 * log(x^2), na.rm=TRUE)
  }
  return(test)
}

cpgram.test <- function(y, p=0.05, taper=0.1)
{
  K <- length(y)
  test <- numeric(K)
  
  for(k in 1:K) {
    x <- y[[k]]
    x <- x[!is.na(x)]
    x <- spec.taper(scale(x, center=TRUE, scale=FALSE), p=taper)
    y <- Mod(fft(x))^2/length(x)
    y[1] <- 0
    n <- length(x)
    x <- (0:(n/2))/n
    if(length(x) %% 2 == 0) {
      n <- length(x) - 1
      y <- y[1:n]
      x <- x[1:n]
    }
    else y <- y[1:length(x)]
    mp <- length(x) - 1
    if(p == 0.05)
      crit <- 1.358/(sqrt(mp) + 0.12 + 0.11/sqrt(mp))
    else {
      if(p == 0.01) crit <- 1.628/(sqrt(mp) + 0.12 + 0.11/sqrt(mp))
      else stop("critical value is not known")
    }
    D <- abs(cumsum(y)/sum(y) -  0:mp/mp)
    if(max(D) < crit) test[k] <- 1
  }
  return(test)
}

portmanteau.test <- function(y, p = 0.05, type = "Box-Pierce")
{
  K <- length(y)
  test <- numeric(K)

  for(k in 1:K) {
  x <- y[[k]]
  x <- x[!is.na(x)]
  n <- length(x)
  h <- trunc(n/2)
  x.acf <- my.acf(x)[1:(h+1)]
  x.acf <- x.acf / x.acf[1];
  if(type == "Box-Pierce")
    test[k] <- ifelse(n * sum((x.acf[-1])^2) > qchisq(1-p, h), 0, 1)
  else 
    test[k] <- ifelse(n*(n+2) * sum((x.acf[-1])^2 / (n - h:1)) > 
                      qchisq(1-p, h), 0, 1)
  }
  return(test)
}

dwpt.sim <- function(N, wf, delta, fG, M=2, adaptive=TRUE, epsilon=0.05) {
  M <- M*N
  J <- log(M, 2)
  jn <- rep(1:J, 2^(1:J))
  jl <- length(jn)

  if( adaptive ) {
    Basis <- find.adaptive.basis(wf, J, fG, epsilon) 
  } else {
    Basis <- numeric(jl)
    a <- min((1:jl)[jn == J])
    b <- max((1:jl)[jn == J])
    Basis[a:b] <- 1
  }

  Index <- (1:jl)[as.logical(Basis)]
  Length <- 2^jn

  variance <- bandpass.var.spp(delta, fG, J, Basis, Length)
  z <- vector("list", jl)
  class(z) <- "dwpt"
  attr(z, "wavelet") <- wf

  for(i in Index)
    z[[i]] <- rnorm(M/Length[i], sd=sqrt(Length[i]*variance[i]))

  x <- idwpt(z, Basis)
  xi <- trunc(runif(1, 1, M-N))
  return(x[xi:(xi+N-1)])
}

find.adaptive.basis <- function(wf, J, fG, eps) {
  H <- function(f, L) {
    H <- 0
    for(l in 0:(L/2-1))
      H <- H + choose(L/2+l-1,l) * cos(pi*f)^(2*l)
    H <- 2 * sin(pi*f)^L * H
    return(H)
  }

  G <- function(f, L) {
    G <- 0
    for(l in 0:(L/2-1))
      G <- G + choose(L/2+l-1,l) * sin(pi*f)^(2*l)
    G <- 2 * cos(pi*f)^L * G
    return(G)
  }
  
  L <- wave.filter(wf)$length
  jn <- rep(1:J, 2^(1:J))
  jl <- length(jn)
  U <- numeric(jl)
  U[1] <- G(fG, L)
  U[2] <- H(fG, L)
  for(j in 2:J) {
    jj <- min((1:jl)[jn == j])
    jp <- (1:jl)[jn == j-1]
    for(n in 0:(2^j/2-1)) {
      if (n%%2 == 0) {
        U[jj + 2 * n + 1] <- U[jp[n+1]] * H(2^(j-1)*fG, L)
        U[jj + 2 * n] <- U[jp[n+1]] * G(2^(j-1)*fG, L)
      } else {
        U[jj + 2 * n] <- U[jp[n+1]] * H(2^(j-1)*fG, L)
        U[jj + 2 * n + 1] <- U[jp[n+1]] * G(2^(j-1)*fG, L)
      }
    }
  }
  return(ortho.basis(U < eps))
}

bandpass.var.spp <- function(delta, fG, J, Basis, Length) {
  a <- unlist(sapply(2^(1:J)-1, seq, from=0, by=1)) / (2*Length)
  b <- unlist(sapply(2^(1:J), seq, from=1, by=1)) / (2*Length)
  bp.var <- rep(0, length(Basis))
  for(jn in (1:length(Basis))[as.logical(Basis)]) {
    if(fG < a[jn] | fG > b[jn])
      bp.var[jn] <- 2*integrate(spp.sdf, a[jn], b[jn], d=delta, fG=fG)$value
    else {
      result1 <- 2*integrate(spp.sdf, a[jn], fG, d=delta, fG=fG)$value
      result2 <- 2*integrate(spp.sdf, fG, b[jn], d=delta, fG=fG)$value
      bp.var[jn] <- result1 + result2
    }
  }
  return(bp.var)
}

dwt <- function(x, wf="la8", n.levels=4, boundary="periodic")
{
  switch(boundary,
    "reflection" =  x <- c(x, rev(x)),
    "periodic" = invisible(),
    stop("Invalid boundary rule in dwt"))
  N <- length(x)
  J <- n.levels
  if(N/2^J != trunc(N/2^J))
    stop("Sample size is not divisible by 2^J")
  if(2^J > N)
    stop("wavelet transform exceeds sample size in dwt")

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  y <- vector("list", J+1)
  names(y) <- c(paste("d", 1:J, sep=""), paste("s", J, sep=""))
  for(j in 1:J) {
    W <- V <- numeric(N/2^j)
    out <- .C("dwt", as.double(x), as.integer(N/2^(j-1)), L, h, g, 
              W=as.double(W), V=as.double(V), PACKAGE="waveslim")[6:7]
    y[[j]] <- out$W
    x <- out$V
  }
  y[[J+1]] <- x
  class(y) <- "dwt"
  attr(y, "wavelet") <- wf
  attr(y, "boundary") <- boundary
  return(y)
}

dwt.nondyadic <- function(x)
{
  M <- length(x)
  N <- 2^(ceiling(log(M, 2)))
  xx <- c(x, rep(0, N - M))
  y <- dwt(xx)
  
  J <- length(y) - 1
  for(j in 1:J)
    y[[j]] <- y[[j]][1:trunc(M/2^j)]
  return(y)
}

idwt <- function(y)
{
  ctmp <- class(y)
  if(is.null(ctmp) || all(ctmp != "dwt"))
    stop("argument `y' is not of class \"dwt\"")

  J <- length(y) - 1

  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  jj <- paste("s", J, sep="")
  X <- y[[jj]]
  for(j in J:1) {
    jj <- paste("d", j, sep="")
    N <- length(X)
    XX <- numeric(2 * length(y[[jj]]))
    X <- .C("idwt", as.double(y[[jj]]), as.double(X), as.integer(N), L, 
            h, g, out=as.double(XX), PACKAGE="waveslim")$out
  }
  return(X)
}

modwt <- function(x, wf="la8", n.levels=4, boundary="periodic")
{
  switch(boundary,
    "reflection" =  x <- c(x, rev(x)),
    "periodic" = invisible(),
    stop("Invalid boundary rule in modwt"))
  N <- length(x)
  storage.mode(N) <- "integer"
  J <- n.levels
  if(2^J > N)
    stop("wavelet transform exceeds sample size in modwt")

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  ht <- dict$hpf / sqrt(2)
  storage.mode(ht) <- "double"
  gt <- dict$lpf / sqrt(2)
  storage.mode(gt) <- "double"

  y <- vector("list", J+1)
  names(y) <- c(paste("d", 1:J, sep=""), paste("s", J, sep=""))
  W <- V <- numeric(N)
  storage.mode(W) <- "double"
  storage.mode(V) <- "double"
  
  for(j in 1:J) {
    out <- .C("modwt", as.double(x), N, as.integer(j), L, ht, gt, 
              W=W, V=V, PACKAGE="waveslim")[7:8]
    y[[j]] <- out$W
    x <- out$V
  }
  y[[J+1]] <- x
  class(y) <- "modwt"
  attr(y, "wavelet") <- wf
  attr(y, "boundary") <- boundary
  return(y)
}

imodwt <- function(y)
{
  ctmp <- class(y)
  if(is.null(ctmp) || all(ctmp != "modwt"))
    stop("argument `y' is not of class \"modwt\"")

  ##if(attributes(y)$boundary != "periodic")
  ##  stop("Invalid boundary rule in imodwt")
  J <- length(y) - 1

  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  ht <- dict$hpf / sqrt(2)
  storage.mode(ht) <- "double"
  gt <- dict$lpf / sqrt(2)
  storage.mode(gt) <- "double"

  jj <- paste("s", J, sep="")
  X <- y[[jj]]
  N <- length(X)
  storage.mode(N) <- "integer"
  XX <- numeric(N)
  storage.mode(XX) <- "double"
  for(j in J:1) {
    jj <- paste("d", j, sep="")
    X <- .C("imodwt", as.double(y[[jj]]), as.double(X), N, as.integer(j), 
            L, ht, gt, out=XX, PACKAGE="waveslim")$out
  }
  return(X)
}

brick.wall <- function(x, wf, method="modwt")
{
  m <- wave.filter(wf)$length
  for(j in 1:(length(x)-1)) {
    if(method == "dwt")
      n <- ceiling((m - 2) * (1 - 1/2^j))
    else
      n <- (2^j - 1) * (m - 1)
    n <- min(n, length(x[[j]]))
    x[[j]][1:n] <- NA
  }
  x[[j+1]][1:n] <- NA
  return(x)
}

phase.shift <- function(z, wf, inv=FALSE)
{
  coe <- function(g)
    sum(0:(length(g)-1) * g^2) / sum(g^2)

  J <- length(z) - 1
  g <- wave.filter(wf)$lpf
  h <- wave.filter(wf)$hpf

  if(!inv) {
    for(j in 1:J) {
      ph <- round(2^(j-1) * (coe(g) + coe(h)) - coe(g), 0)
      Nj <- length(z[[j]])
      z[[j]] <- c(z[[j]][(ph+1):Nj], z[[j]][1:ph])
    }
    ph <- round((2^J-1) * coe(g), 0)
    J <- J + 1
    z[[J]] <- c(z[[J]][(ph+1):Nj], z[[J]][1:ph])
  } else {
    for(j in 1:J) {
      ph <- round(2^(j-1) * (coe(g) + coe(h)) - coe(g), 0)
      Nj <- length(z[[j]])
      z[[j]] <- c(z[[j]][(Nj-ph+1):Nj], z[[j]][1:(Nj-ph)])
    }
    ph <- round((2^J-1) * coe(g), 0)
    J <- J + 1
    z[[J]] <- c(z[[J]][(Nj-ph+1):Nj], z[[J]][1:(Nj-ph)])
  }
  return(z)
}

mra <- function(x, wf="la8", J=4, method="modwt", boundary="periodic")
{
  switch(boundary,
         "reflection" =  x <- c(x, rev(x)),
         "periodic" = invisible(),
         stop("Invalid boundary rule in mra"))
  n <- length(x)
  
  if(method == "modwt")
    x.wt <- modwt(x, wf, J, "periodic")
  else
    x.wt <- dwt(x, wf, J, "periodic")
  x.mra <- vector("list", J+1)
  
  ## Smooth
  zero <- vector("list", J+1)
  names(zero) <- c(paste("d", 1:J, sep=""), paste("s", J, sep=""))
  class(zero) <- method
  attr(zero, "wavelet") <- wf
  attr(zero, "boundary") <- boundary
  zero[[J+1]] <- x.wt[[J+1]]
  if(method == "modwt") {
    for(k in 1:J)
      zero[[k]] <- numeric(n)
    x.mra[[J+1]] <- imodwt(zero)
  } else {
    for(k in 1:J)
      zero[[k]] <- numeric(n/2^k)
    x.mra[[J+1]] <- idwt(zero)
  }

  ## Details
  for(j in J:1) {
    zero <- vector("list", j+1)
    names(zero) <- c(paste("d", 1:j, sep=""), paste("s", j, sep=""))
    class(zero) <- method
    attr(zero, "wavelet") <- wf
    attr(zero, "boundary") <- boundary
    zero[[j]] <- x.wt[[j]]
    if(method == "modwt") {
      if(j != 1) {
        for(k in c(j+1,(j-1):1))
          zero[[k]] <- numeric(n)
      } else {
        zero[[j+1]] <- numeric(n)
      }
      x.mra[[j]] <- imodwt(zero)
    } else {
      zero[[j+1]] <- numeric(n/2^j)
      if(j != 1) {
        for(k in (j-1):1)
          zero[[k]] <- numeric(n/2^k)
      }
      x.mra[[j]] <- idwt(zero)
    }
  }

  names(x.mra) <- c(paste("D", 1:J, sep=""), paste("S", J, sep=""))
  if(boundary == "reflection") { 
    for(j in (J+1):1)
      x.mra[[j]] <- x.mra[[j]][1:(n/2)]
    return(x.mra)
  } else {
    return(x.mra)
  }
}

fdp.mle <- function(y, wf, J=log(length(y),2))
{
  fdpML <- function(d, y) {
    
    y.dwt <- y[[1]]
    n <- y[[2]]
    J <- y[[3]]

    ## Establish the limits of integration for the band-pass variances
    a <- c(1/2^c(1:J+1), 0)
    b <- 1/2^c(0:J+1)

    ## Define some useful parameters for computing the likelihood
    length.j <- n / c(2^(1:J), 2^J)
    scale.j <- c(2^(1:J+1), 2^(J+1))

    ## Initialize various parameters for computing the approximate ML
    bp.var <- numeric(J+1)

    ## Compute the band-pass variances according to d
    omega.diag <- NULL
    for(j in 1:(J+1)) {
      bp.var[j] <- integrate(fdp.sdf, a[j], b[j], d=d)$value
      omega.diag <- c(omega.diag, scale.j[j] * rep(bp.var[j], length.j[j]))
    }
    
    ## Compute approximate maximum likelihood 
    n * log(sum(y.dwt^2 / omega.diag) / n) +
      sum(length.j * log(scale.j * bp.var)) 
  }

  n <- length(y)
  y.dwt <- as.vector(unlist(dwt(y, wf, n.levels=J)))

  ## Compute MLE of d (limited to stationary region)
  result <- optimize(fdpML, int=c(-0.5,0.5), max=FALSE, y=list(y.dwt, n, J))

  ## Compute MLE of sigma_epsilon^2
  a <- c(1/2^c(1:J+1), 0)
  b <- 1/2^c(0:J+1)
  length.j <- n / c(2^(1:J), 2^J)
  scale.j <- c(2^(1:J+1), 2^(J+1))
  bp.var <- numeric(J+1)
  omega.diag <- NULL
  for(j in 1:(J+1)) {
    bp.var[j] <- integrate(fdp.sdf, a[j], b[j], d=result$minimum)$value
    omega.diag <- c(omega.diag, scale.j[j] * rep(bp.var[j], length.j[j]))
  }
  sigma2 <- sum(y.dwt^2 / omega.diag) / n

  list(parameters=c(result$minimum, sigma2), objective=result$objective)
}

########################################################################

dwt.hilbert <- function(x, wf, n.levels=4, boundary="periodic", ...)
{
  switch(boundary,
         "reflection" =  x <- c(x, rev(x)),
         "periodic" = invisible(),
         stop("Invalid boundary rule in dwt.hilbert"))
  N <- length(x)
  J <- n.levels
  if(N/2^J != trunc(N/2^J)) 
    stop("Sample size is not divisible by 2^J")
  if(2^J > N) 
    stop("Wavelet transform exceeds sample size in dwt")

  dict <- hilbert.filter(wf)
  L <- dict$length; storage.mode(L) <- "integer"
  h0 <- dict$lpf[[1]]; storage.mode(h0) <- "double"
  g0 <- dict$lpf[[2]]; storage.mode(g0) <- "double"
  h1 <- dict$hpf[[1]]; storage.mode(h1) <- "double"
  g1 <- dict$hpf[[2]]; storage.mode(g1) <- "double"

  y <- vector("list", J+1)
  names(y) <- c(paste("d", 1:J, sep=""), paste("s", J, sep=""))
  x.h <- x.g <- x
  for(j in 1:J) {
    W <- V <- numeric(N/2^j)
    out.h <- .C("dwt", as.double(x.h), as.integer(N/2^(j-1)), L, h1, h0, 
                W = W, V = V, PACKAGE = "waveslim")[6:7]
    out.g <- .C("dwt", as.double(x.g), as.integer(N/2^(j-1)), L, g1, g0, 
                W = W, V = V, PACKAGE = "waveslim")[6:7]
    y[[j]] <- complex(real = out.h$W, imaginary = out.g$W)
    x.h <- out.h$V
    x.g <- out.g$V
  }
  y[[J+1]] <- complex(real = x.h, imaginary = x.g)
  attr(y, "wavelet") <- wf
  attr(y, "levels") <- n.levels
  attr(y, "boundary") <- boundary
  return(y)
}

########################################################################

dwt.hilbert.nondyadic <- function(x, ...)
{
  M <- length(x)
  N <- 2^(ceiling(log(M, 2)))
  xx <- c(x, rep(0, N - M))
  y <- dwt.hilbert(xx, ...)
  
  J <- length(y) - 1
  for(j in 1:J) {
    y[[j]] <- y[[j]][1:trunc(M/2^j)]
  }

  return(y)
}

########################################################################

idwt.hilbert <- function(y)
{
  switch(attributes(y)$boundary,
    "reflection" =  x <- c(x, rev(x)),
    "periodic" = invisible(),
    stop("Invalid boundary rule in dwt.dbp"))
  J <- attributes(y)$levels
 
  dict <- my.filter(attributes(y)$wavelet)
  L <- dict$length; storage.mode(L) <- "integer"
  h <- dict$hpf; storage.mode(h) <- "double"
  g <- dict$lpf; storage.mode(g) <- "double"

  jj <- paste("s", J, sep="")
  X <- y[[jj]]
  for(j in J:1) {
    jj <- paste("d", j, sep="")
    XX <- numeric(2 * length(y[[jj]]))
    X <- .C("idwt", y[[jj]], as.double(X), as.integer(length(X)),
            L, h, g, XX=XX, PACKAGE="waveslim")$XX
  }
  return(X)
}

########################################################################

modwt.hilbert <- function(x, wf, n.levels=4, boundary="periodic", ...)
{
  switch(boundary,
         "reflection" =  x <- c(x, rev(x)),
         "periodic" = invisible(),
         stop("Invalid boundary rule in modwt"))
  N <- length(x)
  storage.mode(N) <- "integer"
  J <- n.levels
  if(2^J > N) stop("wavelet transform exceeds sample size in modwt")

  dict <- hilbert.filter(wf)
  L <- dict$length; storage.mode(L) <- "integer"
  h0 <- dict$lpf[[1]] / sqrt(2); storage.mode(h0) <- "double"
  g0 <- dict$lpf[[2]] / sqrt(2); storage.mode(g0) <- "double"
  h1 <- dict$hpf[[1]] / sqrt(2); storage.mode(h1) <- "double"
  g1 <- dict$hpf[[2]] / sqrt(2); storage.mode(g1) <- "double"

  y <- vector("list", J+1)
  names(y) <- c(paste("d", 1:J, sep=""), paste("s", J, sep=""))
  W <- V <- numeric(N)

  x.h <- x.g <- x  
  for(j in 1:J) {
    out.h <- .C("modwt", as.double(x.h), N, as.integer(j), L, h1, h0,
                W = W, V = V, PACKAGE="waveslim")[7:8]
    out.g <- .C("modwt", as.double(x.g), N, as.integer(j), L, g1, g0,
                W = W, V = V, PACKAGE="waveslim")[7:8]
    y[[j]] <- complex(real = out.h$W, imaginary = out.g$W)
    x.h <- out.h$V
    x.g <- out.g$V
  }
  y[[J+1]] <- complex(real = x.h, imaginary = x.g)
  attr(y, "wavelet") <- wf
  attr(y, "boundary") <- boundary
  attr(y, "levels") <- n.levels
  return(y)
}

########################################################################

imodwt.hilbert <- function(y)
{
  if(attribute(y)$boundary != "periodic")
    stop("Invalid boundary rule in imodwt")
  J <- length(y) - 1

  dict <- my.filter(attributes(y)$wavelet)
  L <- dict$length
  ht <- dict$hpf / sqrt(2)
  gt <- dict$lpf / sqrt(2)

  jj <- paste("s", J, sep="")
  X <- y[[jj]]; N <- length(X)
  XX <- numeric(N)
  for(j in J:1) {
    jj <- paste("d", j, sep="")
    X <- .C("imodwt", y[[jj]], X, as.integer(N), as.integer(j), 
      as.integer(L), ht, gt, XX, PACKAGE="waveslim")[[8]]
  }
  return(X)
}

########################################################################

hilbert.filter <- function(name)
{
  select.K3L3 <- function() {
    L <- 12
    h0 <- c(1.1594353e-04, -2.2229002e-03, -2.2046914e-03, 4.3427642e-02,
            -3.3189896e-02, -1.5642755e-01, 2.8678636e-01, 7.9972652e-01,
            4.9827824e-01, 2.4829160e-02, -4.2679177e-02, -2.2260892e-03)
    h1 <- qmf(h0)
    g0 <- c(1.6563361e-05, -5.2543406e-05, -6.1909121e-03, 1.9701141e-02,
            3.2369691e-02, -1.2705043e-01, -1.5506397e-02, 6.1333712e-01,
            7.4585008e-01, 2.1675412e-01, -4.9432248e-02, -1.5582624e-02)
    g1 <- qmf(g0)
    return(list(length = L, hpf = list(h1, g1), lpf = list(h0, g0)))
  }
    select.K3L5 <- function() {
    L <- 12
    h0 <- c(5.4258791e-06, -2.1310518e-04, -2.6140914e-03,  1.0212881e-02,
            3.5747880e-02, -4.5576766e-02, 3.9810341e-03, 5.3402475e-01,
            7.8757164e-01, 2.6537457e-01, -1.3008915e-01, -5.9573795e-02,
            1.2733976e-02, 2.8641011e-03, -2.2992683e-04, -5.8541759e-06)
    h1 <- qmf(h0)
    g0 <- c(4.9326174e-07, 3.5727140e-07, -1.1664703e-03, -8.4003116e-04,
            2.8601474e-02, 9.2509748e-03, -7.4562251e-02, 2.2929480e-01,
            7.6509138e-01, 5.8328559e-01, -4.6218010e-03, -1.2336841e-01,
            -6.2826896e-03, 9.5478911e-03, 4.6642226e-05, -6.4395935e-05)
    g1 <- qmf(g0)
    return(list(length = L, hpf = list(h1, g1), lpf = list(h0, g0)))
  }
  select.K4L2 <- function() {
    L <- 12
    h0 <- c(-1.7853301e-03, 1.3358873e-02, 3.6090743e-02, -3.4722190e-02,
            4.1525062e-02, 5.6035837e-01, 7.7458617e-01, 2.2752075e-01,
            -1.6040927e-01, -6.1694251e-02, 1.7099408e-02, 2.2852293e-03)
    h1 <- qmf(h0)
    g0 <- c(-3.5706603e-04, -1.8475351e-04, 3.2591486e-02, 1.3449902e-02,
            -5.8466725e-02, 2.7464308e-01, 7.7956622e-01, 5.4097379e-01,
            -4.0315008e-02, -1.3320138e-01, -5.9121296e-03, 1.1426146e-02)
    g1 <- qmf(g0)
    return(list(length = L, hpf = list(h1, g1), lpf = list(h0, g0)))
  }
  select.K4L4 <- function() {
    L <- 16
    h0 <- c(2.5734665593981519e-05, -6.6909066441298817e-04, 
            -5.5482443985275260e-03, 1.3203474646343588e-02,
            3.8605327384848696e-02, -5.0687259299773510e-02,
            8.1364447220208733e-03, 5.3021727476690994e-01,
            7.8330912249663232e-01, 2.7909546754271131e-01,
            -1.3372674246928601e-01, -6.9759509629953295e-02,
            1.6979390952358446e-02, 5.7323570134311854e-03,
            -6.7425216644469892e-04, -2.5933188060087743e-05)
    h1 <- qmf(h0)
    g0 <- c(2.8594072882201687e-06, 1.9074538622058143e-06,
            -2.9903835439216066e-03, -1.9808995184875909e-03,
            3.3554663884350758e-02, 7.7023844121478988e-03,
            -7.7084571412435535e-02, 2.3298110528093252e-01,
            7.5749376288995063e-01, 5.8834703992067783e-01,
            5.1708789323078770e-03, -1.3520099946241465e-01,
            -9.1961246067629732e-03, 1.5489641793018745e-02,
            1.5569563641876791e-04, -2.3339869254078969e-04)
    g1 <- qmf(g0)
    return(list(length = L, hpf = list(h1, g1), lpf = list(h0, g0)))
  }

  switch(name,
         "k3l3" = select.K3L3(),
         "k3l5" = select.K3L5(),
         "k4l2" = select.K4L2(),
         "k4l4" = select.K4L4(),
         stop("Invalid selection for hilbert.filter"))
}

########################################################################

phase.shift.hilbert <- function(x, wf)
{
  coe <- function(g)
    sum(0:(length(g)-1) * g^2) / sum(g^2)

  J <- length(x) - 1

  h0 <- hilbert.filter(wf)$lpf[[1]]
  h1 <- hilbert.filter(wf)$hpf[[1]]
  
  for(j in 1:J) {
    ph <- round(2^(j-1) * (coe(h0) + coe(h1)) - coe(h0), 0)
    Nj <- length(x[[j]])
    x[[j]] <- c(x[[j]][(ph+1):Nj], x[[j]][1:ph])
  }

  ph <- round((2^J-1) * coe(h0), 0)
  J <- J + 1
  x[[J]] <- c(x[[J]][(ph+1):Nj], x[[J]][1:ph])

  return(x)
}

########################################################################

modwpt.hilbert <- function(x, wf, n.levels=4, boundary="periodic")
{
  N <- length(x)
  storage.mode(N) <- "integer"
  J <- n.levels
  if(2^J > N) stop("wavelet transform exceeds sample size in modwpt")

  dict <- hilbert.filter(wf)
  L <- dict$length; storage.mode(L) <- "integer"
  h0 <- dict$lpf[[1]] / sqrt(2); storage.mode(h0) <- "double"
  g0 <- dict$lpf[[2]] / sqrt(2); storage.mode(g0) <- "double"
  h1 <- dict$hpf[[1]] / sqrt(2); storage.mode(h1) <- "double"
  g1 <- dict$hpf[[2]] / sqrt(2); storage.mode(g1) <- "double"

  y <- vector("list", sum(2^(1:J)))
  yn <- length(y)
  crystals1 <- rep(1:J, 2^(1:J))
  crystals2 <- unlist(apply(as.matrix(2^(1:J) - 1), 1, seq, from=0))
  names(y) <- paste("w", crystals1, ".", crystals2, sep="")

  W <- V <- numeric(N)
  storage.mode(W) <- storage.mode(V) <- "double"
  for(j in 1:J) {
    ## cat(paste("j =", j, fill=T))
    index <- 0
    jj <- min((1:yn)[crystals1 == j])
    for(n in 0:(2^j / 2 - 1)) {
      index <- index + 1
      if(j > 1)
        x <- y[[(1:yn)[crystals1 == j-1][index]]]
      else
        x <- complex(real=x, imaginary=x)

      if(n %% 2 == 0) {
        zr <- .C("modwt", as.double(Re(x)), N, as.integer(j), L, h1, h0, 
                 W = W, V = V, PACKAGE="waveslim")[7:8]
        zc <- .C("modwt", as.double(Im(x)), N, as.integer(j), L, g1, g0, 
                 W = W, V = V, PACKAGE="waveslim")[7:8]
        y[[jj + 2*n + 1]] <- complex(real=zr$W, imaginary=zc$W)
        y[[jj + 2*n]] <- complex(real=zr$V, imaginary=zc$V)
      }
      else {
        zr <- .C("modwt", as.double(Re(x)), N, as.integer(j), L, h1, h0, 
                 W = W, V = V, PACKAGE="waveslim")[7:8]
        zc <- .C("modwt", as.double(Im(x)), N, as.integer(j), L, g1, g0, 
                 W = W, V = V, PACKAGE="waveslim")[7:8]
        y[[jj + 2*n]] <- complex(real=zr$W, imaginary=zc$W)
        y[[jj + 2*n + 1 ]] <- complex(real=zr$V, imaginary=zc$V)
      }
    }
  }
  attr(y, "wavelet") <- wf
  return(y)
}

########################################################################

phase.shift.hilbert.packet <- function(x, wf)
{
  coe <- function(g)
    sum(0:(length(g)-1) * g^2) / sum(g^2)

  dict <- hilbert.filter(wf)
  h0 <- dict$lpf[[1]]; h1 <- dict$hpf[[1]]
  g0 <- dict$lpf[[2]]; g1 <- dict$hpf[[2]]

  xn <- length(x)
  N <- length(x[[1]])
  J <- trunc(log(xn,2))
  jbit <- vector("list", xn)
  jbit[[1]] <- FALSE; jbit[[2]] <- TRUE
  crystals1 <- rep(1:J, 2^(1:J))

  for(j in 1:J) {
    jj <- min((1:xn)[crystals1 == j])
    for(n in 0:(2^j - 1)) {
      if(j > 1) {
        jp <- min((1:xn)[crystals1 == j-1])
        if(n %% 4 == 0 | n %% 4 == 3)
          jbit[[jj + n]] <- c(jbit[[jp + floor(n/2)]], FALSE)
        else
          jbit[[jj + n]] <- c(jbit[[jp + floor(n/2)]], TRUE)
      }
      Sjn0 <- sum((1 - jbit[[jj + n]]) * 2^(0:(j-1)))
      Sjn1 <- sum(jbit[[jj + n]] * 2^(0:(j-1)))
      ph <- round(Sjn0 * coe(h0) + Sjn1 * coe(h1), 0)
      x[[jj + n]] <- c(x[[jj + n]][(ph+1):N], x[[jj + n]][1:ph])
    }
  }
  return(x)
}

########################################################################

modhwt.coh <- function(x, y, f.length = 0)
{
  filt <- rep(1, f.length + 1)
  filt <- filt / length(filt)

  J <- length(x) - 1
  coh <- vector("list", J)
  for(j in 1:J) {
    co.spec <- filter(Re(x[[j]] * Conj(y[[j]])), filt) 
    quad.spec <- filter(-Im(x[[j]] * Conj(y[[j]])), filt)
    x.spec <- filter(Mod(x[[j]])^2, filt)
    y.spec <- filter(Mod(y[[j]])^2, filt)
    coh[[j]] <- (co.spec^2 + quad.spec^2) / x.spec / y.spec
  }
  coh
}

########################################################################

modhwt.phase <- function(x, y, f.length = 0)
{
  filt <- rep(1, f.length + 1)
  filt <- filt / length(filt)

  J <- length(x) - 1
  phase <- vector("list", J)
  for(j in 1:J) {
    co.spec <- filter(Re(x[[j]] * Conj(y[[j]])), filt) 
    quad.spec <- filter(-Im(x[[j]] * Conj(y[[j]])), filt)
    phase[[j]] <- Arg(co.spec - 1i * quad.spec)
  }
  phase
}

########################################################################

modhwt.coh.seasonal <- function(x, y, S=10, season=365)
{
  J <- length(x) - 1
  coh <- shat <- vector("list", J)
  for(j in 1:J) {
    xj <- x[[j]]
    yj <- y[[j]]

    ## Cospectrum
    co <- matrix(Re(xj * Conj(yj)), ncol=season, byrow=TRUE)
    co.spec <- c(apply(co, 2, mean, na.rm=TRUE))
    gamma.c <- my.acf(as.vector(co))
    omega.c <- sum(gamma.c[c(1, rep(seq(season+1, S*season, by=season),
                                    each=2))])
    
    ## Quadrature spectrum
    quad <- matrix(-Im(xj * Conj(yj)), ncol=season, byrow=TRUE)
    quad.spec <- c(apply(quad, 2, mean, na.rm=TRUE))
    gamma.q <- my.acf(as.vector(quad))
    omega.q <- sum(gamma.q[c(1, rep(seq(season+1, S*season, by=season),
                                    each=2))])

    gamma.cq <- my.ccf(as.vector(co), as.vector(quad))
    omega.cq <- sum(gamma.cq[S*season + seq(-S*season+1, S*season, by=season)])

    ## Autospectrum(X)
    autoX <- matrix(Mod(xj)^2, ncol=season, byrow=TRUE)
    x.spec <- c(apply(autoX, 2, mean, na.rm=TRUE))

    ## Autospectrum(Y)
    autoY <- matrix(Mod(yj)^2, ncol=season, byrow=TRUE)
    y.spec <- c(apply(autoY, 2, mean, na.rm=TRUE))

    shat[[j]] <- 4 * (co.spec*omega.c + quad.spec * omega.q +
                      2*co.spec*quad.spec*omega.cq) / x.spec^2 / y.spec^2
    coh[[j]] <- (co.spec^2 + quad.spec^2) / x.spec / y.spec
  }
  list(coh = coh, var = shat)
}

########################################################################

modhwt.phase.seasonal <- function(x, y, season=365)
{
  J <- length(x) - 1
  phase <- vector("list", J)
  for(j in 1:J) {
    co.spec <- Re(x[[j]] * Conj(y[[j]]))
    co.spec <- c(apply(matrix(co.spec, ncol=season, byrow=TRUE), 2,
                       mean, na.rm=TRUE))
    quad.spec <- -Im(x[[j]] * Conj(y[[j]]))
    quad.spec <- c(apply(matrix(quad.spec, ncol=season, byrow=TRUE), 2,
                         mean, na.rm=TRUE))
    phase[[j]] <- Arg(co.spec - 1i * quad.spec)
  }
  phase
}

hosking.sim <- function(n, acvs) {
  .C("hosking", tseries=rnorm(n), as.integer(n), as.double(acvs[1:n]),
     PACKAGE="waveslim")$tseries
}

my.acf <- function(x)
{
  n <- length(x)
  x <- c(x, rep(0, n))
  Re(fft(Mod(fft(x))^2, inv = TRUE)/2/n^2)[1:n]
}

my.ccf <- function(a, b) {
  n <- length(a)
  a <- c(a, rep(0, n))
  b <- c(b, rep(0, n))
  x <- Re(fft(fft(a) * Conj(fft(b)), inv=TRUE))/2/n^2
  x[c((n+2):(2*n), 1:n)]
}
mra.2d <-
  function(x, wf="la8", J=4, method="modwt", boundary="periodic")
{
  m <- dim(x)[1]
  n <- dim(x)[2]

  switch(boundary,
    "periodic" = invisible(),
    stop("Invalid boundary rule in mra"))

  if(method == "modwt") {
    x.wt <- modwt.2d(x, wf, J, "periodic")
  } else {
    x.wt <- dwt.2d(x, wf, J, "periodic")
  }
  
  x.mra <- vector("list", 3*J+1)

  ## Smooth
  zero <- vector("list", 3*J+1)
  names(zero) <-
    c(matrix(rbind(paste("LH", 1:J, sep=""), paste("HL", 1:J, sep=""),
                   paste("HH", 1:J, sep="")), nrow=1), paste("LL", J, sep=""))
  attr(zero, "J") <- J
  attr(zero, "wavelet") <- wf
  attr(zero, "boundary") <- boundary
  zero[[3*J+1]] <- x.wt[[3*J+1]]
  if(method == "modwt") {
    for(k in 1:(3*J))
      zero[[k]] <- matrix(0, m, n)
    x.mra[[3*J+1]] <- imodwt.2d(zero)
  } else {
    for(k in 1:J)
      zero[[3*(k-1)+1]] <- zero[[3*(k-1)+2]] <- zero[[3*k]] <-
        matrix(0, m/2^k, n/2^k)
    x.mra[[3*J+1]] <- idwt.2d(zero)
  }

## Details
  for(j in (3*J):1) {
    Jj <- ceiling(j/3)
    zero <- vector("list", 3*Jj+1)
    names(zero) <-
      c(matrix(rbind(paste("LH", 1:Jj, sep=""), paste("HL", 1:Jj, sep=""),
                     paste("HH", 1:Jj, sep="")), nrow=1),
        paste("LL", Jj, sep=""))
    attr(zero, "J") <- Jj
    attr(zero, "wavelet") <- wf
    attr(zero, "boundary") <- boundary
    zero[[j]] <- x.wt[[j]]
    if(method == "modwt") {
      for(k in names(zero)[-charmatch(names(zero)[j], names(zero))])
        zero[[k]] <- matrix(0, m, n)
      x.mra[[j]] <- imodwt.2d(zero)
    } else {
      for(k in 1:Jj)
        zero[[3*(k-1)+1]] <- zero[[3*(k-1)+2]] <- zero[[3*k]] <-
          matrix(0, m/2^k, n/2^k)
      zero[[3*Jj+1]] <- matrix(0, m/2^Jj, n/2^Jj)
      zero[[j]] <- x.wt[[j]]
      x.mra[[j]] <- idwt.2d(zero)
    }
  }

  names(x.mra) <-
    c(matrix(rbind(paste("LH", 1:J, sep=""), paste("HL", 1:J, sep=""),
                   paste("HH", 1:J, sep="")), nrow=1), paste("LL", Jj, sep=""))
   return(x.mra)
}

rotcumvar <- function(x) {
  x <- x[!is.na(x)]
  n <- length(x)
  plus <- 1:n/(n-1) - cumsum(x^2)/sum(x^2)
  minus <- cumsum(x^2)/sum(x^2) - 0:(n-1)/(n-1)
  pmax(abs(plus), abs(minus))
}

testing.hov <- function(x, wf, J, min.coef=128, debug=FALSE) {
  n <- length(x)
  change.points <- NULL

  x.dwt <- dwt(x, wf, J)
  x.dwt.bw <- brick.wall(x.dwt, wf, method="dwt") 
  x.modwt <- modwt(x, wf, J)
  x.modwt.bw <- brick.wall(x.modwt, wf)

  for(j in 1:J) {
    cat("##### Level ", j, " #####", fill=TRUE)
    Nj <- n/2^j
    dwt.list <- list(dwt = (x.dwt.bw[[j]])[!is.na(x.dwt.bw[[j]])],
                     left = min((1:Nj)[!is.na(x.dwt.bw[[j]])]) + 1,
                     right = sum(!is.na(x.dwt.bw[[j]])))
    modwt.list <- list(modwt = (x.modwt.bw[[j]])[!is.na(x.modwt.bw[[j]])],
                       left = min((1:n)[!is.na(x.modwt.bw[[j]])]) + 1,
                       right = sum(!is.na(x.modwt.bw[[j]])))
    if(debug) cat("Starting recursion; using", dwt.list$left,
                  "to", dwt.list$right - 1, "...  ")
    change.points <-
      rbind(change.points,
            mult.loc(dwt.list, modwt.list, wf, j, min.coef, debug))
  }
  dimnames(change.points) <-
    list(NULL, c("level", "crit.value", "loc.dwt", "loc.modwt"))
  return(change.points)
}

mult.loc <- function(dwt.list, modwt.list, wf, level, min.coef, debug)
{
  Nj <- length(dwt.list$dwt)
  N <- length(modwt.list$modwt)
  crit <- 1.358
  change.points <- NULL
  
  if(Nj > min.coef) {
    ## test statistic using the DWT
    P <- cumsum(dwt.list$dwt^2) / sum(dwt.list$dwt^2)
    test.stat <- pmax((1:Nj) / (Nj-1) - P, P - (1:Nj - 1) / (Nj-1))
    loc.dwt <- (1:Nj)[max(test.stat) == test.stat]
    test.stat <- max(test.stat)

    ## location using the MODWT
    P <- cumsum(modwt.list$modwt^2) / sum(modwt.list$modwt^2)
    loc.stat <- pmax((1:N) / (N-1) - P, P - (1:N - 1) / (N-1))
    loc.modwt <- (1:N)[max(loc.stat) == loc.stat]

    if(test.stat > sqrt(2) * crit / sqrt(Nj)) {
      if(debug) cat("Accepted!", fill=TRUE)
      ## Left
      if(debug) cat("Going left; using", dwt.list$left,
                    "to", loc.dwt + dwt.list$left - 1, "...  ")
      temp.dwt.list <- list(dwt = dwt.list$dwt[1:(loc.dwt-1)],
                            left = dwt.list$left,
                            right = loc.dwt + dwt.list$left - 1)
      temp.modwt.list <- list(modwt = modwt.list$modwt[1:(loc.modwt-1)],
                              left = modwt.list$left,
                              right = loc.modwt + modwt.list$left - 1)
      change.points <-
        rbind(c(level, test.stat, loc.dwt + dwt.list$left,
                loc.modwt + modwt.list$left),
              Recall(temp.dwt.list, temp.modwt.list, wf, level, min.coef, debug))
      ## Right
      if(debug) cat("Going right; using", loc.dwt + dwt.list$left + 1,
                    "to", dwt.list$right, "...  ")
      temp.dwt.list <- list(dwt = dwt.list$dwt[(loc.dwt+1):Nj],
                            left = loc.dwt + dwt.list$left + 1,
                            right = dwt.list$right)
      temp.modwt.list <- list(modwt = modwt.list$modwt[(loc.modwt+1):N],
                              left = loc.modwt + modwt.list$left + 1,
                              right = modwt.list$right)
      change.points <-
        rbind(change.points,
              Recall(temp.dwt.list, temp.modwt.list, wf, level, min.coef, debug))
    }
    else
      if(debug) cat("Rejected!", fill=TRUE)
  }
  else
    if(debug) cat("Sample size does not exceed ", min.coef, "!",
                  sep="", fill=TRUE)

  return(change.points)
}

fdp.sdf <- function(freq, d, sigma2=1)
  sigma2 / ((2*sin(pi * freq))^2)^d

bandpass.fdp <- function(a, b, d)
  2 * integrate(fdp.sdf, lower=a, upper=b, d=d)$value

spp.sdf <- function(freq, d, fG, sigma2=1)
  sigma2 * abs(2 * (cos(2*pi*freq) - cos(2*pi*fG)))^(-2*d)

spp2.sdf <- function(freq, d1, f1, d2, f2, sigma2=1) {
  sigma2 * abs(2 * (cos(2*pi*freq) - cos(2*pi*f1)))^(-2 * d1) * 
    abs(2 * (cos(2*pi*freq) - cos(2*pi*f2)))^(-2 * d2)
}

sfd.sdf <- function(freq, s, d, sigma2=1)
  sigma2 / (2 * (1 - cos(s * 2*pi*freq)))^d

bandpass.spp <- function(a, b, d, fG) {
  if(fG > a && fG < b) {
    result1 <- integrate(spp.sdf, lower=a, upper=fG, d=d, fG=fG)$value
    result2 <- integrate(spp.sdf, lower=fG, upper=b, d=d, fG=fG)$value
  }
  else {
    result1 <- integrate(spp.sdf, lower=a, upper=b, d=d, fG=fG)$value
    result2 <- 0
  }
  return(2*(result1 + result2))
}

bandpass.spp2 <- function(a, b, d1, f1, d2, f2) {
  a1 <- a
  b1 <- b
  if(a1 < f1 && b1 > f2) {
    a2 <- f1
    b2 <- f2
    result1 <- integrate(spp2.sdf, a1, a2, d1=d1, f1=f1, d2=d2, f2=f2)$value
    result2 <- integrate(spp2.sdf, a1, b2, d1=d1, f1=f1, d2=d2, f2=f2)$value
    result3 <- integrate(spp2.sdf, b2, b1, d1=d1, f1=f1, d2=d2, f2=f2)$value
  }
  else {
    if(a1 < f1 && b1 < f2) {
      a2 <- f1
      result1 <- integrate(spp2.sdf, a1, a2, d1=d1, f1=f1, d2=d2, f2=f2)$value
      result2 <- integrate(spp2.sdf, a2, b1, d1=d1, f1=f1, d2=d2, f2=f2)$value
      result3 <- 0
    }
    else {
      if(a1 < f1 && b1 > f1 && b1 < f2) {
        a2 <- f1
        result1 <- integrate(spp2.sdf, a1, a2, d1=d1, f1=f1, d2=d2, f2=f2)$value
        result2 <- integrate(spp2.sdf, a2, b1, d1=d1, f1=f1, d2=d2, f2=f2)$value
        result3 <- 0
      }
      else {
        if(a1 > f1 && a1 < f2 && b1 > f2) {
          a2 <- f2
        result1 <- integrate(spp2.sdf, a1, a2, d1=d1, f1=f1, d2=d2, f2=f2)$value
        result2 <- integrate(spp2.sdf, a2, b1, d1=d1, f1=f1, d2=d2, f2=f2)$value
        result3 <- 0
        }
        else {
          result1 <- integrate(spp2.sdf, a1, b1, d1=d1, f1=f1, d2=d2, f2=f2)$value
          result2 <- 0
          result3 <- 0
        }
      }
    }
  }
  return(2*(result1 + result2 + result3))
}

Hypergeometric <- function(a, b, c, z) {
  ## Recursive implementation taken from Numerical Recipes in C (6.12)
  ## Press, Teukolsky, Vetterling and Flannery (1992)
  fac <- 1
  temp <- fac
  aa <- a
  bb <- b
  cc <- c
  for(n in 1:1000) {
    fac <- fac * (aa * bb) / cc
    fac <- fac * z / n
    series <- temp + fac
    if(series == temp)
      return(series)
    temp <- series
    aa <- aa + 1
    bb <- bb + 1
    cc <- cc + 1
  }
  stop("convergence failure in Hypergeometric")
}

spp.var <- function(d, fG, sigma2=1) {
  ## Hypergeometric series representation of the variance taken from
  ## Lapsa (1997)
  omega <- 2*pi*fG
  A <- sigma2/2/sqrt(pi) * gamma(1-2*d) / gamma(3/2 - 2*d) * sin(omega)^(1-4*d)
  P1 <- Hypergeometric(1-2*d, 1-2*d, 3/2 - 2*d, sin(omega/2)^2)
  P2 <- Hypergeometric(1-2*d, 1-2*d, 3/2 - 2*d, cos(omega/2)^2)
  return(A * (P1 + P2))
}
shift.2d <- function(z, wf, inv=FALSE)
{
  coe <- function(g) sum(0:(length(g)-1) * g^2) / sum(g^2)
  
  h <- wave.filter(wf)$hpf
  g <- wave.filter(wf)$lpf
  
  J <- (length(z) - 1) / 3
  n <- dim(z[[1]])[1]
  
  nu.H <- round(2^(1:J-1) * (coe(g) + coe(h)) - coe(g), 0)
  nu.H <- ifelse(nu.H/n < 1, nu.H, nu.H - trunc(nu.H/n) * n)
  nu.G <- round((2^(1:J) - 1) * coe(g), 0)
  nu.G <- ifelse(nu.G/n < 1, nu.G, nu.G - trunc(nu.G/n) * n)
  
  if(!inv) {
    ## Apply the phase shifts
    for(j in 0:(J-1)) {
      H.order <- c((nu.H[j+1]+1):n, 1:nu.H[j+1])
      G.order <- c((nu.G[j+1]+1):n, 1:nu.G[j+1])
      z[[3*j+1]] <- z[[3*j+1]][G.order, H.order]
      z[[3*j+2]] <- z[[3*j+2]][H.order, G.order]
      z[[3*j+3]] <- z[[3*j+3]][H.order, H.order]
    } 
    z[[3*J+1]] <- z[[3*J+1]][G.order, G.order]
  }
  else {
    ## Apply the phase shifts "reversed"
    for(j in 0:(J-1)) {
      H.order <- c((n-nu.H[j+1]+1):n, 1:(n-nu.H[j+1]))
      G.order <- c((n-nu.G[j+1]+1):n, 1:(n-nu.G[j+1]))
      z[[3*j+1]] <- z[[3*j+1]][G.order, H.order]
      z[[3*j+2]] <- z[[3*j+2]][H.order, G.order]
      z[[3*j+3]] <- z[[3*j+3]][H.order, H.order]
    }
    z[[3*J+1]] <- z[[3*J+1]][G.order, G.order]
  }
  return(z)
}

spp.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01, frac=1)
{
  sppLL <- function(x, y) {
    delta <- x[1]
    fG <- x[2]

    ## cat("Parameters are: d =", delta, ", and f =", fG, fill=TRUE)
    
    y.dwpt <- y[[1]]
    y.basis <- y[[2]]
    n <- y[[3]]
    J <- y[[4]]
    
    ## Establish the limits of integration for the band-pass variances
    a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    
    ## Define some useful parameters for the wavelet packet tree
                                        # n <- length(y)
    length.jn <- n / rep(2^(1:J), 2^(1:J))
    scale.jn <- rep(2^(1:J+1), 2^(1:J))
    
    ## Initialize various parameters for the reduced LL
    Basis <- (1:length(y.basis))[y.basis]
    bp.var <- numeric(length(Basis))
    delta.n <- 100

    ## Compute the band-pass variances according to \delta and f_G
    omega.diag <- NULL
    for(i in 1:sum(y.basis)) {
      jn <- Basis[i]
      bp.var[i] <- bandpass.spp(a[jn], b[jn], delta, fG)
      omega.diag <- c(omega.diag,
                      scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
    }
    
    ## Compute reduced log-likelihood 
    rLL <- n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
      sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var))
    rLL
  }
  
  n <- length(y)
  x0 <- numeric(2)

  ## Perform discrete wavelet packet transform (DWPT) on Y
  y.dwpt <- dwpt(y, wf, n.levels=J)
  n <- length(y)
  if(frac < 1) {
    for(i in 1:length(y.dwpt)) {
      vec <- y.dwpt[[i]]
      ni <- length(vec)
      j <- rep(1:J, 2^(1:J))[i]
      vec[trunc(frac * n/2^j):ni] <- NA
      y.dwpt[[i]] <- vec
    }
  }
  y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p)))
  y.dwpt <- as.matrix(unlist(y.dwpt[y.basis]))

  ## Compute initial estimate of the Gegenbauer frequency
  y.per <- per(y - mean(y))
  x0[2] <- (0:(n/2)/n)[max(y.per) == y.per]

  ## Compute initial estimate of the fractional difference parameter
  muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
          2^(rep(1:J, 2^(1:J))) + 
          unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /  
          2^(rep(1:J, 2^(1:J)))) / 4
  y.modwpt <- modwpt(y, wf=wf, n.levels=J)
  y.varJ <- rep(2^(1:J), 2^(1:J)) *
    unlist(lapply(y.modwpt,
                  FUN=function(x)sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
  x0[1] <- min(-0.5 * lsfit(log(abs(muJ[y.basis] - x0[2])), 
                            log(y.varJ[y.basis]))$coef[2], 0.49)

  cat(paste("Initial parameters are: delta =", round(x0[1],4), 
            "freqG =", round(x0[2],4), "\n"))
  result <- optim(par=x0, fn=sppLL, method="L-BFGS-B",
                  lower=c(0.001,0.001), upper=c(0.499,0.499),
                  control=list(trace=0, fnscale=2),
                  y=list(y.dwpt, y.basis, n, J))
  return(result)
}

spp2.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01,
                     dyadic=TRUE, frac=1)
{

  spp2LL <- function(x, y) {
    d1 <- x[1]
    f1 <- x[2]
    d2 <- x[3]
    f2 <- x[4]
    ## cat("Parameters are: d1 =", round(d1,6), ", and f1 =", round(f1,6),
    ##     ", d2 =", round(d2,6), ", and f2 =", round(f2,6), fill=TRUE)

    y.dwpt <- y[[1]]
    y.basis <- y[[2]]
    n <- y[[3]]
    J <- y[[4]]

    ## Establish the limits of integration for the band-pass variances
    a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) / 
      2^(rep(1:J, 2^(1:J))) / 2
    
    ## Define some useful parameters for the wavelet packet tree
    length.jn <- n / rep(2^(1:J), 2^(1:J))
    scale.jn <- rep(2^(1:J+1), 2^(1:J))
    
    ## Initialize various parameters for the reduced LL
    Basis <- (1:length(y.basis))[y.basis]
    bp.var <- numeric(length(Basis))
    delta.n <- 100

    ## Compute the band-pass variances according to \delta and f_G
    omega.diag <- NULL
    for(i in 1:sum(y.basis)) {
      jn <- Basis[i]
      bp.var[i] <- bandpass.spp2(a[jn], b[jn], d1, f1, d2, f2)
      omega.diag <- c(omega.diag,
                      scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
    }
    
    ## Compute reduced log-likelihood 
    n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
      sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var), na.rm=TRUE)
  }

  n <- length(y)
  x0 <- numeric(4)

  ## Perform discrete wavelet packet transform (DWPT) on Y
  y.dwpt <- dwpt(y, wf, n.levels=J)
  if(!dyadic) {
    for(i in 1:length(y.dwpt)) {
      vec <- y.dwpt[[i]]
      ni <- length(vec)
      j <- rep(1:J, 2^(1:J))[i]
      vec[trunc(frac * n/2^j):ni] <- NA
      y.dwpt[[i]] <- vec
    }
  }
  y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p, type="other")))
  y.dwpt <- as.vector(unlist(y.dwpt[y.basis]))

  ## Compute initial estimate of the Gegenbauer frequencies
  if(dyadic)
    y.per <- per(y - mean(y))
  else
    y.per <- per(y[1:(frac*n)] - mean(y[1:(frac*n)]))
  freq.y <- (0:(frac*n %/% 2))/(frac*n)
  x0[2] <- freq.y[max(y.per) == y.per]
  x0[4] <- freq.y[max(y.per[freq.y > x0[2] + freq.y[10] | 
    freq.y < x0[2] - freq.y[10]]) == y.per]
  if(x0[2] > x0[4]) {
    xx <- x0[2]
    x0[2] <- x0[4]
    x0[4] <- xx
    rm(xx)
  }

  ## Compute initial estimate of the fractional difference parameters
  muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) / 
      2^(rep(1:J, 2^(1:J))) + 
      unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /  
      2^(rep(1:J, 2^(1:J)))) / 4
  y.modwpt <- modwpt(y, wf=wf, n.levels=J)
  y.varJ <- rep(2^(1:J), 2^(1:J)) *
    unlist(lapply(y.modwpt,
                  FUN = function(x) sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
  x0.mid <- (x0[2] + x0[4]) / 2
  muJ <- muJ[y.basis]
  y.varJ <- y.varJ[y.basis]
  x0[1] <- min(-0.5 * lsfit(log(abs(muJ[muJ < x0.mid] - x0[2])), 
                            log(y.varJ[muJ < x0.mid]))$coef[2], 0.49)
  x0[3] <- min(-0.5 * lsfit(log(abs(muJ[muJ > x0.mid] - x0[4])), 
                            log(y.varJ[muJ > x0.mid]))$coef[2], 0.49)

  cat(paste("Initial parameters: d1 = ", round(x0[1],4), 
            ", f1 = ", round(x0[2],4), ", d2 = ", round(x0[3],4), 
            ", f2 = ", round(x0[4],4), sep=""), fill=TRUE)
  result <- optim(par=x0, fn=spp2LL, method="L-BFGS-B",
                  lower=rep(0.001,4), upper=rep(0.499,4),
                  control=list(trace=1, fnscale=2),
                  y=list(y.dwpt, y.basis, n, J))
  return(result)
}
stack.plot <-
  function (x, plot.type = c("multiple", "single"), panel = lines, 
    log = "", col = par("col"), bg = NA, pch = par("pch"), cex = par("cex"),
    lty = par("lty"), lwd = par("lwd"), ann = par("ann"), xlab = "Time",
    main = NULL, oma = c(6, 0, 5, 0), layout = NULL,
    same.scale = 1:dim(x)[2], ...) 
{
  addmain <- function(main, cex.main = par("cex.main"),
                      font.main = par("font.main"), 
                      col.main = par("col.main"), ...) {
    mtext(main, 3, 3, cex = cex.main, font = font.main, col = col.main, 
          ...)
  }
  plot.type <- match.arg(plot.type)
  panel <- match.fun(panel)
  nser <- NCOL(x)
  if (plot.type == "single" || nser == 1) {
    m <- match.call()
    m[[1]] <- as.name("plot.ts")
    m$plot.type <- "single"
    return(eval(m, parent.frame()))
  }
  if (nser > 10) 
    stop("Can't plot more than 10 series")
  if (is.null(main)) 
    main <- deparse(substitute(x))
  nm <- colnames(x)
  if (is.null(nm)) 
    nm <- paste("Series", 1:nser)
  nc <- if (nser > 4) 
    2
  else 1
  oldpar <- par("mar", "oma", "mfcol")
  on.exit(par(oldpar))
  par(mar = c(0, 5.1, 0, 2.1), oma = oma)
  nr <- ceiling(nser/nc)
  ## Begin added code
  if(!is.null(same.scale)) {
    unique.scales <- length(unique(same.scale))
    ylim <- vector("list", unique.scales)
    for (i in 1:unique.scales)
      ylim[[i]] <- range(x[, same.scale==i])
  }
  else
    for (i in 1:dim(x)[2])
      ylim[[i]] <- range(x[,i])
  if(is.null(layout))
    par(mfcol = c(nr, nc))
  else {
    par(mfcol = layout)
    nr <- layout[1]
  }
  ## End added code
  for (i in 1:nser) {
    plot(x[, i], axes = FALSE, xlab = "", ylab = "", log = log, 
         col = col, bg = bg, pch = pch, ann = ann, type = "n", 
         ylim=ylim[[same.scale[i]]], ...)
    panel(x[, i], col = col, bg = bg, pch = pch, ...)
    box()
    axis(2, xpd = NA)
    mtext(nm[i], 2, 3)
    if (i%%nr == 0 || i == nser) 
      axis(1, xpd = NA)
  }
  if (ann) {
    mtext(xlab, 1, 3, ...)
    if (!is.null(main)) {
      par(mfcol = c(1, 1))
      addmain(main, ...)
    }
  }
  invisible()
}
dpss.taper <- function(n, k, nw = 4, nmax = 2^(ceiling(log(n,2)))) {
  if(n > nmax)
    stop("length of taper is greater than nmax")
  w <- nw/n
  if(w > 0.5)
    stop("half-bandwidth parameter (w) is greater than 1/2")
  if(k <= 0)
    stop("positive dpss order (k) required")
  v <- matrix(0, nrow = nmax, ncol = (k + 1))
  storage.mode(v) <- "double"
  out <- .Fortran("dpss",
                  nmax = as.integer(nmax),
                  kmax = as.integer(k),
                  n = as.integer(n),
                  w = as.double(w),
                  v = v,
                  sig = double(k + 1),
                  totit = integer(1),
                  sines = double(n),
                  vold = double(n),
                  u = double(n),
                  scr1 = double(n),
                  ifault = integer(1),
                  PACKAGE="waveslim")
  ##list(v = out$v[1:n, 1:k], eigen = out$sig[-1] + 1, iter = 
  ##     out$totiTRUE, n = out$n, w = out$w, ifault = out$ifault)
  return(out$v[1:n, 1:k])
}

sine.taper <- function(n, k) {
  tapers <- NULL
  for(i in 1:k)
    tapers <- cbind(tapers, sqrt(2/(n+1)) * sin((pi*i*1:n)/(n+1)))
  return(tapers)
}
###########################################################################
###########################################################################
###########################################################################

dwt.3d <- function(x, wf, J=4, boundary="periodic")
{
  nx <- dim(x)[1]
  storage.mode(nx) <- "integer"
  ny <- dim(x)[2]
  storage.mode(ny) <- "integer"
  nz <- dim(x)[3]
  storage.mode(nz) <- "integer"

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  z <- array(0, dim=c(nx,ny,nz)/2)
  storage.mode(z) <- "double"

  x.wt <- vector("list", 7*J+1)
  x.names <- NULL
  for(j in 1:J) {
    out <- .C("three_D_dwt", "cube"=as.double(x), "NX"=nx, "NY"=ny,
              "NZ"=nz, "filter.length"=L, "hpf"=h, "lpf"=g, "LLL"=z,
              "HLL"=z, "LHL"=z, "LLH"=z, "HHL"=z, "HLH"=z, "LHH"=z,
              "HHH"=z, PACKAGE="waveslim")[8:15]
    if(j < J) {
      index <- (7*(j-1)+1):(7*j)
      x.wt[index] <- out[-1]
      x.names <- c(x.names, sapply(names(out)[-1], paste, j, sep=""))
      x <- out[[1]]
      nx <- dim(x)[1]
      storage.mode(nx) <- "integer"
      ny <- dim(x)[2]
      storage.mode(ny) <- "integer"
      nz <- dim(x)[3]
      storage.mode(nz) <- "integer"
      z <- array(0, dim=c(nx,ny,nz)/2)
      storage.mode(z) <- "double"
    }
    else {
      index <- (7*(j-1)+1):(7*j+1)
      x.wt[index] <- out[c(2:8,1)]
      x.names <- c(x.names, sapply(names(out)[c(2:8,1)], paste, j, sep=""))
    }
  }

  names(x.wt) <- x.names
  class(x.wt) <- "dwt.3d"
  attr(x.wt, "J") <- J
  attr(x.wt, "wavelet") <- wf
  attr(x.wt, "boundary") <- boundary
  return(x.wt)
}

###########################################################################
###########################################################################
###########################################################################

idwt.3d <- function(y)
{
  J <- attributes(y)$J
  LLL <- paste("LLL", J, sep="")

  wf <- attributes(y)$wavelet
  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  y.in <- y$LLL

  for(j in J:1) {
    HLL <- paste("HLL", j, sep="")
    LHL <- paste("LHL", j, sep="")
    LLH <- paste("LLH", j, sep="")
    HHL <- paste("HHL", j, sep="")
    HLH <- paste("HLH", j, sep="")
    LHH <- paste("LHH", j, sep="")
    HHH <- paste("HHH", j, sep="")
    
    nx <- dim(y.in)[1]
    storage.mode(nx) <- "integer"
    ny <- dim(y.in)[2]
    storage.mode(ny) <- "integer"
    nz <- dim(y.in)[3]
    storage.mode(nz) <- "integer"

    z <- array(0, dim=2*c(nx, ny, nz))
    storage.mode(z) <- "double"

    out <- .C("three_D_idwt", as.double(y.in), as.double(y[[HLL]]),
              as.double(y[[LHL]]), as.double(y[[LLH]]),
              as.double(y[[HHL]]), as.double(y[[HLH]]),
              as.double(y[[LHH]]), as.double(y[[HHH]]), 
              nx, ny, nz, L, h, g, "Y"=z, PACKAGE="waveslim")

    y.in <- out$Y
  }
  zapsmall(y.in)
}

###########################################################################
###########################################################################
###########################################################################

modwt.3d <- function(x, wf, J=4, boundary="periodic")
{
  nx <- dim(x)[1]
  storage.mode(nx) <- "integer"
  ny <- dim(x)[2]
  storage.mode(ny) <- "integer"
  nz <- dim(x)[3]
  storage.mode(nz) <- "integer"

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf / sqrt(2)
  storage.mode(h) <- "double"
  g <- dict$lpf / sqrt(2)
  storage.mode(g) <- "double"

  z <- array(0, dim=c(nx,ny,nz))
  storage.mode(z) <- "double"

  x.wt <- vector("list", 7*J+1)
  x.names <- NULL
  for(j in 1:J) {
    out <- .C("three_D_modwt", "cube"=as.double(x), "NX"=nx, "NY"=ny,
              "NZ"=nz, "J"=j, "filter.length"=L, "hpf"=h, "lpf"=g,
              "LLL"=z, "HLL"=z, "LHL"=z, "LLH"=z, "HHL"=z, "HLH"=z,
              "LHH"=z, "HHH"=z, PACKAGE="waveslim")[9:16]
    if(j < J) {
      index <- (7*(j-1)+1):(7*j)
      x.wt[index] <- out[-1]
      x.names <- c(x.names, sapply(names(out)[-1], paste, j, sep=""))
      x <- out[[1]]
      nx <- dim(x)[1]
      storage.mode(nx) <- "integer"
      ny <- dim(x)[2]
      storage.mode(ny) <- "integer"
      nz <- dim(x)[3]
      storage.mode(nz) <- "integer"
      z <- array(0, dim=c(nx,ny,nz))
      storage.mode(z) <- "double"
    }
    else {
      index <- (7*(j-1)+1):(7*j+1)
      x.wt[index] <- out[c(2:8,1)]
      x.names <- c(x.names, sapply(names(out)[c(2:8,1)], paste, j, sep=""))
    }
  }

  names(x.wt) <- x.names
  class(x.wt) <- "modwt.3d"
  attr(x.wt, "J") <- J
  attr(x.wt, "wavelet") <- wf
  attr(x.wt, "boundary") <- boundary
  return(x.wt)
}

###########################################################################
###########################################################################
###########################################################################

imodwt.3d <- function(y)
{
  J <- attributes(y)$J
  LLL <- paste("LLL", J, sep="")

  wf <- attributes(y)$wavelet
  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf / sqrt(2)
  storage.mode(h) <- "double"
  g <- dict$lpf / sqrt(2)
  storage.mode(g) <- "double"

  y.in <- y$LLL

  for(j in J:1) {
    HLL <- paste("HLL", j, sep="")
    LHL <- paste("LHL", j, sep="")
    LLH <- paste("LLH", j, sep="")
    HHL <- paste("HHL", j, sep="")
    HLH <- paste("HLH", j, sep="")
    LHH <- paste("LHH", j, sep="")
    HHH <- paste("HHH", j, sep="")
    
    nx <- dim(y.in)[1]
    storage.mode(nx) <- "integer"
    ny <- dim(y.in)[2]
    storage.mode(ny) <- "integer"
    nz <- dim(y.in)[3]
    storage.mode(nz) <- "integer"

    z <- array(0, dim=c(nx, ny, nz))
    storage.mode(z) <- "double"

    out <- .C("three_D_imodwt", as.double(y.in), as.double(y[[HLL]]),
              as.double(y[[LHL]]), as.double(y[[LLH]]),
              as.double(y[[HHL]]), as.double(y[[HLH]]),
              as.double(y[[LHH]]), as.double(y[[HHH]]), 
              nx, ny, nz, j, L, h, g, "Y"=z, PACKAGE="waveslim")

    y.in <- out$Y
  }
  zapsmall(y.in)
}

###########################################################################
###########################################################################
###########################################################################

mra.3d <-
  function(x, wf="la8", J=4, method="modwt", boundary="periodic")
{
  nx <- dim(x)[1]
  ny <- dim(x)[2]
  nz <- dim(x)[3]

  if(method == "modwt") {
    x.wt <- modwt.3d(x, wf, J, "periodic")
  } else {
    x.wt <- dwt.3d(x, wf, J, "periodic")
  }
  
  x.mra <- vector("list", 7*J+1)
  names(x.mra) <-
    c(matrix(rbind(paste("HLL", 1:J, sep=""), paste("LHL", 1:J, sep=""),
                   paste("LLH", 1:J, sep=""), paste("HHL", 1:J, sep=""),
                   paste("HLH", 1:J, sep=""), paste("LHH", 1:J, sep=""),
                   paste("HHH", 1:J, sep="")), nrow=1),
      paste("LLL", J, sep=""))

  ## Smooth
  zero <- vector("list", 7*J+1)
  names(zero) <- names(x.mra)
  attr(zero, "J") <- J
  attr(zero, "wavelet") <- wf
  attr(zero, "boundary") <- boundary
  zero[[7*J+1]] <- x.wt[[7*J+1]]
  if(method == "modwt") {
    class(x.wt) <- "modwt.3d"
    for(k in 1:(7*J))
      zero[[k]] <- array(0, dim=c(nx,ny,nz))
    x.mra[[7*J+1]] <- imodwt.3d(zero)
  } else {
    class(x.wt) <- "dwt.3d"
    for(k in 1:J)
      zero[[7*(k-1)+1]] <- zero[[7*(k-1)+2]] <- zero[[7*(k-1)+3]] <-
        zero[[7*(k-1)+4]] <- zero[[7*(k-1)+5]] <- zero[[7*(k-1)+6]] <-
          zero[[7*k]] <- array(0, dim=c(nx,ny,nz)/2^k)
    x.mra[[7*J+1]] <- idwt.3d(zero)
  }

  ## Details
  for(j in (7*J):1) {
    Jj <- ceiling(j/7)
    zero <- vector("list", 7*Jj+1)
    names(zero) <-
      c(matrix(rbind(paste("HLL", 1:Jj, sep=""), paste("LHL", 1:Jj, sep=""),
                     paste("LLH", 1:Jj, sep=""), paste("HHL", 1:Jj, sep=""),
                     paste("HLH", 1:Jj, sep=""), paste("LHH", 1:Jj, sep=""),
                     paste("HHH", 1:Jj, sep="")), nrow=1),
        paste("LLL", Jj, sep=""))
    attr(zero, "J") <- Jj
    attr(zero, "wavelet") <- wf
    attr(zero, "boundary") <- boundary
    zero[[j]] <- x.wt[[j]]
    if(method == "modwt") {
      for(k in names(zero)[-charmatch(names(zero)[j], names(zero))])
        zero[[k]] <- array(0, dim=c(nx,ny,nz))
      x.mra[[j]] <- imodwt.3d(zero)
    } else {
      for(k in 1:Jj)
        zero[[7*(k-1)+1]] <- zero[[7*(k-1)+2]] <- zero[[7*(k-1)+3]] <-
          zero[[7*(k-1)+4]] <- zero[[7*(k-1)+5]] <- zero[[7*(k-1)+6]] <-
            zero[[7*k]] <- array(0, dim=c(nx,ny,nz)/2^k)
      zero[[7*Jj+1]] <- array(0, dim=c(nx,ny,nz)/2^Jj)
      zero[[j]] <- x.wt[[j]]
      x.mra[[j]] <- idwt.3d(zero)
    }
  }
  return(x.mra)
}
###########################################################################
###########################################################################
###########################################################################

dwt.2d <- function(x, wf, J=4, boundary="periodic")
{
  m <- dim(x)[1]
  storage.mode(m) <- "integer"
  n <- dim(x)[2]
  storage.mode(n) <- "integer"

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  z <- matrix(0, m/2, n/2)
  storage.mode(z) <- "double"

  x.wt <- vector("list", 3*J+1)
  x.names <- NULL
  for(j in 1:J) {
    out <- .C("two_D_dwt", "Image"=as.double(x), "Rows"=m, "Cols"=n, 
                "filter.length"=L, "hpf"=h, "lpf"=g, "LL"=z, "LH"=z,
                "HL"=z, "HH"=z, PACKAGE="waveslim")[7:10]
    if(j < J) {
      index <- (3*j-2):(3*j)
      x.wt[index] <- out[-1]
      x.names <- c(x.names, sapply(names(out)[-1], paste, j, sep=""))
      x <- out[[1]]
      m <- dim(x)[1]
      storage.mode(m) <- "integer"
      n <- dim(x)[2]
      storage.mode(n) <- "integer"
      z <- matrix(0, m/2, n/2)
      storage.mode(z) <- "double"
    }
    else {
      index <- (3*j):(3*(j+1)) - 2
      x.wt[index] <- out[c(2:4,1)]
      x.names <- c(x.names, sapply(names(out)[c(2:4,1)], paste, j, sep=""))
    }
  }

  names(x.wt) <- x.names
  attr(x.wt, "J") <- J
  attr(x.wt, "wavelet") <- wf
  attr(x.wt, "boundary") <- boundary
  x.wt
}

###########################################################################
###########################################################################
###########################################################################

idwt.2d <- function(y)
{
  J <- attributes(y)$J

  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  LL <- paste("LL", J, sep="")
  y.in <- y[[LL]]

  for(j in J:1) {
    LH <- paste("LH", j, sep="")
    HL <- paste("HL", j, sep="")
    HH <- paste("HH", j, sep="")
    
    m <- dim(y.in)[1]
    storage.mode(m) <- "integer"
    n <- dim(y.in)[2]
    storage.mode(n) <- "integer"
    x <- matrix(0, 2*m, 2*n)
    storage.mode(x) <- "double"

    out <- .C("two_D_idwt", as.double(y.in), as.double(y[[LH]]),
              as.double(y[[HL]]), as.double(y[[HH]]), m, n, L, h, g,
              "Y"=x, PACKAGE="waveslim")
    y.in <- out$Y
  }
  zapsmall(y.in)
}

###########################################################################
###########################################################################
###########################################################################

modwt.2d <- function(x, wf, J=4, boundary="periodic")
{
  m <- dim(x)[1]
  storage.mode(m) <- "integer"
  n <- dim(x)[2]
  storage.mode(n) <- "integer"

  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf / sqrt(2)
  storage.mode(h) <- "double"
  g <- dict$lpf / sqrt(2)
  storage.mode(g) <- "double"

  z <- matrix(0, m, n)
  storage.mode(z) <- "double"

  x.wt <- vector("list", 3*J+1)
  x.names <- NULL
  for(j in 1:J) {
    out <- .C("two_D_modwt", "Image"=as.double(x), "Rows"=m, "Cols"=n,
              "Level"=j, "filter.length"=L, "hpf"=h, "lpf"=g, "LL"=z,
              "LH"=z, "HL"=z, "HH"=z, PACKAGE="waveslim")[8:11]
    if(j < J) {
      index <- (3*j-2):(3*j)
      x.wt[index] <- out[-1]
      x.names <- c(x.names, sapply(names(out)[-1], paste, j, sep=""))
      x <- out$LL
    }
    else {
      index <- (3*j):(3*(j+1)) - 2
      x.wt[index] <- out[c(2:4,1)]
      x.names <- c(x.names, sapply(names(out)[c(2:4,1)], paste, j, sep=""))
    }
  }

  names(x.wt) <- x.names
  attr(x.wt, "J") <- J
  attr(x.wt, "wavelet") <- wf
  attr(x.wt, "boundary") <- boundary
  x.wt
}

###########################################################################
###########################################################################
###########################################################################

imodwt.2d <- function(y)
{
  J <- attributes(y)$J

  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf / sqrt(2)
  storage.mode(h) <- "double"
  g <- dict$lpf / sqrt(2)
  storage.mode(g) <- "double"

  LL <- paste("LL", J, sep="")
  y.in <- y[[LL]]

  for(j in J:1) {
    LH <- paste("LH", j, sep="")
    HL <- paste("HL", j, sep="")
    HH <- paste("HH", j, sep="")
    
    m <- dim(y.in)[1]
    storage.mode(m) <- "integer"
    n <- dim(y.in)[2]
    storage.mode(n) <- "integer"
    x <- matrix(0, m, n)
    storage.mode(x) <- "double"

    out <- .C("two_D_imodwt", as.double(y.in), as.double(y[[LH]]),
              as.double(y[[HL]]), as.double(y[[HH]]), m, n, j, L,
              h, g, "Y"=x, PACKAGE="waveslim")
    y.in <- out$Y
  }
  zapsmall(y.in)
}

###########################################################################
###########################################################################
###########################################################################

plot.dwt.2d <- function(x, cex.axis=1, plot=TRUE, ...)
{
  J <- attributes(x)$J
  X <- x[[paste("LL", J, sep="")]]
  for(j in J:1) {
    x.names <- sapply(c("LH","HL","HH"), paste, j, sep="")
    X <- rbind(cbind(X, x[[x.names[2]]]),
               cbind(x[[x.names[1]]], x[[x.names[3]]]))
  }
  M <- dim(X)[1]; N <- dim(X)[2]
  if(plot) {
    image(1:M, 1:N, X, col=rainbow(128), axes=FALSE, xlab="", ylab="", ...)
    x.label <- NULL
    lines(c(0,N,N,0,0) + 0.5, c(0,0,M,M,0) + 0.5)
    for(j in J:1) {
      lines(c(M/2^j,M/2^j) + 0.5, 2*c(0,N/2^j) + 0.5)
      lines(2*c(0,M/2^j) + 0.5, c(N/2^j,N/2^j) + 0.5)
    }
    at <- c((3*N+2)/2^(1:J+1),(N+2)/2^(J+1))
    labs <- c(paste("H",1:J,sep=""), paste("L",J,sep=""))
    axis(side=1, at=at, label=labs, tick=FALSE, cex.axis=cex.axis)
    axis(side=2, at=at, label=labs, tick=FALSE, cex.axis=cex.axis)
  }
  else
    return(X)
  invisible()
}

###########################################################################
###########################################################################
###########################################################################

denoise.dwt.2d <- function(x, wf = "la8", J = 4, method = "universal", 
                           H = 0.5, noise.dir = 3, rule = "hard")
{
  soft <- function(x, delta) sign(x) * pmax(abs(x) - delta, 0)
  hard <- function(x, delta) ifelse(abs(x) > delta, x, 0)

  n <- length(x)
  x.dwt <- dwt.2d(x, wf, J)
  if(noise.dir == 3)
    sigma.mad <- list(HH = mad(x.dwt$HH1), HL = mad(x.dwt$HL1), 
  		      LH = mad(x.dwt$LH1))
  else {
    noise <- x.dwt$jj
    sigma.mad <- list(HH = mad(noise), HL = mad(noise), LH = mad(noise))
  }
    thresh <- list(HH = rep(sqrt(2 * sigma.mad$HH^2 * log(n)), J), 
		   HL = rep(sqrt(2 * sigma.mad$HL^2 * log(n)), J),
		   LH = rep(sqrt(2 * sigma.mad$LH^2 * log(n)), J))

  if(method == "long-memory")
    thresh <- lapply(thresh, function(x,J,H) 2^(0:(J-1)*(H-1/2))*x, J=J, H=H)
  for(j in 1:J) {
    jj <- paste("HL", j, sep = "")
    if(rule == "hard")
      x.dwt[[jj]] <- hard(x.dwt[[jj]], thresh$HL[j])
    else 
      x.dwt[[jj]] <- soft(x.dwt[[jj]], thresh$HL[j])
    jj <- paste("LH", j, sep = "")
    if(rule == "hard")
      x.dwt[[jj]] <- hard(x.dwt[[jj]], thresh$LH[j])
    else 
      x.dwt[[jj]] <- soft(x.dwt[[jj]], thresh$LH[j])
    jj <- paste("HH", j, sep = "")
    if(rule == "hard")
      x.dwt[[jj]] <- hard(x.dwt[[jj]], thresh$HH[j])
    else 
      x.dwt[[jj]] <- soft(x.dwt[[jj]], thresh$HH[j])
  }
  idwt.2d(x.dwt)
}

###########################################################################
###########################################################################
###########################################################################

denoise.modwt.2d <- function(x, wf = "la8", J = 4, method = "universal", 
  H = 0.5, rule = "hard")
{
  soft <- function(x, delta) sign(x) * pmax(abs(x) - delta, 0)
  hard <- function(x, delta) ifelse(abs(x) > delta, x, 0)
  n <- length(x)
  x.modwt <- modwt.2d(x, wf, J)
  sigma.mad <- list(HH = sqrt(2) * mad(x.modwt$HH1),
		    HL = sqrt(2) * mad(x.modwt$HL1),
		    LH = sqrt(2) * mad(x.modwt$LH1))
    thresh <- list(HH = rep(sqrt(2 * sigma.mad$HH^2 * log(n))/2^(1:J), J), 
		   HL = rep(sqrt(2 * sigma.mad$HL^2 * log(n))/2^(1:J), J), 
		   LH = rep(sqrt(2 * sigma.mad$LH^2 * log(n))/2^(1:J), J))
  if(method == "long-memory")
    thresh <- lapply(thresh, function(x,J,H) 2^(0:(J-1)*(H-1/2))*x, J=J, H=H)
  for(j in 1:J) {
    jj <- paste("HL", j, sep = "")
    if(rule == "hard")
      x.modwt[[jj]] <- hard(x.modwt[[jj]], thresh$HL[j])
    else 
      x.modwt[[jj]] <- soft(x.modwt[[jj]], thresh$HL[j])
    jj <- paste("LH", j, sep = "")
    if(rule == "hard")
      x.modwt[[jj]] <- hard(x.modwt[[jj]], thresh$LH[j])
    else 
      x.modwt[[jj]] <- soft(x.modwt[[jj]], thresh$LH[j])
    jj <- paste("HH", j, sep = "")
    if(rule == "hard")
      x.modwt[[jj]] <- hard(x.modwt[[jj]], thresh$HH[j])
    else 
     x.modwt[[jj]] <- soft(x.modwt[[jj]], thresh$HH[j])
  }
  imodwt.2d(x.modwt)
}

###########################################################################
###########################################################################
###########################################################################

dwpt.2d <- function(x, wf="la8", J=4, boundary="periodic")
{
  ## x <- xbox
  ## Define image dimensions (assign mode for C) and perform simple
  ## diagnostics.
  m <- dim(x)[1]
  storage.mode(m) <- "integer"
  n <- dim(x)[2]
  storage.mode(n) <- "integer"
  if(log(m, 2) != trunc(log(m, 2)) | log(n, 2) != trunc(log(n, 2)))
    stop("One dimension is not a power of 2")
  if(2^J > m | 2^J > n)
    stop("Wavelet transform exceeds sample size in one dimension of DWPT")

  ## Extract wavelet and scaling filter coefficients, along with filter
  ## length, from the filter name provided.
  dict <- wave.filter(wf)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"

  ## Create names for wavelet packet nodes (quad-tree structure).
  N <- sum(4^(1:J))
  level <- rep(1:J, 4^(1:J))
  x.wpt <- vector("list", N)
  c1 <- rep(1:J, 2^(1:J))
  c2 <- unlist(apply(as.matrix(2^(1:J) - 1), 1, seq, from=0))
  cry <- paste("w", c1, ".", c2, sep="")
  x.wpt.names <- NULL
  for(j in 1:J) {
    xx <- matrix(cry[c1 == j], 2^j, 2^j)
    yy <- matrix(cry[c1 == j], 2^j, 2^j, byrow=TRUE)
    x.wpt.names <- c(x.wpt.names, as.matrix(paste(xx, "-", yy, sep="")))
  }
  names(x.wpt) <- x.wpt.names
  rm(j,xx,yy,c1,c2,cry)
  ## Define initial zero matrix to store wavelet sub-images.
  z <- matrix(0, m/2, n/2)
  storage.mode(z) <- "double"

  ## Implement the 2D DWPT in a nested loop structure.
  for(j in 1:J) {
    ## cat("j =", j, fill=TRUE)
    for(k in 0:(4^j/4-1)) {
      if(j > 1) {
        ## if j > 1, grab wavelet coefficient image and also its name.
        index <- min((1:N)[level == j-1]) + k
        parent <- x.wpt.names[index]
        ## cat("parent =", parent, fill=TRUE)
        x <- x.wpt[[parent]]
        tmp <- unlist(strsplit(parent, "\\-"))
      }
      else
        tmp <- c("w0.0", "w0.0")
      ## Deconstruct name into nodes for the x and y dimensions.
      node <- unlist(strsplit(tmp, "\\."))
      node <- as.integer(node[-c(1,3)])
      ## Preliminary assignments in order to keep wavelet coefficient
      ## sub-images in sequency order.
      if(node[1] %% 2 == 0) {
        Xlow <- paste("w", j, ".", 2 * node[1], sep="")
        Xhigh <- paste("w", j, ".", 2 * node[1] + 1, sep="")
      }
      else {
        Xlow <- paste("w", j, ".", 2 * node[1] + 1, sep="")
        Xhigh <- paste("w", j, ".", 2 * node[1], sep="")
      }
      if(node[2] %% 2 == 0) {
        Ylow <- paste("w", j, ".", 2 * node[2], sep="")
        Yhigh <- paste("w", j, ".", 2 * node[2] + 1, sep="")
      }
      else {
        Ylow <- paste("w", j, ".", 2 * node[2] + 1, sep="")
        Yhigh <- paste("w", j, ".", 2 * node[2], sep="")
      }
      ## Create names for the new wavelet coefficient images.
      LL <- paste(Xlow, "-", Ylow, sep="")
      LH <- paste(Xlow, "-", Yhigh, sep="")
      HL <- paste(Xhigh, "-", Ylow, sep="")
      HH <- paste(Xhigh, "-", Yhigh, sep="")
      ## cat(matrix(c(LH,LL,HH,HL), 2, 2), fill=TRUE)
      ## Perform the DWPT
      out <- .C("two_D_dwt", "Image"=as.double(x), "Rows"=m, "Cols"=n, 
                "filter.length"=L, "hpf"=h, "lpf"=g, "LL"=z, "LH"=z,
                "HL"=z, "HH"=z, PACKAGE="waveslim")[7:10]
      ## Pass wavelet coefficient images into the DWPT object.
      x.wpt[[LL]] <- out[["LL"]]
      x.wpt[[LH]] <- out[["LH"]]
      x.wpt[[HL]] <- out[["HL"]]
      x.wpt[[HH]] <- out[["HH"]]
    }
    ## Redefine zero matrix to its new (decimated) size.
    m <- dim(out[["LL"]])[1]
    storage.mode(m) <- "integer"
    n <- dim(out[["LL"]])[2]
    storage.mode(n) <- "integer"
    z <- matrix(0, m/2, n/2)
    storage.mode(z) <- "double"
  }
  attr(x.wpt, "J") <- J
  attr(x.wpt, "wavelet") <- wf
  attr(x.wpt, "boundary") <- boundary
  return(x.wpt)
}

###########################################################################
###########################################################################
###########################################################################

idwpt.2d <- function(y, y.basis)
{
  ## Error checking
  if(length(y) != length(y.basis))
    stop("DWPT object and basis selection must be the same length")
  ## Number of wavelet scales
  J <- attributes(y)$J
  ## Define wavelet/scaling filter coefficients and length
  dict <- wave.filter(attributes(y)$wavelet)
  L <- dict$length
  storage.mode(L) <- "integer"
  h <- dict$hpf
  storage.mode(h) <- "double"
  g <- dict$lpf
  storage.mode(g) <- "double"
  ## Nested for loops
  names(y.basis) <- names(y)
  for(j in J:1) {
    for(nx in seq(0, 2^j - 1, by = 2)) {
      for(ny in seq(0, 2^j - 1, by = 2)) {
        ## Name the four wavelet coefficients sub-images
        LL <- paste("w", j, ".", nx, "-", "w", j, ".", ny, sep="")
        LH <- paste("w", j, ".", nx, "-", "w", j, ".", ny+1, sep="")
        HL <- paste("w", j, ".", nx+1, "-", "w", j, ".", ny, sep="")
        HH <- paste("w", j, ".", nx+1, "-", "w", j, ".", ny+1, sep="")
        if(any(y.basis[LL], y.basis[LH], y.basis[HL], y.basis[HH])) {
          m <- nrow(y[[LL]])
          storage.mode(m) <- "integer"
          n <- ncol(y[[LL]])
          storage.mode(n) <- "integer"
          XX <- matrix(0, 2*m, 2*n)
          storage.mode(XX) <- "double"
          ## parent indices to construct string
          pnx <- floor(nx / 2)
          pny <- floor(ny / 2)
          if((pnx %% 2 != 0) & (pny %% 2 != 0))
            ## Upper right-hand corner
            out <- .C("two_D_idwt", as.double(y[[HH]]),
                      as.double(y[[HL]]), as.double(y[[LH]]),
                      as.double(y[[LL]]), m, n, L, h, g, "Y"=XX,
                      PACKAGE="waveslim")$Y
          else {
            ## Upper left-hand corner
            if((pnx %% 2 == 0) & (pny %% 2 != 0))
              out <- .C("two_D_idwt", as.double(y[[LH]]),
                        as.double(y[[LL]]), as.double(y[[HH]]),
                        as.double(y[[HL]]), m, n, L, h, g, "Y"=XX,
                        PACKAGE="waveslim")$Y
            else {
              ## Lower right-hand corner
              if((pnx %% 2 != 0) & (pny %% 2 == 0))
                out <- .C("two_D_idwt", as.double(y[[HL]]),
                          as.double(y[[HH]]), as.double(y[[LL]]),
                          as.double(y[[LH]]), m, n, L, h, g, "Y"=XX,
                          PACKAGE="waveslim")$Y
              else {
                ## Lower left-hand corner
                if((pnx %% 2 == 0) & (pny %% 2 == 0))
                  out <- .C("two_D_idwt", as.double(y[[LL]]),
                            as.double(y[[LH]]), as.double(y[[HL]]),
                            as.double(y[[HH]]), m, n, L, h, g, "Y"=XX,
                            PACKAGE="waveslim")$Y
                else
                  stop("Ouch!")
              }
            }
          }
          if(j > 1) {
            pname <- paste("w", j-1, ".", pnx, "-", "w", j-1, ".", pny, sep="")
            y[[pname]] <- out
            y.basis[pname] <- 1
          }
        }
      }
    }
  }
  return(out)
}

up.sample <- function(x, f, y=NA) {
  n <- length(x)
  as.vector(rbind(x, matrix(rep(y, (f-1)*n), nrow=f-1)))
}
wave.filter <- function(name)
{
  select.haar <- function() {
    L <- 2
    g <- c(0.7071067811865475, 0.7071067811865475)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.d4 <- function() {
    L <- 4
    g <- c(0.4829629131445341, 0.8365163037378077, 0.2241438680420134, 
	  -0.1294095225512603)
    h <- qmf(g)

    return(list(length = L, hpf = h, lpf = g))
  }
  select.mb4 <- function() {
    L <- 4
    g <- c(4.801755e-01, 8.372545e-01, 2.269312e-01, -1.301477e-01)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.bs3.1 <- function() {
    L <- 4
    g <- c(0.1767767, 0.5303301, 0.5303301, 0.1767767)
    h <- qmf(g)
    gd <- c(0.3535534, 1.06066, -1.06066, -0.3535534)
    hd <- qmf(g)
    return(list(length = L, hpf = h, lpf = g, dhpf = hd, dlpf = gd))
  }
  select.w4 <- function() {
    L <- 4
    g <- c(-1, 3, 3, -1) / 8
    h <- c(-1, 3, -3, 1) / 8
    return(list(length = L, hpf = h, lpf = g))
  }
  select.fk4 <- function() {
    L <- 4
    g <- c(.6539275555697651, .7532724928394872, .5317922877905981e-1,
          -.4616571481521770e-1)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.d6 <- function() {
    L <- 6
    g <- c(0.3326705529500827, 0.8068915093110928, 0.4598775021184915,
          -0.1350110200102546, -0.0854412738820267, 0.0352262918857096)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.fk6 <- function() {
    L <- 6
    g <- c(.4279150324223103, .8129196431369074, .3563695110701871,
          -.1464386812725773, -.7717775740697006e-1, .4062581442323794e-1)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.d8 <- function() {
    L <- 8
    g <- c(0.2303778133074431, 0.7148465705484058, 0.6308807679358788,
          -0.0279837694166834, -0.1870348117179132, 0.0308413818353661,
          0.0328830116666778, -0.0105974017850021)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.fk8 <- function() {
    L <- 8
    g <- c(.3492381118637999, .7826836203840648, .4752651350794712,
          -.9968332845057319e-1, -.1599780974340301, .4310666810651625e-1,
          .4258163167758178e-1, -.1900017885373592e-1)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.la8 <- function() {
    L <- 8
    g <- c(-0.07576571478935668, -0.02963552764596039, 0.49761866763256290, 
	  0.80373875180538600, 0.29785779560560505, -0.09921954357695636, 
	  -0.01260396726226383, 0.03222310060407815)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.mb8 <- function() {
    L <- 8
    g <- rev(c(-1.673619e-01, 1.847751e-02, 5.725771e-01, 7.351331e-01,
	   2.947855e-01, -1.108673e-01, 7.106015e-03, 6.436345e-02))
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.bl14 <- function() {
    L <- 14
    g <- c( 0.0120154192834842, 0.0172133762994439, -0.0649080035533744,
	  -0.0641312898189170, 0.3602184608985549, 0.7819215932965554,
	   0.4836109156937821, -0.0568044768822707, -0.1010109208664125,
	   0.0447423494687405, 0.0204642075778225, -0.0181266051311065,
	  -0.0032832978473081, 0.0022918339541009)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.fk14 <- function() {
    L <- 14
    g <- c(.2603717692913964, .6868914772395985, .6115546539595115,
          .5142165414211914e-1, -.2456139281621916, -.4857533908585527e-1,
          .1242825609215128, .2222673962246313e-1, -.6399737303914167e-1,
          -.5074372549972850e-2, .2977971159037902e-1, -.3297479152708717e-2,
          -.9270613374448239e-2, .3514100970435962e-2)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.d16 <- function() {
    L <- 16
    g <- c(0.0544158422431049, 0.3128715909143031, 0.6756307362972904,
	  0.5853546836541907, -0.0158291052563816, -0.2840155429615702,
	  0.0004724845739124, 0.1287474266204837, -0.0173693010018083,
	 -0.0440882539307952, 0.0139810279173995, 0.0087460940474061,
	 -0.0048703529934518, -0.0003917403733770, 0.0006754494064506,
	 -0.0001174767841248)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.la16 <- function() {
    L <- 16
    g <- c(-0.0033824159513594, -0.0005421323316355, 0.0316950878103452, 
	   0.0076074873252848, -0.1432942383510542, -0.0612733590679088, 
	   0.4813596512592012, 0.7771857516997478, 0.3644418948359564, 
	  -0.0519458381078751, -0.0272190299168137, 0.0491371796734768, 
	   0.0038087520140601, -0.0149522583367926, -0.0003029205145516, 
	   0.0018899503329007)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.mb16 <- function() {
    L <- 16
    g <- rev(c(-1.302770e-02, 2.173677e-02, 1.136116e-01, -5.776570e-02, 
	  -2.278359e-01, 1.188725e-01, 6.349228e-01, 6.701646e-01, 
	   2.345342e-01, -5.656657e-02, -1.987986e-02, 5.474628e-02, 
	  -2.483876e-02, -4.984698e-02, 9.620427e-03, 5.765899e-03))
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.la20 <- function() {
    L <- 20
    g <- c(0.0007701598091030, 0.0000956326707837, -0.0086412992759401,
	 -0.0014653825833465, 0.0459272392237649, 0.0116098939129724,
	 -0.1594942788575307, -0.0708805358108615, 0.4716906668426588,
	  0.7695100370143388, 0.3838267612253823, -0.0355367403054689,
	 -0.0319900568281631, 0.0499949720791560, 0.0057649120455518,
	 -0.0203549398039460, -0.0008043589345370, 0.0045931735836703,
	  0.0000570360843390, -0.0004593294205481)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.bl20 <- function() {
    L <- 20
    g <- c(0.0008625782242896, 0.0007154205305517, -0.0070567640909701,
	  0.0005956827305406, 0.0496861265075979, 0.0262403647054251,
	 -0.1215521061578162, -0.0150192395413644, 0.5137098728334054,
	  0.7669548365010849, 0.3402160135110789, -0.0878787107378667,
	 -0.0670899071680668, 0.0338423550064691, -0.0008687519578684,
	 -0.0230054612862905, -0.0011404297773324, 0.0050716491945793,
	  0.0003401492622332, -0.0004101159165852)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.fk22 <- function() {
    L <- 22
    g <- c(.1938961077599566, .5894521909294277, .6700849629420265,
          .2156298491347700, -.2280288557715772, -.1644657152688429,
          .1115491437220700, .1101552649340661, -.6608451679377920e-1,
          -.7184168192312605e-1, .4354236762555708e-1, .4477521218440976e-1,
          -.2974288074927414e-1, -.2597087308902119e-1, .2028448606667798e-1,
          .1296424941108978e-1, -.1288599056244363e-1, -.4838432636440189e-2,
          .7173803165271690e-2, .3612855622194901e-3, -.2676991638581043e-2,
          .8805773686384639e-3)
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }
  select.mb24 <- function() {
    L <- 24
    g <- rev(c(-2.132706e-05, 4.745736e-04, 7.456041e-04, -4.879053e-03,
	  -1.482995e-03, 4.199576e-02, -2.658282e-03, -6.559513e-03,
	   1.019512e-01, 1.689456e-01, 1.243531e-01, 1.949147e-01,
	   4.581101e-01, 6.176385e-01, 2.556731e-01, -3.091111e-01,
	  -3.622424e-01, -4.575448e-03, 1.479342e-01, 1.027154e-02,
	  -1.644859e-02, -2.062335e-03, 1.193006e-03, 5.361301e-05))
    h <- qmf(g)
    return(list(length = L, hpf = h, lpf = g))
  }

  switch(name,
    "haar" = select.haar(),
    "d4" = select.d4(),
    "mb4" = select.mb4(),
    "w4" = select.w4(),
    "bs3.1" = select.bs3.1(),
    "fk4" = select.fk4(),
    "d6" = select.d6(),
    "fk6" = select.fk6(),
    "d8" = select.d8(),
    "fk8" = select.fk8(),
    "la8" = select.la8(),
    "mb8" = select.mb8(),
    "bl14" = select.bl14(),
    "fk14" = select.fk14(),
    "d16" = select.d16(),
    "la16" = select.la16(),
    "mb16" = select.mb16(),
    "la20" = select.la20(),
    "bl20" = select.bl20(),
    "fk22" = select.fk22(),
    "mb24" = select.mb24(),
    stop("Invalid selection for wave.filter"))
}

qmf <- function(g, low2high = TRUE) {
  L <- length(g)
  if(low2high)
    h <- (-1)^(0:(L - 1)) * rev(g)
  else
    h <- (-1)^(1:L) * rev(g)
  return(h)
}
.First.lib <- function(lib, pkg) library.dynam("waveslim", pkg, lib)
