The discrete Fourier transform

Motivation

Let’s assume we have an \(L\)-periodic function \(f(x)\), which is defined on the interval \([0, L)\). We define the \(N\) equidistant points on this interval as \(x_k = k \frac{L}{N}\), where \(k = 0, 1, \ldots, N-1\) with corresponding function values \(f_k = f(x_k)\).

Based on the points \(x_k\) and sampled function values \(f_k\) we now want to compute/approximate the Fourier series for the function \(f(x)\). As we have \(N\) sampling points, it seems natural to compute \(N\) Fourier coefficients \(c_k(f)\) for the Fourier series expansion of \(f(x)\). To compute the integrals for the Fourier coefficients, we recall the definition of the composite trapezoidal rule to approximate integrals. For a given \(L\)-periodic \(g(x)\), the integral of \(g(x)\) over the interval \([0, L)\) can be approximated by

(37)\[\begin{align} \int_0^L g(x) \, dx &\approx \frac{L}{N} \left( \frac{g_0}{2} + g_1 + g_2 + \ldots + g_{N-1} + \frac{g_N}{2} \right) \\ &= \frac{L}{N} \left( {g_0} + g_1 + g_2 + \ldots + g_{N-1} \right) \end{align}\]

where we set \(g_k = g(x_k)\) and used the fact that \(g_0 = g_N\) for \(L\)-periodic functions.

We apply this formula to approximate the Fourier coefficients \(c_k(f)\) of the function \(f(x)\):

(38)\[\begin{align} c_k(f) = \widehat{f}(n) &= \frac{1}{L} \int_0^L f(x) e^{-i 2 \pi k x/L} \, dx \\ &\approx \frac{1}{N} \sum_{l=0}^{N-1} f_l e^{-i 2 \pi k x_l / L} \\ &= \frac{1}{N} \sum_{l=0}^{N-1} f_l e^{-i 2 \pi k l / N} \\ &=\frac{1}{N}\sum_{l=0}^{N-1} f_k \omega_N^{-k l} \end{align}\]

Definition 12 (Discrete Fourier transform)

The discrete Fourier transform (DFT) of a sequence \(\boldsymbol{f} = \{f_0, f_1, \ldots, f_{N-1}\} \in \mathbb{C}^N\) is itself a sequence \(\widehat{\boldsymbol{f}} = \{\widehat{f}_0, \widehat{f}_1, \ldots, \widehat{f}_{N-1}\} \in \mathbb{C}^N\) defined by

(39)\[ \widehat{f}_k = \frac{1}{N} \sum_{l=0}^{N-1} f_l \omega_N^{-l k}, \]

where \(\omega_N = e^{-i 2 \pi / N}\).

In matrix notation, the DFT can be written as

\[ \widehat{\boldsymbol{f}} = \mathcal{F}_N \boldsymbol{f} \]

where \(\mathcal{F}_N\) is the (symmetric!) Fourier matrix with elements \(F_{k,l} = \omega_N^{-k l}\), i.e.

