# [R] Layout of Fourier frequencies

David Brahm brahm at alum.mit.edu
Thu Apr 11 18:10:30 CEST 2002

```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==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 <- 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 <- 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), 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```