Chapter 4: Multitaper Spectral Analysis of High-frequency Seismograms
ABSTRACT
Spectral estimation procedures which employ several prolate spheroidal sequences as tapers have been shown to yield better results than standard single-taper spectral analysis when used on a variety of engineering data. We apply the adaptive multitaper spectral estimation method of Thomson (1982) to a number of high-resolution digital seismic records and compare the results to those obtained using standard single-taper spectral estimates. Single-taper smoothed-spectrum estimates are plagued by a trade-off between the variance of the estimate and the bias caused by spectral leakage. Applying a taper to reduce bias discards data, increasing the variance of the estimate. Using a taper also unevenly samples the record. Throwing out data from the ends of the record can result in a spectral estimate which does not adequately represent the character of the spectrum of nonstationary processes like seismic waveforms. For example, a discrete Fourier transform of an untapered record (i.e., using a boxcar taper) produces a reasonable spectral estimate of the large-amplitude portion of the seismic source spectrum but cannot be trusted to provide a good estimate of the high-frequency roll-off. A discrete Fourier transform of the record multiplied by a more severe taper (like the Hann taper) which is resistant to spectral leakage leads to a reliable estimate of high-frequency spectral roll-off, but this estimate weights the analyzed data unequally. Therefore single-taper estimators which are less affected by leakage not only have increased variance but also can misrepresent the spectra of nonstationary data. The adaptive multitaper algorithm automatically adjusts between these extremes. We demonstrate its advantages using 16-bit seismic data recorded by instruments in the Anza Telemetered Seismic Network. We also present an analysis demonstrating the superiority of the multitaper algorithm in providing low-variance spectral estimates with good leakage resistance which do not overemphasize the central portion of the record.
1. Introduction
Spectral estimation is a powerful method of data analysis which is often used to study geophysical processes. The estimation of the spectra of background noise, line components, and transient signals is central to the analysis of electric, magnetic, and seismic time series. There have been many techniques developed which are effective for the analysis of long records of stationary processes. Unfortunately, these techniques are not universally applicable to seismic data sets. In many studies it is necessary to estimate a spectrum from a short time series. This situation can occur if some of the data are missing or if the data of interest (e.g., a seismic phase) are contained in a short segment of a longer record.
A new approach for estimating the spectra of short time series, known as multitaper spectral analysis, has been developed recently. We have applied this technique, which was first presented by Thomson [1982], to several dozen seismograms; in this paper we analyze two representative records. The spectra estimated using the multitaper technique are compared with several direct spectral estimates employing commonly used single tapers. We show that the multitaper approach can yield superior results when applied to high-frequency seismic data.
Spectral analysis of specific phases within a seismogram, particularly those at regional or local distances, can be difficult. It is often impossible to isolate a particular phase. If one isolates a major phase by discarding the rest of the record and then makes a direct estimate of the waveform's spectrum without first tapering the data (i.e., using a boxcar taper), the high-frequency roll-off of the estimated spectrum can be severely biased by spectral leakage. Therefore it is standard practice to multiply the time series by a taper before performing a discrete Fourier transform (DFT) to reduce spectral leakage (an extensive review of tapering is provided by Harris [1978]).
The cosine or Hann taper is popular in seismic analysis, being both effective and easy to calculate. The utility of the Hann taper is bought dearly, however. If one views each data point in a time series as a constraint on the estimated frequency content of the record, the Hann taper discards 5/8 of the statistical information in a given time series. This can be easily seen from the graph of the Hann taper in Figure 1. The data points at the extremes of the record are weighted weakly, while the center of the time series is emphasized. This unequal weighting causes the statistical variance of a direct spectral estimate using a Hann taper to be greater than the variance of a periodogram spectral estimate.
Blackman-Tukey tapers are designed to circumvent this loss of information somewhat by applying the cosine weighting to only the extremes of the record. For instance, the 20% cosine taper (Figure 1) discards only 12.5% of the available data variance constraints. However, Harris [1978] shows that this taper has less resistance to spectral leakage than a Hann taper. As long as only a single data taper is used, there will be a trade-off between the resistance to spectral leakage and the variance of a spectral estimate.
Thomson [1982] introduced the multitaper spectral analysis technique. First, the data are multiplied by not one, but several leakage-resistant tapers. This yields several tapered time series from one record. Taking the DFTs of each of these time series, several "eigenspectra" are produced which are combined to form a single spectral estimate.
The tapers are constructed so that each taper samples the time series in a different manner while optimizing resistance to spectral leakage. The statistical information discarded by the first taper is partially recovered by the second taper, the information discarded by the first two tapers is partially retrieved by the third taper, and so on. Only a few low-order tapers are employed, as the higher-order tapers allow an unacceptable level of spectral leakage. One can use these tapers to produce an estimate that is not hampered by the trade-off between leakage and variance that plagues single-taper estimates, as we will demonstrate.
Single-taper spectral estimates have relatively large variance (increasing as a larger fraction of the data is discarded and the bias of the estimate is reduced) and are inconsistent estimates (i.e., the variance of the estimate does not drop as one increases the number of data). To counteract this, it is conventional to smooth the single-taper spectral estimate by applying a moving average to the estimate. This reduces the variance of the estimate but results in a short-range loss of frequency resolution and therefore an increase in the bias of the estimate.
The multitaper spectral estimates are formed as a weighted sum of the eigenspectra. Therefore the multitaper spectral estimate is already a smooth estimate; it has less variance than single taper spectral estimates which have been designed to reduce bias, and it is also a consistent estimator. The comparison between the bias and variance properties of the single taper and multitaper estimates is discussed further in sections 3 and 4.
Another difficulty with seismic data is that the records are nonstationary; that is, the statistical character of the data changes with position in the record. Therefore a spectral estimator which weights the data in the center of the time series more heavily than data at the ends can overemphasize the signal energy in the middle of the record. This can result in a misrepresentation of the spectrum, as we demonstrate in section 3. The multitaper estimate, which discards very little data from the record and weights the data relatively evenly, is not subject to this problem.
Section 2 presents an outline of the basic algorithm. This outline contains sufficient detail to allow the reader to implement the algorithm but avoids derivations that can be found elsewhere. Section 3 describes the seismic data used in this study and presents comparisons of the spectra of seismic time series generated by both regional and local events. We demonstrate the trade-off between spectral leakage resistance and variance of the spectral estimates produced using the boxcar, 20% cosine, and cosine tapers. We compare the bias and variance of these conventional single-taper direct spectral estimates with the bias and variance of the multitaper spectral estimates in section 4. A numerical method for calculating the prolate eigentapers is given in the appendix.
2. The Multitaper Algorithm
The multitaper method is based on a family of tapers which are resistant to spectral leakage. We outline the multitaper method here, and note that more detailed treatments can be found in the works by Thomson [1982] and Lindberg [1986].
Suppose that we are given the finite time series { xt ; t = 0, 1, ..., N -1 } which is a set of discrete samples of a continuous time process { xt ; 0 ≤ t ≤ N -1 } (we assume a unit sample interval &Tau = 1, without loss of generality). If xt is harmonizable (a wide class of processes including stationary and several forms of non-stationarity) then it has the Cramer spectral representation [Doob, 1953]:
x sub t ~=~ int from {- half} to {{~^}half} ~ e sup { i2 pi ft} X(f)df
We wish to estimate the power spectrum S (f ) = E { | X (f ) |2 } (where E { } denotes expected value) of the continuous time process { xt ; 0 ≤ t ≤ N -1 } from the time series { xt }t =0N -1. A conventional direct spectral estimate | Xa (f )|2 of S (f ) is found by multiplying the data { xt }t =0N -1 by a sequence { at }t =0N -1 called a taper, applying a DFT,
X hihat sub a (f) ~=~ sum from t=0 to N-1 ~ a sub t ^ x sub t ^ e sup {-i2 pi ft}
and finally taking the squared modulus of the resulting function X a (f ). Although t is discrete, f is continuous, with | f | ≤ 1/2 (as the Nyquist frequency f Nyquist = 1/2). We normalize the taper so that
sum from t=0 to N-1 ~ a sub t sup 2 ~=~ 1
The spectral leakage properties of the data taper at , t = 0, 1, 2, ..., N -1 can be deduced from its DFT:
A (f) ~=~ sum from t=0 to N-1 ~ a sub t e sup {-i2 pi ft}
For conventional tapers the function | A | has a broad main lobe and a succession of smaller side lobes. For example, for the boxcar taper, at ≡ 1/ sqrt N and
A(f) ~=~ 1 over {sqrt N}^ {1-e sup {-i2 pi fN}} over {1-e sup {-i2 pi f}} ~=~ {e sup {-i(N-1)^ pi f}} over {sqrt N} ^ left ( { sin^ N pi f} over {sin^ pi f } right )
In this case the function | A (f ) | is readily observed to have a central lobe flanked by smaller side peaks. (The phase factor e -i (N -1) π f results from choosing the time origin t = 0 to coincide with the first data point. It does not affect the leakage resistance of the taper.) The larger the side lobes, the more spectral leakage is encountered which biases the estimate X a away from its desired value. This can be seen by observing that
X hihat sub a (f) ~=~ int from {- half} to {{~^}half} ~ A (f-f prime ) X( f prime ) df prime
which follows from substituting the Cramer spectral representation of xt in the definition of X a , and therefore
| X hihat sub a (f) | sup 2 ~=~ int from {- half} to {{~^} half} ~ |A(f-f prime ) | sup 2 S(f prime ) df prime
A good data taper should have a spectral window A ( f - f′ ) whose amplitude is large in the central lobe region where | f - f′ | is small and has low side lobes at more distant frequencies. This reduces the bias in the estimate by preventing the energy in X at distant frequencies from leaking over to affect the estimate | X a |2 at frequency f .
Suppose we wish to minimize the bias at frequency f due to spectral leakage from outside the frequency band | f′ - f | ≤ W, where 2W is some chosen bandwidth. We maximize the fraction of energy of A within the chosen band:
λ (N, W ) ~=~ {int from {-W} to {{~^}W} ~ |A(f)| sup 2 df } over {int from {- half} to {{~^}half} ~ |A(f)| sup 2 df}
Since no finite time series can be completely band-limited, λ (N, W ) < 1 for finite N and nontrivial W . The functional λ can be interpreted as follows: in a single-taper direct estimate of the spectrum of a white noise process at frequency f , λ is the fraction of spectral energy in that estimate that derives from the frequency interval | f - f′ | ≤ W ; 1 - λ is the fraction of spectral energy that leaks in from outside that band.
It is convenient from a computational viewpoint to substitute (4.1) into (4.3) to express λ in terms of the data tapers themselves rather than their transforms. If we seek a taper for an N point time series, the sequence { at }t =0N -1 can be represented as an N vector a . This notation allows us to express (4.3) in matrix form (following the derivation of equation (2.5) of Park et al. [1987], letting the decay rate α = 0)
lambda (N,W ) ~=~ {bold a ~ cdot ~ bold C ~ cdot ~ bold a} over {bold a ~ cdot bold a}
where the matrix C has components
C sub {tt prime} ~=~ {sin^ [2 pi W (t-t prime )]} over {pi (t-t prime )} ^; ~~~~~~~ t, t prime ~=~ 0,^1,^2,^...,^ N-1
We seek those values of a for which the functional λ is stationary. This leads to the matrix eigenvalue problem
C . a - λ (N , W )a = 0
which has as its solutions the ordered eigenvalues 1 >~λ0 > λ1 > λ 2 > ... > λ N -1 ~>~ 0 and associated eigenvectors v(k )( N , W ) ; k = 0, 1, 2, ..., N - 1 (which have components v t(k ); t = 0, 1, 2, ..., N -1 ). The v(k )( N , W ) are discrete prolate spheroidal sequences [Slepian, 1978], which we also refer to as prolate eigentapers. We will suppress the explicit dependence of v(k ) on N and W in the following. A prolate eigentaper with a time-bandwidth product of P = N W is called a P π prolate taper; it concentrates spectral energy in frequency bands of width 2W = 2 P / N . As the Rayleigh frequency 1/N is the fast Fourier transform (FFT) frequency bin spacing, a P π prolate taper will have a main lobe which is 2 P "frequency bins" wide. For instance, tapers for which NW = 4 minimize the spectral leakage at frequency f from outside the frequency band defined by | f′ - f | ≤ 4/N. For large N (> 100) one can construct a set of the v(k ) for any value of the time-bandwidth product N W . As noted in the appendix, this allows the user to calculate one set of eigentapers v(k ) for a fixed value of N and to interpolate this set to construct tapers for time series of various lengths. We have restricted the following discussion to 4π prolate tapers, but similar behavior is found for other choices of the time-bandwidth product.
The five lowest-order eigentapers v(k ), k = 0, 1, 2, 3, 4 shown in Figure 2 have been made for a time series of length N = 128 and time-bandwidth product N W = 4. The lowest-order taper (k = 0) is the familiar 4π prolate taper advocated by Thomson [1971, 1977a,b] and Eberhard [1973] and has a shape similar to conventional tapers such as the cosine taper (Figure 1). The higher-order eigentapers are markedly different from ordinary data tapers. For even values of k , the v(k ) are symmetric about the midpoint of the time series. For odd values of k , the v(k ) are antisymmetric about the midpoint. All the tapers, except the lowest-order one, have regions of positive and negative data weighting. We normalize the tapers so that
sum from t=0 to N-1 ~ (v sub t sup (k) ) sup 2 ~=~ 1
As the eigentapers v(k ) are solutions to (4.5), they are orthogonal:
v(k ) . v(k ') = δ k k′
(This can be clearly seen in Figure 2.) This relation shows that each v(k ) can be used to provide an orthogonal sample of the data { xt }t =0N -1 .
Taking discrete Fourier transforms of the prolate eigentapers produces the spectral windows
U sub k (N,W ; f) = {epsilon sub k} ~ sum from t=0 to N-1 ~ v sub t sup (k) (N,W ) e sup {-i2 pi f(t-(N-1)/2)}
where we have used the time-centered transform to eliminate spurious phase factors in the definition. The function Ε k = 1 if k is even; Ε k = i if k is odd. The use of Ε k = is a notational convention so that Uk is real-valued. Plots of the Uk for N =128 and N W = 4 appear in Figure 3 for k = 0, 1, ..., 4. Most of the energy of the Uk is concentrated within the specified frequency band as was required by maximizing (4.3). The spectral windows corresponding to the lowest-order eigentapers have impressively small side lobes, but spectral leakage resistance becomes progressively poorer as the order of the taper increases. The lowest-order 2NW eigentapers (e.g., the eight lowest-order 4π prolate tapers) have eigenvalues λ k close enough to unity that they are useful for minimizing spectral leakage.
The eigenvalues λk of the eight lowest-order eigentapers with time-bandwidth products 4, 3, and 2 are given in Table 1 for N = 128. For reference, the value of the functional (4.3) is given for a boxcar taper which concentrates spectral energy within frequency bands of the same width.
To construct a multitaper spectral estimate, one first calculates the complex "eigencoefficients" yk (f ) by taking a DFT of the product of the data with each {vt(k)}t=0N-1
y sub k (f) ~=~ sum from t=0 to N-1 ~ v sub t sup (k) x sub t e sup {-^i2 pi ft}
An estimate of the spectrum can be made from weighted sums of the eigenspectra |yk|2. Thomson [1982] formulates the problem of estimating the spectrum of a record as an integral equation. The solution of the integral equation is averaged over (f - W, f + W) to produce the smoothed high-resolution spectral estimate:
S hibar (f) ~=~ K sup -1 ~ sum from k=0 to K-1 ~ ( lambda sub k ) sup -1 | y sub k (f) | sup 2
where K is the number of tapers used. If K is not large, the smoothed high-resolution estimate (4.9) differs little from an arithmetic average of the eigenspectra as λ k ~ 1 for the lowest-order eigentapers.
Although straightforward, (4.9) is not the best multitaper spectral estimate to use. An adaptive spectral estimate
S hihat (f) ~=~ {sum from k=0 to K-1 ~ |d sub k (f) y sub k (f) | sup 2} over {sum from k=0 to K-1 ~ |d sub k (f) | sup 2}
can be devised which has frequency-dependent weights dk (f ) chosen to reduce bias from spectral leakage [Thomson, 1982]. This technique proves extremely useful in the analysis of highly-colored spectral processes. At frequencies f where the spectrum is reasonably flat, the weights dk (f ) ~ 1, reducing the variance of the spectral estimates. At frequencies f where spectrum has a steep slope, the contribution from the higher-order eigentapers, which have poorer leakage resistance, is reduced. The trade-off between spectral leakage and variance of the spectral estimate is balanced at each frequency.
The optimal weights dk can be found by minimizing the misfit of the estimated spectrum to the true spectrum S (f ). This misfit, although unknown, can be estimated statistically. The resulting equation for the weight dk (f ) is
d sub k (f) ~=~ {sqrt {lambda sub k} ^ S (f)} over {lambda sub k S(f) + E "{" B sub k (f) "}" }
where S (f ) is the true value of the spectrum at frequency f and Bk(f ) is the spectral energy at frequency f that leaks in from outside the frequency band (f - W, f + W ). We replace the unknown value S (f ) by its estimate S (f ). Thomson [1982] found it adequate to approximate E { Bk(f)} ~ σ 2 (1- λ k), i.e., as a constant fraction of the total variance of the time series:
sigma sup 2 ~=~ sum from t=0 to N-1 ~ x sub t sup 2
We find the estimate S (f ) by iteration. We take the arithmetic average of |y0(f )|2 and |y1(f )|2 as an initial estimate of S (f ), then substitute this value into (4.11) to produce first guesses of the weights dk (f ). These weights are then used in (4.10) to generate a new spectral estimate S (f ), and the process is repeated. Convergence is usually satisfactory within a few cycles.
Careful examination of the adaptive spectral estimate shows that Parseval's theorem is not explicitly satisfied, i.e., there is no requirement that the energy of the spectrum estimate, integrated over frequency, equal the total variance of the time series. This arises from the way that the multitaper algorithm attempts to compensate for the effects of spectral leakage. If the expected broadband bias E{Bk(f)} were to vanish, then (4.11) would become dk (f ) = λ k-1/2, and the adaptive estimate (4.10) would reduce to the smoothed high-resolution estimate (4.9) (except for a small multiplicative factor due to the departure of the eigenvalues λ k from unity). This would occur if the true spectrum were zero outside the frequency band | f -f′ | ≤ W. As (1- λ k) of the process variance within the frequency interval | f -f′ | < W is leaked outside the band, the limiting case dk (f ) = λ k-1/2 represents an attempt by the estimator to compensate for this spectral energy lost to leakage by boosting the coefficients of the higher-order eigenspectra in the weighted sum. When the spectrum has a steep slope, the higher-order eigenspectra are downweighted and the adaptive spectral estimate tends toward the least biased eigenspectral estimate | y0(f) | 2.
Thomson [1982] analyzed two synthetic time series using multitaper methods. Both series had fewer than 100 data points and a numerical precision of roughly 20 bits. In the first example, it was demonstrated that a multitaper approach could accurately estimate a spectrum with a dynamic range of more than seven decades and accurately infer the existence of harmonic lines (i.e., coherent sinusoids) in the data. Thomson also analyzed a 64-point time series used by Kay and Marple [1981] in a spectrum analysis "shootout" comparing 11 spectral estimation techniques, including the maximum entropy method as well as a single-taper direct spectral estimate and several other popular spectral estimates. Unlike any of the techniques tested by Kay and Marple [1981], a multitaper technique was able to produce a spectral estimate which was similar to the true spectrum of the synthetic time series.
3. Spectral Comparisons Using Siesmic Data
We compare a number of single-taper direct spectral estimates with the adaptive multitaper spectral estimate techniques on wide dynamic range, high-resolution seismic data. The advent of digital arrays with 16-bit data loggers and the proposed 22- or 24-bit precision instruments demand an improved sophistication in data analysis techniques. We may soon have seismic data which are recorded to the same precision as the synthetic examples of Thomson [1982].
The data used in this paper were recorded on seismometers in the Anza Seismic Telemetered Array. The Anza array was designed to record high-frequency seismic signals from local earthquakes. The instruments in this array measure surface velocity, and the data are recorded as 16-bit numbers (this allows a dynamic range of 96 dB). See Chapter II for a more detailed description of the Anza array.
The multitaper spectral estimate has a smaller variance at each frequency than a single-taper direct spectral estimate. To make a fair comparison between the various direct spectral estimates and the adaptive multitaper method, we will smooth each single-taper estimate using a moving average so that each estimate averages information over roughly the same frequency band as a multitaper estimate using seven 4 π prolate eigentapers.
The effect of smoothing single-taper direct spectral estimates in this way is shown in Figure 4. The section of the seismogram which is analyzed is shown between the vertical dashed lines at the top of Figure 4. The unsmoothed spectral estimates are shown below on the left, and the smoothed estimates are displayed on the right below the record. The upper traces are direct estimates using the Hann taper, the middle traces are spectral estimates made with a 20% cosine, and the lower traces are spectrum estimates which employ a boxcar taper. The amplitude is plotted on a logarithmic scale on the vertical axis, and frequency is plotted on a linear scale on the horizontal axis. Each trace is offset by a multiplicative factor of 50 from the adjacent traces. Notice that if one studies the unsmoothed spectral estimates, it is difficult to distinguish any specific features common to each of the estimates except for a general linear trend. In comparison, the smoothed spectral estimates have many of the same features. Each major peak or trough appears at the same frequency in each of the smoothed estimates.
Unfortunately, since we are using real data, it is impossible to know the true spectrum for any of the examples. However, the work of Thomson [1982] demonstrates that the multitaper method provides a reasonable spectral estimate. This is confirmed by a study comparing the multitaper estimate with the smoothed direct estimates on a synthetic seismic wave train with a known spectrum (C. Lindberg et al., unpublished manuscript, 1989).
To study how tapering affects the spectra of body wave pulses, we isolate a phase in the middle of a seismogram, produce spectral estimates using each of the four methods, and compare the results. The upper graph in Figure 5 shows the transverse horizontal seismogram of an earthquake which had an epicentral distance of 100 km from the recording station PFO (in Pinyon Flat, California). We extract that section of the seismogram corresponding to the shear wave arrival (between vertical dashed lines) and estimate its amplitude spectrum by each method. The spectral estimates are plotted on a linear-linear scale in the lower portion of Figure 5 and for clarity are plotted in dimensionless velocity units on the vertical axis. Each of the four spectral estimates have two main peaks in the frequency band from 0 to 20 Hz, near 4 and 14 Hz.
These estimates are interesting to compare. Three of the estimated spectra (those plotted using solid and dashed lines) have almost identical features (except for the offset between them). In these spectral estimates the amplitude of the peak at 14 Hz is about 20% less than the amplitude of the peak at 4 Hz. The other estimated spectrum (curve d, plotted with asterisks) does not resemble the other three estimates closely. The peak at 14 Hz is 10% higher than any other peak in this estimate. This change in the relative amplitude of the two spectral peaks would influence the choice of a corner frequency if these spectra were converted from velocity to displacement or acceleration.
The three spectral estimates which exhibit similar characteristics are the multitaper estimate (curve a, plotted as a solid line), the 20% cosine direct estimate (curve b, the upper dashed line), and the boxcar direct estimate (curve c, the lower dashed line). The spectrum showing a different distribution of spectral energy was estimated using the Hann taper (curve d). The Hann direct spectral estimate is unlike the other three estimates because it imposes a different emphasis on the time series. Referring back to Figure 1, it is easy to see that the boxcar applies equal weighting to the entire time series and the 20% cosine taper weights 80% of the series equally. Not surprisingly, using either of these two tapers produces essentially the same result. However, the multitaper spectral estimate also gives essentially equal importance to every data point, like the boxcar and 20% cosine estimates (see Figure 2). The Hann taper puts over 80% of its emphasis on the middle 50% of the time series and gives the data in the first and last 25% of the series less weight. This rejection of data near the ends of the series causes the apparent misrepresentation of the distribution of spectra energy shown in Figure 5.
We also compared estimates of the spectrum of a vertical recording of a nuclear explosion. This event had an epicentral range of 412 km and also was recorded at PFO. The section of data which was analyzed is bounded by the vertical dashed lines in the upper trace in Figure 6. The analysis procedures were identical to those used in the previous example except that the log amplitudes of the spectra were plotted on the vertical axis.
The spectrum of the nuclear test has a large dynamic range and has most of its energy concentrated below 20 Hz. By examining the estimated spectra, one can see that some estimates are more effected by spectral leakage than others. The two estimates which are less subject to spectral leakage, the Hann direct estimate (curve d, plotted with asterisks) and the adaptive multitaper estimate (curve a, the solid line), are very similar. Both of these estimated spectra clearly show the spectrum of the signal from 0 to 20 Hz; from 20 to 60 Hz the spectrum of the ground noise is visible. The antialias filters of the recording system are 6 pole Butterworth filters which have a corner frequency of 62.5 Hz. The effect of the filters is visible in the 60-80 Hz band. In the band from 80 to 125 Hz the ground noise is less than the instrument noise. The variance of the adaptive multitaper spectrum is larger in the low-amplitude portion of the spectrum and hence appears unsmoothed. This is because the downweighting of the higher-order eigenspectra minimizes spectral leakage at the cost of reducing the effective number of degrees of freedom of the estimate at each frequency. If smaller variance is desired in the low-amplitude portion of the adaptive multitaper spectrum, then prolate tapers with a larger time-bandwidth product could be used.
The spectra shown in Figure 6 which were obtained using the 20% cosine and boxcar tapers suffer from the effects of spectral leakage. The spectrum estimate employing the 20% cosine (curve c, the lower dashed line) suffers less from spectral leakage than the estimate utilizing the boxcar, as expected. The leakage of spectrum estimated using a 20% cosine taper hides nearly all the features in the ground noise between 20 and 60 Hz. The effect of the antialias filters is completely obscured. The apparent energy in the 20% cosine spectrum estimate is larger than the instrument noise in the 80-125 Hz band by a factor of 10. The performance of the spectral estimate obtained using a boxcar taper (curve b, the upper dashed line) is even worse, since it does not exhibit any of the features of the true spectrum between 20 and 125 Hz.
These examples show that each of the spectral estimates has different advantages. The smoothed spectrum estimate employing a boxcar taper produces a good estimate of the large-amplitude portions of the spectrum but has very poor spectral leakage properties and is not of much use for spectra which have a large dynamic range. The smoothed spectrum estimate using a Hann taper is less affected by spectral leakage, but this estimate can misrepresent the large-amplitude portion of the spectrum. A smoothed spectral estimate incorporating the 20% cosine taper combines the best properties of the spectral estimates which use the boxcar and the Hann tapers. It retrieves the large-amplitude features almost as well as the boxcar estimate and has spectral leakage properties which are sufficient for many geophysical applications. The adaptive multitaper estimate has even better performance, representing the large-amplitude spectral components as accurately as the boxcar estimate and having excellent spectral leakage properties.
4. Statistical Comparisons
We compare the broadband bias and variance of the smoothed single-taper direct spectral estimates with the smoothed high-resolution and adaptive multitaper estimates. We consider smoothed single-taper and multitaper spectral estimates whose values at some frequency f are formed by averaging seven direct spectral estimates which concentrate the spectral energy at frequency f mainly within the frequency band (f - W, f + W), where W = 4/N (4 times the Rayleigh frequency 1/N ). Therefore we use 4π prolate eigentapers for the multitaper estimates; the seven lowest-order 4π eigentapers have good resistance to spectral leakage (see Table 1), but we do not use the seventh-order 4π prolate eigentaper; it allows excessive spectral leakage, as λ 7 = 0.721233936. We compare the multitaper estimates with a smoothed single-taper estimate Sa(f ) which is formed by averaging the seven direct spectral estimates |Xa(f+j/N)|2; j = -3, -2, ..., 2, 3 obtained using a taper { a(t) }t =0N -1, i.e.,
S hibar sub a (f) ~=~ 1/7 ~ sum from j=-3 to 3 ~ | X hihat sub a (f+j/N)| sup 2
This estimate is mostly an average of spectral energy from the band (f - 4/N, f + 4/N). (The main lobes of tapers other than the boxcar taper are wider than 2/N, so this is not strictly correct, but we will make this approximation.)
4.1. Variance
To gain some idea of how smooth the estimators are, we compare the variances of each spectral estimate for a time series composed of Gaussian white noise. For single-taper estimates we define the covariance matrix
matrix { ccol { LAMBDA sub ij above {~} } ccol { = above = } lcol { 1 over {sigma sup 2} E "{" X hihat sub a sup star (f^+^ i/N ) X hihat sub a (f^+^ j/N ) "}" above int from {- half} to {{~^}half} ~ A sup star (f^+^ i/N ) A (f^+^ j/N ) ^ df } }
for i , j = -3, -2, ..., 2, 3 [see Thomson, 1982, equation 4.1], where σ 2 is the process variance defined in (4.12) and A (f ) is the spectral window introduced in (4.1). If the single-taper estimates Xa (f + i/N ) and Xa (f + j/N ) are uncorrelated for i ≠ j, then Λ is a diagonal matrix. If the amount of correlation of the estimates Xa (f + i/N ) is such that one or more of the Xa (f + i/N ) can be expressed as a linear combination of the others, then Λ is a singular matrix, with at least one zero eigenvalue. In practice, Λ has a behavior which lies somewhere between these extremes.
For white noise data the expected value of Sa(f ) is
E "{" S hibar sub a (f) "}" ~=~ 1 over 7 ~ sigma sup 2 ~ sum from i=-3 to 3 ~ LAMBDA sub ii ~=~ {sigma sup 2} over 7 ~ tr ^ "{" fat LAMBDA "}"
where tr denotes the trace of a matrix, but as
LAMBDA sub ii ~=~ sum from t=0 to N-1 ~ a sup 2 (t) ~=~ 1
for all i ,
sum from i=-3 to 3 ~ LAMBDA sub ii ~=~ 7
and E{Sa(f)} = σ2. The variance of Sa(f) is
E "{" ( S hibar sub a (f)) sup 2 "}" ~-~ (E "{" S hibar sub a (f) "}" ) sup 2
The first term of (4.15) is
E "{" ( S hibar sub a (f)) sup 2 "}" ~=~ 1 over 49 ~ sum from {i,^j^=-3} to 3 ~ E "{" | X hihat sub a (f ^+^ i over N ) | sup 2 | X hihat sub a (f ^+^ j over N ) | sup 2 "}"
As each function | X hihat sub a (f^+^ i/N)| sup 2$ is the sum of squares of two Gaussian random variables, we can show that
E "{" ( S hibar sub a (f) ) sup 2 "}" ~=~ {sigma sup 4} over 49 ~ sum from i=-3 to 3 ~ sum from j=-3 to 3 ~ [ LAMBDA sub ii LAMBDA sub jj ^+^ | LAMBDA sub ij | sup 2 ]
[Papoulis, 1977, chapter 11], so
Var ^ "{" S hibar sub a (f) "}" ~=~ {sigma sup 4} over 49 ~ sum from i=-3 to 3 ~ sum from j=-3 to 3 ~ | LAMBDA sub ij | sup 2 ~=~ {sigma sup 4} over 49 ~ || LAMBDA || sub F sup 2
where || || F2 denotes the squared Frobenius norm.
For the smoothed periodogram estimate (i.e., direct estimate using a boxcar taper), $LAMBDA sub ij ^=^ delta sub ij$, and $roman Var "{" S hibar sub a (f) "}"$ $=$ $( sigma sup 4 / 7)$. Values of $(7 / {sigma sup 4}) ^ roman Var ^ "{" S hibar sub a (f) "}"$ for the smoothed periodogram, 20% cosine taper and Hann taper direct spectral estimates are tabulated in Table 2. Notice that as more data are discarded by the taper, the variance of the spectral estimate increases.
For the smoothed high-resolution multitaper spectral estimate (4.9),
matrix { ccol { E "{" S hibar (f) "}" above {~} } ccol {= above =} lcol { 1 over K ~ sum from k=0 to K-1 ~ ( lambda sub k sup -1 ) E ^ "{" |y sub k (f) | sup 2 "}" above {sigma sup 2} over K ~ sum from k=0 to K-1 ~ ( lambda sub k ) sup -1 } }
When the K = 7 lowest-order 4π prolate eigentapers are used, then E {S(f )} = (1.0095) σ2, so the estimate S (f ) is mildly biased for white noise data. Also,
Var ^ "{" S hibar (f) "}" ~=~ E ^ "{" [ S hibar (f) ] sup 2 "}" ~-~ (E ^ "{" S hibar (f) "}" ) sup 2
But E ^ "{" |y sub k (f) | sup 2 |y sub {k prime} (f) | sup 2 "}" ~=~ sum from k=0 to K-1 ~ sum from {k prime =0} to K-1 ~ [ E "{" |y sub k (f) | sup 2 "}" ~cdot ~E "{" | y sub {k prime} (f) | sup 2 "}" ~+~ sigma sup 4 delta sub {kk prime}]
as the eigentapers are orthonormal (equation (4.6)) so
matrix { ccol { roman Var ^ "{" S hibar (f) "}" above {~} } ccol {= above =} lcol { {sigma sup 4} over {K sup 2} ~ sum from k=0 to K-1 ~ sum from {k prime =0} to K-1 ~ {delta sub {kk prime}} over {lambda sub k lambda sub {k prime}} above {sigma sup 4} over {K sup 2} ~ sum from k=0 to K-1 ~ 1 over {lambda sub k sup 2} } }
When the K = 7 lowest-order 4π prolate eigentapers are used, then Var{S(f )} = σ4(1.0196)} / 7. Therefore the smoothed high-resolution multitaper estimate has only slightly more variance than the smoothed periodogram estimate for Gaussian white noise data.
For the adaptive spectral estimator, dk = sqrt (λk) for white noise data (equation (4.5) of Thomson [1982]). Therefore
matrix { ccol { E "{" S hihat (f) "}" above {~} } ccol {= above =} lcol { {sigma sup 2 ~ sum from k=0 to K-1 ~ lambda sub k } over {sum from k=0 to K-1 ~ lambda sub k } above sigma sup 2 } }
and
matrix { ccol { roman Var ^ "{" S hihat (f) "}" above {~} } ccol {= above =} lcol { {sigma sup 4 ~ sum from k=0 to K-1 ~ sum from k-0 to K-1 ~ lambda sub k lambda sub {k prime} delta sub {kk prime}} over {left ( sum from k=0 to K-1 ~ lambda sub k right ) sup 2 } above {sigma sup 4 ~ sum from k=0 to K-1 ~ ( lambda sub k ) sup 2} over {left ( sum from k=0 to K-1 ~ lambda sub k right ) sup 2 } } }
For the K = 7 lowest-order 4π prolate eigentapers, Var{S(f )} = σ4(1.00038) / 7. Therefore the adaptive multitaper estimate has slightly less variance than the smoothed high-resolution estimate and slightly more variance than the periodogram estimate.
4.2. Bias
It is not useful to compare the bias performance of these spectral estimates for white noise data. One is most interested in a measure of broadband bias. Broadband bias is caused when spectral energy at one frequency leaks away to affect the spectral estimate at a distant frequency and is an important factor to consider in the estimation of spectra of colored processes.
We take as our measure of broadband bias the fraction of energy (1- λ) in the frequency band | f - f′ | < W that leaks out to affect the estimated spectrum at other frequencies. Suppose that the record consists of a single sinusoid, so that the spectrum is highly colored. The k th eigenspectrum retains λ =λk of the spectral energy of the sinusoid within a frequency band of width 2W centered on the sinusoid frequency. The fractional leakage of the smoothed high-resolution spectral estimate is
matrix { ccol { 1~-~ lambda above {~} } ccol { = above = } lcol { 1~-~ {sum from k=0 to K-1 ~ 1 over {lambda sub k} ~ int from -W to {{^}W} ~ |U sub k (f) | sup 2 ^ df } over {sum from k=0 to K-1 ~ 1 over {lambda sub k} ~ int from {- half} to {{~^} half} ~ |U sub k (f) | sup 2 ^ df } above 1-~ K over {sum from k=0 to K-1 ~ 1 over {lambda sub k} } } }
If we use the seven lowest-order 4π prolate eigentapers in the estimate, λ ~ 0.99057, so 1- λ ~ 0.00943. For the adaptive multitaper spectral estimate, a numerical calculation shows that 1 - λ ~ 0.00256.
The smoothed single-taper direct spectral estimates are also biased when the process has a colored spectrum. A single periodogram estimate allows
1~-~ lambda ~=~ 1~-~ int from {-1/N} to {{~~}1/N} ~ |A(f)| sup 2 ^ df ~approx~ 1~-~ 0.903 ~=~ 0.097
of the energy of a single sinusoid to leak outside its main lobe. The smoothed periodogram estimate allows
1~-~ lambda ~=~ 1~-~ 1 over 7 ~ sum from j=-3 to 3 ~ int from -4/N to {{~~}4/N} ~ |A(f^+^j over N )| sup 2 df ~approx~ 0.0367
of the energy in | f - f ′ | ≤ 4/ N to leak out. For the smoothed Hann taper estimate, we find 1- λ ~ 0.00934, while the smoothed 20% cosine taper estimate allows 1 - λ ~ 0.0192 of the sinusoid's energy to leak out of | f - f ′ | ≤ W.
The Hann and 20% cosine tapers do not permit as much spectral leakage as the boxcar taper, but only the smoothed Hann taper estimates exhibit broadband bias characteristics which are as good as the multitaper estimates. Numerical experiments using the ω-square and ω-cube source spectrum models of Aki [1967] demonstrate that spectral estimates employing a boxcar taper are inadequate for representing the source spectrum roll-off. The other tapers are sufficient unless the spectrum rolls off more steeply than f -4.
Clearly, for smoothed single-taper spectral estimates there is a trade-off. The more severe the taper, the less spectral leakage contaminates the estimate but also the larger the variance of the estimate. The multitaper estimates manage to defeat this trade-off by using several orthogonal leakage-resistant tapers in a single estimate. The relative variances and fractional spectral leakage that are associated with each spectral estimate are listed in Table 2 and are plotted in Figure 7 for comparison.
5. Conclusions
Multitaper spectral analysis techniques offer the seismologist formal and practical advantages over single-taper techniques. Adaptive reweighting of eigenspectra according to the predicted level of spectral leakage enables well-constrained smoothed spectral estimates in portions of the spectrum that have large amplitude, while retaining excellent resistance to spectral leakage in the region where earthquake spectra exhibit a steep roll-off. Comparisons between direct spectral estimates produced using boxcar, Hann, and 20% cosine tapers show that the boxcar taper estimate is contaminated by spectral leakage, that the Hann taper estimates can be misleading in the high-amplitude portion of the spectrum, and that the 20% cosine taper offers a compromise between these two extremes. Therefore a 20% cosine taper may be adequate in many cases but would not be suitable for the analysis of either an unusually dispersive or unusually band-limited seismic signal. However, these pathological situations present no difficulty for the adaptive multitaper estimate.
There are drawbacks to using the multitaper method. The adaptive multitaper algorithm consumes more computer time, since several FFTs must be computed for each time series and one needs to calculate a set of prolate tapers for each time series length. The computational burden is becoming a less serious problem as computer speeds increase and the cost of computation drops. Also, we have found it adequate to calculate the prolate eigentapers once for a time series of length 128 and generate tapers of other lengths N > 128 by interpolating the 128-point tapers using cubic splines (see the appendix).
Acknowledgments. This chapter has appeared in the Journal of Geophysical Research, Vol. 92, pp. 12,675-12,684, 1987 under the title "Multitaper spectral analysis of high-frequency seismograms," J. Park, C.R. Lindberg, and F.L. Vernon III. Frank Vernon participated in all phases of the development of this research except for the programming of the computer algorithms.
Appendix: Calculating P π Prolate Eigentapers
We follow a procedure similar to that outlined in the appendix of Thomson [1982] to calculate eigentapers which have a given time-bandwidth product NW and length N. A set of standard tapers of length N′ > 100 is constructed. We use N′ = 128, as a series whose length is a power of two is convenient for calculating DFTs of the taper. The matrix eigenvalue problem (4.5) is then solved for the largest 2 N′ W eigenvalues and the associated eigenvectors using EISPACK routines TRED1, BISECT, TINVIT, and TRBAK1 [see Garbow et al, 1977]. This procedure determines only the largest eigenvalues and their eigenvectors of a matrix, avoiding the numerical burden of fully decomposing the matrix. In this manner one calculates the prolate eigentapers for a time series of length N′ . Using the algorithm described by Thomson [1982], one approximates the discrete time tapers with the continuous time prolate spheroidal wave functions in order to set up an eigenvalue problem based on Gaussian quadrature. One obtains discrete tapers at nonuniform sample points that can be interpolated to produce tapers with even sampling and of a given length. In our applications we have chosen to use spline interpolation routines to interpolate the evenly spaced tapers of length N′ to produce tapers of length N >N′ .
References
Aki, K. Scaling law of seismic spectrum, J. Geophys. Res., 72, 1217-1231, 1967.
Berger, J., L. M. Baker, J. N. Brune, J. B. Fletcher, T. C. Hanks, and F. L. Vernon III, The Anza array: A high-dynamic range, broadband, digitally radio telemetered seismic array, Bull. Seismol. Soc. Am., 74, 1469-1481, 1984.
Doob, J. L., Stochastic Processes, John Wiley, New York, 1953.
Eberhard, A. An optimum discrete window for the calculation of power spectra, IEEE Trans. Audio Electroacoust., AU-21, 37-43, 1973.
Garbow, B., J. Boyle, J. Dongarra, and C. Molar, Matrix Eigensystem Routines -- EISPACK Guide Extension, Lecture Notes in Computer Science, Springer-Verlag, New York, 1977.
Harris, F., On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. Instr. Electr. Eng., 66, 51-83, 1978.
Kay, S. M., and S. L. Marple, Jr., Spectrum analysis--A modern perspective, Proc. IEEE, 69, 1380-1419, 1981.
Lindberg, C. R., Multiple taper spectral analysis of terrestrial free oscillations, Ph.D. thesis, 182 pp., Univ. of Calif., San Diego, La Jolla, 1986.
Papoulis, A., Signal Analysis, McGraw-Hill, New York, 1977.
Park, J., C. R. Lindberg, and D. J. Thomson, Multiple taper spectral analysis of terrestrial free oscillations: Part I, Geophys. J. R. Astron. Soc., 91, 755-794, 1987.
Park, J., F. L. Vernon III, and C. R. Lindberg, Frequency dependent polarization analysis of high-frequency seismograms, J. Geophys. Res.. 92, 12,664-12,674, 1987.
Slepian, D., Prolate spheroidal wave functions, Fourier analysis, and uncertainty, V, The discrete case, Bell Syst. Tech. J., 57, 1371-1429, 1978.
Thomson, D. J., Spectral analysis of short series, Ph.D. dissertation, Polytech. Inst. of Brooklyn, New York, 1971.
Thomson, D. J., Spectrum estimation techniques for characterization and development of WT4 waveguide, I, Bell Syst. Tech. J., 56, 1769-1815, 1977a.
Thomson, D. J., Spectrum estimation techniques for characterization and development of WT4 waveguide, II, Bell Syst. Tech. J., 56, 1983-2005, 1977b.
Thomson, D. J., Spectrum estimation and harmonic analysis, IEEE Proc., 70, 1055-1096, 1982.
