R/hilbert.R

In neuroconductor/waveslim: Basic Wavelet Routines for One-, Two-, and Three-Dimensional Signal Processing

#' Discrete Hilbert Wavelet Transforms
#' 
#' The discrete Hilbert wavelet transforms (DHWTs) for seasonal and
#' time-varying time series analysis.  Transforms include the usual orthogonal
#' (decimated), maximal-overlap (non-decimated) and maximal-overlap packet
#' transforms.
#' 
#' @usage dwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...)
#' @usage dwt.hilbert.nondyadic(x, ...) 
#' @usage idwt.hilbert(y)
#' @usage modwt.hilbert(x, wf, n.levels = 4, boundary = "periodic", ...)
#' @usage imodwt.hilbert(y)
#' @usage modwpt.hilbert(x, wf, n.levels = 4, boundary = "periodic")
#' @aliases dwt.hilbert dwt.hilbert.nondyadic idwt.hilbert modwt.hilbert
#' imodwt.hilbert modwpt.hilbert
#' @param x Real-valued time series or vector of observations.
#' @param wf Hilbert wavelet pair
#' @param n.levels Number of levels (depth) of the wavelet transform.
#' @param boundary Boundary treatment, currently only \code{periodic} and
#' \code{reflection}.
#' @param \ldots Additional parametes to be passed on.
#' @param y An object of S3 class \code{dwt.hilbert}.
#' @return Hilbert wavelet transform object (list).
#' @author B. Whitcher
#' @seealso \code{\link{hilbert.filter}}
#' @references Selesnick, I. (200X). \emph{IEEE Signal Processing Magazine}
#' 
#' Selesnick, I. (200X). \emph{IEEE Transactions in Signal Processing}
#' 
#' Whither, B. and P.F. Craigmile (2004). Multivariate Spectral Analysis Using
#' Hilbert Wavelet Pairs, \emph{International Journal of Wavelets,
#' Multiresolution and Information Processing}, \bold{2}(4), 567--587.
#' @keywords ts
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(C_dwt, as.double(x.h), as.integer(N/2^(j-1)), L, h1, h0, 
                W = W, V = V)[6:7]
    out.g <- .C(C_dwt, as.double(x.g), as.integer(N/2^(j-1)), L, g1, g0, 
                W = W, V = V)[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" =  y <- c(y, rev(y)),
    "periodic" = invisible(),
    stop("Invalid boundary rule in dwt.dbp"))
  J <- attributes(y)$levels
 
  dict <- hilbert.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(C_idwt, y[[jj]], as.double(X), as.integer(length(X)),
            L, h, g, XX=XX)$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(C_modwt, as.double(x.h), N, as.integer(j), L, h1, h0,
                W = W, V = V)[7:8]
    out.g <- .C(C_modwt, as.double(x.g), N, as.integer(j), L, g1, g0,
                W = W, V = V)[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(attributes(y)$boundary != "periodic")
    stop("Invalid boundary rule in imodwt")
  J <- length(y) - 1

  dict <- hilbert.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(C_imodwt, y[[jj]], X, as.integer(N), as.integer(j), 
      as.integer(L), ht, gt, XX)[[8]]
  }
  return(X)
}

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



