R/wavelet_components.R

In HydrosystemsGroup/Weather-Generator: Daily weather generator for performing a climate stress test.

#' Reconstruct Wavelet Components from Wavelet Transform
#'
#' @param x annual timeseries of observed climate variable as numeric vector
#' @param wt wavelet transform generated by wavelet_analysis()
#' @param n.periods number of final composite periods
#' @param sig.periods integer array of significance periods
#' @param n.comp.periods number of significance periods for each composite period
#' @export
#' @examples
#' wavelet_components(x=x, n.periods=1, sig.periods=seq(8,15), n.comp.periods=8)
wavelet_components <- function(x, wt, n.periods, sig.periods, n.comp.periods) {
  freq.comps <- array(0,c(length(x),n.periods))
  for (i in 1:n.periods) {
    cur.periods <- sig.periods[1:n.comp.periods[i]]
    if (i>1) {
      cur.periods <- sig.periods[(1 + (i-1)*n.comp.periods[i-1]):(n.comp.periods[i] + (i-1)*n.comp.periods[i-1])]
    }
    sj <- wt$scale[cur.periods]
    #for Morlet Wavelet with freq = 6
    Cdelta <- .776
    w0_0 <- pi^(-1/4)
    components <- (wt$dj*sqrt(wt$dt)/(Cdelta*w0_0))*Re(wt$wave)[cur.periods,]/sqrt(sj)
    if (length(cur.periods)>1) {
      components <- apply(components, 2, sum)
    }
    freq.comps[,i] <- components
  }

  noise <- x - apply(freq.comps, 1, sum)

  colnames(freq.comps) <- paste("COMPONENT", 1:n.periods, sep=' ')

  return(cbind(NOISE=noise, freq.comps))
}