pwelch {oce} | R Documentation |
Compute periodogram using the Welch (1967) method.
First, x
is broken up into chunks,
overlapping as specified by noverlap
. These chunks are then
detrended with detrend
, multiplied by the window, and then
passed to spectrum
. The resulting spectra are then averaged,
with the results being stored in spec
of the return value. Other
entries of the return value mimic those returned by spectrum
.
pwelch(x, window, noverlap, nfft, fs, spectrumtype, esttype, plot = TRUE, debug = getOption("oceDebug"), ...)
x |
a vector or timeseries to be analyzed. If a timeseries, then there
is no need to specify |
window |
window specification, either a single value giving the number of windows to use, or a vector of window coefficients. If not specified, then 8 windows are used, each with a Hamming (raised half-cosine) window. |
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 |
fs |
frequency of time-series. If |
spectrumtype |
not used (yet) |
esttype |
not used (yet) |
plot |
logical, set to |
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
|
List mimicking the return value from spectrum
,
containing frequency freq
, spectral power spec
, degrees of
freedom df
, bandwidth bandwidth
, etc.
Both bandwidth and degrees of freedom are just copied from the values for one of the chunk spectra, and are thus incorrect. That means the cross indicated on the graph is also incorrect.
Dan Kelley
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.
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, plot=FALSE) lines(w$freq, w$spec, col="red") w2 <- pwelch(xts, nfft=75, plot=FALSE) lines(w2$freq, w2$spec, col='green') abline(v=200, col="blue", lty="dotted") cat("Checking spectral levels with Parseval's theorem:\n") cat("var(x) = ", var(x), "\n") cat("2 * sum(s$spec) * diff(s$freq[1:2]) = ", 2 * sum(s$spec) * diff(s$freq[1:2]), "\n") cat("sum(w$spec) * diff(s$freq[1:2]) = ", sum(w$spec) * diff(w$freq[1:2]), "\n") cat("sum(w2$spec) * diff(s$freq[1:2]) = ", sum(w2$spec) * diff(w2$freq[1:2]), "\n") ## co2 par(mar=c(3,3,2,1), mgp=c(2,0.7,0)) s <- spectrum(co2, plot=FALSE) plot(log10(s$freq), s$spec * s$freq, xlab=expression(log[10]*Frequency), ylab="Power*Frequency", type='l') title("Variance-preserving spectrum") pw <- pwelch(co2, nfft=256, plot=FALSE) lines(log10(pw$freq), pw$spec * pw$freq, col='red')