Fast Fourier Transform (FFT)

What is the Fast Fourier Transform?

The Fast Fourier Transform (FFT) is an algorithm that computes the Discrete Fourier Transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa.

We can take advantage of the symmetry and periodicity of the complex exponential kernel in the DFT to compute the DFT in $\mathcal{O}[N\log N]$ time, where $N$ is the number of points in the sequence. This is a significant improvement over the $\mathcal{O}[N^2]$ time complexity of the naive DFT algorithm.

How does the FFT work?

As discussed in class, the FFT is a fast, $\mathcal{O}[N\log N]$ algorithm to compute the DFT of a sequence, or its inverse. The DFT is defined as:

$$F(k) = \sum_{n=0}^{N-1} f(n) e^{-i2\pi kn/N}$$

where $f(n)$ is the input sequence, $F(k)$ is the transformed sequence, and $N$ is the number of points in the sequence.

FFT Algorithm (Radix-2 DIT)

Overview

The FFT algorithm computes the DFT in $\mathcal{O}[N\log N]$ time by using a divide-and-conquer approach. The algorithm recursively divides the input sequence into two smaller sequences, computes their DFTs, and then combines them to form the DFT of the original sequence. The algorithm is based on the observation that the DFT of a sequence of even length $N$ can be expressed in terms of the DFTs of its subsequences of length $N/2$.

Fun Fact

The algorithm, along with its recursive application, was invented by Carl Friedrich Gauss. Cooley and Tukey independently rediscovered and popularized it 160 years later.

Data Flow Diagram for $N=8$ (Radix-2 DIT)

The data flow diagram above¹ shows the computation of the DFT of a sequence of length $N=8$ using the Radix-2 Decimation-In-Time (DIT) FFT algorithm. The diagram illustrates the divide-and-conquer approach of the algorithm, where the input sequence is divided into two smaller sequences, and their DFTs are computed recursively. The diagram also shows the "butterfly" operations, which combine the DFTs of the smaller sequences to form the DFT of the original sequence.

Understanding the $8^{th}$ Root ($W$)

The $8^{th}$ root of unity, denoted as $W$, is fundamental to the FFT algorithm. It's defined as $W = e^{-i2\pi/8}$, where $i$ is the imaginary unit. This root of unity represents the complex numbers on the unit circle in the complex plane, equidistantly spaced at angles corresponding to their fractional part of a full circle (in this case, $2\pi/8$, or $\pi/4$, radians apart).

For $N=8$, we have:

• $W^0 = e^{0} = 1$
• $W^1 = e^{-i\pi/4} = \frac{1 - i}{\sqrt{2}}$
• $W^2 = e^{-i\pi/2} = -i$
• $W^3 = e^{-3i\pi/4} = \frac{-1 - i}{\sqrt{2}}$
• $W^4 = e^{-i\pi} = -1$
• $W^5 = e^{-5i\pi/4} = \frac{-1 + i}{\sqrt{2}}$
• $W^6 = e^{-3i\pi/2} = i$
• $W^7 = e^{-7i\pi/4} = \frac{1 + i}{\sqrt{2}}$

These values of $W^k$ illustrate the symmetry and periodicity inherent in the DFT and are critical for the "butterfly" calculations in the FFT algorithm.

Matrix Representation

The DFT matrix for $N=8$, incorporating the simplified $W^k$ values, provides a clear visualization of the FFT's operation:

$$\begin{align*} \begin{bmatrix} A_0 \\ A_1 \\ A_2 \\ A_3 \\ A_4 \\ A_5 \\ A_6 \\ A_7 \\ \end{bmatrix} &= \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & W & W^2 & W^3 & W^4 & W^5 & W^6 & W^7 \\ 1 & W^2 & W^4 & W^6 & 1 & W^2 & W^4 & W^6 \\ 1 & W^3 & W^6 & W^1 & W^4 & W^7 & W^2 & W^5 \\ 1 & W^4 & 1 & W^4 & 1 & W^4 & 1 & W^4 \\ 1 & W^5 & W^2 & W^7 & W^4 & W & W^6 & W^3 \\ 1 & W^6 & W^4 & W^2 & 1 & W^6 & W^4 & W^2 \\ 1 & W^7 & W^6 & W^5 & W^4 & W^3 & W^2 & W \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7 \\ \end{bmatrix} \end{align*}$$

