R/func_pilot.R

Defines functions pilot_spec.default pilot_spec.matrix pilot_spec.ts pilot_spec

Documented in pilot_spec pilot_spec.default pilot_spec.matrix pilot_spec.ts

#' Calculate initial power spectral density estimates
#'
#' @description
#' This PSD is used as the starting point -- the pilot spectrum -- for
#' the adaptive estimation routine.
#'
#' @details
#' A fixed number
#' of tapers is applied across all frequencies using \code{\link{psdcore}}, and
#' subsequent taper-refinements are based on the spectral derivatives
#' of this spectrum; hence, changes in the number of tapers can affect
#' how many adaptive stages may be needed (though there are no formal convergence
#' criteria to speak of).
#'
#' The taper series of the returned spectrum is constrained using
#' \code{as.tapers(..., minspan=TRUE)}.
#'
#' The default behavior (\code{remove.AR <= 0}) is to remove the standard linear 
#' model \eqn{[f(x) = \alpha x + \beta]} from the data; however,
#' the user can model the effect of an autoregressive process by specifying
#' \code{remove.AR}.
#'
#' @section Removing an AR effect from the spectrum:
#' If \code{remove.AR > 0} the argument is used as \code{AR.max} in 
#' \code{\link{prewhiten}}, from which an AR-response spectrum is calculated using
#' the best fitting model.
#'
#' If the value of \code{remove.AR} is too low the spectrum 
#' could become distorted,
#' so use with care.
#' \emph{Note, however, that the 
#' value of \code{remove.AR} will be restricted to within the 
#' range \eqn{[1,100]}.}
#' If the AR order is much larger than this, it's unclear how \code{\link{prewhiten}}
#' will perform and whether the AR model is appropriate.
#'
#' \emph{Note that this function does not produce a parametric spectrum estimation; rather,
#' it will return the amplitude response of the best-fitting AR model as \code{\link{spec.ar}}
#' would. \strong{Interpret these results with caution, as an AR response spectrum
#' can be misleading.}}
#'
#' @name pilot_spec
#' @aliases pilot_spectrum spec.pilot
#' @export
#' @author A.J. Barbour
#' @seealso \code{\link{psdcore}}, \code{\link{prewhiten}}, \code{\link{spec.ar}}
#'
#' @param x  vector; the data series to find a pilot spectrum for
#' @param x.frequency  scalar; the sampling frequency (e.g. Hz) of the series
#' @param ntap  scalar; the number of tapers to apply during spectrum estimation
#' @param remove.AR  scalar; the max AR model to be removed from the data.
#' @param plot  logical; should a plot be created?
#' @param verbose  logical; should messages be given?
#' @param fast logical; use fast method in \code{\link{psdcore}}?
#' @param ...  additional parameters passed to \code{\link{psdcore}}
#' @return Invisibly, an object with class \code{'spec'}, and 
#' \code{"pilot_psd"} in the working environment.
#'
#' @example inst/Examples/rdex_pilotspec.R
pilot_spec <- function(x, ...) UseMethod("pilot_spec")

#' @rdname pilot_spec
#' @export
pilot_spec.ts <- function(x, ...){
  stopifnot(is.ts(x))
  frq <- stats::frequency(x)
  args <- list(...)
  args[['x']] <- x
  args[['x.frequency']] <- frq
  do.call('pilot_spec.default', args)
}

#' @rdname pilot_spec
#' @aliases pilot_spec.matrix
#' @export
pilot_spec.matrix <- function(x, x.frequency, ...){
  frq <- if (missing(x.frequency)){
    stats::frequency(x)
  } else {
    x.frequency
  }
  pilot_spec(stats::ts(x, frequency=frq), ...)
}

