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 which is defined on the interval . Due to memory limitations, a computer can only store values at a finite number of sample points, i.e. . For our purposes these points will be equally spaced, for example , and so we can write
where are the sample points, is the number of sample points and
It is convenient to use the standard interval, for which . Rewriting in terms of standard interval yields
Notice how is omitted; periodicity implies that the value of the function at is the same as the value of the function at , 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 at the sample points as a finite dimensional vector
The DFT decomposes the sampled function into a linear combination of complex exponentials, where is an index. Since
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 points, the highest frequency of oscillation that can be resolved will have oscillations. Any frequencies higher than in the original function are not adequately resolved and cause an aliasing error (see, for example, Boyd or Uecker 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 can be approximated by an interpolating trigonometric polynomial so that
The symbol means that and agree on each sample point, i.e. for each , but the interpolated polynomial is only an approximation of the true solution away from the sample points.. The are called discrete Fourier coefficients and are what we will be looking to solve for. represents the values of interpolating trigonometric polynomial of degree , 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
where . The interpolation conditions, , can also be rewritten in vectorial form
Here is a vector evaluated at the sample points, which is decomposed into vectors , much as a vector in three dimensional space can be decomposed into the components in the , and directions. The discrete Fourier transform allows us to compute the coefficients 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, , than it is to work with the original function. In order to solve for , we use the orthonormality of the basis elements . We now explain how this is done (For a more detailed explanation see Olver and Shakiban.).
Define . We observe that
For this reason is known as the primitive root of unity. Note also that for , we have that , so all other roots of unity when taken to the power can be obtained from the primitive root of unity. We will use this to perform the DFT algorithm to calculate the coefficients in eq. 1 . The main idea behind the DFT algorithm is to use orthogonality of the vectors . To show the orthogonality between the vectors and , we let denote the complex conjugate of , and then take the inner product of and and find that
if and otherwise.
To deduce the last part, if then , and if , 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
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, using the DFT from the definition can be very slow for large values of . Computing the Fourier coefficients requires complex multiplications and complex additions. In 1960, Cooley and Tukey 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 . For large values of , 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 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.
Cooley, J.W.; Tukey, J.W. (1965). "An algorithm for the machine calculation of complex Fourier series". Mathematics of Computation 19: 297-301.
Shakiban, C. (2006). Applied Linear Algebra. Prentice Hall.