This is the second in a series of three posts about the Finite Fourier Transform. This post is about the fast FFT algorithm itself. A recursive divide and conquer algorithm is implemented in an elegant MATLAB function named ffttx.

Contents

FFFT

The acronym FFT is ambiguous. The first F stands for both "fast" and "finite". A more accurate abbreviation would be FFFT, but nobody wants to use that. In MATLAB, the expression fft(x) computes the finite Fourier transform of any vector x. The computation is fastest if the integer n = length(x) is the product of powers of small primes.

FFT

Direct application of the definition

$$ Y_k = \sum_{j=0}^{n-1} \omega^{jk} y_j, \ k = 0, \ldots, n-1 $$

requires $n$ multiplications and $n$ additions for each of the $n$ components of $Y$ for a total of $2 n^2$ floating-point operations. This does not include the generation of the powers of $\omega$. A computer capable of doing one multiplication and addition every microsecond would require a million seconds, or about 11.5 days, to do a million-point FFT.

Several people discovered fast FFT algorithms independently and many people have since contributed to their development, but it was a 1965 paper by John Tukey of Princeton University and John Cooley of IBM Research that is generally credited as the starting point for the modern usage of the FFT.

Modern fast FFT algorithms have computational complexity $O(n \mbox{ log}_2 n)$ instead of $O(n^2)$. If $n$ is a power of 2, a one-dimensional FFT of length $n$ requires fewer than $3 n \mbox{ log}_2 n$ floating-point operations. For $n = 2^{20}$, that's a factor of almost 35,000 faster than $2 n^2$. Even if $n = 2^{10}$, the factor is about 70.

With MATLAB 8.3 and a 2.7 GMHz Intel Core i7 laptop, the time required for fft(x) if length(x) is 2^23 = 8388608 is about 0.5 second. The built-in fft function is based on FFTW, "The Fastest Fourier Transform in the West," developed at M.I.T. by Matteo Frigo and Steven G. Johnson in 1998.

Fast Algorithm

The key to the fast FFT algorithms is the fact that the square of the $2n$ th root of unity is the $n$ th root of unity. Using complex notation,

$$ \omega = \omega_n = e^{-2 \pi i/n} $$

we have

$$ \omega_{2n}^2 = \omega_n $$

The derivation of the fast algorithm starts with the definition of the finite Fourier transform:

$$ Y_k = \sum_{j=0}^{n-1} \omega^{jk} y_j, \ k = 0, \ldots, n-1 $$

Divide the sum into terms with even subscripts and terms with odd subscripts.

$$ Y_k = \sum_{even \ j} \omega^{jk} y_j \ + \ \sum_{ odd \ j} \omega^{jk} y_j $$

Assume that $n$ is even, that $k \leq n/2-1$ and repurpose the summation subscripts.

$$ Y_k = \sum_{j=0}^{n/2-1} \omega^{2jk} y_{2j} \ + \ \omega^k \sum_{j=0}^{n/2-1} \omega^{2jk} y_{2j+1} $$

The two sums are components of the FFTs of length $n/2$ of the portions of $y$ with even and odd subscripts. In order to get the entire FFT of length $n$, we have to do two FFTs of length $n/2$, multiply one of these by powers of $\omega$, and concatenate the results.

The relationship between an FFT of length $n$ and two FFTs of length $n/2$ can be expressed compactly in MATLAB.

  omega = exp(-2*pi*i/n);
  k = (0:n/2-1)';
  w = omega .^ k;
  u = fft(y(1:2:n-1));
  v = w.*fft(y(2:2:n));

then

  fft(y) = [u+v; u-v];

Now, if $n$ is not only even but actually a power of 2, the process can be repeated. The FFT of length $n$ is expressed in terms of two FFTs of length $n/2$, then four FFTs of length $n/4$, then eight FFTs of length $n/8$, and so on until we reach $n$ FFTs of length one. An FFT of length one is just the number itself. If $n = 2^p$, the number of steps in the recursion is $p$. There is $O(n)$ work at each step, so the total amount of work is

$$ O(np) = O(n \mbox{ log}_2 n) $$

If $n$ is not a power of two, it is still possible to express the FFT of length $n$ in terms of several shorter FFTs. An FFT of length 100 is two FFTs of length 50, or four FFTs of length 25. An FFT of length 25 can be expressed in terms of five FFTs of length 5. If $n$ is not a prime number, an FFT of length $n$ can be expressed in terms of FFTs whose lengths divide $n$. Even if $n$ is prime, it is possible to embed the FFT in another whose length can be factored. We do not go into the details of these algorithms here.

ffttx

The textbook function ffttx combines two basic ideas. If $n$ is a power of 2, it uses the $O(n \mbox{ log}_2 n)$ fast algorithm. If $n$ has an odd factor, it uses the fast recursion until it reaches an odd length, then sets up the discrete Fourier matrix and uses matrix-vector multiplication.

ffttx is one of my favorite MATLAB functions. Here is the entire code.

  function y = ffttx(x)
  %FFTTX Textbook Fast Finite Fourier Transform.
  % FFTTX(X) computes the same finite Fourier transform
  % as FFT(X).  The code uses a recursive divide and
  % conquer algorithm for even order and matrix-vector
  % multiplication for odd order.  If length(X) is m*p
  % where m is odd and p is a power of 2, the computational
  % complexity of this approach is O(m^2)*O(p*log2(p)).

   x = x(:);
   n = length(x);
   omega = exp(-2*pi*i/n);

   if rem(n,2) == 0
      % Recursive divide and conquer
      k = (0:n/2-1)';
      w = omega .^ k;
      u = ffttx(x(1:2:n-1));
      v = w.*ffttx(x(2:2:n));
      y = [u+v; u-v];
   else
      % The Fourier matrix
      j = 0:n-1;
      k = j';
      F = omega .^ (k*j);
      y = F*x;
   end

Fourier Matrix

Here is a plot of the Fourier matrix of order 16. Plots like this are the subject of the next post in this series. For now, just notice that this is not a complete graph with 16 nodes.

   plot(fft(eye(16)))
   axis equal
   axis off

References

Steve Eddins has written extensively about the FFT in his blog. Select "Categories more" and then click on "Fourier transforms". <https://blogs.mathworks.com/steve>

James W. Cooley and John W. Tukey, "An algorithm for the machine calculation of complex Fourier series", Math. Comput. 19: 297-301, 1965. <http://dx.doi.org/10.2307%2F2003354>

Matteo Frigo and Steven G. Johnson, "FFTW, Fastest Fourier Transform in the West", <http://www.fftw.org>.

Get the MATLAB code

Published with MATLAB® R2014a