Parallel Spectral Numerical Methods/One-Dimensional Discrete Fourier Transforms

The discrete Fourier transform (DFT) takes a function sampled at a finite number of points and finds the coefficients for the linear combination of trigonometric polynomials that best approximates the function; the number of trigonometric polynomials used is equal to the number of sample points. Suppose we have a function $f (x)$ which is defined on the interval $a \leq x \leq b$ . Due to memory limitations, a computer can only store values at a finite number of sample points, i.e. $a \leq x_{0} < x_{1} < . . . < x_{n} \leq b$ . For our purposes these points will be equally spaced, for example $x_{1} - x_{0} = x_{3} - x_{2}$ , and so we can write

$x_{j} = a + j h, j = 0, 1, 2, . . ., n$

where $x_{j}$ are the sample points, $n$ is the number of sample points and

$h = \frac{b - a}{n} .$

It is convenient to use the standard interval, for which $0 \leq x \leq 2 π$ . Rewriting $x$ in terms of standard interval yields

$x_{0} = 0, x_{1} = \frac{2 π}{n}, x_{2} = \frac{4 π}{n}, x_{j} = \frac{2 j π}{n}, . . ., x_{n - 1} = \frac{2 (n - 1) π}{n}$

Notice how $x_{n} = 2 π$ is omitted; periodicity implies that the value of the function at $2 π$ is the same as the value of the function at $0$ , so it need not be included. We will introduce the DFT using the language of linear algebra. Much of this formalism carries over to continuous functions that are being approximated. It also makes it easier to understand the computer implementation of the algorithms. Many computer packages and programs are optimized to perform calculations through matrix operations, so the formalism is also useful when actually calculating transforms. We write the approximation to $f (x)$ at the sample points as a finite dimensional vector

$f = (f_{0}, f_{1}, . . ., f_{n - 1})^{T} = (f (x_{0}), f (x_{1}), . . ., f (x_{n - 1}))$

where

$f_{j} = f (x_{j}) = f (\frac{2 j π}{n}) .$

The DFT decomposes the sampled function $f (x)$ into a linear combination of complex exponentials, $\exp (i k x)$ where $k$ is an index. Since

$\exp (i k x) = \cos (k x) + i \sin (k x),$

we also obtain an expansion in trigonometric functions, which may be more familiar from courses in calculus and differential equations. Since the function is sampled at $n$ points, the highest frequency of oscillation that can be resolved will have $n$ oscillations. Any frequencies higher than $n$ in the original function are not adequately resolved and cause an aliasing error (see, for example, Boyd^[1] or Uecker^[2] for more on this). This error can be reduced by sampling at a greater number of points so that the number of approximating exponentials functions can also be increased. There is a tradeoff between increasing the accuracy of the simulation and the time required for the simulation to complete. For many cases of scientific and practical interest, simulations with up to thousands of grid points can be computed relatively quickly. Below we explain how a function $f (x)$ can be approximated by an interpolating trigonometric polynomial $p (x)$ so that

$f (x) \approx p (x) = c_{0} + c_{1} e^{2 i x} + c_{2} e^{2 i x} + . . . + c_{n - 1} e^{(n - 1) i x} = \sum_{k = 0}^{n - 1} c_{k} e^{i k x}$

The $\approx$ symbol means that $f (x)$ and $p (x)$ agree on each sample point, i.e. $f (x_{j}) = p (x_{j})$ for each $j = 0, 1, . . . n - 1$ , but the interpolated polynomial $p (x)$ is only an approximation of the true solution $f (x)$ away from the sample points.. The $c_{n}$ are called discrete Fourier coefficients and are what we will be looking to solve for. $p (x)$ represents the values of interpolating trigonometric polynomial of degree $\leq n - 1$ , so if we have the values of these coefficients then we have a function we can use. Since we are working in a finite-dimensional vector space, a useful approach is to rewrite the discrete Fourier series as a vector. We let

$ω_{k} = (e^{i k x_{0}}, e^{i k x_{1}}, e^{i k x_{2}}, . . ., e^{i k x_{n}})^{T}$

$= (1, e^{2 k π i / n}, e^{4 k π i / n}, . . ., e^{2 (n - 1) k π i / n})^{T},$

