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.
Maintainer: Andrew J. Barbour andy.barbour@gmail.com (ORCID)
Authors:
Jonathan Kennel (ORCID)
Robert L. Parker
Andrew J. Barbour <andy.barbour@gmail.com>, Jonathan Kennel, and Robert L. Parker
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.
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.
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.
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.
Useful links:
tools:::Rd_expr_doi("10.1016/j.cageo.2013.09.015")
Report bugs at https://github.com/abarbour/psd/issues
pspectrum (main function); psdcore and riedsid