The Multi-Taper 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 (Ghil and Vautard, 1991; Kuo et al, 1990; Mann and Park, 1993;1994;1996b; Lall and Mann, 1995; Mann et al, 1995a; Thomson 1995; Mann and Park, 1996a) paleoclimate data (Berger et al, 1991; Chappellaz et al, 1990; Mann et al, 1995b; Mann and Lees, 1996; Mommersteeg et al, 1995; Park and Maasch, 1993; Thomson, 1990ab; Yiou et al, 1991;1994;1995;1995) geochemical tracer data (Koch and Mann, 1995), and seismological data (Park et al, 1987; Lees, 1995). Time-frequency ``evolutive'' analyses based on moving window adaptations of MTM have also been applied (e.g., Yiou et al, 1991; Birchfield and Ghil, 1993; Mann et al, 1995b; Mann and Park, 1996b).

MTM, like the method of Blackman and Tukey (1958),
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 (Thomson, 1982, Percival and Walden,1993)
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) and defined as the
eigenvectors of a suitable Rayleigh-Ritz minimization problem, were extensively studied
by Slepian (1978).
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 (Thomson, 1990b).
The tapers are the discrete set of eigenfunctions which solve
the variational problem of minimizing leakage outside
of a frequency band of
half bandwidth where is
the Rayleigh frequency, and *p* is an integer. Note that the
case *p*=1 reduces to the trivial Blakman-Tukey case of a single
tapered Discrete Fourier Transform (DFT).
Because the windowing
functions or eigentapers are the specific solution to an
appropriate variational problem, this method is less
heuristic than traditional nonparametric techniques
(see Box and Jenkins, 1970; Jenkins and Watts, 1968).

Detailed algorithms for the calculation of the eigentapers
are readily available (e.g., Thomson, 1990b; Percival
and Walden, 1993)
In practice, only
the first 2*p* - 1 tapers provide usefully small spectral-leakage
(Slepian, 1978; Thomson, 1982; Park et al, 1987).
Thus, the number of tapers used *K* should
be less than 2*p*-1 in any application of MTM.
The choice of the bandwidth
and number of tapers *K* thus represents the classical
tradeoff between spectral resolution and the stability or ``variance''
properties of the spectral estimate (Thomson, 1982).
The trivial case *K*=1 and *p*=1 is, again, the single-tapered
DFT of Blackman and Tukey.
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
(see e.g., Mann and Park, 1993).
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.

We show in Figure 1 the *K*=3 usefully leakage-resistant
Slepian tapers for
bandwidth parameter *p*=2. In this example,
the effective half bandwidth spectral
resolution is for the spectral estimate centered
at any particular frequency , twice that of an unsmoothed
DFT (a spectral estimate at *f*=0.2 cycle/yr
for a *N*=588 month (49 year) series for example,
averages over the frequency
interval *f*=0.16 to *f*=0.24 cycle/yr).
However, the variance of the spectral
estimate at is decreased threefold by averaging the
independent eigenspectra-a clearly superior variance/resolution
tradeoff to conventional periodogram smoothing. The process
of MTM spectral estimation is described in more detail in the
subsequent section.