where $k = 0, 1, . . ., n - 1$ . The interpolation conditions, $f (x_{j}) = p (x_{j})$ , can also be rewritten in vectorial form

Template:NumBlk

Here $f$ is a vector evaluated at the sample points, which is decomposed into vectors $ω_{k}$ , much as a vector in three dimensional space can be decomposed into the components in the $x$ , $y$ and $z$ directions. The discrete Fourier transform allows us to compute the coefficients $c_{i}$ given the value of the function at the sample points. This may at first seem unmotivated, but in many applications, such as solving differential equations, it is easier to manipulate a linear combination of trigonometric polynomials, $ω_{0}, . . ., ω_{n - 1}$ , than it is to work with the original function. In order to solve for $c_{k}$ , we use the orthonormality of the basis elements $ω_{0}, . . ., ω_{n - 1}$ . We now explain how this is done (For a more detailed explanation see Olver and Shakiban^[3].).

Define $ξ_{n} = e^{2 π i / n}$ . We observe that

${(ξ_{n})}^{n} = \exp (\frac{2 π i n}{n}) = \cos (2 π) + i \sin (2 π) = 1$

For this reason $ξ_{n}$ is known as the primitive $n^{th}$ root of unity. Note also that for $0 \leq k < n$ , we have that $(ξ_{n}^{k})^{n} = 1$ , so all other roots of unity when taken to the power $n$ can be obtained from the primitive $n^{th}$ root of unity. We will use this to perform the DFT algorithm to calculate the coefficients $c_{0}, . . ., c_{k - 1}$ in eq. Template:EquationNote. The main idea behind the DFT algorithm is to use orthogonality of the vectors $ω_{k}$ . To show the orthogonality between the vectors $ω_{k}$ and $ω_{l}$ , we let $ω_{l}^{*}$ denote the complex conjugate of $ω_{l}$ , and then take the inner product of $ω_{k}$ and $ω_{l}$ and find that

$< ω_{k}, ω_{l} >$

$= \frac{1}{n} \sum_{m = 0}^{n - 1} \exp (\frac{2 π i k m}{n}) {[\exp (\frac{2 π i l m}{n})]}^{*}$

$= \frac{1}{n} \sum_{m = 0}^{n - 1} \exp (\frac{2 π i (k - l) m}{n})$

$= \frac{1}{n} \sum_{m = 0}^{n - 1} \cos (\frac{π (k - l) m}{n}) + i \sin (\frac{π (k - l) m}{n})$

$= 1$ if $k = l$ and $0$ otherwise.

To deduce the last part, if $k = l$ then $\exp (0) = 1$ , and if $k \neq l$ , then we are sampling the sine and cosine functions at equally spaced points on over an integral number of wavelengths. Since these functions have equal magnitude positive and negative parts, they sum to zero, much as the integral of a sine or cosine over an integral number of wavelengths is zero. This implies that we can compute the Fourier coefficients in the discrete Fourier sum by taking inner products

$c_{k} = < f, ω_{k} > = \frac{1}{n} \sum_{m = 0}^{n - 1} ξ_{n}^{- m k} f_{j} .$

We note the close connection between the continuous and discrete settings, where an integral is replaced by a sum.

Fast Fourier Transform

Computing the Fourier coefficients, $c_{0}, . . ., c_{n - 1}$ using the DFT from the definition can be very slow for large values of $n$ . Computing the Fourier coefficients $c_{0}, . . . c_{n - 1}$ requires $n^{2} - n$ complex multiplications and $n^{2} - n$ complex additions. In 1960, Cooley and Tukey^[4] rediscovered a much more efficient way of computing DFT by developing an algorithm known as the Fast Fourier Transforms (FFT). The FFT cuts the number of arithmetic operations down to $O (n \log n)$ . For large values of $n$ , this can make a huge difference in computation time compared to the standard DFT. The reason why the FFT is so important is that it is heavily used spectral methods. The basic FFT algorithm used by Cooley and Tukey^[5] is well documented in many places, however, there are other implementations of the algorithm and the best version of the algorithm to use depends heavily on computer architecture. We therefore do not give further descriptions here.