Fast Fourier transform

These notes derive the radix-2 Fast Fourier Transform (FFT) by viewing the discrete Fourier transform (DFT) as a linear operator on an $N$-dimensional Hilbert space. The key steps are: an even–odd decomposition of the Hilbert space, expressing the $N$-point DFT in terms of two $M$-point DFTs ($N = 2M$), and writing an explicit operator factorization in Dirac notation, followed by a block-matrix form and explicit $N=4$ and $N=8$ examples.

II. Even–odd decomposition of the Hilbert space

Assume $N = 2M$. We identify

\[ H_N \simeq H_2 \otimes H_M, \]

by

\[ \lvert 2m\rangle \leftrightarrow \lvert 0\rangle \otimes \lvert m\rangle, \qquad \lvert 2m+1\rangle \leftrightarrow \lvert 1\rangle \otimes \lvert m\rangle, \qquad m = 0,\dots,M-1. \]

A general input vector is

\[ \lvert x\rangle = \sum_{n=0}^{N-1} x_n \lvert n\rangle = \sum_{m=0}^{M-1} x_{2m}\lvert 2m\rangle + \sum_{m=0}^{M-1} x_{2m+1}\lvert 2m+1\rangle. \]

Using the tensor-product structure

\[ \lvert 2m\rangle = \lvert 0\rangle \otimes \lvert m\rangle, \qquad \lvert 2m+1\rangle = \lvert 1\rangle \otimes \lvert m\rangle, \]

we obtain

\[ \lvert x\rangle = \lvert 0\rangle \otimes \left( \sum_{m=0}^{M-1} x_{2m}\lvert m\rangle \right) + \lvert 1\rangle \otimes \left( \sum_{m=0}^{M-1} x_{2m+1}\lvert m\rangle \right). \]

Define the even and odd subvectors

\[ \lvert x^{(e)}\rangle := \sum_{m=0}^{M-1} x_{2m}\lvert m\rangle, \qquad \lvert x^{(o)}\rangle := \sum_{m=0}^{M-1} x_{2m+1}\lvert m\rangle. \]

Then

\[ \lvert x\rangle = \lvert 0\rangle \otimes \lvert x^{(e)}\rangle + \lvert 1\rangle \otimes \lvert x^{(o)}\rangle. \]

IV. Operator factorization in Dirac notation

We now express the above procedure as a factorization of the operator $\hat{F}_N$ on $H_N \simeq H_2 \otimes H_M$.

A. Two smaller Fourier transforms

Define

\[ \hat{G} := \hat{I}_2 \otimes \hat{F}_M. \]

Acting on the decomposed state,

\[ \begin{aligned} \hat{G}\lvert x\rangle &= (\hat{I}_2 \otimes \hat{F}_M) \left( \lvert 0\rangle \otimes \lvert x^{(e)}\rangle + \lvert 1\rangle \otimes \lvert x^{(o)}\rangle \right) \\ &= \lvert 0\rangle \otimes \hat{F}_M\lvert x^{(e)}\rangle + \lvert 1\rangle \otimes \hat{F}_M\lvert x^{(o)}\rangle \\ &= \lvert 0\rangle \otimes \lvert E\rangle + \lvert 1\rangle \otimes \lvert O\rangle. \end{aligned} \]

B. Twiddle-factor operator

Define

\[ \hat{D} := \sum_{k=0}^{M-1} \omega_N^k \lvert k\rangle\langle k\rvert \]

on $H_M$, and the operator

\[ \hat{T} := \lvert 0\rangle\langle 0\rvert \otimes \hat{I}_M + \lvert 1\rangle\langle 1\rvert \otimes \hat{D}. \]

Then

\[ \hat{T} \left( \lvert 0\rangle \otimes \lvert E\rangle + \lvert 1\rangle \otimes \lvert O\rangle \right) = \lvert 0\rangle \otimes \lvert E\rangle + \lvert 1\rangle \otimes \hat{D}\lvert O\rangle. \]

Let