**Figure 1:** *The first K=3 eigentapers or `DPSS' for
the case p=2,
computed for the N=588 available monthly values of the SOI series.
Note that the 3rd eigentaper in the sequence has the least
optimal leakage-resistance properties,
with substantially finite
amplitude at the interval boundaries. The eigentapers shown
have been normalized
for unit power.*

MTM can provide estimates of both the singular components (ie,
the ``lines'') and the continuous component (ie, the
background) of the spectrum.
Once the tapers are computed for a chosen frequency
bandwidth, the total power spectrum can
be estimated by averaging the individual spectra given by each
tapered version of the data set. We
call the
*k*th eigenspectrum, where where
is the DFT of . The
*high-resolution* multitaper spectrum is a simple weighted
sum of the *K* eigenspectra,

Its resolution in frequency is , which means that ``line'' components will
be detected as peaks or bumps of width . For a white
noise or ``locally white noise'' process (see Thomson, 1982)
the high-resolution
spectrum is chi-squared distributed
with 2*K* degrees of freedom.

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 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*K* of
the high-resolution multitaper spectrum.

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. Other approaches are
described by Macdonald (1989).

Assume the time series *X*(*t*) is the sum of a sinusoid of
frequency
and amplitude , plus a ``noise'' which is
the sum of other sinusoids and white noise. One can then write

If are the first *K*
eigentapers and the
DFT of , a least-square fit in the frequency domain yields an
estimate of the amplitude :

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* (Kendall and Stuart, 1979).
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
-- 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
(see Thomson, 1982).

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 does not depend -- to first order -- on the magnitude of
. This feature is an important advantage of MTM over
standard methods where
error bars are essentially
proportional to the amplitude of a peak (e.g., Jenkins and
Wattes, 1968).
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 (c.f. the discussion by Park, 1992) 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 (Mann and Lees, 1996). 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 (e.g., Hasselmann, 1976). 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 (see the discussion in Mann and Lees, 1996).

The power spectrum of the AR(1) process is given by (e.g., Bartlett, 1966),

for frequency *f* where is the average value of the power spectrum,
related
to the white-noise variance by,

, the Nyquist frequency, is the highest frequency that can be resolved for sampling rate . The characteristic noise decay time scale can be estimated

For periodicities much larger than , 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 ) 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 as a compromise between describing the full variation of the background spectrum over the Nyquist interval, and insensitivity to narrow spectral features.

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 (see Mann and Lees, 1996). A reshaped spectrum is determined in which the contributions from harmonic signals are removed (see Thomson, 1982), 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.

To illustrate the revised MTM procedure, we apply the
approach to the SOI series discussed earlier.
Consistent with the results from Singular Spectrum Analysis
of this series discussed earlier,
the MTM analysis (Figure 2)
recognizes
two highly signficant interannual peaks, one centered at
*f*=0.2 cycle/yr (roughly 5 year period) and another centered
at *f*=0.39 cycle/year (roughly 2.5 year period) which we associate
with the low-frequency (LF) band and high-frequency (HF) band
ENSO signals, respectively. The signals are significant at
well above the 99% level, and are unlikely to have arisen from
chance coincidence; for *p*=2 and *N*=588 monthly samples,
there are fewer than 50 statistically independent bandwidths in
the spectral estimate within
the subannual (*f*;*SPMlt*;0.5 cycle/yr) band, and not even 1 coincident
interannual or lower frequency peak is expected at the 99% level.
A low frequency variation inseparable
from trend is also isolated as significant at the 99% level
relative to the estimated red noise backround. It is interesting
to note that the AR(1) red noise model does *not* provide
a good description of the noise background
in a neighborhood of the spectrum surrounding *f*=1.0 cycle/yr-more
than anything else, this discrepancy appears to be associated with the
severe impact of the deseasonalization procedure on
the spectral character of the series (see also Thomson, 1995)
implicit in the construction
of the standardized SOI index.
Away from *f*=1.0 cycle/yr, the robustly
estimated red noise background provides a good visual fit to the
background spectrum. It is assuring that the robust noise background
estimation procedure
is entirely insensitive to the clearly anomalous behavior
of the spectrum near the annual cycle.
The prominent high-frequency peak
near *f*=5.4 *cycles*/*yr* (approximately 65 day period),
which may be associated with intraseasonal
oscillatory behavior of the tropical atmosphere.

**Figure 2:** *Adaptively weighted MTM spectrum of the SOI time series.
The robust red noise background fit and associated 90significance levels are show by 4 smooth curves. Three significant
low-frequency signals-the trend, and LF and HF band interannual signals-are
detected. Several higher frequency peaks are also significant.
The bandwidth parameter is p=2 and K=3 tapers were used.*

To better understand the distinction between the inferences from
the conventional MTM procedure and the
revised procedure of Mann and Lees (1996) we compare the
results of the harmonic peak detection approach and the
revised approach (Figure 3). The F-test criterion for
harmonic signals yields 7 peaks at the 99% confidence
level, and 31 peaks at the 95% confidence level. Many of
those peaks are associated with weak power in the spectrum.
Note that the harmonic peak test has the higher Raleigh
frequency resolution , and not the broader bandwidth
of the
adaptive multitaper spectral estimate.
In a series
of *N*=588 samples, we would expect about 6 and 30 peaks at
the 99% and 95% confidence levels respectively from
chance coincidence alone, which is indistinguishable from
the results of the F-test for the SOI above.
Thus, it is difficult to
distinguish any significant features from
random noise, based on the
harmonic peak detection approach alone. However,
the results of the harmonic peak detection are nonetheless
useful in the reshaping procedure used by Mann and Lees (1996).
Figure 3 (bottom) also
shows the reshaped and adaptively weighted multitaper
spectra along with the signficance levels relative to the robustly
estimated red noise background. Dashed peaks indicate peaks
that satisfy the harmonic detection test at level
and are significant relative to red noise at level.
There are 6 such peaks, including both an LF and HF ENSO band
peak. The inference in this case is that the LF and HF ENSO
band signals are significant well above the 99noise, and are furthermore approximated as harmonic, phase-coherent
oscillations at a 95% level of confidence.
If a 99% F-test threshold were required in the reshaping
procedure, neither of those peaks are indicated
as passing the harmonic test. Nonetheless, they are still highly significant
narrowband features of the spectrum.
The reader is thus warned that, while the silmultaneous fullfullment
of both the harmonic peak and red noise significance tests may
seem to provide an even more robust significance criterion,
the harmonic peak test
is an imperfect detection tool for many of the kinds of signals
encountered in geophysical systems. For the SOI, we have
shown that only the revised test for narrowband features
significant relative to a properly isolated red noise
background yields
robust evidence at the 99% level for interannual timescales signals.
Such signals
are known to exist for physical, as well as statistical reasons,
and provide a good benchmark for comparison of spectral analysis
methods.
The true ENSO signal is almost certainly associated with amplitude and phase
modulation over time-features not appropriately modeled
in the harmonic peak test alone.

**Figure 3:** *Harmonic peak (F variance-ratio) test with
median, 90reshaped vs. unreshaped adaptively weighted MTM spectrum (below)
based on p=2 and K=3, and a 95% F-test signficance criterion
for reshaping. The solid spectrum indicates the reshaped spectrum,
while the dashed spectrum represents the discrepancy between
the reshaped and unreshaped spectra, providing a measure of the
portion of the spectrum associated with harmonic features in the
spectrum. Signficance levels for the spectrum are shown as in
Figure 2.*

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 is written

or, for the discrete case at hand,

where the envelope function *A*(*t*) is determined from
a time-domain inversion of the spectral domain information
contained in the *K* complex eigenspectra (see Park 1992;
Park and Maasch, 1993; Mann and Park, 1994; 1996b).
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
is determined from a discrete
inverse problem using the complex amplitudes
of each of the the *K* eigenspectra
and appropriate boundary conditions (see Park, 1992;
Park and Maasch, 1993).
The three lowest order boundary contraints in this
inversion involve minimization of the envelope
itself near the boundaries, numerical
mimimization of the slope of the envelope
near the boundaries, or optimization of the smoothness
of the envelope (ie, numerical minimization of the 2nd derivative
near the boundaries).
Any one of these choices might be favoured
if a priori information regarding the characteristics of
the signal is available (see the discussion of Park, 1992).
The amplitude of the seasonal cycle
in surface temperature, for example, is nearly constant in time,
and a minimum slope constraint is most
appropriate for its reconstruction (see Mann and Park, 1996a).
Generally, however, a nearly optimal reconstruction
can always be obtained
through seeking
the weighted linear combination of these 3 contraints that minimizes
the mean-square misfit of the reconstructed signal with
respect to the raw data series (Mann and Park, 1996b).
We favour this method of time-domain inversion.

In Figure 4 we show the reconstructed components (RCs) of the
southern oscillation index (SOI) corresponding to the three
low-frequency signals that are detected as significant at
level in Figure 2-the
trend and low-frequency (LF) and high-frequency (HF)
interannual components.
The LF and HF interannual signals each exhibit roughly
2 unit peak-to-peak variation. The amplitude modulations
of both oscillatory components are similar, exhibiting
weak amplitude in the middle of the 50 year data
interval.
The long-term trend in the SOI is rarely discussed
and, while weak, nonetheless explains an
important
roughly 0.5 unit peak-to-peak
variation, about 25% of the amplitude of the peak-to-peak
variations associated with the LF and HF band signals.
A low-frequency peak near *f*=0.6 cycle/yr significant at the 95%
level was not reconstructed, but may explain a significant
quasibiennial component of ENSO.
It is apparent that the model of two
amplitude-modulated quasi-oscillatory components and
a low-frequency trend provided by the MTM decomposition
provides a parsimonious description of a relatively
large (24the series.
Figure 5 shows the sum over the 3 RCs shown in Figure 4
against the raw monthly SOI series, providing a similar
description of the low-frequency signal to that of the
SSA reconstruction described
earlier.

References

Table of Contents

This file last modified: Mar. 4, 1998