Discrete Fourier Transform

Section 25.3 Discrete Fourier Transform

There is a discrete version of the Fourier series theory outlined in Section 25.2. In it we fix an integer \(n\) and replace a function \(f\colon [0, 1]\to \mathbb C\) by its sampled values \(z_k = f(k/n)\text{,}\) \(k=0, 1, \dots, n-1\text{.}\) The discrete version of inner product (25.2.1) is

\begin{equation} \langle f, g \rangle_d = \sum_{\ell=0}^{n-1}f(\ell/n)\overline{g(\ell/n)} = \sum_{\ell=0}^{n-1} z_\ell \overline{w_\ell}\label{eq-inner-product-discrete}\tag{25.3.1} \end{equation}

The functions \(f_k = \exp(2\pi i k)\text{,}\) \(k=0, \dots, n-1\text{,}\) are orthogonal with respect to this inner product. (Not orthonormal, as \(\langle f_k, f_k\rangle = n\text{.}\)) Note we only need \(n\) of these functions (why?).

The Fourier coefficients of a given discretized function \((z_0, z_1, \dots, z_{n-1})\) are computed as:

\begin{equation} c_k = \sum_{\ell=0}^{n-1} z_\ell \exp(-2\pi i k\ell/n) \label{eq-fourier-coefficients-discrete}\tag{25.3.2} \end{equation}

and from there, the function is precisely reconstructed as

\begin{equation} z_\ell = \frac{1}{n} \sum_{k=0}^{n-1} c_k \exp(2\pi i k\ell/n)\label{eq-fourier-sum-discrete}\tag{25.3.3} \end{equation}

Some people put \(1/n\) in (25.3.2) instead of (25.3.3), and some put \(1/\sqrt{n}\) in both, for symmetry.

The computation of Fourier coefficients (25.3.2) amounts to matrix-vector multiplication \(c = Wz\) where \(W\) is the Fourier matrix defined as follows. Let \(\omega = \exp(-2\pi i /n)\text{.}\) The matrix \(W\) has size \(n\times n\text{,}\) with entries given by \(W_{k\ell} = \omega^{(k-1)(\ell-1)}\text{.}\) Explicitly,

\begin{equation} W = \begin{pmatrix} 1 & 1 & 1 & 1 & \dots & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \dots & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \dots & \omega^{2(n-1)} \\ 1 & \omega^3 & \omega^6 & \omega^9 & \dots & \omega^{3(n-1)} \\ \cdots \\ 1 & \omega^{n-1} & \omega^{2(n-1)} & \omega^{3(n-1)} & \dots & \omega^{(n-1)(n-1)} \\ \end{pmatrix}\label{eq-fourier-matrix}\tag{25.3.4} \end{equation}

Normally, matrix-vector multiplication takes on the order of \(n^2\) operations. The Fast Fourier Transform (FFT) is an algorithm that uses the specific structure of \(W\) to reduce the number of computations to something of the order \(n\log n\text{.}\) The reconstruction of function values from coefficients (25.3.3) is very similar. Besides the factor of \(1/n\text{,}\) the only other difference is replacing matrix \(W\) with matrix \(W^*\) of the same form (25.3.4) but with \(\omega = \exp(2\pi i /n)\text{.}\)

In Matlab, the computation (25.3.2) is implemented by command fft while (25.3.3) is implemented by ifft (“i” for “inverse”). Try applying them to standard basis vectors and compare the output to the columns of \(W\) or \(\frac{1}{n} W^*\text{.}\) (The commands also allow matrix input: fft(eye(n)) and ifft(eye(n)) return the matrices \(W\) and \(\frac{1}{n} W^*\text{.}\))