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 2p - 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 2p-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 kth 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 2K 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 2K 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 X(t), using K eigentapers, to the residual variance. By expanding the variance of the model one finds that it is the sum of two terms,
that are respectively the ``explained'' and ``unexplained'' contributions to the variance.
The random variable
would obey a Fisher-Snedecor law with 2 and 2K-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 f0 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 f0. 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.