Discrete Fourier Transforms

Consider first two one-dimensional arrays \(\boldsymbol{u} = \{u_j\}_{j=0}^{N-1}\) and \(\boldsymbol{\hat{u}} =\{\hat{u}_k\}_{k=0}^{N-1}\). We define the forward and backward Discrete Fourier transforms (DFT), respectively, as

(2)\[\begin{split}\hat{u}_k &= \frac{1}{N}\sum_{j=0}^{N-1}u_j e^{-2\pi i j k / N}, \quad \forall \, k\in \textbf{k}=0, 1, \ldots, N-1, \\ u_j &= \sum_{k=0}^{N-1}\hat{u}_k e^{2\pi i j k / N}, \quad \forall \, j\in\textbf{j}=0, 1, \ldots, N-1,\end{split}\]

where \(i=\sqrt{-1}\). Discrete Fourier transforms are computed efficiently using algorithms termed Fast Fourier Transforms, known in short as FFTs.

Note

The index set for wavenumbers \(\textbf{k}\) is usually not chosen as \([0, 1, \ldots, N-1]\), but \(\textbf{k}=[-N/2, -N/2-1, \ldots, N/2-1]\) for even \(N\) and \(\textbf{k}=[-(N-1)/2, -(N-1)/2+1, \ldots, (N-1)/2]\) for odd \(N\). See numpy.fft.fftfreq. Also note that it is possible to tweak the default normalization used above when calling either forward or backward transforms.

A more compact notation is commonly used for the DFTs, where the 1D forward and backward transforms are written as

\[\begin{split}\boldsymbol{\hat{u}} &= \mathcal{F}(\boldsymbol{u}), \\ \boldsymbol{u} &= \mathcal{F}^{-1}(\boldsymbol{\hat{u}}).\end{split}\]

Numpy, Scipy, and many other scientific softwares contain implementations that make working with Fourier series simple and straight forward. These 1D Fourier transforms can be implemented easily with just Numpy as, e.g.:

import numpy as np
N = 16
u = np.random.random(N)
u_hat = np.fft.fft(u)
uc = np.fft.ifft(u_hat)
assert np.allclose(u, uc)

However, there is a minor difference. Numpy performs by default the \(1/N\) scaling with the backward transform (ifft) and not the forward as shown in (2). These are merely different conventions and not important as long as one is aware of them. We use the scaling on the forward transform simply because this follows naturally when using the harmonic functions \(e^{i k x}\) as basis functions when solving PDEs with the spectral Galerkin method or the spectral collocation method (see chap. 3).

With mpi4py-fft the same operations take just a few more steps, because instead of executing ffts directly, like in the calls for np.fft.fft and np.fft.ifft, we need to create the objects that are to do the transforms first. We need to plan the transforms:

from mpi4py_fft import fftw
u = fftw.aligned(N, dtype=np.complex)
u_hat = fftw.aligned_like(u)
fft = fftw.fftn(u, flags=(fftw.FFTW_MEASURE,))        # plan fft
ifft = fftw.ifftn(u_hat, flags=(fftw.FFTW_ESTIMATE,)) # plan ifft
u[:] = np.random.random(N)
# Now execute the transforms
u_hat = fft(u, u_hat, normalize=True)
uc = ifft(u_hat)
assert np.allclose(uc, u)

The planning of transforms makes an effort to find the fastest possible transform of the given kind. See more in The fftw module.

Multidimensional transforms

It is for multidimensional arrays that it starts to become interesting for the current software. Multidimensional arrays are a bit tedious with notation, though, especially when the number of dimensions grow. We will stick with the index notation because it is most straightforward in comparison with implementation.

We denote the entries of a two-dimensional array as \(u_{j_0, j_1}\), which corresponds to a row-major matrix \(\boldsymbol{u}=\{u_{j_0, j_1}\}_{(j_0, j_1) \in \textbf{j}_0 \times \textbf{j}_1}\) of size \(N_0\cdot N_1\). Denoting also \(\omega_m=j_m k_m / N_m\), a two-dimensional forward and backward DFT can be defined as

(3)\[\begin{split}\hat{u}_{k_0,k_1} &= \frac{1}{N_0}\sum_{j_0 \in \textbf{j}_0}\Big( e^{-2\pi i \omega_0} \frac{1}{N_1} \sum_{j_1\in \textbf{j}_1} \Big( e^{-2\pi i \omega_1} u_{j_0,j_1}\Big) \Big), \quad \forall \, (k_0, k_1) \in \textbf{k}_0 \times \textbf{k}_1, \\ u_{j_0, j_1} &= \sum_{k_1\in \textbf{k}_1} \Big( e^{2\pi i \omega_1} \sum_{k_0\in\textbf{k}_0} \Big( e^{2\pi i \omega_0} \hat{u}_{k_0, k_1} \Big) \Big), \quad \forall \, (j_0, j_1) \in \textbf{j}_0 \times \textbf{j}_1.\end{split}\]

Note that the forward transform corresponds to taking the 1D Fourier transform first along axis 1, once for each of the indices in \(\textbf{j}_0\). Afterwords the transform is executed along axis 0. The two steps are more easily understood if we break things up a little bit and write the forward transform in (3) in two steps as

