Estimate the process (cross) spectral density function via nonparametric models.
1 2 3 4 
x 
a vector or matrix containing uniformlysampled realvalued time series.
If a 
blocksize 
an integer representing the number
of points (width) of each block in the WOSA estimator scheme.
Default: 
center 
a logical value. If 
method 
a character string denoting the method to use in estimating the SDF.
Choices are 
n.taper 
an integer defining the number of tapers
to use in a multitaper scheme. This value is overwritten if
the 
npad 
an integer representing the total length of each
time series to analyze after padding with zeros. This argument
allows the user to control the spectral resolution of the SDF
estimates: the normalized frequency interval is
deltaf=1/npad.
This argument must be set such that
npad > 2.
Default: 
overlap 
a numeric value on [0,1] denoting the fraction
of window overlap for the WOSA estimator. Default: 
recenter 
a logical value. If 
sampling.interval 
a numeric value representing the interval
between samples in the input time series 
single.sided 
a logical value. If 
taper. 
an object of class

window 
an object of class 
Let x(t) be a uniformly sampled realvalued time series of length N, Let an estimate of the process spectral density function be denoted as S(f) where f are frequencies on the interval 1/(2*deltat),1/(2*deltat) where deltat is the sampling interval. The supported SDF estimators are:
The direct SDF estimator is defined as
S(f)=sum[t=0,...,N1]{h(t)*x(t)*exp(i*2*pi*f*t)}^2,
where h(t) is a data taper normalized such that
sum[t=0,...,N1]{h(t)^2} = 1. If
h(t)=1/sqrt(N) then we obtain the definition
of the periodogram
S(f)=(1/N) * sum[t=0,...,N1]{x(t)*exp(i*2*pi*f*t)}^2.
See the taper
function for more details on supported window types.
The lag window SDF estimator is defined as
S(f)=sum[k=(N1),...,(N1)]{w(k)*s(k)*exp(i*2*pi*f*k)}^2,
where s(k) is the autocovariance
sequence estimator corresponding to some
direct spectral estimator (often the periodogram)
and w(k) is a lag window (popular choices
are the Parzen, Papoulis, and Daniell windows). See the
taper
function for more details.
Welch's Overlapped Segment Averaging SDF estimator is defined as
S(f)=(1/Nb)*sum[j=0,...,Nb1]{S(j*No,f)}
where
S(l,f) =sum[t=0,...,Ns1]{h(t)*x(t+l)*exp(i*2*pi*f*t)}^2
Here, No is a positive integer that controls how much overlap there is between segments and that must satisfy both No <= Ns and No * (Nb  1) = N  Ns, while h(t) is a data taper appropriate for a series of length Ns (i.e., sum[t=0,...,Ns1]{h_t^2} = 1).
A multitaper spectral estimator is given by
S(f) = (1/K) * sum[k=0,...,K1] S(k,f)
where S(k,f) = sum[t=0,...,N1]{h(k,t) * X(t) * exp(i*2*pi*f*t)}^2 and h(k,t) for k=0,...,K1, is a set of K orthonormal data tapers.
See reference(s) for further details.
Popular choices for multitapers include sinusoidal tapers and
discrete prolate spheroidal sequences (DPSS). See the
taper
function for more details.
Cross spectral density function estimation:
If the input x
is a matrix, where each column contains
a different time series, then the results are returned in
a matrix whose columns correspond to all possible unique combinations
of crossSDF estimates. For example, if x
has three columns,
then the output will be a matrix whose columns are
{S11, S12, S13, S22, S23, S33}
where Sij is the crossSDF estimate of the i
th
and j
th column of x
. All crossspectral density
function estimates are returned as complexvalued series to maintain
the phase relationships between components. For all Sij
where i=j, however, the imaginary portions will be zero (up to a
numerical noise limit).
an object of class SDF
.
converts the (cross)SDF estimate(s) as a matrix. Optional arguments
are passed directly to the matrix
function during the conversion.
plots the (cross)SDF estimate(s). Optional arguments are:
a character string defining the scaling to perform on the (common) frequency vector
of the SDF estimates. See the scaleData
function for supported choices. Default: "linear"
.
a character string defining the scaling to perform on the SDF estimates.
See the scaleData
function for supported choices. Default: "linear"
.
a single character defining the plot type (ala the par
function) of the
SDF plots. Default: ifelse(numRows(x) > 100, "l", "h")
.
a character string representing the xaxis label. Default: "FREQUENCY (Hz)"
.
a (vector of) character string(s), one per (cross)SDF estimate,
representing the yaxis label(s). Default: in the multivariate case, the strings
"Sij"
are used for the yaxis labels, where i and j are the indices of the
different variables. For example, if the user supplies a 2column matrix for x
,
the labels "S11"
, "S12"
, and "S22"
are used to label the yaxes of the corresponding
(cross)SDF plots. In the univariate case, the default string "SDF"
prepended with a string
describing the type of SDF performed (such as "Multitaper"
) is used to label the yaxis.
a logical value. If TRUE
, the SDF value at normalized frequency f=0
is plotted for each SDF. This frequency is associated with the sample mean
of the corresponding time series. A relatively large mean value dominates
the spectral patterns in a plot and thus the corresponding frequency is typically not plotted.
Default: !attr(x,"center")
.
an integer defining the maximum number of SDF plots to place onto a single graph.
Default: 3
.
a post processing function to apply to the SDF values prior to plotting. Supported
functions are Mod
, Im
, Re
and Arg
. See each of these functions for details.
If the SDF is purely real (no crossSDF is calculated), this argument is coerced to the Mod
function.
Default: Mod
.
A logical value. If TRUE
, the plot is added using the current
par()
layout. Otherwise a new plot is produced. Default: FALSE
.
additional plot parameters passed directly to the genPlot
function used
to plot the SDF estimates.
prints the object. Available options are:
text justification ala prettPrintList
. Default: "left"
.
header separator ala prettyPrintList
. Default: ":"
.
Additional print arguments sent directly to the prettyPrintList
function.
Percival, Donald B. and Constantine, William L. B. (2005) “Exact Simulation of Gaussian Time Series from Nonparametric Spectral Estimates with Application to Bootstrapping", Journal of Computational and Graphical Statistics, accepted for publication.
D.B. Percival and A. Walden (1993), Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, Cambridge, UK.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22  ## calculate various SDF estimates for the
## sunspots series. remove mean component for a
## better comparison.
data < as.numeric(sunspots)
methods < c("direct","wosa","multitaper",
"lag window")
S < lapply(methods, function(x, data) SDF(data, method=x), data)
x < attr(S[[1]], "frequency")[1]
y < lapply(S,function(x) decibel(as.vector(x)[1]))
names(y) < methods
## create a stack plot of the data
stackPlot(x, y, col=1:4)
## calculate the crossspectrum of the same
## series: all spectra should be the same in
## this case
SDF(cbind(data,data), method="lag")
## calculate the SDF using npad=31
SDF(data, npad=31, method="multitaper")

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.