Nonparametric (cross) spectral density function estimation
Description
Estimate the process (cross) spectral density function via nonparametric models.
Usage
1 2 3 4 
Arguments
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 
Details
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:
 direct
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. lag window
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. wosa
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).
 multitaper
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).
Value
an object of class SDF
.
S3 METHODS
 as.matrix
converts the (cross)SDF estimate(s) as a matrix. Optional arguments are passed directly to the
matrix
function during the conversion. plot
plots the (cross)SDF estimate(s). Optional arguments are:
 xscale
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"
. yscale
a character string defining the scaling to perform on the SDF estimates. See the
scaleData
function for supported choices. Default:"linear"
. type
a single character defining the plot type (ala the
par
function) of the SDF plots. Default:ifelse(numRows(x) > 100, "l", "h")
. xlab
a character string representing the xaxis label. Default:
"FREQUENCY (Hz)"
. ylab
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 forx
, 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. plot.mean
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")
. n.plot
an integer defining the maximum number of SDF plots to place onto a single graph. Default:
3
. FUN
a post processing function to apply to the SDF values prior to plotting. Supported functions are
Mod
,Im
,Re
andArg
. See each of these functions for details. If the SDF is purely real (no crossSDF is calculated), this argument is coerced to theMod
function. Default:Mod
. add
A logical value. If
TRUE
, the plot is added using the currentpar()
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:
 justify
text justification ala
prettPrintList
. Default:"left"
. sep
header separator ala
prettyPrintList
. Default:":"
. ...
Additional print arguments sent directly to the
prettyPrintList
function.
References
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.
See Also
taper
, ACVS
.
Examples
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")
