Timothy H. Keitt <Timothy.Keitt at stonybrook.edu>
wrote:> I'm doing convolutions in the frequency domain and need to know the
> layout of the Fourier modes returned by fft...
> I think in 1D the pattern is:
> 0 1 2 3 -2 1 (even)
> 0 1 2 3 -3 2 1 (odd)
I believe that's correct. Perhaps you'll find my fft cheat-sheet
(below)
handy, though it's only relevant to a univariate real series "z".
Define:
N <- length(z)
w <- fft(z)/N # Division by N here is a little
unconventional
jmax <- floor(N/2+1)
tv <- 0:(N-1) # Time values, starting at
0
w[1]==mean(z) is the constant term.
For j=2:jmax, w[j] is the Fourier coefficient for frequency (j-1)/N, and
w[N+2-j] = Conj(w[j]). Note that for even N, this implies w[jmax] is real.
You could derive the Fourier transform yourself (more slowly) with:
w <- rep(NA, N)
w[1] <- mean(z)
for (j in 2:jmax) {
w[j] <- sum(z * exp(-2i*pi*tv*(j-1)/N)) / N
w[N+2-j] <- Conj(w[j])
}
or even:
w <- rep(NA, N)
w[1] <- mean(z)
for (j in 2:jmax) {
phi <- 2*pi*tv*(j-1)/N
w[j] <- sum(z * (cos(phi) - 1i*sin(phi))) / N
w[N+2-j] <- sum(z * (cos(phi) + 1i*sin(phi))) / N
}
You can reconstruct the original series "z" either by:
zz <- Re(fft(w, inverse=T))
or (more slowly) by:
zz <- rep(Re(w[1]), N)
for (j in 2:jmax) {
dup <- if (N%%2==0 && j==jmax) 1 else 2 # Don't double
jmax if N even
zz <- zz + dup * Mod(w[j]) * cos(Arg(w[j]) + 2*pi*tv*(j-1)/N)
}
The power spectrum ("periodogram"):
s <- N * abs(w[2:jmax])^2
plot((2:jmax-1)/N, s, type="l", log="y",
xlab="frequency", ylab="spectrum")
can also be obtained from the "spectrum" function:
s1 <- spectrum(z, fast=F, taper=0, detrend=F) # Plots
automatically
plot(s1$freq, s1$spec, type="l", log="y",
xlab="frequency", ylab="spectrum")
Note I had to override the default to detrend, taper, and pad z. Smoother
(more meaningful?) spectra are obtained with, e.g.:
s2 <- spectrum(z, spans=c(2*round(N/20)+1, 2*round(N/10)+1))
s3 <- spec.ar(z)
Among other things, the power spectrum shows the distribution in frequency
space of the "sum-of-squares" errors. Again you must watch that
tricky jmax:
sp <- s; if (N%%2==0) sp[jmax-1] <- sp[jmax-1]/2
but then (N-1)*var(z) == 2*sum(sp).
--
-- David Brahm (brahm at alum.mit.edu)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._