welchPSD {bspec}

Power spectral density estimation using Welch's method.

Description

Estimates a time series' power spectral density using Welch's method, i.e., by subdividing the data into segments, computing spectra for each, and averaging over the results.

Usage

welchPSD(x, seglength, two.sided = FALSE, windowfun = tukeywindow,
         method = c("mean", "median"), windowingPsdCorrection = TRUE, ...)

Arguments

x

a time series (ts object).

seglength

the length of the subsegments to be used (in units of time relative to x).

two.sided

a logical flag indicating whether the result is supposed to be one-sided (default) or two-sided.

windowfun

The windowing function to be used.

method

The "averaging" method to be used – either "mean" or "median".

windowingPsdCorrection

a logical flag indicating whether an overall correction for the windowing is supposed to be applied to the result – this essentially specifies whether the result is supposed to reflect the PSD of the windowed or of the "un-windowed" time series.

...

other parameters passed to the windowing function.

Details

The time series will be divided into overlapping sub-segments, each segment is windowed and its "empirical" spectrum is computed. The average of these spectra is the resulting PSD estimate. For robustness, the median may also be used instead of the mean.

Value

A list containing the following elements:

frequency

the Fourier frequencies.

power

the estimated spectral power.

kappa

the number of (by definition) non-zero imaginary components of the Fourier series.

two.sided

a logical flag indicating one- or two-sidedness.

segments

the number of (overlapping) segments used.

Author(s)

Christian Roever, [e-mail]

References

Welch, P. D. The use of Fast Fourier Transform for the estimation of Power Spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, AU-15(2):70–73, 1967. doi: 10.1109/TAU.1967.1161901.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. Numerical recipes in C. Cambridge University Press, 1992.

See Also

empiricalSpectrum, tukeywindow, spectrum

Examples

# load example data:
data(sunspots)
# compute and plot the "plain" spectrum:
spec1 <- empiricalSpectrum(sunspots)
plot(spec1$frequency, spec1$power, log="y", type="l")

# plot Welch spectrum using segments of length 10 years:
spec2 <- welchPSD(sunspots, seglength=10)
lines(spec2$frequency, spec2$power, col="red")

# use 20-year segments and a flatter Tukey window:
spec3 <- welchPSD(sunspots, seglength=20, r=0.2)
lines(spec3$frequency, spec3$power, col="blue")