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 halfcosine) 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 timeseries. 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.
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, AU15, 7073.
library(oce) Fs < 1000 t < seq(0, 0.296, 1/Fs) x < cos(2 * pi * t * 200) + rnorm(n=length(t)) X < ts(x, frequency=Fs) s < spectrum(X, spans=c(3,2), main="random + 200 Hz", log='no')cat("Checking spectral levels with Parseval's theorem:\n")#> Checking spectral levels with Parseval's theorem:#> var(x) = 1.588179#> 2 * sum(s$spec) * diff(s$freq[1:2]) = 1.541786#> sum(w$spec) * diff(s$freq[1:2]) = 0.8662452#> sum(w2$spec) * diff(s$freq[1:2]) = 0.8741546## 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("Variancepreserving spectrum")