psd-package {psd}

Adaptive power spectral density estimation using optimal sine multitapers

Description

Estimate the power spectral density (PSD) of a timeseries using the sine multitapers, adaptively; the number of tapers (and hence the resolution and uncertainty) vary according to spectral shape. The main function to be used is pspectrum.

Details

In frequency ranges where the spectrum (S) is relatively flat, more tapers are taken and so a higher accuracy is attained at the expense of lower frequency resolution. The program makes a pilot estimate of the spectrum, then uses Riedel and Sidorenko's (1995) estimate of the MSE (minimum square error), which is based on an estimate of the second derivative of the PSD (S''). The process is repeated niter times; further iteration may be necessary to reach convergence, or an acceptably low spectral variance. In this context the term "acceptable" is rather subjective: one can usually detect an unconverged state by a rather jagged appearance of the spectrum, but this is uncommon in our experience.

Adaptive estimation

The adaptive process used is as follows. A quadratic fit to the logarithm of the PSD within an adaptively determined frequency band is used to find an estimate of the local second derivative of the spectrum. This is used in an equation like R-S equation (13) for the MSE taper number, with the difference that a parabolic weighting is applied with increasing taper order. Because the FFTs of the tapered series can be found by resampling the FFT of the original time series (doubled in length and padded with zeros) only one FFT is required per series, no matter how many tapers are used. The spectra associated with the sine tapers are weighted before averaging with a parabolically varying weight. The expression for the optimal number of tapers given by R-S must be modified since it gives an unbounded result near points where S'' vanishes, which happens at many points in most spectra. This program restricts the rate of growth of the number of tapers so that a neighboring covering interval estimate is never completely contained in the next such interval.

Resolution and uncertainty

The sine multitaper adaptive process introduces a variable resolution and error in the frequency domain. See documentation for spectral_properties details on how these are computed.

Author(s)

Andrew J. Barbour <[e-mail]>, Jonathan Kennel, and Robert L. Parker

References

Barbour, A. J. and R. L. Parker, (2014), psd: Adaptive, sine multitaper power spectral density estimation for R, Computers and Geosciences, 63, 1–8, doi: 10.1016/j.cageo.2013.09.015

Percival, D. B., and A.T. Walden (1993), Spectral analysis for physical applications, Cambridge University Press

Prieto, G. A., R. L. Parker, D. J. Thomson, F. L. Vernon, and R. L. Graham (2007), Reducing the bias of multitaper spectrum estimates, Geophysical Journal International, 171, 1269–1281, doi: 10.1111/j.1365-246X.2007.03592.x

Riedel, K. S., & Sidorenko, A. (1995), Minimum bias multiple taper spectral estimation, Signal Processing, IEEE Transactions on, 43(1), 188–195.

See Also

pspectrum (main function); psdcore and riedsid