**The MultiTaper method (MTM)** of spectral analysis provides a
novel means for spectral estimation (Thomson, 1982;
Percival and Walden, 1993) and signal reconstruction
(e.g., Park, 1992)
of a time series which is believed to exhibit a spectrum
containing both continuous and singular components.
This method has been widely applied
to problems in geophysical signal analysis, including
analyses of atmospheric and oceanic data, paleoclimate data, geochemical tracer data, and
seismological data.

MTM, like the method of Blackman and Tukey,
offers the appeal of being nonparametric,
in that it does not prescribe an a priori (e.g., autoregressive)
model for the process generating the time
series under analysis. MTM
attempts to reduce the variance of spectral estimates by
using a small set of tapers
rather than the unique data taper or spectral window used by
Blackman-Tukey methods.
A set of independent estimates of the power spectrum is computed, by
pre-multiplying the data by orthogonal tapers which are
constructed to minimize the spectral leakage due
to the finite length of the data set. The optimal tapers or
``eigentapers'' belong to a family of functions known as
discrete prolate spheroidal sequences (DPSS) which are defined as the
eigenvectors of a suitable Rayleigh-Ritz minimization problem.
Averaging over this (small) ensemble of spectra yields a better and more
stable estimate -- i.e., one with lower variance -- than do single-taper
methods.
The tapers are the discrete set of eigenfunctions which solve
the variational problem of minimizing leakage outside
of a frequency band of
half bandwidth *pf _{n}* where

In practice, only
the first 2*p* - 1 tapers provide usefully small spectral-leakage,
and thus the number of tapers used *K* should
be less than 2*p*-1 in any application of MTM.
The choice of the bandwidth 2*pf _{n}* and number of tapers

For typical length instrumental climate records, the
choice *K*=3, *p*=2 offers a good compromise between
the required frequency resolution for resolving distinct
climate signals (e.g., ENSO and decadal-scale variability)
and the benefit of multiple spectral degrees of freedom.
Longer datasets
can admit the use of a greater number *K* of tapers while
maintaining a desired frequency resolution, and the optimal
choice of *K* and *p* is in general most decidedly
application specific.

MTM can provide estimates of both the singular components (i.e,
the `** `lines''**) and the

and the weights *&mu _{k}* are based on the eigenvalues from the DPSS solution. The

By adjusting the relative weights on the contributions
from each of the *K* eigenspectra, a more leakage-resistant
spectral estimate termed the ** adaptively weighted multitaper
spectrum** is obtained,

where* b _{k}(f)* is a weighting function that further guards
against broadband leakage for a non-white (``colored'') but
locally-white process. The adaptive
spectrum estimate has an effective degrees of freedom that
generally departs only slightly from the nominal value 2

The purpose of *harmonic analysis* is to determine the line
components-ie the spikes in the spectrum corresponding to
a
periodic or quasi-periodic signal-in terms of their frequency, amplitude,
and phase. The Fourier
transform of a clean periodic signal of infinite length yields a *Dirac function*
at the frequency of the signal, viz., a line (or peak of zero width)
with infinite magnitude. A
spectral estimate based on the methods discussed in earlier
sections,
gives indirect information on the amplitude of the
signal at a given frequency, through the area under the peak centered at
that frequency and whose width is, roughly speaking, inversely
proportional to the length *N* of the time series; this area is
nearly constant, since the height of the peak is also proportional to *N*. Harmonic analysis attempts, instead, to
determine directly the (finite) amplitude of a (pure) line in the
spectrum of a time series of finite length. We
explain next how this is done within MTM.

Assume the time series *x*(*t*) is the sum of a sinusoid of
frequency
*f _{0}* and amplitude

