8.2. Discrete Fourier Transform

If we want to know \(\tilde{f}(\omega)\) just for a particular value of \(\omega\), the transformation (\ref{eq:fourier_fwd}) is simply an improper integral and thus the numerical methods we studied in Section 5.3.1 is sufficient. In general that is not what we want. We want to know the function \(\tilde{f}(\omega)\) for \(-\infty < \omega <\infty\). That is a big challenge for digital computers. First of all, the variables \(t\) and \(\omega\) are continuous, which digital computer cannot handle and thus they must be digitized as we did in all previous chapters, If \(t\) is discretized with \(N\) points and \(\omega\) with \(M\) points, the transformation needs the order of \(N M\) operations. In a three-dimensional space, \(N\) can be easily \(1 \times 10^{6}\). \(M\) is also in the similar order. Hence, the number of operations can be huge. Fortunately, there is an ingenious algorithm called fast fourier tranform or FFT. But before going to FFT, we need to express Fourier transform in a discrete form known as discrete Fourier transform. The algorithm of FFT will be explained in next section.

First, we need to replace the infinite integration interval \((-\infty, +\infty)\) in Eq. () with \([-T/2,+T/2]\) where \(T>0\) is assumed to be sufficiently large so that \(f(t)\approx 0\) for \(|t|\ge T/2\). Second, we introduce discrete time \(t_n = n \Delta t\) (\(n=-N/2, \cdots, N/2\)) where \(\Delta t=T/N\). Then, the transform () is replaced with a summation

()\[ \tilde{f}(\omega) = \sum_{n=-N/2}^{N/2-1} f(t_n)\, e^{i n \omega \Delta t} \, \Delta t\,. \]

Here we use the trapezoidal rule for numerical integration (see Chapter Section 5).\footnote{It looks like the rectangular rule but recall that when \(f(-T/2)=f(T/2)=0\), the rectangular rule is identical to trapezoidal rule.} Since we have chosen \(T\) so that \(f(-T/2)=f(T/2) \approx 0\), the rectangular rule is OK provided that \(\Delta t\) is small enough. Note that there are \(N+1\) sampling points of \(t\) but only \(N\) rectangles to sum up. That is why the upper limit of the summation is \(N/2-1\). We need to do the same for \(\omega\). Introducing the integration interval \([-\Omega/2, +\Omega/2]\) and discrete frequency \(\omega_m = m \Delta \omega\) (\(m=-M/2, \cdots, M/2\)) where \(\Delta \omega = \Omega/M\), the inverse transform (\ref{eq:fourier_inv}) is expressed as a numerical integral

()\[ f(t) = \frac{1}{2\pi}\sum_{m=-M/2}^{M/2-1} \tilde{f}(\omega_m) \, e^{-im \Delta \omega t}\, \Delta \omega\,. \]

Now, the discrete forms of transformations must consistent with the Fourier theorem (). Substituting Eq. () to Eq. (), we find the discrete version of the Fourier theorem

()\[ f(t_k) = \frac{1}{2\pi} \sum_{m=-M/2}^{M/2}\quad \sum_{n=-N/2}^{N/2} f(t_n) e^{i (n-k) m \Delta \omega \Delta t} \Delta \omega \Delta t. \]

\[ \tilde{f}(\omega_k) = \frac{1}{2\pi} \sum_{n=-N/2}^{N/2}\quad \sum_{m=-M/2}^{M/2} \tilde{f}(\omega_m) e^{-i (m-k) n \Delta \omega \Delta t} \Delta \omega \Delta t. \]

These equations hold simultaneously if \(N=M\) and \(\Delta \omega \Delta t = 2\pi/N\). Commonly, \(\Delta t = T/N\) and \(\Delta \omega = 2\pi /T\) are used. The bound of \(\omega\) is now \(\Omega=2\pi N/T\).

In practice, the choice of \(T\) is sometime tricky. Normally, we choose \(T\) such that \(f(T/2) \approx 0\) and \(N\) such that the resolution \(\Delta t=T/N\) is small enough. However, we also need a reasonable resolution of frequency \(\Delta\omega=2\pi/T\). If \(T\) is too small, the resolution of the frequency becomes poor. Therefore, a larger \(T\) is better for the frequency. On the other hand, if \(T\) is large, \(N\) has to be large so that the resolution of time is small enough. We will see an example below.

Considering the periodicity of the exponential function, the bound of the summation \((-N/2,N/2-1)\) may be shifted to \((0,N-1)\). In summary, the discrete version of Fourier transforms (DFT) are defined by

()\[\begin{split} \begin{eqnarray} \tilde{f}_m &=& T \left ( \frac{1}{N}\sum_{n=0}^{N-1} f_n\, e^{2 \pi i m n / N} \right )\\ f_n &=& \frac{1}{T} \left ( \sum_{m=0}^{N-1} \tilde{f}_m\, e^{- 2 \pi i n m /N} \right ) \end{eqnarray} \end{split}\]

where short-handed notations \(f_n \equiv f(t_n)\) and \(\tilde{f}_m\equiv \tilde{f}(\omega_m)\) are used.

An important consequence of the discretization is that the discretized functions are periodic even when the original functions are not. Explicitly writing it,

\[\begin{split} \begin{eqnarray} f_{n+N} = f_n \qquad &\text{or}& \qquad f(t+T)=f(t) \\ \tilde{f}_{n+N}=\tilde{f}_n \qquad &\text{or}& \qquad \tilde{f}(\omega+\Omega) = \tilde{f}(\omega). \end{eqnarray} \end{split}\]

DFT () can be expressed in a matrix form,

()\[ \tilde{\mathbf{f}} = \mathcal{F}\, \mathbf{f}, \qquad \mathbf{f} = \mathcal{F}^{-1}\, \tilde{\mathbf{f}} \]

where the matrix elements are defined by

()\[ \mathcal{F}_{mn} = \frac{T}{N}\, e^{2\pi i m n/N}, \quad n,m = 0, \cdots, N-1 \]

and the functions are expressed as vectors

\[\begin{split} \mathbf{f} = \begin{bmatrix} f_1 \\ f_2 \\ \vdots \\ f_N \end{bmatrix}, \qquad \tilde{\mathbf{f}} = \begin{bmatrix} \tilde{f}_1 \\ \tilde{f}_2 \\ \vdots \\ \tilde{f}_N \end{bmatrix} \end{split}\]

The multiplication of a matrix and a vector involves \(N^2\) of multiplications and \(N^2\) of additions, which can be too large if we need to perform Fourier transform many times.

---

Written on 9/23/2024 by Ryoichi Kawai.