\[ \lvert O'\rangle := \hat{D}\lvert O\rangle. \]

C. Hadamard on the first qubit

Define the Hadamard operator

\[ \hat{H} = \frac{1}{\sqrt{2}} \big( \lvert 0\rangle\langle 0\rvert + \lvert 0\rangle\langle 1\rvert + \lvert 1\rangle\langle 0\rvert - \lvert 1\rangle\langle 1\rvert \big) = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}. \]

Consider $\hat{U} := \hat{H} \otimes \hat{I}_M$. Acting on

\[ \lvert \psi\rangle := \lvert 0\rangle \otimes \lvert E\rangle + \lvert 1\rangle \otimes \lvert O'\rangle, \]

\[ \begin{aligned} \hat{U}\lvert \psi\rangle &= (\hat{H} \otimes \hat{I}_M) \left( \lvert 0\rangle \otimes \lvert E\rangle + \lvert 1\rangle \otimes \lvert O'\rangle \right) \\ &= (\hat{H}\lvert 0\rangle) \otimes \lvert E\rangle + (\hat{H}\lvert 1\rangle) \otimes \lvert O'\rangle \\ &= \frac{1}{\sqrt{2}} (\lvert 0\rangle + \lvert 1\rangle) \otimes \lvert E\rangle + \frac{1}{\sqrt{2}} (\lvert 0\rangle - \lvert 1\rangle) \otimes \lvert O'\rangle \\ &= \frac{1}{\sqrt{2}} \left[ \lvert 0\rangle \otimes (\lvert E\rangle + \lvert O'\rangle) + \lvert 1\rangle \otimes (\lvert E\rangle - \lvert O'\rangle) \right]. \end{aligned} \]

Writing

\[ \lvert E\rangle = \sum_{k=0}^{M-1} E_k \lvert k\rangle, \qquad \lvert O'\rangle = \sum_{k=0}^{M-1} \omega_N^k O_k \lvert k\rangle, \]

the coefficient of $\lvert 0\rangle \otimes \lvert k\rangle$ is

\[ \frac{1}{\sqrt{2}} (E_k + \omega_N^k O_k) = X_k, \]

and the coefficient of $\lvert 1\rangle \otimes \lvert k\rangle$ is

\[ \frac{1}{\sqrt{2}} (E_k - \omega_N^k O_k) = X_{k+M}. \]

Identifying $\lvert 0\rangle \otimes \lvert k\rangle$ with $\lvert k\rangle$ and $\lvert 1\rangle \otimes \lvert k\rangle$ with $\lvert k+M\rangle$, we see that $\hat{U}\hat{T}(\hat{I}_2 \otimes \hat{F}_M)$ reproduces the DFT amplitudes.

Up to a permutation $\hat{P}$ that maps the standard basis $\{\lvert n\rangle\}$ to the even–odd tensor basis, we obtain

\[ \hat{F}_N = (\hat{H} \otimes \hat{I}_M)\,\hat{T}\, (\hat{I}_2 \otimes \hat{F}_M)\,\hat{P}. \]

VI. Explicit examples: $N=4$ and $N=8$

A. The case $N = 4$

For $N = 4$ we have $M = N/2 = 2$ and

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

The $2 \times 2$ Fourier matrix $F_2$ is

\[ F_2 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}. \]

The twiddle-factor matrix $D$ (of size $M \times M$) for $N = 4$ is

\[ D_2 = \mathrm{diag}(1,\omega_4^1) = \mathrm{diag}(1,-i) = \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix}. \]

We can then compute

\[ D_2 F_2 = \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ -i & i \end{pmatrix}. \]

In the even–odd basis, where the input is ordered as $(x_0,x_2,x_1,x_3)^{T}$, the factorization

\[ F_4^{(\mathrm{even\text{-}odd})} = \frac{1}{\sqrt{2}} \begin{pmatrix} F_2 & D_2 F_2 \\ F_2 & -D_2 F_2 \end{pmatrix} \]

yields, after simplifying, the matrix often written as

\[ F_4^{(\mathrm{even\text{-}odd})} = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & -i & i \\ 1 & 1 & -1 & -1 \\ 1 & -1 & i & -i \end{pmatrix}. \]

To relate this to the standard 4-point DFT in the natural ordering $(x_0,x_1,x_2,x_3)^{T}$, introduce the permutation matrix

\[ P = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}, \qquad P \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} x_0 \\ x_2 \\ x_1 \\ x_3 \end{pmatrix}. \]