#' Select a Hilbert Wavelet Pair
#' 
#' Converts name of Hilbert wavelet pair to filter coefficients.
#' 
#' Simple \code{switch} statement selects the appropriate HWP.  There are two
#' parameters that define a Hilbert wavelet pair using the notation of
#' Selesnick (2001,2002), \eqn{K} and \eqn{L}.  Currently, the only implemented
#' combinations \eqn{(K,L)} are (3,3), (3,5), (4,2) and (4,4).
#' 
#' @param name Character string of Hilbert wavelet pair, see acceptable names
#' below (e.g., \code{"k3l3"}).
#' @return List containing the following items: \item{L}{length of the wavelet
#' filter} \item{h0,g0}{low-pass filter coefficients} \item{h1,g1}{high-pass
#' filter coefficients}
#' @author B. Whitcher
#' @seealso \code{\link{wave.filter}}
#' @references Selesnick, I.W. (2001). Hilbert transform pairs of wavelet
#' bases. \emph{IEEE Signal Processing Letters} \bold{8}(6), 170--173.
#' 
#' Selesnick, I.W. (2002). The design of approximate Hilbert transform pairs
#' of wavelet bases. \emph{IEEE Transactions on Signal Processing} 
#' \bold{50}(5), 1144--1152.
#' @keywords ts
#' @examples
#' 
#' hilbert.filter("k3l3")
#' hilbert.filter("k3l5")
#' hilbert.filter("k4l2")
#' hilbert.filter("k4l4")
#' 
#' @export hilbert.filter
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)))
  }
  select.K5L7 <- function() {
    L <- 24
    h0 <- c(-2.5841959496364648e-10, 6.0231243121018760e-10,
            2.1451486802217960e-06,  -4.9989222844980982e-06,
            -2.2613489535132104e-04, 5.1967501391358343e-04,
            3.4011963595840899e-03,  -7.1996997688061597e-03,
            -1.7721433874932836e-02, 3.5491112173858148e-02,
            3.0580617312936355e-02,  -1.3452365188777773e-01,
            2.1741748603083836e-03, 5.8046856094922639e-01,
            7.4964083145768690e-01, 2.6775497264154541e-01,
            -7.9593287728224230e-02,  -4.3942149960221458e-02,
            1.9574969406037097e-02, 8.8554643330725387e-03,
            -7.2770446614145033e-04,  -3.1310992841759443e-04,
            1.4045333283124608e-06, 6.0260907100656169e-07)
    h1 <- qmf(h0)
    g0 <- c(-3.8762939244546978e-09, 2.9846463282743695e-07,
            5.6276030758515370e-06, -7.7697066311187957e-05,
            -2.1442686434841905e-04, 2.1948612668324223e-03,
            9.5408758453423542e-04, -1.7149735951945008e-02,
            1.5212479104581677e-03, 5.6600564413983846e-02,
            -4.8900162376504831e-02, -1.3993440493611778e-01,
            2.7793346796113222e-01, 7.6735603850281364e-01,
            5.4681951651005178e-01, 3.6275855872448776e-02,
            -8.8224410289407154e-02, 3.2821708368951431e-05,
            1.7994969189524142e-02, 1.8662128501760204e-03,
            -7.8622878632753014e-04, -5.8077443328549205e-05,
            3.0932895975646042e-06, 4.0173938067104100e-08)
    g1 <- qmf(g0)
    return(list(length = L, hpf = list(h1, g1), lpf = list(h0, g0)))
  }
  select.K6L6 <- function() {
    L <- 24
    h0 <- c(1.4491207137947255e-09  -3.4673992369566253e-09,
            -6.7544152844875963e-06, 1.6157040144070828e-05,
            4.0416340595645441e-04, -9.4536696039781878e-04,
            -4.2924086033924620e-03, 9.0688042722858742e-03,
            1.8690864167884680e-02, -3.7883945370993717e-02,
            -2.7337592282061701e-02, 1.3185812419468312e-01,
            -2.1034481553730465e-02, -5.9035515013747486e-01,
            -7.4361804647499452e-01, -2.5752016951708306e-01,
            9.2725410672739983e-02, 4.9100676534870831e-02,
            -2.4411085480175867e-02, -1.1190458223944993e-02,
            1.7793885751382626e-03, 7.4715940333597059e-04,
            -6.2392430013359510e-06, -2.6075498267775052e-06)
    h1 <- qmf(h0)
    g0 <- c(1.8838569279331431e-08, -1.1000360697229965e-06,
            -1.4600820117782769e-05, 1.6936567299204319e-04,
            2.6967189953984829e-04, -3.1633669438102655e-03,
            -7.2081460313487946e-04, 1.9638595542490079e-02,
            -3.0968325940269846e-03, -5.6722348677476261e-02,
            5.2260784738219289e-02, 1.2763836788794369e-01,
            -2.9566169882112192e-01, -7.6771793937333599e-01,
            -5.3818432160802543e-01, -2.4023872575927138e-02,
            9.9019132161496132e-02, -1.2059411664071501e-03,
            -2.2693488886969308e-02, -1.8724943382560243e-03,
            1.7270823778712107e-03, 1.5415480681200776e-04,
            -1.1712464100067407e-05, -2.0058075590596196e-07)
    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(),
         "k5l7" = select.K5L7(),
         "k6l6" = select.K6L6(),
         stop("Invalid selection for hilbert.filter"))
}

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



