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) {
    my.acf <- function(x)
      {
        n <- sum(!is.na(x))
        x <- c(x[!is.na(x)], rep(0, n))
        Re(fft(Mod(fft(x))^2, inv = TRUE)/2/n^2)[1:n]
      }
    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)
    ## cat(paste("  max{M_j/2^j, 1} =", eta3, "\n\n"))
    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))
}

##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

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

## 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=T)
  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)
  }

  if(hard) {
    # Hard thresholding
    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))
      s <- cumsum(sxi) + ((ni - 1):0) * sxi
      risk <- (ni - (2 * (1:ni)) + s) / ni
      surethresh <- sqrt(sxi[order(risk)[1]])
      wc.shrink[[i]] <- wci * (abs(wci) > surethresh)
    }
  }
  else {
    # Soft thresholding
    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))
      s <- cumsum(sxi) + ((ni - 1):0) * sxi
      risk <- (ni - (2 * (1:ni)) + s) / ni
      surethresh <- sqrt(sxi[order(risk)[1]])
      wc.shrink[[i]] <- sign(wci) * (abs(wci) - surethresh) *
        (abs(wci) > surethresh)
    }
  }
  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)
}
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))
        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))
        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))$out
        else
          X <- .C("idwt", as.double(y[[n]]), as.double(y[[n+1]]), 
                  as.integer(m), L, h, g, out = as.double(XX))$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=F)
{
  ## 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)[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)[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)
}

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]
}

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)
}
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))[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))$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)[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)$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)
}
hosking.sim  <-  function(n, acvs) {
  .C("hosking", tseries=rnorm(n), as.integer(n), as.double(acvs[1:n]))$tseries
}
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))
  ##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)
}
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)
