# 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. 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

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

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 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.