#' Phase Shift for Hilbert Wavelet Coefficients
#' 
#' Wavelet coefficients are circularly shifted by the amount of phase shift
#' induced by the discrete Hilbert wavelet transform.
#' 
#' The "center-of-energy" argument of Hess-Nielsen and Wickerhauser (1996) is
#' used to provide a flexible way to circularly shift wavelet coefficients
#' regardless of the wavelet filter used.
#' 
#' @aliases phase.shift.hilbert phase.shift.hilbert.packet
#' @param x Discete Hilbert wavelet transform (DHWT) object.
#' @param wf character string; Hilbert wavelet pair used in DHWT
#' @return DHWT (DHWPT) object with coefficients circularly shifted.
#' @author B. Whitcher
#' @seealso \code{\link{phase.shift}}
#' @references Hess-Nielsen, N. and M. V. Wickerhauser (1996) Wavelets and
#' time-frequency analysis, \emph{Proceedings of the IEEE}, \bold{84}, No. 4,
#' 523-540.
#' @keywords ts
#' @export phase.shift.hilbert
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(C_modwt, as.double(Re(x)), N, as.integer(j), L, h1, h0, 
                 W = W, V = V)[7:8]
        zc <- .C(C_modwt, as.double(Im(x)), N, as.integer(j), L, g1, g0, 
                 W = W, V = V)[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(C_modwt, as.double(Re(x)), N, as.integer(j), L, h1, h0, 
                 W = W, V = V)[7:8]
        zc <- .C(C_modwt, as.double(Im(x)), N, as.integer(j), L, g1, g0, 
                 W = W, V = V)[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)
}

#' Time-varying and Seasonal Analysis Using Hilbert Wavelet Pairs
#' 
#' Performs time-varying or seasonal coherence and phase anlaysis between two
#' time seris using the maximal-overlap discrete Hilbert wavelet transform
#' (MODHWT).
#' 
#' The idea of seasonally-varying spectral analysis (SVSA, Madden 1986) is
#' generalized using the MODWT and Hilbert wavelet pairs.  For the seasonal
#' case, \eqn{S} seasons are used to produce a consistent estimate of the
#' coherence and phase.  For the non-seasonal case, a simple rectangular
#' (moving-average) filter is applied to the MODHWT coefficients in order to
#' produce consistent estimates.
#' 
#' @usage modhwt.coh(x, y, f.length = 0)
#' @usage modhwt.phase(x, y, f.length = 0)
#' @usage modhwt.coh.seasonal(x, y, S = 10, season = 365)
#' @usage modhwt.phase.seasonal(x, y, season = 365)
#' @aliases modhwt.coh modhwt.phase modhwt.coh.seasonal modhwt.phase.seasonal
#' @param x MODHWT object.
#' @param y MODHWT object.
#' @param f.length Length of the rectangular filter.
#' @param S Number of "seasons".
#' @param season Length of the "season".
#' @return Time-varying or seasonal coherence and phase between two time
#' series.  The coherence estimates are between zero and one, while the phase
#' estimates are between \eqn{-\pi}{-pi} and \eqn{\pi}{pi}.
#' @author B. Whitcher
#' @seealso \code{\link{hilbert.filter}}
#' @references Madden, R.A. (1986). Seasonal variation of the 40--50 day
#' oscillation in the tropics. \emph{Journal of the Atmospheric Sciences}
#' \bold{43}(24), 3138--3158.
#' 
#' Whither, B. and P.F. Craigmile (2004). Multivariate Spectral Analysis Using
#' Hilbert Wavelet Pairs, \emph{International Journal of Wavelets,
#' Multiresolution and Information Processing}, \bold{2}(4), 567--587.
#' @keywords ts
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
}