Adaptive, Sine-Multitaper Power Spectral Density Estimation

Produces power spectral density estimates through iterative refinement of the optimal number of sine-tapers at each frequency. This optimization procedure is based on the method of Riedel and Sidorenko (1995), which minimizes the Mean Square Error (sum of variance and bias) at each frequency, but modified for computational stability.

Adaptive, sine multitaper power spectral density estimation for R

by Andrew J Barbour and Robert L Parker

###Description### This is an R package for computing univariate power spectral density estimates with little or no tuning effort. We employ sine multitapers, allowing the number to vary with frequency in order to reduce mean square error, the sum of squared bias and variance, at each point. The approximate criterion of Riedel and Sidorenko (1995) is modified to prevent runaway averaging that otherwise occurs when the curvature of the spectrum goes to zero. An iterative procedure refines the number of tapers employed at each frequency. The resultant power spectra possess significantly lower variances than those of traditional, non-adaptive estimators. The sine tapers also provide useful spectral leakage suppression. Resolution and uncertainty can be estimated from the number of degrees of freedom (twice the number of tapers).

This technique is particularly suited to long time series, because it demands only one numerical Fourier transform, and requires no costly additional computation of taper functions, like the Slepian functions. It also avoids the degradation of the low-frequency performance associated with record segmentation in Welch's method. Above all, the adaptive process relieves the user of the need to set a tuning parameter, such as time-bandwidth product or segment length, that fixes frequency resolution for the entire frequency interval; instead it provides frequency-dependent spectral resolution tailored to the shape of the spectrum itself.

psd elegantly handles spectra with large dynamic range and mixed-bandwidth features|features typically found in geophysical datasets.

##How to Cite##

Bob and I have published a paper in Computers & Geosciences to accompany this software (download a pdf, 1MB); it describes the theory behind the estimation process, and how we apply it in practice. If you find psd useful in your research, we kindly request you cite our paper.


##Getting Started##

Firstly you'll need to install the package and it's dependencies from CRAN (from within the R environment):

install.packages("psd", dependencies=TRUE)

then load the package library


We have included a dataset to play with, named Tohoku, which represents recordings of high-frequency borehole strainmeter data during teleseismic waves from the 2011 Mw 9.0 Tohoku earthquake (source). Access and inspect these data with:


The 'preseismic' data has interesting spectral features, so we subset it, and use the areal strain (the change in borehole diameter):

Dat <- subset(Tohoku, epoch=="preseismic")
Areal <- ts(Dat$areal)

For the purposes of spectral estimation, we remove a linear trend:

Dat <- prewhiten(Areal, plot=FALSE)

Now we can calculate the adaptive PSD:

mtpsd <- pspectrum(Dat$prew_lm, plot=TRUE)

We can visualize the spectrum with builtin methods:

plot(mtpsd, log="dB")

and the spectral uncertainty:

sprop <- spectral_properties(mtpsd)
Ntap <- sprop$taper/max(sprop$taper)
plot(Ntap, type="h", ylim=c(0,2), col="dark grey") 

##Installing the Development Version##

Should you wish to install the development version of this software, the devtools library will be useful:

install.packages("devtools", dependencies=TRUE)


Reference manual

It appears you don't have a PDF plugin for this browser. You can click here to download the reference manual.


1.0-1 by Andrew J. Barbour, 2 years ago, Barbour and Parker (2014):, Riedel and Sidorenko (1995):

Report a bug at

Browse source code at

Authors: Andrew J. Barbour and Robert L. Parker, with contributions from David Myer

Documentation:   PDF Manual  

Task views: Time Series Analysis

GPL (>= 2) license

Imports Rcpp, RColorBrewer, signal, zoo

Depends on utils, stats, graphics, grDevices

Suggests bspec, fftw, ggplot2, knitr, multitaper, plyr, RSEIS, rbenchmark, reshape2, testthat

Linking to Rcpp, RcppArmadillo

Suggested by multitaper.

See at CRAN