An example of Discrete Fourier Transform using FFT Algorithm

First, we will explore the Theory behind the FFT (Fast Fourier Transform) Algorithm, applied to a 1D array. We will show that it is faster than the traditional method using an easy example.

Consider a sequence of numbers: \(x_0, x_1, \ldots, x_{N-1}\). Recall that the DFT (Discrete Fourier Transform) is defined by:

\[X_k=\sum_{m=1}^{N-1} x_m e^{-\frac{2 \pi i}{N} m k}\]

for \(k\) ranging from \(0\) to \(N-1\).

We can decompose the array into two, according to even-indexed elements and odd-indexed elements. The idea is that DFTs of even-indexed inputs and odd-indexed inputs are computed separated then combined to provide the DFT of the whole sequence.

The decomposition can be written as:

\[X_k=\sum_{m=0}^{N / 2-1} x_{2 m} e^{-\frac{2 \pi i }{N}(2 m) k}+\sum_{m=1}^{N / 2-1} x_{2 m+1} e^{-\frac{2 \pi i}{N}(2 m+1) k}.\]

Hence,

\[X_k=\underbrace{\sum_{m=0}^{N / 2-1} x_{2 m } e^{-\frac{2 \pi i }{N / 2} m k}}_{\text {Even Part} }+e^{-\frac{2 \pi i}{N} k} \underbrace{\sum_{m=0}^{N / 2-1} x_{2 m+1} e^{-\frac{2 \pi i}{N / 2} m k}}_{\text {Odd Part}}=E_k+e^{-\frac{2 \pi i}{N} k} O_k \quad \text { for } k=0, \ldots, \frac{N}{2}-1.\]

Similarly:

\[\begin{aligned} X_{k+\frac{N}{2}} & =\sum_{m=0}^{N / 2-1} x_{2 m} e^{-\frac{2 \pi i}{N / 2} m \left(k+\frac{N}{2}\right)}+e^{-\frac{2 \pi i}{N}\left(k+\frac{N}{2}\right)} \sum_{m=0}^{N / 2-1} x_{2 m+1} e^{-\frac{2 \pi i}{N / 2} m\left(k+\frac{N}{2}\right)} \\ & =\sum_{m=0}^{N / 2-1} x_{2 m} e^{-\frac{2 \pi i}{N / 2} m k} e^{-2 \pi m i}+e^{-\frac{2 m i}{N} k} e^{-\pi i} \sum_{m=0}^{N / 2-1} x_{2 m+1} e^{-\frac{2 \pi i}{N / 2} m k} e^{-2 \pi m i} \\ & =\sum_{m=0}^{N / 2-1} x_{2 m} e^{-\frac{2 m i}{N / 2} m k}-e^{-\frac{2 m i}{N} k} \sum_{m=0}^{N / 2-1} x_{2 m+1} e^{-\frac{2 m i}{N / 2} m k} \\ & =E_k-e^{-\frac{2 m}{N} k} O_k. \end{aligned}\]

To conclude, we have:

\[\begin{aligned} X_k & =E_k+e^{-\frac{2 N i}{N} k} O_k \\ X_{k+\frac{N}{2}} & =E_k-e^{-\frac{2 \pi i}{N} k} O_k. \end{aligned}\]

Example

We want to Fourier Transform this array:

\([1,2,3,4]\).

How can we do it?

Calculation using Traditional Method

We can calculate the DFT using:

\[X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-2\pi i \frac{kn}{N}} \quad \text{for} \; k = 0, 1, 2, 3\]

Where the answer can be expressed in:

\([X_0, X_1, X_2, X_3]\).

For \(X_0\):

\[X_0 = \sum_{n=0}^{3} x_n \cdot e^{-2\pi i \frac{0 \cdot n}{4}} = 1 \cdot e^{0} + 2 \cdot e^{0} + 3 \cdot e^{0} + 4 \cdot e^{0} = 1 + 2 + 3 + 4 = 10.\]

For \(X_1\):

\[\begin{align*} X_1 &= 1 \cdot e^{-2\pi i \frac{1 \cdot 0}{4}} + 2 \cdot e^{-2\pi i \frac{1 \cdot 1}{4}} + 3 \cdot e^{-2\pi i \frac{1 \cdot 2}{4}} + 4 \cdot e^{-2\pi i \frac{1 \cdot 3}{4}} \\ &= 1 \cdot 1 + 2 \cdot e^{-i \frac{\pi}{2}} + 3 \cdot e^{-i \pi} + 4 \cdot e^{-i \frac{3\pi}{2}} \\ &= 1 + 2 \cdot (-i) + 3 \cdot (-1) + 4 \cdot i \\ &= 1 - 2i - 3 + 4i \\ &= -2 + 2i. \end{align*}\]