#' @rdname pilot_spec
#' @export
pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, plot=FALSE, verbose=FALSE, fast = FALSE, ...){
  if (is.null(ntap)) ntap <- 7
  if (is.null(remove.AR)) remove.AR <- 0
  if (is.null(x.frequency)) x.frequency <- 1
  stopifnot(length(ntap)==1 | length(ntap)==nrow(x) %/% 2)
  stopifnot(length(remove.AR)==1)
  stopifnot(length(x.frequency)==1)

  # AR spectrum or no?
  REMAR <- ifelse(remove.AR > 0, TRUE, FALSE)
  
  #restrict maximum ar orders to within [1,100]
  if (REMAR) remove.AR <- max(1, min(100, abs(remove.AR)))
  
  xprew <- prewhiten(x, x.fsamp=x.frequency, AR.max=remove.AR, detrend=TRUE, impute=TRUE, plot=FALSE, verbose=verbose)
  # (result is ts)

  ## Remove and AR model
  if (REMAR){
    # AR fit
    ordAR <- xprew[['ardfit']][['order']]
    if (ordAR==0){
      warning("AR(0) was the highest model found!\n\t\tConsider fitting a linear model instead ( remove.AR = 0 ).")
    } else {
      if (verbose) message(sprintf("removed AR(%s) effects from the spectrum", ordAR))
    }
    # calculate PSD of the AR fit
    xar <- xprew[['prew_ar']] # ts object
    Pspec_ar <- psdcore(xar, ntaper=ntap, AR=TRUE, preproc=FALSE, refresh=TRUE, verbose=FALSE, fast = fast)
    arvar <- stats::var(Pspec_ar[['spec']])
    mARs <- mean(Pspec_ar[['spec']])
    Pspec_ar[['spec']] <- Pspec_ar[['spec']] / mARs
  }
  
  ## Calculate spectrum of non-AR (lm) model:
  xlm <- xprew[['prew_lm']] # ts object
  Pspec <- psdcore(xlm, ntaper=ntap, AR=FALSE, preproc=FALSE, refresh=TRUE, verbose=FALSE, fast = fast)
  Ptap <- Pspec[['taper']]
  num_tap <- length(Ptap)
  num_frq <- Pspec[['numfreq']]
  
  ## generate a series, if necessary
  if (num_tap < num_frq) Ptap <- rep.int(Ptap[1], num_frq)
  
  ## return tapers object
  Pspec[['taper']] <- as.tapers(Ptap, setspan=TRUE)
  
  ## remove the spectrum of the AR process
  if (REMAR){
    stopifnot(exists("Pspec") & exists("Pspec_ar"))
    if (verbose) message(sprintf("Removing AR(%s) effects from spectrum", ordAR))
    # reup the spectrum
    Ospec <- psd_envAssignGet("pre_AR_psd", Pspec)
    psd_envAssign("AR_psd", Pspec_ar)
    p.lin <- Pspec[['spec']]
    p.ar <- Pspec_ar[['spec']]
    Pspec[['spec']] <-  p.lin / p.ar
  }
  
  if (plot){
    ttl <- "Pilot spectrum estimation"
    llog <- 'dB'
    try({
      if (REMAR){
        if (verbose) message('Plotting,', tolower(ttl))
        graphics::par(las=1)
        plot(Ospec, log=llog, col="red", main=ttl)
        graphics::mtext(sprintf("(with AR(%s) response)", ordAR), line=0.4)
        # rescale
        Pspec_ar[['spec']] <- Pspec_ar[['spec']] * mARs
        plot(Pspec_ar, log=llog, col="blue", add=TRUE)
        plot(Pspec, log=llog, add=TRUE, lwd=2)
        legend("bottomleft", 
               c("original PSD",
                 sprintf("AR-innovations PSD\n(mean %.01f +- %.01f dB)", dB(mARs), dB(sqrt(arvar))/4),
                 "AR-filter response"), 
               lwd=2, col=c("red","blue","black"))
      } else {
        plot(Pspec, log=llog, main=ttl)
      }
    })
  }
  return(invisible(psd_envAssignGet("pilot_psd", Pspec)))
}
abarbour/psd documentation built on Aug. 15, 2023, 8:56 a.m.