Contents

multitaper()

Usage:

multitaper(y, W, K [,degree:D,nfreq:Nfreq],deltat:delta_t,wts:wts]),
  REAL vector y, scalar W > 0, integer K >= 1, integer D, integer Nfreq
  >= length(y), REAL scalar delta_t > 0, REAL vector wts with wts[i] > 0

Keywords: spectrum analysis, frequency domain

multitaper(y, W, K, degree:D, deltat:delta_t, nfreq:Nfreq, wts:Wts)
computes multitaper spectrum estimates using K discrete prolate
spheroidal sequences (DPSS) as tapers, with total bandwidth
approximately 2*W cycles per unit time.

The keyword phrases are optional.  See below for defaults.

y is a REAL matrix with no MISSING values whose columns are discrete
parameter time series observed at equally spaced times delta_t units
apart.

delta_t > 0 is a scalar; see below for default.

REAL scalar W specifies the half bandwidth W; 0 < W <.5/delta_t is
required.

Integer K > 0 is the number of DPSS tapers used.  Sensible values for K
should satisfy K < 2*nrows(y)*W*delta_t, but this is not checked.  At a
frequency where the spectrum is smooth, the approximate equivalent
degrees of freedom EDF (see topic 'bandwidth') of an estimate is
2*sum(Wts)^2/sum(Wts^2).  For equal weights, EDF = 2*K.

Wts is a REAL vector of length K with Wts[i] > 0.  Default is Wts =
rep(1/K, K), that is, equal weights.

The estimated spectrum is the weighted average, using weights in Wts,
of K periodograms of detrended columns of y, tapered using DPSS tapers.

See Spectral Analysis for Physical Applications by D. B. Percival and
A. T. Walden (Cambridge Univ. Press, 1993) for details on multitaper
spectrum estimation.

Keyword phrase deltat:delta_t specifies the assumed interval between
observations.  It affects only the interpretation of W which is
interpreted as cycles per unit time.  Because of this, W must satisfy 0
< W < .5/delta_t.

For example, with hourly data and W in units of cycles per day, you
would use deltat:1/24 and W must satisfy 0 < W < 12.  See below for the
default value of delta_t.

Keyword phrase degree:D specifies that the columns of y are detrended
using a polynomial of degree D fit by least squares.  When D < 0 no
detrending is done.  When D = 0 column means are subtracted.

When nfreq:Nfreq is an argument, estimated spectra will be computed at
Nfreq frequencies.  It is an error if Nfreq < nrows(y) or if Nfreq has
a prime factor > 29.

When nfreq:Nfreq is not an argument and S is a positive integer
variable, Nfreq = S.  It is an error if S has a prime factor > 29.
When such an S does not exist, Nfreq = smallest integer >= 2*nrows(y)
with no prime factor > 29, that is Nfreq = goodfactors(@*nrows(y))

The default values of omitted keywords are as follows:
  Keyword    Default value
  degree     D = 0 (subtract only the mean).
  deltat     delta_t = variable DELTAT if it exists; otherwise 1
  nfreq      Nfreq = variable S if it exists or the smallest integer >=
             2*nrows(y) otherwise; it is an error if S has a prime
             factor > 29.
  wts        Wts = rep(1/K,K)

multitaper() uses macros dpss() and detrend().  If they are not defined,
an attempt is made to read them from tser.mac using getmacros.