This representation, with the explicit inclusion of $W^k$ values, elucidates the FFT's efficiency in decomposing the DFT into smaller, manageable components.

Butterfly Operations Explained

Butterfly operations are the core computational elements of the FFT algorithm. They represent the pairwise combination of DFT outputs from smaller sequences. A "butterfly" operation takes two complex numbers, multiplies one of them by a twiddle factor (a power of $W$), and then performs addition and subtraction to produce two new complex numbers.

For a single stage in an $N=8$ FFT, a butterfly operation on inputs $a$ and $b$ with a twiddle factor $W^k$ can be represented as:

• Output 1: $a + W^k \cdot b$
• Output 2: $a - W^k \cdot b$

These operations halve the computation required at each stage by reusing results, demonstrating the FFT's divide-and-conquer strategy. The recursive application of butterfly operations across the stages of the FFT leads to its $\mathcal{O}[N\log N]$ complexity, a significant improvement over the $\mathcal{O}[N^2]$ complexity of the direct DFT computation.

Pseudocode

The following pseudocode describes the FFT algorithm for a sequence of length $N = 2^m$:

function FFT(f, N, s):
    if N == 1:
        return f
    else:
        f_even = [f[2*n] for n in range(N//2)]
        f_odd = [f[2*n + 1] for n in range(N//2)]
        F_even = FFT(f_even, N//2, 2*s)
        F_odd = FFT(f_odd, N//2, 2*s)
        F = [0]*N
        for k in range(N//2):
            F[k] = F_even[k] + exp(-2j*pi*k/N)*F_odd[k]
            F[k + N//2] = F_even[k] - exp(-2j*pi*k/N)*F_odd[k]
        return F

The function FFT takes as input a sequence f of length N and a stride s (which is used to compute the twiddle factors). If the length of the sequence is 1, the function returns the sequence itself. Otherwise, it divides the sequence into two subsequences of even and odd indices, computes their DFTs recursively, and then combines them to form the DFT of the original sequence².

Rough Implementation

The following example demonstrates how to use the FFT algorithm to compute the DFT of a sequence:

def FFT(f, N, s):
    if N == 1:
        return f
    else:
        f_even = [f[2*n] for n in range(N//2)]
        f_odd = [f[2*n + 1] for n in range(N//2)]
        F_even = FFT(f_even, N//2, 2*s)
        F_odd = FFT(f_odd, N//2, 2*s)
        F = [0]*N
        for k in range(N//2):
            F[k] = F_even[k] + exp(-2j*pi*k/N)*F_odd[k]
            F[k + N//2] = F_even[k] - exp(-2j*pi*k/N)*F_odd[k]
        return F

# Example usage
f = [1, 2, 3, 4]
N = len(f)
s = 1
F = FFT(f, N, s)
print(F)

How to use the FFT in Python

In your homework, you are asked to implement the FFT. I recommend, as a sanity check, to verify your implementation using the built-in, or in this case, the numpy function. The FFT can be used in Python by using the numpy.fft module. The following code snippet demonstrates how to use the FFT to compute the DFT of a sequence:

import numpy as np
import matplotlib.pyplot as plt

# Number of sample points
N = 600

# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N, endpoint=False)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = np.fft.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)

plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.grid()
plt.show()

Fourier Transforms of Images

As we discussed before, the two-dimensional discrete Fourier transform is a simple generalization of the standard 1-D DFT:

$$F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-i2\pi(\frac{ux}{M} + \frac{vy}{N})}$$

where $f(x, y)$ is the input image, $F(u, v)$ is the transformed image, and $M$ and $N$ are the dimensions of the image. This is good for rectangular images whose dimensions are not powers of 2. If you use the DFT with a rectangular image, the cost is $\mathcal{O}[N^4]$.

For an $N \times N$ image, $N$ even, the cost of a 2D FFT is proportional to $N^2 \log N$. Let's prove this.

Proof: The 2D FFT is computed by performing a 1D FFT on each row of the image, followed by a 1D FFT on each column of the result. The cost of a 1D FFT is $\mathcal{O}[N\log N]$, so the cost of a 2D FFT is $\mathcal{O}[N^2\log N]$.

Problem 1

Compute the DFT of the following sequence using the FFT algorithm:

$$x = [-j, 0, j, 1]$$

Just a Heads Up

This is slightly different than the homework problem, so don't go copying this exactly. Nonetheless, this should give you a solid start!

Solution

To compute the Discrete Fourier Transform (DFT) of the sequence $x = [-j, 0, j, 1]$ using the Fast Fourier Transform (FFT) algorithm, we'll follow the Cooley-Tukey FFT algorithm's radix-2 decimation-in-time approach, assuming $N=4$, which is a power of 2.

The FFT algorithm decomposes the DFT of a sequence into smaller DFTs to reduce the computational complexity. The decomposition follows:

1. Divide the sequence into two sequences: one containing the even-indexed elements and the other containing the odd-indexed elements.
2. Compute the DFT of these smaller sequences.
3. Combine these DFTs to get the DFT of the original sequence.

Given $x = [-j, 0, j, 1]$, let's denote the DFT of $x$ as $X[k]$, where $k = 0, 1, 2, 3$.

Step 1: Divide

• Even-indexed elements: $x_{even} = [-j, j]$
• Odd-indexed elements: $x_{odd} = [0, 1]$

Step 2: Compute DFT of Smaller Sequences

• For $N=2$, the DFT can be computed directly:
  • DFT of $x_{even}$: $X_{even}[0] = -j + j = 0$, $X_{even}[1] = -j - j = -2j$
  • DFT of $x_{odd}$: $X_{odd}[0] = 0 + 1 = 1$, $X_{odd}[1] = 0 - 1 = -1$

Step 3: Combine

• Using the FFT combination step, where the twiddle factors for $N=4$ are $W_N^k = e^{-j\frac{2\pi}{N}k}$ for $k=0,1,2,3$, we get:
  • $X[0] = X_{even}[0] + W_4^0 \cdot X_{odd}[0] = 0 + 1 \cdot 1 = 1$
  • $X[1] = X_{even}[1] + W_4^1 \cdot X_{odd}[1] = -2j + e^{-j\frac{\pi}{2}} \cdot (-1) = -2j + j = -j$
  • $X[2] = X_{even}[0] - W_4^0 \cdot X_{odd}[0] = 0 - 1 \cdot 1 = -1$
  • $X[3] = X_{even}[1] - W_4^1 \cdot X_{odd}[1] = -2j - e^{-j\frac{\pi}{2}} \cdot (-1) = -2j - j = -3j$

Therefore, the DFT $X[k]$ of the sequence $x = [-j, 0, j, 1]$ using the FFT algorithm is:

$$X = [1, -j, -1, -3j]$$

This solution demonstrates the basic principles behind the FFT algorithm's efficiency: by breaking down the DFT computation into smaller parts and reusing results, we significantly reduce the computational effort compared to directly computing the DFT. Once again, if I made any mistakes, please let me know!

Problem 2

In image processing, the FFT is used to transform an image from the spatial domain to the frequency domain. The following code snippet demonstrates how to use the FFT to transform an image:

import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack

# Load the Lena image
lena = plt.imread('/path/to/lena.png')


# Compute the 2D FFT of the image
lena_fft = 


# Plot the magnitude of the FFT
plt.imshow(np.abs(lena_fft), cmap='gray')
plt.colorbar()
plt.show()

Using this example code, compute the 2D FFT of the Lena image and plot the magnitude of the FFT. What do you observe?

Footnotes

1. https://en.wikipedia.org/wiki/Cooley--Tukey_FFT_algorithm ↩

2. E. Oran Brigham, The Fast Fourier Transform, Prentice Hall, Englewood Cliffs, NJ, 1974. ↩