Fast Fourier Transformation (FFT) algorithm

We have seen that the brute force implementation of the discrete Fourier transform scales with the square of the number of data (\(N^2\)) which makes Fourier transforming large data sets extremely time-consuming. However, there is a very clever trick first discovered by Gauss that allows to perform the discrete Fourier transformation much more efficiently and leads to the \( \mathcal{O}(N\times \log N)\) scaling.

The resulting algorithm is called Fast Fourier Transform (FFT) and is a basis for a tremendous number of technological and scientific approximations.

The FFT-algorithm is considered to be the most far-reaching algorithm in computing and has revolutionalized signal processing. Although originally introduced by Gauss who tried to calculate orbits of asteroid Juno, it was rediscovered in the 20th century by the American mathematicians Cooley and Tukey.

In the following, we develop and implement one variant of the Cooley-Tukey's algorithm.

We recall once more the definition of the discrete Fourier transform of a vector \(\mathbf{f}\) and rewrite it in a slightly modified form:

$$NF_{k}=\sum_{j=0}^{N-1}f_{j}W_{N}^{kj}\quad k=0,\cdots N-1$$

with

$$ W_{N}=\exp\left[-\frac{2\pi \mathrm{i}}{N}\right] $$

Most variants of the FFT-algorithm are formulated in the simplest way if we assume that the number of data points N is a power of two (\(N=2^{n}\)).

In this case, the above sum can be also rewritten by summing only over the first half of the summation indices such that:

$$ NF_{k}=\sum_{j=0}^{N/2-1}\left(f_{j}W_{N}^{kj}+f_{j+\frac{N}{2}}W_{N}^{\left(j+\frac{N}{2}\right)k}\right) $$

By utilizing the properties of the complex exponentials we can show that:

$$ W_{N}^{\left(k+\frac{N}{2}\right)j}=W_{N}^{kj}W_{N}^{\frac{N}{2}k} $$

Furthermore, $$ W_{N}^{\frac{N}{2}k}=\exp\left[-\frac{2\pi \mathrm{i}}{N}\times\frac{N}{2}\right]^{k}=\left(\exp(-\mathrm{i}\pi)\right)^{k}=(-1)^{k} $$

With this, we can rewrite the discrete Fourier transform as:

$$NF_{k}=\sum_{j=0}^{\frac{N}{2}-1}\left[f_{j}+(-1)^{k}f_{j+\frac{N}{2}}\right]W_{N}^{kj}$$

Now, it is natural to separate the Fourier coefficients into the ones with odd and even indices. With this, we can rewrite the discrete Fourier transform as:

$$NF_{l}^{e}=NF_{2l}=\sum_{j=0}^{\frac{N}{2}-1}\left[f_{j}+(-1)^{2l}f_{j+\frac{N}{2}}\right]W_{N}^{2lj}\quad l=1,\cdots\frac{N}{2}$$

$$NF_{l}^{o}=NF_{2l+1}=\sum_{j=0}^{\frac{N}{2}-1}\left[f_{j}+(-1)^{2l+1}f_{j+\frac{N}{2}}\right]W_{N}^{\left(2l+1\right)j}\quad l=1,\cdots\frac{N}{2}$$

By noting that: \(W_{N}^{2lj}=W_{N/2}^{lj}\) and \(W_{N}^{\left(2l+1\right)j}=W_{N/2}^{lj}W_{N}^{j}\)

we can further simplify the above expressions by introducing wo auxiliary arrays:

$$g_{j}=f_{j}+f_{j+\frac{N}{2}}\quad j=1,\cdots,\frac{N}{2}$$

$$h_{j}=\left[f_{j}-f_{j+\frac{N}{2}}\right]W_{N}^{j}\quad j=1,\cdots,\frac{N}{2}.$$

which gives us for the even and odd components of the Fourier coefficients the following expression:

$$NF_{l}^{e}= \sum_{j=0}^{\frac{N}{2}-1}g_{j}W_{N/2}^{lj}$$ $$NF_{l}^{o}= \sum_{j=0}^{\frac{N}{2}-1}h_{j}W_{N/2}^{lj}$$

The central observation now is to realize that these are just two new discrete Fourier transforms with half as many of the data points. This leads to a conceptually very simple recursive FFT-algorithm!