matlab pmtm,Multitaper power spectral density estimate

The periodogram is not a consistent estimator of the true power spectral

density (PSD) of a wide-sense stationary process. To reduce the variability in the periodogram —

and thus produce a consistent estimate of the PSD — the multitaper method averages modified

periodograms obtained using a family of mutually orthogonal windows or

tapers. In addition to mutual orthogonality, the tapers also have

optimal time-frequency concentration properties. Both the orthogonality and time-frequency

concentration of the tapers are critical to the success of the multitaper technique. See Discrete Prolate Spheroidal (Slepian) Sequences for a brief description of the Slepian

sequences used in Thomson’s multitaper method.

The multitaper method uses K modified periodograms, each one obtained

using a different Slepian sequence as the window. Let

Sk(f)=Δt|∑n=0N−1gk(n)x(n)e−j2πfnΔt|2

denote the modified periodogram obtained with the kth

Slepian sequence, gk(n). In its

simplest form, the multitaper method simply averages the K modified

periodograms to produce the multitaper PSD estimate:

S(MT)(f)=1K∑k=0K−1Sk(f).

Thomson's multitaper approach, introduced in [4], resembles

Welch’s overlapped segment averaging method, in that both average over approximately

uncorrelated estimates of the PSD. However, the two approaches differ in how they produce these

uncorrelated PSD estimates. The multitaper method uses the entire signal in each modified

periodogram. The orthogonality of the Slepian tapers decorrelates the different modified

periodograms. Welch’s approach uses segments of the signal in each modified periodogram, and the

segmenting decorrelates the different modified periodograms.

The equation for S(MT)(f) corresponds to the 'unity' option in

pmtm. However, as explained in Discrete Prolate Spheroidal (Slepian) Sequences, the Slepian sequences do not possess

equal energy concentration in the frequency band of interest. The higher the order of the

Slepian sequence, the less concentrated the sequence energy is in the band [–W,W] with the concentration given by the eigenvalue. Consequently, it can be

beneficial to use the eigenvalues to weight the K modified periodograms prior

to averaging. This corresponds to the 'eigen' option in

pmtm.

Using the sequence eigenvalues to produce a weighted average of modified periodograms

accounts for the frequency concentration properties of the Slepian sequences. However, it does

not account for the interaction between the power spectral density of the random process and the

frequency concentration of the Slepian sequences. Specifically, frequency regions where the

random process has little power are less reliably estimated in the modified periodograms using

higher-order Slepian sequences. This argues for a frequency-dependent adaptive process, which

accounts not only for the frequency concentration of the Slepian sequence but also for the power

distribution in the time series. This adaptive weighting corresponds to the

'adapt' option in pmtm and is the default for

computing the multitaper estimate.