Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT License. (c)2015 David I. Ketcheson

Version 0.1 - March 2016

The FFT, aliasing, and filtering

Welcome to lesson 3. Here we'll learn about the Fast Fourier Transform, which we've been using all along. We'll also learn about a numerical pathology of pseudospectral methods (known as aliasing) and one way to avoid it (known as filtering or dealiasing).

The fast Fourier transform

We won't go into great detail regarding the FFT algorithm, since there is already an excellent explanation of the Fast Fourier Transform in Jupyter Notebook form available on the web:

• Understanding the FFT Algorithm by Jake Vanderplas

Suffice it to say that the FFT is a fast algorithm for computing the discrete Fourier transform (DFT):

$$ \hat{u}_\xi = \sum_{j=0}^{m-1} u_j e^{-2\pi i \xi j/m} $$

or its inverse. The DFT, as we know, is linear and can be computed by multiplying $u$ by a certain $m\times m$ dense matrix $F$. Multiplication by a dense matrix requeires ${\mathcal O}(m^2)$ operations. The FFT is a shortcut to compute that matrix-vector product in just ${\mathcal O}(m \log m)$ operations by taking advantage of the special structure of $F$.

This is very important for pseudospectral methods, since most of the computational work occurs in computing the Fourier transform and its inverse. It's also important that we make use of a compiled version of the FFT, since native Python code is relatively slow. The np.fft function is an interface to a compiled FFT library.

Ordering of wavenumbers in FFT output

The vector returned by np.fft contains the Fourier coefficients of $u$ corresponding to the wavenumbers

$$ \frac{2\pi}{L} \{-m/2, -m/2 + 1, \dots, m/2 - 1\}. $$

However, for computational efficiency the output vector does not use the natural ordering above. The ordering it uses can be obtained with the following command.

As you can see, the return vector starts with the nonnegative wavenumbers, followed by the negative wavenumbers. It may seem strange to you that the range of wavenumbers returned is not symmetric; in the case above, it includes $-8$ but not $+8$. This apparent asymmetry can be explained once one understands the phenomenon known as aliasing.

Aliasing

A numerical grid has a limited resolution. If you try to represent a rapidly-oscillating function with relatively few grid points, you will observe an effect known as aliasing. This naturally limits the range of frequencies that can be modelled on a given grid. It can also lead to instabilities in pseudospectral simulations, when generation of high frequencies leads to buildup of lower-frequency energy due to aliasing.

The code below plots a sine wave of a given frequency, along with its representation on a grid with $m$ points. Try changing $p$ and notice how for $m<2p$ the function looks like a lower-frequency mode.

Exercise

Try to answer the questions below with pencil and paper; then check them by modifying the code above.

1. For a given number of grid points $m$, which wavenumbers $p$ will be aliased to the $p=0$ mode? Which will be aliased to $p=1$? Can you explain why?
2. What is the highest frequency mode that can be represented on a given grid?

After completing the exercise, explain why the sequence of wavenumbers given by np.fft.fftfreq is not, in fact, asymmetric after all.

Aliasing as a source of numerical instability

As we have seen, aliasing means that wavenumbers of magnitude greater than $\pi m/L$ are incorrectly represented as lower wavenumbers on a grid with $m$ points. This suggests that we shouldn't allow larger wavenumbers in our numerical solution. For linear problems, this simply means that we should represent the initial condition by a truncated Fourier series containing modes with wavenumbers less than $\pi m/L$. This happens naturally when we sample the function at the grid points. As we evolve in time, higher frequencies are not generated due to the linearity of the problem.

Nonlinear problems are a different story. Let's consider what happens when we have a quadratic term like $u^2$, as in Burgers' equation. In general, if the grid function $u$ contains wavenumbers up to $\pi m/L$, then $u^2$ contains frequencies up to $2 \pi m/L$. So each time we compute this term, we generate high frequencies that get aliased back to lower frequencies on our grid. Clearly this has nothing to do with the correct mathematical solution and will lead to errors. Even worse, this aliasing effect can, as it is repeated at every step, lead to an instability that causes the numerical solution to blow up.

