SpecMTM: Multitaper spectral estimate

View source: R/spec-estimation.R

SpecMTMR Documentation

Multitaper spectral estimate

Description

This is a wrapper for the spec.mtm function from package multitaper, which sets convenient defaults for the spectral estimate and makes the degrees of freedom accessible as a direct list element of the output.

Usage

SpecMTM(
  timeSeries,
  k = 3,
  nw = 2,
  nFFT = "default",
  centre = c("Slepian"),
  dpssIN = NULL,
  returnZeroFreq = FALSE,
  Ftest = FALSE,
  jackknife = FALSE,
  jkCIProb = 0.95,
  maxAdaptiveIterations = 100,
  plot = FALSE,
  na.action = stats::na.fail,
  returnInternals = FALSE,
  detrend = TRUE,
  bPad = FALSE,
  ...
)

Arguments

timeSeries

A time series of equally spaced data, this can be created by the ts() function where deltat is specified.

k

k a positive integer, the number of tapers, often 2*nw.

nw

nw a positive double precision number, the time-bandwidth parameter.

nFFT

This function pads the data before computing the fft. nFFT indicates the total length of the data after padding.

centre

The time series is centred using one of three methods: expansion onto discrete prolate spheroidal sequences ('Slepian'), arithmetic mean ('arithMean'), trimmed mean ('trimMean'), or not at all ('none').

dpssIN

Allows the user to enter a dpss object which has already been created. This can save computation time when Slepians with the same bandwidth parameter and same number of tapers are used repeatedly.

returnZeroFreq

Boolean variable indicating if the zeroth frequency (DC component) should be returned for all applicable arrays.

Ftest

Boolean variable indicating if the Ftest result should be computed and returned.

jackknife

Boolean variable indicating if jackknifed confidence intervals should be computed and returned.

jkCIProb

Decimal value indicating the jackknife probability for calculating jackknife confidence intervals. The default returns a 95% confidence interval.

maxAdaptiveIterations

Maximum number of iterations in the adaptive multitaper calculation. Generally convergence is quick, and should require less than 100 iterations.

plot

Boolean variable indicating if the spectrum should be plotted.

na.action

Action to take if NAs exist in the data, the default is to fail.

returnInternals

Return the weighted eigencoefficients, complex mean values, and so on. These are necessary for extensions to the multitaper, including magnitude-squared coherence (function mtm.coh in this package). Note: The internal ($mtm) variables eigenCoefs and eigenCoefWt correspond to the multitaper eigencoefficients. The eigencoefficients correspond to equation (3.4) and weights, eigenCoefWt, correspond to sqrt(|d_k(f)|^2) from equation (5.4) in Thomson's 1982 paper. This is because the square root values contained in eigenCoefWt are commonly used in additional calculations (example: eigenCoefWt * eigenCoefs). The values returned in mtm$cmv correspond to the the estimate of the coefficients hat(mu)(f) in equation (13.5) in Thomson (1982), or to the estimate of hat(C)_1 at frequency 1 in equation (499) form Percival and Walden (1993)

detrend

Shall the time series data be linearly detrended before computing the spectrum? Defaults to TRUE.

bPad

Shall the time series data be padded before computing the spectrum? Defaults to FALSE.

...

further parameters passed to spec.mtm.

Value

an object of class "spec", with the additional list component dof: a numeric vector of the same length as the spectrum with the degrees of freedom of the spectral estimate (copied from the spec.mtm default output).

Author(s)

Thomas Laepple

References

Thomson, D.J (1982) Spectrum estimation and harmonic analysis. _Proceedings of the IEEE_ Volume *70*, Number 9, pp. 1055-1096.

Percival, D.B. and Walden, A.T. (1993) _Spectral analysis for physical applications_ Cambridge University Press.

See Also

spec.mtm

Examples

x <- ts(arima.sim(list(ar = 0.9), 1000))
spec <- proxysnr:::SpecMTM(x)
proxysnr:::LPlot(spec, col = 'grey', ylim = c(0.01, 100))
proxysnr:::LLines(proxysnr:::LogSmooth(spec), lwd = 2)

EarthSystemDiagnostics/proxysnr documentation built on Sept. 15, 2024, 7:47 a.m.