For \(X_2\):

\[\begin{align*} X_2 &= 1 \cdot e^{-2\pi i \frac{2 \cdot 0}{4}} + 2 \cdot e^{-2\pi i \frac{2 \cdot 1}{4}} + 3 \cdot e^{-2\pi i \frac{2 \cdot 2}{4}} + 4 \cdot e^{-2\pi i \frac{2 \cdot 3}{4}} \\ &= 1 \cdot 1 + 2 \cdot e^{-i \pi} + 3 \cdot e^{-i 2\pi} + 4 \cdot e^{-i 3\pi} \\ &= 1 + 2 \cdot (-1) + 3 \cdot 1 + 4 \cdot (-1) \\ &= 1 - 2 + 3 - 4 \\ &= -2. \end{align*}\]

For \(X_3\):

\[\begin{align*} X_3 &= 1 \cdot e^{-2\pi i \frac{3 \cdot 0}{4}} + 2 \cdot e^{-2\pi i \frac{3 \cdot 1}{4}} + 3 \cdot e^{-2\pi i \frac{3 \cdot 2}{4}} + 4 \cdot e^{-2\pi i \frac{3 \cdot 3}{4}} \\ &= 1 \cdot 1 + 2 \cdot e^{-i \frac{3\pi}{2}} + 3 \cdot e^{-i 3\pi} + 4 \cdot e^{-i \frac{9\pi}{2}} \\ &= 1 + 2 \cdot i + 3 \cdot (-1) + 4 \cdot (-i) \\ &= 1 + 2i - 3 - 4i \\ &= -2 - 2i. \end{align*}\]

So the answer is:

\([10, -2+2i, -2, -2-2i]\).

Calculation using FFT

Even-indexed elements: \(E_k = [x_0, x_2] = [1, 3]\)

Odd-indexed elements: \(O_k = [x_1, x_3] = [2, 4]\)

Even Part

DFT for \(k = 0\):

\[E_0 = 1 \cdot e^{-2\pi i \cdot 0 \cdot 0 / 2} + 3 \cdot e^{-2\pi i \cdot 0 \cdot 1 / 2} = 1 + 3 = 4.\]

DFT for \(k = 1\):

\[E_1 = 1 \cdot e^{-2\pi i \cdot 1 \cdot 0 / 2} + 3 \cdot e^{-2\pi i \cdot 1 \cdot 1 / 2} = 1 + 3 \cdot e^{-i\pi} = 1 - 3 = -2.\]

Thus:

\[E_k = [4, -2]\]

Odd Part

DFT for \(k = 0\):

\[O_0 = 2 \cdot e^{-2\pi i \cdot 0 \cdot 0 / 2} + 4 \cdot e^{-2\pi i \cdot 0 \cdot 1 / 2} = 2 + 4 = 6.\]

DFT for \(k = 1\):

\[O_1 = 2 \cdot e^{-2\pi i \cdot 1 \cdot 0 / 2} + 4 \cdot e^{-2\pi i \cdot 1 \cdot 1 / 2} = 2 + 4 \cdot e^{-i\pi} = 2 - 4 = -2.\]

Thus:

\[O_k = [6, -2]\]

Combining Results

We use the results of the smaller DFTs to compute the DFT of the original sequence.

Before that, we shall use \(W_n=e^{-2 \pi i/ n}\) for convenience.

Then, \(W^k_n=e^{-2 \pi i k/ n}\).

This implies that:

\[X_k = E_k + W_N^k \cdot O_k.\] \[X_{k+N/2} = E_k - W_N^k \cdot O_k.\]

Now Recall:

\[W_4 = e^{-2\pi i / 4} = i.\]

For \(k = 0\):

\[X_0 = E_0 + W_4^0 \cdot O_0 = 4 + 1 \cdot 6 = 10.\]

For \(k = 1\):

\[X_1 = E_1 + W_4^1 \cdot O_1 = -2 + i \cdot (-2) = -2 - 2i.\]

For \(k = 2\):

\[X_2 = E_0 - W_4^0 \cdot O_0 = 4 - 1 \cdot 6 = -2.\]

For \(k = 3\):

\[X_3 = E_1 - W_4^1 \cdot O_1 = -2 - i \cdot (-2) = -2 + 2i.\]

So, the DFT of the sequence \([1, 2, 3, 4]\) using the FFT algorithm is:

\[[10, -2 - 2i, -2, -2 + 2i].\]

Which agrees with the calculation above!

It can be seen that the Fast Fourier Transform is indeed much faster than the traditional method. This method is very useful in image processing, where we convolve the “pixels” in an image with a filter to blur the image or sharpen the edges. More can be seen in this post.

Credits

The first part is from the Wikipedia Page on Cooley-Tukey FFT Algorithm.