oce (version 0.2-1)

pwelch: Welch periodogram

Description

Compute periodogram using the Welch (1967) method

Usage

pwelch(x, window, noverlap, nfft, fs, spectrumtype, esttype,
     debug=getOption("oceDebug"), ...)

Arguments

x

a vector or timeseries to be analyzed. If a timeseries, then there is no need to specify fs.

window

window specification, either a single value giving the number of windows to use, or a vector of values, specifying the window coefficients. If not specified, then 8 windows are used, each with a Hamming (raised half-cosine) window. BUG: i

noverlap

number of points to overlap between windows. If not specified, this will be set to half the window length.

nfft

length of FFT. This cannot be given if window is given, and the latter is a single integer.

fs

frequency of time-series. If x is a time-series, and if fs is supplied, then time-series is altered to have frequency fs.

spectrumtype

not used (yet)

esttype

not used (yet)

debug

a flag that turns on debugging. Set to 1 to get a moderate amount of debugging information, or to 2 to get more.

...

optional extra arguments to be passed to spectrum. Unless specified in this list, spectrum is called with plot=FALSE to prevent pl

Details

The Welch (1967) method is used.

References

Welch, P. D., 1967. 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 Electroacoustics, AU-15, 70--73.

Examples

library(oce)
Fs <- 1000
t <- seq(0, 0.296, 1/Fs)
x <- cos(2 * pi * t * 200) + rnorm(n=length(t))
xts <- ts(x, frequency=Fs)
s <- spectrum(xts, spans=c(3,2), main="random + 200 Hz", log='no')
w <- pwelch(xts)
lines(w$freq, w$spec, col="red")
w2 <- pwelch(xts, nfft=75)
lines(w2$freq, w2$spec, col='green')
abline(v=200, col="blue", lty="dotted")
cat("Checking spectral levels with Parseval's theorem:
")
cat("var(x)                              = ", var(x), "")
cat("2 * sum(s$spec) * diff(s$freq[1:2]) = ", 2 * sum(s$spec) * diff(s$freq[1:2]), "")
cat("sum(w$spec) * diff(s$freq[1:2])     = ", sum(w$spec) * diff(w$freq[1:2]), "")
cat("sum(w2$spec) * diff(s$freq[1:2])    = ", sum(w2$spec) * diff(w2$freq[1:2]), "")