[R] Layout of Fourier frequencies

[David Brahm]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._