R/windowSmooth.R

In adept: Adaptive Empirical Pattern Transformation

#' @title
#' Fast Computation of Moving Window Average
#'
#' @description
#' Compute moving window average of a time-series \code{x}.
#'
#' @param x A numeric vector. A time-series for which a moving window average is computed.
#' @param W A numeric scalar. A length of a moving window, expressed in time (seconds).
#' @param x.fs Frequency of a time-series \code{x}, expressed in a number of observations
#' per second. Defaults to \code{1}.
#'
#' @details
#' Time-series frequency \code{x.fs} and a length  of a moving window (expressed in time) \code{W}
#' together determine
#' \code{W.vl = round(W * x.fs)}, a length of a moving window expressed in a length of
#' \code{x} vector object.
#' Note: \code{W.vl} must be equal or greater than \code{3}.
#' \itemize{
#'   \item If \code{W.vl < 3} then an error is thrown.
#'   \item If \code{W.vl} is an even number then \code{(W.vl-1)} value is silently
#'   used instead as a length of a moving window expressed in \code{x} vector length.
#' }
#'
#' @return A numeric vector. Moving window average of a time-series \code{x}. Note:
#' head and tail of the
#' output vector where the moving window is undefined are filled with \code{NA}.
#'
#' @export
#' @import dvmisc
#'
#' @examples
#' ## Time-series defined as a function f(x) = x
#' N <- 100
#' W  <- 20
#' x <- 1:N
#' x.smoothed <- windowSmooth(x, W)
#' plot(x, type = "l")
#' points(x.smoothed, col = "red")
#'
#' ## Time-series defined as a function f(x) = sin(x) + noise
#' N <-  1000
#' W  <- 100
#' x <- sin(seq(0, 4 * pi, length.out = N)) + rnorm(N, sd = 0.1)
#' x.smoothed <- windowSmooth(x, W)
#' plot(x, type = "l")
#' points(x.smoothed, col = "red")
#'
windowSmooth <- function(x, W, x.fs = 1){

  ## Check function arguments for correctness
  if (!(is.numeric(x) & is.atomic(x))) stop("x must be a numeric (atomic) vector.")
  if (!(is.numeric(W) & is.atomic(W) & length(W) == 1))  stop("W must be a numeric scalar.")
  if (!(is.numeric(x.fs) & is.atomic(x.fs) & length(x.fs) == 1))  stop("x.fs must be a numeric scalar.")

  W.vl <- round(W * x.fs)
  if (W.vl < 3) stop("W.vl (refer to function's details description) must not be smaller than 3 vector indices. Define wider W averaging window length")
  if (W.vl >= length(x)) stop("W.vl (refer to function's details description) must be less than x vector length. Define narrower W averaging window length")

  ## Replace W with closest odd integer no larger than W
  W.vl <-  W.vl + (W.vl %% 2) - 1

  ## Replace head and tail of a signal with NA's
  W.wing <- floor(W.vl/2)
  x.out  <- c(rep(NA, W.wing), moving_mean(x, W.vl), rep(NA, W.wing))

  return(x.out)
}




#' Wrapper function for \code{windowSmooth}
#'
#' Wrapper function for \code{windowSmooth}. Replaces \code{NA} values with
#' which appear in head and tail of smoothed signal as a result of MA not defined
#' for first and last \code{k} elements of a signal. The \code{NA} values are replaced
#' with sample means of subsequent/procceeding non-\code{NA} values that appear
#' in smoothed signal.
#'
#'
#' @param x Numeric vector to be smoothed.
#' @param W A width of a moving window given in time (seconds).
#' @param x.fs Frequency of \code{x} expressed in number of observations collected
#' per second. Defaults to \code{1}.
#' @param NA.repl.source.k Scalar for number of subsequent/procceeding non-\code{NA} values that appear
#' in smoothed signal and are used to provide sample mean replacement for \code{NA} values.
#'
#' @return Smoothed signal with  \code{NA} values replaced with sample means of
#' subsequent/procceeding non-\code{NA} values that appear
#' in smoothed signal.
#'
#' @noRd
#'
get.x.smoothed <- function(x, W, x.fs = 1, NA.repl.source.k = 4,
                           x.cut.vl = 100 * 60 * 60 * 24){

  W.vl <- W * x.fs

  ## Nested fuction responsible for filling the NAs
  x.smoothed.fill.NA <- function(x.smoothed){
    ## Replace NA's in head/tail of smoothed signal with some neutral average flat line
    ## Vector length of replacement NA's area
    NA.vl <- floor((W.vl + (W.vl %% 2) - 1)/2)
    ## Replace NAs in vector's head
    repl.head.idx        <- 1:NA.vl
    repl.head.source.idx <- 1:(NA.repl.source.k * NA.vl)
    x.smoothed[repl.head.idx] <- mean(x.smoothed[repl.head.source.idx], na.rm = TRUE)
    ## Replace NAs in vector's tail
    repl.tail.idx        <- length(x.smoothed) - ((NA.vl - 1):0)
    repl.tail.source.idx <- length(x.smoothed) - ((NA.repl.source.k * NA.vl - 1):0)
    x.smoothed[repl.tail.idx] <- mean(x.smoothed[repl.tail.source.idx], na.rm = TRUE)
    return(x.smoothed)
  }

  ## Decide whether to smooth on x parts or not
  x.vl <- length(x)
  if (x.vl > x.cut.vl){
    rn <- ceiling(x.vl / x.cut.vl)
    x.split.idx <- unlist(lapply(1:rn, function(i) rep(i, x.cut.vl)))[1:x.vl]
    if ((sum(x.split.idx == rn) < 100 * 60) & (rn > 1)){
      x.split.idx[x.split.idx == rn] <- rn - 1
    }
    x.split <- split(x, x.split.idx)
    ## Apply to each part of x split
    x.smoothed.list <- lapply(x.split, function(x.split.i){
      x.smoothed.i <- windowSmooth(x = x.split.i, W = W, x.fs = x.fs)
      x.smoothed.i <- x.smoothed.fill.NA(x.smoothed.i)
      return(x.smoothed.i)
    })
    x.smoothed <- unlist(x.smoothed.list)

  } else {
    x.smoothed <- windowSmooth(x = x, W = W, x.fs = x.fs)
    x.smoothed <- x.smoothed.fill.NA(x.smoothed)
  }

  return(x.smoothed)
}