Then

\[ F_4 = F_4^{(\mathrm{even\text{-}odd})} P = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix}, \]

which is precisely the 4-point DFT matrix with the normalization $1/\sqrt{N} = 1/2$.

B. The case $N = 8$

For $N = 8$ we have $M = N/2 = 4$ and

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

The $4 \times 4$ Fourier matrix $F_4$ (in the standard ordering) is

\[ F_4 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & \omega_4 & \omega_4^2 & \omega_4^3 \\ 1 & \omega_4^2 & \omega_4^4 & \omega_4^6 \\ 1 & \omega_4^3 & \omega_4^6 & \omega_4^9 \end{pmatrix}, \qquad \omega_4 = e^{-2\pi i/4} = -i. \]

Using $\omega_4^2 = -1$ and $\omega_4^3 = i$, this can be written explicitly as

\[ F_4 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix}. \]

For $N = 8$, the twiddle-factor matrix $D_4$ is

\[ D_4 = \mathrm{diag} \big( 1,\omega_8^1,\omega_8^2,\omega_8^3 \big) = \mathrm{diag} \big( 1,e^{-i\pi/4},e^{-i\pi/2},e^{-3i\pi/4} \big). \]

In the even–odd basis for $N = 8$, the block factorization reads

\[ F_8^{(\mathrm{even\text{-}odd})} = \frac{1}{\sqrt{2}} \begin{pmatrix} F_4 & D_4 F_4 \\ F_4 & -D_4 F_4 \end{pmatrix}. \]

This $8 \times 8$ matrix acts on the input ordered as

\[ (x_0,x_2,x_4,x_6,x_1,x_3,x_5,x_7)^{T}, \]

and its block structure directly encodes the top-level radix-2 FFT step:

• $F_4$ applied independently to the even and odd subsequences,
• multiplication of the “odd” 4-point transform by $D_4$,
• a final “butterfly” layer given by the block matrix \[ \frac{1}{\sqrt{2}} \begin{pmatrix} I_4 & I_4 \\ I_4 & -I_4 \end{pmatrix}. \]

To obtain the standard 8-point DFT matrix $F_8$ in the natural ordering $(x_0,x_1,\dots,x_7)^{T}$, one composes $F_8^{(\mathrm{even\text{-}odd})}$ with the corresponding even–odd permutation matrix $P_8$:

\[ F_8 = F_8^{(\mathrm{even\text{-}odd})} P_8. \]

The structure of $F_8^{(\mathrm{even\text{-}odd})}$ shows the recursive nature of the FFT, as it is entirely built from smaller Fourier transforms ($F_4$) and diagonal twiddle-factor matrices ($D_4$), combined by a simple $2 \times 2$ block “Hadamard”.

import numpy as np

N = 8
n = np.arange(N)
k = n[:, None]
F = np.exp(-2j * np.pi * k * n / N)

x = np.sin(2*np.pi*2*n/N)
X_matrix = F @ x
X_fft = np.fft.fft(x)
print(np.max(np.abs(X_matrix - X_fft)))

def radix2_stages(x):
    """Return the bit-reversed input and every butterfly stage."""
    N = len(x)
    bits = int(np.log2(N))
    reverse = lambda i: int(f"{i:0{bits}b}"[::-1], 2)
    y = np.asarray([x[reverse(i)] for i in range(N)], complex)
    stages = [y.copy()]
    size = 2
    while size <= N:
        for start in range(0, N, size):
            for j in range(size // 2):
                w = np.exp(-2j*np.pi*j/size)
                even = y[start + j]
                odd = w * y[start + j + size//2]
                y[start + j], y[start + j + size//2] = even + odd, even - odd
        stages.append(y.copy())
        size *= 2
    return stages

← Back to Notes