R/decz.R

In tagtools: Work with Data from High-Resolution Biologging Tags

#' Decimate sampling rate recursively.
#'
#' Recursive sampling rate decimator. This function can be run iteratively over a long data set, e.g., to decimate an entire recording that is too large to be read into memory.
#'
#' The first time decz is called, use the following format: y = decz(x,df). The subsequent calls to decz for contiguous input data are: decz(x,Z). The final call when there is no more input data is: decz(x = NULL, Z = Z). Each output y in the above contains a segment of the decimated signal and so these need to be concatenated.
#' Decimation is performed in the same way as for \code{\link{decdc}}. The group delay of the filter is removed. For large decimation factors (e.g., df much greater than 50), it is better to perform several nested decimations with lower factors.
#' @param x A vector, matrix, or tag data list containing the signal(s) to be decimated. If x is a matrix, each column is decimated separately. If inputs \code{df} and \code{Z} are both provided, then the value of \code{df} stored in \code{Z} will override the user-provided \code{df}.
#' @param df The decimation factor. The output sampling rate is the input sampling rate divided by df. df must be an integer greater than 1. df can also be a three element vector in which case: df(1) is the decimation factor; df(2) is the number of output samples spanned by the filter (default value is 12). A larger value makes the filter steeper; df(3) is the fractional bandwidth of the filter (default value is 0.8) relative to the output Nyquist frequency. If df(2) is greater than 12, df(3) can be closer to 1.
#' @param nf The number of output samples spanned by the filter (default value is 12). A larger value makes the filter steeper.
#' @param frbw The fractional bandwidth of the filter (default value is 0.8) relative to the output Nyquist frequency. If \code{nf} is greater than 12, \code{frbw} can be closer to 1.
#' @param Z The 'state' list that is generated by a previous call to decz. This is how the function keeps track of filter internal values (i.e., memory) from call-to-call.
#' @export
#' @seealso \code{\link{decdc}}
#' @return A list with elements:
#' \itemize{
#'   \item{\strong{y: }} The decimated signal vector or matrix. It has the same number of columns as x but has, on average, 1/df of the rows.
#'   \item{\strong{Z: }} The state list (for internal tracking of filter internal values). Contains elements df (the decimation factor), nf (used to compute the filter length), frbw (the bandwidth of the filter relative to the new Nyquist frequency), h (the FIR filter coefficients), n (the filter length), z (padded signal used for filtering), and ov ("overflow" samples to be passed to future iterations).
#' }
#' @examples
#' plott_base(list(Accel = beaked_whale$A)) # acceleration data before decimation
#' a_rows <- nrow(beaked_whale$A$data)
#' a_ind <- data.frame(start = c(1, floor(a_rows / 3), floor(2 * a_rows / 3)))
#' a_ind$end <- c(a_ind$start[2:3] - 1, a_rows)
#' df <- 10
#' Z <- NULL
#' y <- NULL
#' for (k in 1:nrow(a_ind)) {
#'   decz_out <- decz(
#'     x = beaked_whale$A$data[c(a_ind[k, 1]:a_ind[k, 2]), ],
#'     df = df, Z = Z
#'   )
#'   df <- NULL
#'   Z <- decz_out$Z
#'   y <- rbind(y, decz_out$y)
#' }
#'
decz <- function(x, df = NULL, Z = NULL, nf = 12, frbw = 0.8) {
  # input checks
  # ================================================
  if (is.null(df) & is.null(Z)) {
    stop("df is required for an initial call to decz, 
         and Z is required for subsequent calls for contiguous data.")
  }

  first_call <- ifelse(is.null(Z), TRUE, FALSE)
  last_call <- ifelse(is.null(x), TRUE, FALSE)


  # make sure x is a one- or multi-column matrix
  if (!is.null(x)) {
    if (is.list(x)) {
      # if input is a sensor data structure, extract the data
      x <- x$data
    }
    if (!is.matrix(x)) {
      # if x is not a matrix, try to make it one
      x <- as.matrix(x)
      # if x is now a one row matrix, make it one column
      if (dim(x)[1] == 1) {
        x <- matrix(x, ncol = 1)
      }
    }
  }

  # filter settings
  # ================================================
  if (is.null(Z)) {
    Z <- list(df = df, frbw = frbw, nf = nf)
    Z$h <- signal::fir1((Z$df * Z$nf), (Z$frbw / Z$df))
    Z$n <- length(Z$h)
    # below number carefully (must be)
    # chosen to precisely match to decdc.
    # do not change!
    npre <- floor(Z$df * (Z$nf - 1) / 2 - 1)
    # =====================================
    Z$z <- rbind(
      2 * x[1, ] - x[(1 + seq(from = (Z$n - Z$df - npre), to = 1)), ],
      x[1:npre, ]
    )
    Z$ov <- NULL
  }

  # if it's the first call - start of data
  # =============================================
  if (first_call) {
    x <- x[c((npre + 1):nrow(x)), ]
  }

  # if it's the last call - no more data
  # =============================================
  if (last_call) {
    # reuse the last few inputs to squeeze some more output
    # from the filter.
    x <- rbind(Z$z, Z$ov)
    x <- rbind(Z$ov, (2 * x[nrow(x), ] - x[(nrow(x) - (1:(npre - 1))), ]))
  } else {
    # if it's neither the first nor the last call
    # =============================================
    if (length(Z$ov) > 0) {
      # combine overflow data from last call with current input data
      x <- rbind(Z$ov, x)
    }
  }

  Z$ov <- NULL
  for (k in 1:ncol(x)) {
    buff_out <- buffer(
      x = x[, k],
      n = Z$n,
      overlap = Z$n - Z$df,
      opt = Z$z[, k]
    )
    if (k == 1) {
      y <- matrix(0, ncol(buff_out$X), ncol(x))
    }
    y[, k] <- t((Z$h %*% buff_out$X))
    Z$z[, k] <- buff_out$opt
    if (is.null(Z$ov)) {
      Z$ov <- matrix(0, length(buff_out$z), ncol(x))
    }
    Z$ov[, k] <- buff_out$z
  }
  return(list(y = y, Z = Z))
}