MIT 6.300 — Lecture 14

Fast Fourier Transform

See Lectures 22 and 23 of 6.1220 for an alternative perspective on the content from this lecture.

Remark. Note that 6.1220 adopts slightly different conventions than 6.300.

In 6.300, the DFT is $X[k] := \frac{1}{N} \sum_{n = 0}^{N - 1} x[n] e^{-j\frac{2\pi k}{N}n}$ ; refer to Lecture 12.
In 6.1220, the DFT is $X[k] := \sum_{n = 0}^{N - 1} x[n] e^{j\frac{2\pi k}{N}n}$ ; simpler for abstract purposes, but not as meaningful in physical contexts.

§ Recitation: DFT = Polynomial Sampling

We first provide an alternative perspective on the DFT, as suggested by 6.1220.

Theorem. (DFT = Polynomial Sampling) Consider a DT signal $x[n]$ of length $N$ along with the corresponding polynomial $p(z) := \sum_{n = 0}^{N - 1} x[n] z^n$ . Let $\omega = e^{2\pi j / N}$ be a primitive $N^{\text{th}}$ root of unity. Then:

\mathrm{DFT}(x)[k] = \frac{1}{N}p(\omega^{-k}) ~~~ \text{ and } ~~~ \mathrm{DFT}^{-1}(x)[n] = p(\omega^{n}).

Proof: Follows immediately from the DFT analysis and synthesis equations from Lecture 12. $\blacksquare$

The above reinterpretation makes it much easier to argue that the DFT of a circular convolution equals the product of the DFTs.

Theorem. (Convolution = Product) Given DT signals $a[n]$ and $b[n]$ of length $N$ with DFTs $A[k]$ and $B[k]$ , respectively, we have:

\mathrm{DFT}(a \circledast b) = N \cdot A \cdot B ~~~ \text{ and } ~~~ \mathrm{DFT}^{-1}(A \circledast B) = a \cdot b.

Proof: Say $a[n]$ and $b[n]$ correspond to polynomials $p(z)$ and $q(z)$ . The point is that the polynomial corresponding to $(a \circledast b)[n]$ is simply $p(z) \cdot q(z)$ in $\mathbb{C}[z] / (z^N - 1)$ : convolution is just polynomial multiplication. Thus,

\mathrm{DFT}(a \circledast b) = \left[\frac{1}{N} \cdot p(\omega^{-k}) \cdot q(\omega^{-k})\right]_{k = 0}^{N - 1} = N \cdot \left[\frac{1}{N} \cdot p(\omega^{-k})\right]_{k = 0}^{N - 1} \cdot \left[\frac{1}{N} \cdot q(\omega^{-k})\right]_{k = 0}^{N - 1} = N \cdot \mathrm{DFT}(a) \cdot \mathrm{DFT}(b).

The argument for the inverse DFT works the same way. $\blacksquare$

A corollary of this reinterpretation is the following.

Theorem. (DFT = Inverse DFT) If $X[k]$ is the DFT of $x[n]$ , then the DFT of $X[k]$ is $x'[n] := \frac{1}{N} \cdot x[(-n)\%N]$ . Equivalently, writing $\mathrm{Reverse}$ for circular reversal,

\mathrm{DFT}(\mathrm{DFT}(x)) = \frac{1}{N} \cdot \mathrm{Reverse}(x).

Proof: From 6.1220, the Lectures 22 and 23 notes provide a linear-algebra-based argument from the fundamentals. But we don't need that: this follows immediately from the relationship between the DFT and inverse DFT in the polynomial setting. $\blacksquare$

Remark. All of the mathematical content of the above proof has, in some sense, been deferred to the mathematical content of the consistency of the synthesis and analysis equations from Lecture 12.

Thus, any result we discover about the DFT has a natural “dual” in the context of the inverse DFT.

§ Lecture: The FFT Algorithm

Via a divide-and-conquer approach, we may compute the DFT of a DT signal of length $N$ in $O(N \log N)$ time. This algorithm is known as the Fast Fourier Transform (FFT).

from cmath import exp, pi

def FFT(x):
    """Returns the DFT of x, an array whose length is a power of two."""

    # Base Case
    N = len(x)
    if N == 1:
        return x

    # Divide and Conquer
    even_x = x[::2]
    odd_x = x[1::2]
    FFT_even_x = FFT(even_x)
    FFT_odd_x = FFT(odd_x)

    # Combine Step
    FFT_x = []
    for k in range(N // 2):
        FFT_x.append((FFT_even_x[k] + exp(-2j * pi * k/N) * FFT_odd_x[k]) / 2)
    for k in range(N // 2):
        FFT_x.append((FFT_even_x[k] - exp(-2j * pi * k/N) * FFT_odd_x[k]) / 2)

    return FFT_x

Refer to Lectures 22 and 23 for details on how this algorithm works; keep in mind the difference in convention.

An implementation of the inverse DFT would look very similar, except with the exponents negated and the factor of $\frac{1}{2}$ removed. Here's how that looks in code.

# DFT Implementation:
FFT_x.append((FFT_even_x[k] + exp(-2j * pi * k/N) * FFT_odd_x[k]) / 2) # first for loop
FFT_x.append((FFT_even_x[k] - exp(-2j * pi * k/N) * FFT_odd_x[k]) / 2) # second for loop

# Inverse DFT Implementation:
FFT_x.append((FFT_even_x[k] + exp(2j * pi * k/N) * FFT_odd_x[k])) # first for loop
FFT_x.append((FFT_even_x[k] - exp(2j * pi * k/N) * FFT_odd_x[k])) # second for loop