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")
w <- pwelch(X, plot = FALSE)
lines(w$freq, w$spec, col = "red")
w2 <- pwelch(X, 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")
Run the code above in your browser using DataLab