SDF: Nonparametric (cross) spectral density function estimation
In sapa: Spectral Analysis for Physical Applications

Description Usage Arguments Details Value S3 METHODS References See Also Examples

Estimate the process (cross) spectral density function via nonparametric models.

SDF(x, method="direct", taper.=NULL, window=NULL,
    n.taper=5, overlap=0.5, blocksize=NULL,
    single.sided=TRUE, sampling.interval=NULL,
    center=TRUE, recenter=FALSE, npad=2*numRows(x))

`x`	a vector or matrix containing uniformly-sampled real-valued time series. If a `matrix`, each column should contain a different time series.
`blocksize`	an integer representing the number of points (width) of each block in the WOSA estimator scheme. Default: `floor(N/4)` where `N` is the number of samples in each series.
`center`	a logical value. If `TRUE`, the mean of each time series is recentered prior to estimating the SDF. Default: `TRUE`.
`method`	a character string denoting the method to use in estimating the SDF. Choices are `"direct"`, `"lag window"`, `"wosa"` (Welch's Overlapped Segment Averaging), `"multitaper"`. See DETAILS for more information. Default: `"direct"`.
`n.taper`	an integer defining the number of tapers to use in a multitaper scheme. This value is overwritten if the `taper` input is of class `taper`. Default: `5`.
`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: `2*numRows(x)`.
`overlap`	a numeric value on [0,1] denoting the fraction of window overlap for the WOSA estimator. Default: `0.5`.
`recenter`	a logical value. If `TRUE`, the mean of each time series is recentered after (posssibly) tapering the series prior to estimating the SDF. Default: `FALSE`.
`sampling.interval`	a numeric value representing the interval between samples in the input time series `x`. Default: `NULL`, which serves as a flag to obtain the sampling interval via the `deltat` function. If `x` is a list, the default sampling interval is `deltat(x[[1]])`. If `x` is an atomic vector (ala `isVectorAtomic`), then the default samplign interval is established ala `deltat(x)`. Finally, if the input series is a matrix, the sampling interval of the first series (assumed to be in the first column) is obtained ala `deltat(x[,1])`.
`single.sided`	a logical value. If `TRUE`, a single-sided SDF estimate is returned corresponding to the normalized frequency range of [0,1/2]. Otherwise, a double-sided SDF estimate corresponding to the normalized frequency interval [-1/2,1/2] is returned. Default: `TRUE`.
`taper.`	an object of class `taper` or a character string denoting the primary taper. If an object of class `taper`, the length of the taper is checked to ensure compatitbility with the input `x`. See DETAILS for more information. The default values are a function of the `method` as follows: direct normalized rectangular taper lag window normalized Parzen window with a cutoff at N/2 where N is the length of the time series. wosa normalized Hanning taper multitaper normalized Hanning taper
`window`	an object of class `taper` or a character string denoting the (secondary) window for the lag window estimator. If an object of class `taper`, the length of the taper is checked to ensure compatitbility with the input `x`. See DETAILS for more information. Default: Normalized Hanning window.

Let x(t) be a uniformly sampled real-valued 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,...,N-1]{h(t)*x(t)*exp(-i*2*pi*f*t)}|^2, where h(t) is a data taper normalized such that sum[t=0,...,N-1]{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,...,N-1]{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=-(N-1),...,(N-1)]{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,...,Nb-1]{S(j*No,f)}

where

S(l,f) =|sum[t=0,...,Ns-1]{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,...,Ns-1]{h_t^2} = 1).

multitaper

A multitaper spectral estimator is given by

S(f) = (1/K) * sum[k=0,...,K-1] S(k,f)

where S(k,f) = |sum[t=0,...,N-1]{h(k,t) * X(t) * exp(-i*2*pi*f*t)}|^2 and h(k,t) for k=0,...,K-1, 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 cross-SDF 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 cross-SDF estimate of the ith and jth column of x. All cross-spectral density function estimates are returned as complex-valued 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.

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 x-axis label. Default: "FREQUENCY (Hz)".
ylab: a (vector of) character string(s), one per (cross-)SDF estimate, representing the y-axis label(s). Default: in the multivariate case, the strings "Sij" are used for the y-axis labels, where i and j are the indices of the different variables. For example, if the user supplies a 2-column matrix for x, the labels "S11", "S12", and "S22" are used to label the y-axes 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 y-axis.
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 and Arg. See each of these functions for details. If the SDF is purely real (no cross-SDF is calculated), this argument is coerced to the Mod function. Default: Mod.
add: 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.

print

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.

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.

taper, ACVS.

## calculate various SDF estimates for the 
## sunspots series. remove mean component for a 
## better comparison. 

require(ifultools)
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 cross-spectrum 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")