R28. Fast Fourier Transform in R

The Discrete Fourier Transform (DFT) is the pair of transforms for the transform and its inverse:

The Fast Fourier Transform (FFT) is the solution of DFT using an algorithm based on symmetry of equations.

We use the fft function to solve for the transform of a 64-length vector. It returns a complex vector, which is again sent to fft function but with inverse set to TRUE, so the sign in exponent is plus. We compare it to the the result from dft and inverse dft. The result are the same, with errors on order of double-precision tolerance.

# ex28.R

dft <- function(x, inverse = FALSE) {
  s <- if(inverse) 1 else -1
  N <- length(x)
  X <- complex(N)
  for (k in 0:(N-1)) {
    for (n in 0:(N-1)) {
      X[k+1] = X[k+1] + x[n+1]*exp(s*1i*2*pi*k*n/N)
    }
  }
  X
}

x <- c(0:31,31:0)
y <- fft(x)
y1 <- dft(x)
z <- fft(y,TRUE)/length(y)
z1 <- dft(y1,TRUE)/length(y1)
xR <- Re(z)
x1R <- Re(z1)
error.fft <- sd(x-xR)
error.dft <- sd(x-x1R)
out <- sprintf("error.fft = %e\nerror.dft = %e",
               error.fft, error.dft)
cat(out)
# error.fft = 2.976939e-15
# error.dft = 1.067366e-13