\[\begin{split} \mathcal{F}_N = \frac{1}{N} \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_N^{-1} & \omega_N^{-2} & \cdots & \omega_N^{-(N-1)} \\ 1 & \omega_N^{-2} & \omega_N^{-4} & \cdots & \omega_N^{-2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_N^{-(N-1)} & \omega_N^{-2(N-1)} & \cdots & \omega_N^{-(N-1)(N-1)} \end{pmatrix} \end{split}\]

Discrete inner products and orthogonality systems

The approximation of the Fourier coefficients \(c_k(f)\) by the DFT can be interpreted as a discrete inner product of the function values \( f(x_k) = f_k\) with the complex exponentials \(e^{-i 2 \pi k l / N} = \omega_N^{-kl}\). To facilitate the analysis of the DFT, we introduce the following discrete inner product.

Definition 13 (Discrete inner product)

For two complex sequences \(\boldsymbol{f} = \{f_0, f_1, \ldots, f_{N-1}\} \in \mathbb{C}\) and \(\boldsymbol{g} = \{g_0, g_1, \ldots, g_{N-1}\} \in \mathbb{C}\), the discrete inner product is defined as

(40)\[ \langle \boldsymbol{f}, \boldsymbol{g} \rangle_N = \dfrac{1}{N}\sum_{l=0}^{N-1} f_l \overline{g_l} \]

where \(\overline{g_l}\) denotes the complex conjugate of \(g_l\).

As before, we assume that we have a sequence of \(N\) equidistant points \(x_k = k \frac{L}{N}\) on the interval \(I = [0, L)\).

For the given interval \(I\), let us now again consider the complex exponential functions \(\omega^l(x) := e^{i 2 \pi l x/L}\), \(l \in \mathbb{Z}\). Previously, we have seen that these functions form an orthogonal system with respect to the continuous inner product \(\langle f, g \rangle = \int_0^L f(x) \overline{g(x)} \, dx\). Funnily enough, these functions satisfy a very similar orthogonality property with respect to the discrete inner product:

Theorem 12 (Orthogonality of complex exponentials)

\[\begin{split} \langle \omega^l, \omega^m \rangle_N = \begin{cases} 1 & \text{if } (l-m)/N \in \mathbb{Z}, \\ 0 & \text{else.} \end{cases} \end{split}\]

Before we turn to the proof of this theorem, we make the very important observation that the evaluation of \(\omega^l(x)\) at the points \(x_k\) is given by the \(k\)-th power of the \(l\)-th \(N\)-th root of unity \(\omega_N^l = e^{i 2 \pi l / N}\), i.e. \(\omega^l(x_k) = e^{i 2 \pi l k / N} = \omega_N^{lk}\), or in other words

\[ (\omega^l(x_k))_{k=0}^{N-1} = (\omega_N^{lk})_{k=0}^{N-1} = (1, \omega_N^l, \omega_N^{2l}, \ldots, \omega_N^{(N-1)l}). \]

(four:eq:nth_root_vectors)

Proof. First we recall a fundamental identity for geometric sums, namely that for any given \(q\neq 1\), we have

\[ \sum_{k=0}^{N-1} q^{N-1} = \frac{1-q^N}{1-q}. \]

For two complex exponentials \(\omega^l\) and \(\omega^m\), we compute the discrete inner product

(41)\[\begin{align} N \langle \omega^l, \omega^m \rangle_N &= \sum_{k=0}^{N-1} \omega^{l}(x_k) \overline{\omega^{m}(x_k)} \\ &= \sum_{k=0}^{N-1} \omega_N^{lk} \omega_N^{-mk} \\ &= \sum_{k=0}^{N-1} (\omega_N^{l-m})^k \end{align}\]

If \(l-m = p N\) for some \(p \in \mathbb{Z}\), then \(\omega_N^{l-m} = e^{i 2\pi (l-m)/N} = e^{i 2 \pi p} = 1\) and the sum evaluates to \(N\). Otherwise we use the geometric sum identity to obtain

\[ \sum_{k=0}^{N-1} (\omega_N^{l-m})^k = \frac{1 - \omega_N^{(l-m)N}}{1 - \omega_N^{l-m}} = \frac{1 - e^{i 2\pi(l-m)}}{1 - \omega_N^{l-m}} = 0. \]

Let us record a number of important consequences of this orthogonality property.

Corollary 1

1. For \(0\leqslant l, m \leqslant N-1\), the complex exponentials \(\omega^l\) and \(\omega^m\) are orthonormal with respect to the discrete inner product

   \[ \langle \omega^l, \omega^m \rangle_N = \delta_{l m}. \]
   Equivalently, the \(N\) vectors
   \[ (\omega^l(x_k))_{k=0}^{N-1} = (\omega_N^{lk})_{k=0}^{N-1} = (1, \omega_N^l, \omega_N^{2l}, \ldots, \omega_N^{(N-1)l}) \quad {l=0, 1, \ldots, N-1} \]
   are orthonormal with respect to the discrete inner product. In particular, they forma a orthornormal basis of the vector space \(\mathbb{C}^N\).

2. The numerical quadrature rule

   \[ \int_0^L f(x) \, dx \approx \frac{L}{N} \sum_{k=0}^{N-1} f(x_k) \]
   is in fact exact for the integration of the complex exponential functions \(\omega^l(x)\) for \(0\leqslant l \leqslant N-1\). Moreover, the quadrature rule is exact for the integration for \(\sin(\pm 2\pi N/L x)\) but not for \(\cos(\pm 2\pi N/L x)\).

Exercise 27

Prove the last statement using the the discrete orthogonality of the complex exponentials.

Important

The definition of the discrete Fourier transform in Equation (39) can be neatly rewritten in terms of the discrete inner product as

(42)\[\begin{align} \widehat{f}_k &= \langle f(x), \omega^k(x) \rangle_N \\ &= \langle (f_l)_{l=0}^{N-1}, (\omega^k(x_l))_{l=0}^{N-1} \rangle_N \\ &= \dfrac{1}{N} (f_0, f_1, \ldots, f_{N-1}) \cdot (1, \omega_N^l, \omega_N^{2l}, \ldots, \omega_N^{(N-1)l}). \end{align}\]

where \(\cdot\) denotes the usual (complex) scalar product of two vectors.

This strongly resembles the definition of the usual Fourier coefficients \(c_k(f) = \widehat{f}(k)\):

(43)\[\begin{align} \widehat{f}(k) &= \dfrac{1}{L} \int_0^L f(x) \overline{\omega^k(x)} \, dx \\ &= \dfrac{1}{L} \langle f, \omega^k \rangle \end{align}\]

Or in other words: The discrete Fourier coefficient \(\widehat{f}_k\) resulting from the DFT is a discrete inner product of the function values \(f_k\) with the complex exponentials \(\omega^l(x_k)\), while the usual Fourier coefficient \(\widehat{f}(k) = c_k(f)\) is the continuous inner product of the function \(f(x)\) with the complex exponentials \(\omega^l(x)\).

We can use the previous observations to show the the matrix associated with the discrete Fourier transform in invertible and to compute its inverse.

Theorem 13

The matrix of the DFT is invertible and its inverse \(\mathcal{F}^{-1}_N\) is given by

\[\begin{split} \mathcal{F}^{-1}_N = \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_N^{1} & \omega_N^{2} & \cdots & \omega_N^{(N-1)} \\ 1 & \omega_N^{2} & \omega_N^{4} & \cdots & \omega_N^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_N^{(N-1)} & \omega_N^{2(N-1)} & \cdots & \omega_N^{(N-1)(N-1)} \end{pmatrix} \end{split}\]

Proof. We simply compute the product of the DFT matrix \(\mathcal{F}_N\) and its inverse \(\mathcal{F}^{-1}_N\):

(44)\[\begin{align} (\mathcal{F}^{-1}_N \mathcal{F}^{}_N)_{l, m} = \sum_{k=0}^{N-1}F^{-1}_{l, k} F_{k, m} = \dfrac{1}{N}\sum_{k=0}^{N-1} \omega_N^{lk} \omega_N^{-km} = \langle \omega^l, \omega^m \rangle_N = \delta_{l m} \end{align}\]

This gives rise to the following definition of the inverse DFT.

Definition 14 (Inverse discrete Fourier transform)

The inverse discrete Fourier transform (IDFT) of a sequence \( \boldsymbol{c} = \{c_0, c_1, \ldots, c_{N-1}\} \in \mathbb{C}^N\) is itself a sequence \({\boldsymbol{f}} = \{{f}_0, {f}_1, \ldots, {f}_{N-1}\} \in \mathbb{C}^N\) defined by

(45)\[ {f}_l = \sum_{k=0}^{N-1} c_k \omega_N^{lk} \]

Important note!

There a lot of different conventions for the normalization of the DFT and the IDFT. Another common one is to use the normalization factor \(1/\sqrt{N}\) for both the DFT and the IDFT, since the columns (respectively rows) of the resulting matrices are then orthonormal with respect to the inner product (fou:eq:disc_inner_product). and hence they unitary, i.e. \(\mathcal{F}_N^* \mathcal{F}_N = \overline{\mathcal{F}_N}^{\top} \mathcal{F}_N = \mathcal{I}\).

Another possible convention is to use the normalization factor \(1/N\) for the IDFT and \(1\) for the DFT. And sometimes, even the sign in complex exponential is changed!

As a consequence, you should always check the normalization conventions both in the used in the in the literature or in software packages! In scipy.fft module which we will use later, the function which computes the DFT has an optional flag to switch between the different conventions.