(4)\[\begin{split}\tilde{u}_{j_0,k_1} &= \frac{1}{N_1}\sum_{j_1 \in \textbf{j}_1} u_{j_0,j_1} e^{-2\pi i \omega_1}, \quad \forall \, k_1 \in \textbf{k}_1, \\ \hat{u}_{k_0,k_1} &= \frac{1}{N_0}\sum_{j_0 \in \textbf{j}_0} \tilde{u}_{j_0,k_1} e^{-2\pi i \omega_0}, \quad \forall \, k_0 \in \textbf{k}_0.\end{split}\]

The backward (inverse) transform if performed in the opposite order, axis 0 first and then 1. The order is actually arbitrary, but this is how is is usually computed. With mpi4py-fft the order of the directional transforms can easily be configured.

We can write the complete transform on compact notation as

(5)\[\begin{split}\boldsymbol{\hat{u}} &= \mathcal{F}(\boldsymbol{u}), \\ \boldsymbol{u} &= \mathcal{F}^{-1}(\boldsymbol{\hat{u}}).\end{split}\]

But if we denote the two partial transforms along each axis as \(\mathcal{F}_0\) and \(\mathcal{F}_1\), we can also write it as

(6)\[\begin{split}\boldsymbol{\hat{u}} &= \mathcal{F}_0(\mathcal{F}_1(\boldsymbol{u})), \\ \boldsymbol{u} &= \mathcal{F}_1^{-1}(\mathcal{F}_0^{-1}(\boldsymbol{\hat{u}})).\end{split}\]

Extension to multiple dimensions is straight forward. We denote a \(d\)-dimensional array as \(u_{j_0, j_1, \ldots, j_{d-1}}\) and a partial transform of \(u\) along axis \(i\) is denoted as

(7)\[\tilde{u}_{j_0, \ldots, k_i, \ldots, j_{d-1}} = \mathcal{F}_i(u_{j_0, \ldots, j_i, \ldots, j_{d-1}})\]

We get the complete multidimensional transforms on short form still as (5), and with partial transforms as

(8)\[\begin{split}\boldsymbol{\hat{u}} &= \mathcal{F}_0(\mathcal{F}_1( \ldots \mathcal{F}_{d-1}(\boldsymbol{u})), \\ \boldsymbol{u} &= \mathcal{F}_{d-1}^{-1}( \mathcal{F}_{d-2}^{-1}( \ldots \mathcal{F}_0^{-1}(\boldsymbol{\hat{u}}))).\end{split}\]

Multidimensional transforms are straightforward to implement in Numpy

import numpy as np
M, N = 16, 16
u = np.random.random((M, N))
u_hat = np.fft.rfftn(u)
uc = np.fft.irfftn(u_hat)
assert np.allclose(u, uc)

The fftw module

The fftw module provides an interface to most of the FFTW library. In the fftw.xfftn submodule there are planner functions for:

  • fftn() - complex-to-complex forward Fast Fourier Transforms

  • ifftn() - complex-to-complex backward Fast Fourier Transforms

  • rfftn() - real-to-complex forward FFT

  • irfftn() - complex-to-real backward FFT

  • dctn() - real-to-real Discrete Cosine Transform (DCT)

  • idctn() - real-to-real inverse DCT

  • dstn() - real-to-real Discrete Sine Transform (DST)

  • idstn() - real-to-real inverse DST

  • hfftn() - complex-to-real forward FFT with Hermitian symmetry

  • ihfftn() - real-to-complex backward FFT with Hermitian symmetry

All these transform functions return instances of one of the classes fftwf_xfftn.FFT, fftw_xfftn.FFT or fftwl_xfftn.FFT, depending on the requested precision being single, double or long double, respectively. Except from precision, the tree classes are identical. All transforms are non-normalized by default. Note that all these functions are planners. They do not execute the transforms, they simply return an instance of a class that can do it (see docstrings of each function for usage). For quick reference, the 2D transform shown for Numpy can be done using fftw as:

from mpi4py_fft.fftw import rfftn as plan_rfftn, irfftn as plan_irfftn
from mpi4py_fft.fftw import FFTW_ESTIMATE
rfftn = plan_rfftn(u.copy(), flags=(FFTW_ESTIMATE,))
irfftn = plan_irfftn(u_hat.copy(), flags=(FFTW_ESTIMATE,))
u_hat = rfftn(uc, normalize=True)
uu = irfftn(u_hat)
assert np.allclose(uu, uc)

Note that since all the functions in the above list are planners, an extra step is required in comparison with Numpy. Also note that we are using copies of the u and u_hat arrays in creating the plans. This is done because the provided arrays will be used under the hood as work arrays for the rfftn() and irfftn() functions, and the work arrays may be destroyed upon creation.

The real-to-real transforms are by FFTW defined as one of (see definitions and extended definitions)

  • FFTW_REDFT00

  • FFTW_REDFT01

  • FFTW_REDFT10

  • FFTW_REDFT11

  • FFTW_RODFT00

  • FFTW_RODFT01

  • FFTW_RODFT10

  • FFTW_RODFT11

Different real-to-real cosine and sine transforms may be combined into one object using factory.get_planned_FFT() with a list of different transform kinds. However, it is not possible to combine, in one single object, real-to-real transforms with real-to-complex. For such transforms more than one object is required.