If (*w _{k}*(

where the asterisk denotes complex conjugation. A statistical
confidence interval
can be given for the least-squares fit by a
Fisher-Snedecor test, or * F-test*.
This test is
roughly based on the ratio of the variance captured by the filtered
portion of the time series

and

that are respectively the ``explained'' and ``unexplained'' contributions to the variance.

The random variable

would obey a Fisher-Snedecor law with 2 and 2*K*-2 degrees of freedom if
the time series *x*(*t*) were a pure white-noise realization.
One can interpret its numerical value for given data by assuming that B=0 -- i.e., that* x*(*t*) is white- and trying to reject the white noise null hypothesis.
In practice, the spectrum need only be ``locally white'' in the
sense that the *K* eigenspectra which describe the local characteristics
of the spectrum should be distributed as they would be for white noise.

This harmonic-analysis application of MTM is able to detect low-amplitude
harmonic oscillations in a relatively short time series with a high degree of
statistical significance or to
reject
a large amplitude if it failed the *F*-test, because the *F*-value* F*(*f*) does not depend -- to first order -- on the magnitude of *B*(*f*). This feature is an important advantage of MTM over
standard methods where
error bars are essentially
proportional to the amplitude of a peak.
However, the implicit assumption in the harmonic-analysis appproach is that
the time series is produced by a process that consists of a
superposition of separate, purely periodic fixed-amplitude components.
If not, a continuous
spectrum (in the case of a colored noise or a chaotic system) will be broken
down into spurious lines with arbitrary frequencies and possibly high *F*-values. In essence, the above procedure assumes that ``signals''
are represented by lines in the spectrum corresponding to phase-coherent
harmonic signals, while the ``noise'' is represented by the continuous
component of the spectrum. In geophysical applications, signals are
frequently associated with narrowband, but not strictly harmonic
variability and truly harmonic
signals are rarely detected. In such cases, the simple
noise null hypothesis and signal assumptions implicit in the
traditional MTM approach described above loses much of its utility.

The more subtle nature of geophysical signals and noise has led to a more recent generalization of the conventional MTM approach which combines the harmonic signal detection procedure described above, with a criterion for detecting significant narrowband, ``quasi-oscillatory'' signals which may exhibit phase and amplitude modulation, and intermittent oscillatory behavior, but are nonetheless significant relative to some suitably defined null hypothesis. We favour this approach, whichs provides for the detection and significance estimation of both harmonic and anharmonic, narrowband signals in the MTM spectra of time series, making use of a ``robust'' estimate of the background noise.

Information from the harmonic peak test of the traditional MTM procedure is retained, but peaks, whether indicated as harmonic or anharmonic, are tested for significance relative to the null hypothesis of a globally red (or, trivially, white) noise background estimated empirically from the data. This is particularly important in climate studies, where the intrinsic inertia of the system leads to greater power (and greater liklihood of prominent peaks in the spectrum) at lower frequencies, even in the absence of any signals. To accomodate the red noise background assumption requisite in geophysical applications, an AR(1) noise process is assumed, although more complex noise models can be motiviated.

The power spectrum of the AR(1) process is given by:

,

where *P _{0}* is the average value of the power spectrum,
related
to the variance

while *r *is the lag-one autocorrelation, and *f _{N}*=0.5/

For periodicities much larger than *&tau*, the spectrum
behavies like a white spectrum.

The approach of Mann and Lees (1996) performs a ``robust''
estimate of the red noise background
by minimizing (as a function of *r*)
the misfit between an analytical AR(1) red noise
spectrum and the adaptively weighted multitaper spectrum
convoluted with a median smoother. The median smoothing
operation in the fit insures that the estimated noise
background is
insensitive to outliers (most obviously,
peaks associated with signals) that, properly, should
not influence
the estimation of the global noise background. This
guards against the typical problem of inflated estimates
of noise variance and noise
autocorrelation that arise in a
conventional Box-Jenkins approach due to the contamination
of noise parameters by non-noise contributions (e.g., a
significant trend or oscillatory component of the series).
Mann and Lees (1996)
motivate the choice of a median smoothing width of *&Delta* f =min(*f _{N}*/4,

Significance levels for harmonic or narrowband spectral features relative to the estimated noise background can be determined from the appropriate quantiles of the chi-squared distribution, approximating the spectrum as being distributed with 2K degrees of freedom. A reshaped spectrum is determined in which the contributions from harmonic signals are removed, depending on a significance threshold for the F variance-ratio test described above. In this way, noise background, harmonic, and anharmonic narrowband signals are each independently isolated. The harmonic peak detection procedure provides information as to whether the signals are best approximated as harmonic (phase-coherent sinusoidal oscillations) or narrowband (amplitude and phase modulated, perhaps intermittant oscillations). In either case, they must be significant relative to a specified (e.g., red) noise hypothesis to be isolated as significant.

Once significant peaks have been isolated in the spectrum, relative to the specified null hypothesis, the associated signals can be reconstructed in the time domain using the information from the multitaper decomposition. The signals or ``reconstructed components'' are analagous to those of SSA described in previous sections, except that information from a frequency domain eigendecomposition, rather than a correlation-domain eigendecomposition, is used to reconstruct the detected signal. As in the correlation-domain case, a priori assumptions regarding the domain boundaries must be invoked.

The reconstructed signal corresponding to a peak centered
at frequency *f _{0}* can be obtained as

,

where the envelope function *A*(*t*) is determined from
a time-domain inversion of the spectral domain information
contained in the *K* complex eigenspectra.
The envelope *A*(*t*) has *K* complex degrees of freedom,
and allows for phase and amplitude variations in the
time reconstruction of the signal centered
at frequency *f _{0}*.
The discrete time sequence describing
the complex envelope

Magnitude-squared coherence* C _{XY }*is a function of frequency with values between 0 and 1, and is a fraction of a common variance between two time series x and y through a linear relation. It is defined by using the individual power spectra estimates (

In practice,*C _{XY} *has to be computed by using some kind of time-frequency ensemble averaging. Multitaper method with

The confidence limits of estimated coherence at a given level *&alpha** (95% is &alpha*=0.05) can be calculated as:

Cross-power spectrum* P _{XY}* can be obtained with Blackman-Tukey method. Multi-channel SSA is an advanced, data-adaptive method to analyze oscillatory spatio-temporal modes in multivariate time series. In addition to identifying oscillatory peaks in the cross-spectrum, MSSA allows reconstruction of the multivariate oscillatory modes.

1. Mann, M.E. and Lees, J.M., 1996: *Robust estimation of background
noise and signal detection in climatic time series*,
Clim. Change., **33**, 409-445.

2. Percival, D.B., and Walden, A.T.,
1993: *Spectral analysis for physical applications--Multitaper and
conventional univariate techniques*. Cambridge University, 580 pp.

3. Park, J., 1992: *Envelope estimation for quasi-periodic geophysical
signals in noise: A multitaper approach, in Statistics in
the Environmental and Earth Sciences*, edited by A.T.
Walden and P. Guttorp, 189-219, Edward Arnold, London.

4. Thomson, D.J., 1982: *Spectrum estimation and harmonic
analysis*, Proc. IEEE, **70**, 1055-1096.