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
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.
A more compact notation is commonly used for the DFTs, where the 1D forward and backward transforms are written as
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 stright forward 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
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
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
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
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
We get the complete multidimensional transforms on short form still as (5), and with partial transforms as
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 Transformsifftn()
- complex-to-complex backward Fast Fourier Transformsrfftn()
- real-to-complex forward FFTirfftn()
- complex-to-real backward FFTdctn()
- real-to-real Discrete Cosine Transform (DCT)idctn()
- real-to-real inverse DCTdstn()
- real-to-real Discrete Sine Transform (DST)idstn()
- real-to-real inverse DSThfftn()
- complex-to-real forward FFT with Hermitian symmetryihfftn()
- 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.