In the discrete Fourier transform, we want to project a vector v∈Cnv \in \mathbb{C}^n to the orthonormal bases

e0=1n(e0⋅2⋅0⋅πin,e12⋅0⋅πin,...,e(n−1)⋅2⋅0⋅πin),...ek−1=1n(e0⋅2⋅(k−1)⋅πin,e12⋅(k−1)⋅πin,...,e(n−1)⋅2⋅(k−1)⋅πin),...en−1=1n(e0⋅2⋅(n−1)⋅πin,e12⋅(n−1)⋅πin,...,e(n−1)⋅2⋅(n−1)⋅πin). \begin{aligned} &e_0 = \frac{1}{\sqrt n} \left( e^{ 0 \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ 1 \frac{2 \cdot 0 \cdot \pi i}{ n} }, ..., e^{ (n - 1) \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} } \right), \\ &... \\ &e_{k - 1} = \frac{1}{\sqrt n} \left( e^{ 0 \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, e^{ 1 \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, ..., e^{ (n - 1) \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} } \right), \\ &... \\ &e_{n - 1} = \frac{1}{\sqrt n} \left( e^{ 0 \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, e^{ 1 \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, ..., e^{ (n - 1) \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} } \right). \end{aligned}

We use S\mathcal{S} to denote the linear span of e0,e1,...,en−1e_0, e_1, ..., e_{n - 1}.

Orthonormal Bases

It is easy to see that eke_k is a unit vector, ∀k∈Zn\forall k \in \mathbb{Z}_n. The vectors in the bases are also mutually orthogonal, as for k,ℓ∈Zn,k≠ℓk, \ell \in \mathbb{Z}_n, k \neq \ell,

<ek,eℓ>H=∑j∈Znej⋅2(ℓ−k)πin=1−en⋅2(ℓ−k)πin1−e2(ℓ−k)πin=0. \left< e_k, e_\ell \right>_H = \sum_{j \in \mathbb{Z}_n } e^{ j \cdot \frac{ 2 (\ell - k) \pi i }{ n } } = \frac{ 1 - e^{ n \cdot \frac{ 2 (\ell - k) \pi i }{ n } } }{ 1 - e^{ \frac{ 2 (\ell - k) \pi i }{ n } } } = 0.

Projection

The projection of vv to the bases, is defined as the vector in the linear span of the bases that is closest to vv. It is given by

P(v)≐∑k∈Znakek, P(v) \doteq \sum_{k \in \mathbb{Z}_n} a_k e_k,

where ak=<ek,v>H≐eˉk⋅va_k = \left< e_k, v\right>_H \doteq \bar e_k \cdot v.

Lemma. For any vector uu in S\mathcal{S}, it holds that

<u,v>H=<u,P(v)>H. \left< u, v \right>_H = \left< u, P(v) \right>_H.

Or, equivalently,

<u,v−P(v)>H=<v−P(v),u>H=0. \left< u , v - P(v) \right>_H = \left< v - P(v), u \right>_H = 0.

Proof. Let u=∑k∈Znckek:ck∈C,∀k∈Znu = \sum_{k \in \mathbb{Z}_n} c_k e_k : c_k \in \mathbb{C}, \forall k \in \mathbb{Z}_n be a vector in S\mathcal{S}. As e0,e1,...,en−1e_0, e_1, ..., e_{n - 1} are mutually orthogonal,

<u,v>H=∑k∈Zncˉk⋅ak=<u,P(v)>H. \left< u, v \right>_H = \sum_{k \in \mathbb{Z}_n} \bar c_k \cdot a_k =\left< u, P(v) \right>_H.

□\square

Next, we verify that P(v)P(v) minimizes the distance to vv.

Theorem.

P(v)=arg min⁡w∈S∣∣u−w∣∣2. P(v) = \argmin_{w \in \mathcal{S} } ||u - w||^2.

Proof. Let w=∑k∈Zndkek:dk∈C,∀k∈Znw = \sum_{k \in \mathbb{Z}_n} d_k e_k : d_k \in \mathbb{C}, \forall k \in \mathbb{Z}_n and define u=P(v)−wu = P(v) - w. Observe that u∈Su \in \mathcal{S}. By the previous lemma,

∣∣v−w∣∣2=∣∣v−P(v)+P(v)−w∣∣2=∣∣v−P(v)∣∣2+∣∣P(v)−w∣∣2+<v−P(v),u>H+<u,v−P(v)>H=∣∣v−P(v)∣∣2+∣∣P(v)−w∣∣2 \begin{aligned} || v - w ||^2 &= || v - P(v) + P(v) - w ||^2 \\ &= || v - P(v) ||^2 + ||P(v) - w||^2 + \left< v - P(v), u \right>_H + \left< u, v - P(v) \right>_H \\ &= || v - P(v) ||^2 + ||P(v) - w||^2 \end{aligned}

Hence, ∣∣v−w∣∣2|| v - w ||^2 is minimized if and only if w=P(v)w = P(v).

□\square

Matrix Form

The coefficients of a0,a1,...,an−1a_0, a_1, ..., a_{n - 1} can be given in matrix form:

[a0a1...an−1]=[−eˉ0−...−eˉk−...−eˉn−1−]v=1n[e−0⋅2⋅0⋅πin,e−12⋅0⋅πin,...,e−(n−1)⋅2⋅0⋅πin...e−0⋅2⋅(k−1)⋅πin,e−12⋅(k−1)⋅πin,...,e−(n−1)⋅2⋅(k−1)⋅πin...e−0⋅2⋅(n−1)⋅πin,e−12⋅(n−1)⋅πin,...,e−(n−1)⋅2⋅(n−1)⋅πin]v. \begin{bmatrix} a_0 \\ a_1 \\ ... \\ a_{n - 1} \end{bmatrix} = \begin{bmatrix} - \bar e_0 - \\ ...\\ - \bar e_k - \\ ... \\ - \bar e_{n - 1} - \\ \end{bmatrix} v = \frac{1}{ \sqrt n } \begin{bmatrix} e^{ -0 \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot 0 \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} } \\ ...\\ e^{ -0 \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} } \\ ... \\ e^{ -0 \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} } \\ \end{bmatrix} v.