$$ \newcommand{\F}{\mathcal F} \newcommand{\Finv}{{\mathcal F}^{-1}} $$

An illustration of aliasing instability: the Korteweg-de Vries equation

To see aliasing in practice, we'll consider the KdV equation, which describes certain kinds of water waves:

$$ u_t = -u u_x - u_{xxx} $$

A natural pseudospectral discretization is obtained if we compute the spatial derivatives via

\begin{align} u_x & = \Finv(i\xi \F(u)) \\ u_{xxx} & = \Finv(-i\xi^3 \F(u)). \end{align}

This gives $$ U'(t) = -D[U] \Finv(i\xi \F(U)) - \Finv(-i\xi^3 \F(U)). $$

This is identical to our discretization of Burgers' equation, except that now we have a third-derivative term. explain dispersion here

The largest-magnitude eigenvalues of the Jacobian for this semi-discretization are related to the 3rd-derivative term. If we consider only that term, the eigenvalues are

$$-i \xi^3$$

where $\xi$ lies in the range $(-m/2, m/2)$. So we need the time step to satisfy $\Delta t (m/2)^3 \in S$, where $S$ is the region of absolute stability of a given time integrator.

For this example we'll use a 3rd-order Runge-Kutta method:

Let's check the size of the imaginary axis interval contained in this method's absolute stability region:

Now we'll go ahead and implement our solution, making sure to set the time step according to the condition above.

In the output, we're plotting the solution (top plot) and its Fourier transform (bottom plot). There are a lot of interesting things to say about the solution, but for now let's focus on the Fourier transform. Notice how the wavenumbers present in the solution remain in the lower half of those representable on the grid (this region is delimited by the dashed red lines). Because of this, no aliasing occurs.

Now change the code above to use only $m=128$ grid points. What happens?

Explanation

Here we will give a somewhat simplified explanation of the blow-up just observed. First, this blowup has nothing to do with the absoute stability condition -- when we change $m$, the time step is automatically changed in a way that will ensure absolute stability. If you're not convinced, try taking the time step even smaller; you will still observe the blowup.

By taking $m=128$, we cut by half the wavenumbers that can be represented on the grid. As you can see from the plots, this means that some of the wavenumbers present in the initial data are in the upper half of the representable range (i.e., outside the dashed red lines). That means that the highest wavenumbers generated by the quadratic term will be aliased -- and they will be aliased back into that upper-half range. This leads to a gradual accumulation of high-wavenumber modes, easily visible in both plots. Eventually the high-wavenumber modes dominate the numerical solution and lead to blowup.

For a detailed discussion of aliasing instabilities, see Chapter 11 of John Boyd's "Chebyshev and Fourier Spectral Methods".

Filtering

How can we avoid aliasing instability? The proper approach is to ensure that the solution is well resolved, so that the instability never appears. However, this may entail a very substantial computational cost. One way to ensure stability even if the solution is underresolved is by filtering, which is also known as dealiasing. In general it is unwise to rely on filtering, since it can mask the fact that the solution is not resolved (and hence not accurate). But understanding filtering can give a bit more insight into aliasing instability itself.

At the most basic level, filtering means removing the modes that lead to aliasing. This can be done by damping the high wavenumbers or simply zeroing them when computing the $(u^2)_x$ term. The obvious approach would be to filter the upper half of all wavenumbers, but this is overkill. In fact, it is sufficient to filter only the uppermost third. To see why, notice that the aliased modes resulting from the lower two-thirds will appear in the uppermost third of the range of modes, and so will be filtered at the next step.

A simple 2/3 filter is implemented in the code below.

Notice how the solution remains stable, but small wiggles appear throughout the domain. These are a hint that something is not sufficiently resolved.