Define the matrix MnM_n as

Mn≐[e−0⋅2⋅0⋅πin,e−12⋅0⋅πin,...,e−(n−1)⋅2⋅0⋅πin...e−0⋅2⋅(k−1)⋅πin,e−12⋅(k−1)⋅πin,...,e−(n−1)⋅2⋅(k−1)⋅πin...e−0⋅2⋅(n−1)⋅πin,e−12⋅(n−1)⋅πin,...,e−(n−1)⋅2⋅(n−1)⋅πin]. M_n \doteq \begin{bmatrix} e^{ -0 \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot 0 \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot 0 \cdot \pi i}{ n} } \\ ...\\ e^{ -0 \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot (k - 1) \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot (k - 1) \cdot \pi i}{ n} } \\ ... \\ e^{ -0 \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, e^{ -1 \frac{2 \cdot (n - 1) \cdot \pi i}{ n} }, ..., e^{ -(n - 1) \cdot \frac{2 \cdot (n - 1) \cdot \pi i}{ n} } \\ \end{bmatrix}.

The key result of this section is the following theorem.

Theorem. An nn-dimension discrete Fourier transform can be performed in O(nlog⁡n)O(n \log n) time. Equivalently, given v∈Cnv \in \mathbb{C}^n, MnvM_n v can be computed in O(nlog⁡n)O(n \log n) time.

The underlying idea of the theorem is divide-and-conquer. Let T(n)T(n) be the time of performing an nn-dimension discrete Fourier transform. We will prove the following recursion

T(n)=2⋅T(n/2)+O(n)T(n) = 2 \cdot T(n / 2) + O(n).

We begin by a special case where n=8n = 8. As each entry in M8M_8 is complex number with unit length, we can represent it by a unit length vector, as the figure below shows.

There are some patterns in M8M_8 and it can be decomposed into 4 matrices related to M4M_4.

Remark: I get this visualization idea from [1].

First, look the the first half of the rows of MM, denoted as MUM_U. The even columns constitute exactly the matrix M4M_4. The odd columns forms a 4×44 \times 4 matrix that can be obtained by proper rotations of M4M_4.

Let vev_e, vov_o be the even and odd dimensions of vv, respectively. Then

MUv=M4ve+[e−2⋅0⋅πine−2⋅1⋅πine−2⋅2⋅πine−2⋅3⋅πin]×M4v0, M_U v = M_4 v_e + \begin{bmatrix} e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} } \\ e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} } \\ e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} } \\ e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \end{bmatrix} \times M_4 v_0,

where ×\times is element-wise multiplication between two vectors. This can also be replaced by matrix multiplication.

MUv=M4ve+diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)M4v0, M_U v = M_4 v_e + \text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) M_4 v_0,

where diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)\text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) denotes a 4×44 \times 4-diagonal matrix.

First, look the the second half of the rows of MM, denoted as MDM_D. The even columns constitute exactly the matrix M4M_4, which is the same as the even columns of MUM_U.

For the odd columns of MDM_D, they are in the opposite directions of the odds columns of MUM_U. This can also be verified algebraically. An entry Mk,ℓM_{k, \ell} in the odd columns MUM_U (k∈[0,n/2)k \in [0, n / 2) and ℓ\ell is odd) equals to

Mk,ℓ=e−ℓ⋅(k−1)⋅2⋅π⋅in. M_{k, \ell} = e^{ - \ell \cdot (k - 1) \cdot \frac{2 \cdot \pi \cdot i}{n} }.

Its corresponding entry in MDM_D is

Mk+n/2,ℓ=e−ℓ⋅(k+n/2−1)⋅2⋅π⋅in=e−ℓ⋅(k−1)⋅2⋅π⋅in⋅e−ℓ⋅π⋅i=Mk,ℓ⋅(−1). M_{k + n / 2, \ell} = e^{ - \ell \cdot (k + n / 2- 1) \cdot \frac{2 \cdot \pi \cdot i}{n} } = e^{ - \ell \cdot (k - 1) \cdot \frac{2 \cdot \pi \cdot i}{n} } \cdot e^{ - \ell \cdot \pi \cdot i } = M_{k, \ell} \cdot (-1).

The last equality holds since ℓ\ell is odd.

Remark: This also explains why Mk+n/2,ℓ=Mk,ℓM_{k + n / 2, \ell} = M_{k, \ell} when ℓ\ell is even.

In summary, we have

Mv=[MUMD]v=[Mn/2ve+diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)Mn/2v0Mn/2ve−diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)Mn/2v0]. M v = \begin{bmatrix} M_U \\ M_D \end{bmatrix} v = \begin{bmatrix} M_{n / 2} v_e + \text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) M_{n / 2} v_0 \\ M_{n / 2} v_e - \text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) M_{n / 2} v_0 \end{bmatrix}.

And the Fast Fourier Transform algorithm is shown as below.

Algorithm FFT (Fast Fourier Transform).
Input. v∈Cnv \in \mathbb{C}^n.

  1. If n=1n = 1, return vv.
  2. Compute re←Mn/2ver_e \leftarrow M_{n / 2} v_e by FFT(vev_e).
  3. Compute ro←Mn/2vor_o \leftarrow M_{n / 2} v_o by FFT(vov_o).
  4. Return
r←[re+diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)rore−diag(e−2⋅0⋅πin,e−2⋅1⋅πin,e−2⋅2⋅πin,e−2⋅3⋅πin)ro]. r \leftarrow \begin{bmatrix} r_e + \text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) r_o \\ r_e - \text{diag} \left( e^{ -\frac{2 \cdot 0 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 1 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 2 \cdot \pi i}{ n} }, e^{ -\frac{2 \cdot 3 \cdot \pi i}{ n} } \right) r_o \end{bmatrix}.

Both step 2 and 3 have time T(n/2)T(n / 2). Step 4 aggregates the results and has time O(n)O(n), which implies that

T(n)=2⋅T(n/2)+O(n). T(n) = 2 \cdot T(n / 2) + O(n).

Reference

[1] T. Roughgarden and G. Valiant, “CS168: The Modern Algorithmic Toolbox Lecture #15 and #16: The Fourier Transform and Convolution,” p. 12.