Discrete Fourier Transforms

The Discrete Fourier Transform (DFT) is a discretized version of the Fourier transform, which is widely used in numerical simulation and analysis. Given a set of $N$ numbers $\{f_0, f_1, \dots, f_{N-1}\}$ , the DFT produces another set of $N$ numbers $N$ numbers $\{F_0, F_1, \dots, F_{N-1}\}$ , defined as follows:

: $\mathrm{DFT}\Big\{f_0, f_1, \dots, f_{N-1}\Big\} = \Big\{F_0, F_1, \dots, F_{N-1}\Big\} \qquad\mathrm{where}\quad F_n = \sum_{m=0}^{N-1} e^{-2\pi i \frac{mn}{N}}\, f_m.$

The inverse of this transformation is the Inverse Discrete Fourier Transform (IDFT):

: $\mathrm{IDFT}\Big\{F_0, F_1, \dots, F_{N-1}\Big\} = \Big\{f_0, f_1, \dots, f_{N-1}\Big\} \qquad\mathrm{where}\quad f_m = \frac{1}{N} \sum_{n=0}^{N-1} e^{2\pi i \frac{mn}{N}}\, F_n.$

The inverse relationship between the DFT and the IDFT is straightforward to prove, by using the identity

: $\sum_{m=0}^{N-1} e^{\pm 2\pi i \frac{m(n-n')}{N}} = N \delta_{nn'},$

where $\delta_{nn'}$ denotes the Kronecker delta. This identity is derived from the geometric series formula.

Conversion of continuous Fourier transform to DFT

The DFT is commonly encountered when discretizing formulas involving Fourier integrals. Recall the definition of the Fourier transform: given a function $f(x)$ , where $x \in (-\infty, \infty)$ , the Fourier transform is a function $F(k)$ , where $k \in (-\infty, \infty)$ , and these two functions are related by a pair of integral formulas:

: $\begin{align} F(k) &= \displaystyle \;\; \int_{-\infty}^\infty dx\; e^{-ikx}\, f(x) \\ f(x) &=\; \displaystyle \int_{-\infty}^\infty \frac{dk}{2\pi}\; e^{ikx}\, F(k).\end{align}$

Typically, a computer simulation or experimental measurement will produce values of $f(x)$ at certain values of $x$ that are discrete and evenly-spaced. Suppose these points are $\{x_0, x_1, \dots, x_{N-1}\}$ , where the spacing is $\Delta x = x_{m+1}-x_m$ ; the corresponding data points are $\{f(x_0), \dots, f(x_{N-1})\}$ . We are then interested in finding the Fourier spectrum, i.e. plotting either $|F(k)|$ or $|F(k)|^2$ versus $k$ . To do this, we can approximate the Fourier integral by using the mid-point rule:

: $F(k) \approx \Delta x \sum_{m=0}^{N-1} e^{-ik x_m} \, f(x_m).$

Note that this necessitates a truncation of the Fourier integral. The Fourier integral ran over $-\infty < x < \infty$ , but our numerical integral runs over a finite range $x_0 \lesssim x \lesssim x_{N-1}$ . This truncation will have important consequences later. Now, we have to decide the values of $k$ at which to find $F(k)$ . Let us choose a set of $N$ equally-spaced points,

: $k_n \equiv \frac{2\pi n}{N\Delta x}.$

At these points, the discretized Fourier integral takes the form

: $\begin{align}F(k_n) &\approx \Delta x\, \sum_{m=0}^{N-1} \exp\left[- \frac{2 \pi i n (x_0 + m \Delta x)}{N \Delta x}\right] \, f(x_m) \\ &= \Delta x\, \exp\left[- \frac{2 \pi i n x_0}{N \Delta x}\right] \sum_{m=0}^{N-1} e^{-i 2 \pi n m / N} \, f(x_m) \\ &= \Delta x\, \exp\left[- \frac{2 \pi i n x_0}{N \Delta x}\right] \mathrm{DFT}\{f(x_m)\}_n.\end{align}$

Here $\mathrm{DFT}\{f(x_m)\}_n$ denotes the $n$ -th element of the Discrete Fourier Transform (DFT). The $m$ index inside the curly brackets is a dummy index, indicating that the DFT involves an internal sum over this index (we're slightly abusing mathematical notation here). The phase factor, $\exp[- 2 \pi i n x_0/N \Delta x]$ , is determined by the choice of "origin" for the spatial coordinates; it does not affect $|F(k_n)|^2$ (which is what's used to plot the Fourier spectrum).

The DFT and IDFT can be computed very efficiently, in $O(N\log N)$ time, using an algorithm called the Fast Fourier Transform (FFT). We will not discuss the FFT algorithm in this article, but many good explanations can be found elsewhere online. In Python, you can perform an FFT (fast DFT) by calling fft, and an inverse FFT (fast IDFT) by calling ifft

Spectral resolution and range

In the previous section, we showed how a continuous Fourier integral is converted into a DFT. This process involved two distinct approximations. Firstly, the Fourier integral is truncated from its original range, $x \in (-\infty,\infty)$ , to a finite interval of length $N \Delta x$ . Secondly, the integral is discretized by reducing the continuous variable $x$ to a set of discrete points $\{x_0,\dots,x_{N-1}\}$ . Both of these approximations have important consequences for the accuracy of our numerical Fourier spectrum, which we will examine in turn.

Spectral resolution

Fig. 1: Fourier power spectrum from a truncated Fourier transform of $f(x) = \exp(ik_0 x)$ , with $k_0 = 1$ and sampling interval $k \in [0,300]$ (about 48 periods).

The truncation of the Fourier integral limits the spectral resolution of the Fourier spectrum. To see this, suppose we perform truncation without discretization, by taking a continuous Fourier integral and truncating it to a finite range $x \in [0, X]$ :

: $F(k) \;\;\approx\;\; \int_0^X dx\; e^{-ikx} \, f(x).$

Consider a harmonic function $f(x) = e^{ik_0 x}$ . The exact Fourier transform can be shown to be a delta function, $F(k) = 2\pi \, \delta(k - k_0)$ , i.e. an infinitely sharp peak centered at $k = k_0$ . With the above truncation, however, the resulting integral is

: $F(k) \;\;\approx\;\; \int_0^X dx \; e^{-i(k - k_0)x} \;\; =\;\; \frac{2\sin[(k-k_0)X/2]}{k-k_0} \cdot e^{-i(k-k_0)X/2}.$

For $X \rightarrow \infty$ , the above formula approaches a delta function (an infinitesimally-thin peak) centered at $k = k_0$ . But for finite $X$ , the plot of $|F(k)|^2$ versus $\omega$ behaves as shown in Fig 1. Evidently, truncating the Fourier integral has "smeared out" the Fourier spectrum, broadening the infinitesimally-thin delta function peak into a finite-width peak. The peak width, $\Delta k \sim 1/X$ , limits the "resolution" of our Fourier analysis.

In the discretized Fourier transform, the truncation of the Fourier integral has essentially the same effect. As discussed in the previous section, the DFT is defined at $k_n \equiv 2\pi n/X$ ; hence, the resolution of the Fourier spectrum is $\Delta k = 2\pi/X$ .

Spectral range

The other approximation which we made in going from a continuous Fourier transform to the DFT involved sampling $f(x)$ at discrete values of $x$ . This discretization has the effect of limiting the spectral range. To see this, let us look again at the DFT formula, which is dimensionless:

: $F_n = \sum_{m=0}^{N-1} e^{-2\pi i \frac{mn}{N}}\, f_m.$

Normally, we consider only the indices $n = 0, 1, \dots, N-1.$ However, if we replace $n$ with $n+N$ in the right-hand side, the result would be the same:

: $F_{n+N} = \sum_{m=0}^{N-1} e^{-2\pi i \frac{m(n+N)}{N}}\, f_m = \sum_{m=0}^{N-1} e^{-2\pi i \frac{mn)}{N} - 2\pi i m}\, f_m = F_n.$

We can hence regard $F$ as a periodic discrete function of $n$ , with period $N$ . Next, consider how the DFT is related to the physical $x$ and $k$ variables. Taking $x_0 = 0$ for simplicity,

: $F(k_n) = \sum_{m=0}^{N-1} e^{-2\pi i \frac{mn}{N}}\, f(x_m) = \sum_{m=0}^{N-1} e^{-i k_n \,\cdot\, m \Delta x}\, f(x_m).$

If we perform the replacement

: $k_n \rightarrow k_n + K, \quad \mathrm{where}\;\; K \equiv\frac{2\pi}{\Delta x},$

then evidently $F(k_n)$ is left unchanged. Indeed, we could add any integer multiple of $K$ without altering the result. This means that the DFT spectrum is only defined under $k$ modulo $K$ , by contrast with the continuous Fourier transform which is defined over the entire interval $-\infty < k < \infty$ .

The default definition of the DFT gives the integer indices $n = 0, 1, \dots, N-1$ , which corresponds to $0 \le k \lesssim K$ . However, when plotting the DFT spectrum, we usually adjust the range of $k$ to $-K/2 \lesssim k \lesssim K/2$ . This is done by taking the "upper half" of the DFT spectrum, $K/2 \lesssim k \lesssim K$ , and translating it via the replacement $k \rightarrow k - K$ . Due to the periodicity of the DFT, the upper half of the DFT spectrum becomes the negative $k$ part of the spectrum. In terms of the integer indices $n$ , the process is depicted in the figure below:

Fig. 2: A DFT spectrum, $F_n$ , is periodic with period $N$ . By default, the DFT is reported in the spectral range $n \in [0, N-1]$ (red curve in the upper plot). To relate this to the continuous Fourier transform, we re-center the spectrum at $n = 0$ , which is equivalent to translating the upper half of the spectrum to negative values (red curve in the lower plot).

The reason for this adjustment is that, intuitively, the discretized Fourier spectrum contains information about the "low-frequency" part of the spectrum, $|k| < K/2$ , including both positive and negative values of $k$ . On the other hand, the discretized Fourier spectrum lacks information about the "high frequency" part of the spectrum, which correspond to harmonics with periods shorter than the discretization step $\Delta x$ . Hence, it makes sense to "center" our Fourier spectrum around the origin. It can then be shown that as the discretization step approaches zero (and hence $K = 2\pi/\Delta x \rightarrow \infty$ ), the $|k| \ll K$ part of the adjusted DFT spectrum converges to the exact (continuous) Fourier spectrum.

The corollary to the above discussion is that if we have a function which has no frequency components larger than $k_\mathrm{max}$ , then it is sufficient to use a sampling interval $\Delta x = \pi / k_\mathrm{max}$ . This is called the Nyquist-Shannon sampling theorem.

Summary of spectral relations

The results of the previous sections can be summarized in this way:

The total range of $x$ , which is denoted by $X$ , limits the resolution of the spectrum to $\Delta k = 2\pi/X$ .
The resolution of $x$ , which is the discretization step $\Delta x$ , limits the range of the spectrum to $K = 2\pi/\Delta x$ .

These relations are easy to remember, because the "interval length" in the one domain places a limit on the "discretization step" in the other domain. It is very important to keep these relations in mind when working with discrete Fourier transforms! For example, a common mistake that people make is to try to improve the resolution of a Fourier spectrum by increasing the number of discretization steps, $N$ , while keeping the total interval $X$ fixed. This doesn't work; it leaves the spectral resolution unchanged! In order to improve the spectral resolution, one has to increase the total interval instead.

The split-step Fourier method

As an example of the usefulness of the DFT, let us discuss a DFT-based method for performing numerical integration of a partial differential equation, known as the split-step Fourier method. Here, the method will be presented in the context of the time-dependent Schrödinger equation in 1D space:

: $i \frac{d\psi(x,t)}{dt} = \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t)\right]\psi(x,t)$

We have taken $\hbar = m = 1$ for simplicity. At each time, the wavefunction is a continuous function of $x$ . Let us truncate and discretize this spatial coordinate, by defining a computational domain of length $L$ containing $N$ discretization points:

: $\psi_n(t) = \psi(x_n,t), \;\;\;\textrm{where}\;\; x_n = -\frac{L}{2} + n \Delta x, \;\; \Delta x = \frac{L}{N}.$

Thus, the wavefunction at each time is represented by a complex vector, which we call a "state vector":

: $\vec{\psi}(t) = \begin{pmatrix} \psi_0(t) \\ \psi_1(t) \\ \vdots \\ \psi_{N-1}(t)\end{pmatrix}.$

Given an initial state vector $\vec{\psi}(t_a)$ , the problem is to compute $\vec{\psi}(t_b)$ at a later time $t_b$ . Note that this differs from previously-studied numerical ODE problems in one important respect: evolving in time involves taking second-order spatial derivatives. We won't go into the details, but it turns out that standard methods for time-stepping and discretizing space don't work very well here, because the errors from time-stepping and spatial discretization interact badly with one another. The split-step Fourier method provides a better way to solve the problem.

Factorizing the time-evolution operator

The split-step Fourier method is based on the concept of the time-evolution operator. Given a wavefunction $\psi(x,t_j)$ , the wavefunction after a small time step $\tau$ is

: $\psi(x,t_j + \tau) = U(t_j + \tau|t_j)\, \psi(x,t_j), \quad\mathrm{where}\;\;\; U(t_j + \tau|t_j) \approx \exp\left\{-i \tau \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t_j+\tau/2)\right]\right\}\;.$

Here, the $\exp(\cdots)$ refers to the exponential of an operator (one involving spatial derivatives). We call $U(t_b|t_a)$ the time-evolution operator, which evolves the system from time $t_a$ to $t_b$ . The exponential of any operator $\mathcal{A}$ is defined as the infinite series

: $\exp(\mathcal{A}) = I + \mathcal{A} + \frac{1}{2} \mathcal{A}^2 + \frac{1}{6} \mathcal{A}^3 + \cdots$

In this case, the exponential contains the Hamiltonian, which consists of a kinetic energy term and a potential energy term that do not generally commute. Due to this non-commutivity, the exponential cannot be simplified by factorization:

: $\exp\left\{-i \tau \left[-\frac{1}{2}\frac{\partial^2}{\partial x^2} + V(x,t_j+\tau/2)\right]\right\} \;\; \ne\;\; \exp\Bigg\{\frac{i \tau}{2} \left[\frac{\partial^2}{\partial x^2}\right]\Bigg\} \; \exp\Bigg\{-i \tau V(x,t_j+\tau/2)\Bigg\}.$

However, we can obtain an approximate factorization by making use of the series definition of the exponential of an operator. One can show that

: $\exp(\mathcal{A} + \mathcal{B}) = \exp(\mathcal{A}/2) \exp(\mathcal{B}) \exp(\mathcal{A}/2) + O\Big((\mathcal{A},\mathcal{B})^3\Big),$

which is a variant of an important formula known as the Baker–Campbell–Hausdorff formula. On the right-hand side, note that $\exp(\mathcal{B})$ is sandwiched "symmetrically" between two copies of $\exp(\mathcal{A}/2)$ . This symmetric arrangement reduces the approximation error to third order, by the cancellation of lower-order errors (in a manner similar to the mid-point formula for the discretized derivative). Applying this factorization to the time-evolution operator gives

$U(t_j + \tau|t_j) \;\; \approx\;\; U_\mathcal{K} \,\cdot\, U_\mathcal{V}(t_j+\tau/2) \,\cdot\, U_\mathcal{K}, \quad \mathrm{where}\;\; \begin{cases} \;U_\mathcal{K} &= \exp\left[\frac{i\tau}{4} \frac{d^2}{dx^2}\right] \\ \; U_\mathcal{V}(t) &= \exp\big[-i\tau\, V(x,\,t)\big]\end{cases}$

In other words, the time-evolution operator decomposes into three pieces. That's why we call this a "split-step" algorithm: each time step from $t_j$ to $t_j+\tau$ consists of applying a kinetic step, then applying a potential step, then applying another kinetic step, in sequence. As previously noted, we'll be working with state vectors (complex $N$ -component vectors), defined through spatial discretization of the wavefunction. So we need to figure out how the above stepping operators act on these state vectors:

: $U_\mathcal{K}\; \vec{\psi} \;= \;\;??, \qquad U_\mathcal{V}\, \vec{\psi} \;=\;\; ??$

The potential stepping operator is simple to deal with. Since the state vector represents the wavefunction at different points in space, the potential operator is represented by a diagonal matrix, and its exponential is also diagonal:

: $U_\mathcal{V} \begin{pmatrix}\psi_0\\ \vdots \\ \psi_{N-1}\end{pmatrix} \;=\; \begin{pmatrix}\exp\left[-i\tau V(x_0,t+\frac{\tau}{2})\right] & & \\ &\ddots&\\ && \exp\left[-i\tau V(x_{N-1},t+\frac{\tau}{2})\right] \end{pmatrix} \begin{pmatrix}\psi_0\\ \vdots \\ \psi_{N-1}\end{pmatrix} \;=\; \begin{pmatrix}\exp\left[-i\tau V(x_0,t+\frac{\tau}{2})\right] \psi_0\\ \vdots \\ \exp\left[-i\tau V(x_{N-1}, t+\frac{\tau}{2})\right]\psi_{N-1}\end{pmatrix}$

Kinetic step

The kinetic stepping operator, $U_\mathcal{K}$ , is less obvious. It contains spatial derivatives and is thus not diagonal in the current basis. The key thing to realize, however, is that this operator is diagonal in wavenumber space. Let us return to the continuous wavefunction, and write its Fourier representation:

: $\psi(x) = \int_{-\infty}^\infty \frac{dk}{2\pi} \; e^{ikx}\; \Psi_k.$

Then

: $U_\mathcal{K} \; \psi(x) = \int_{-\infty}^\infty \frac{dk}{2\pi} \; \exp\left(-\frac{i\tau}{4}k^2\right) \; e^{ikx}\; \Psi_k.$

Let us discretize space in steps of $\Delta x$ , as discussed earlier, and also discretize the Fourier integrals by steps of $\Delta k$ :

: $\begin{align} x_n &= -\frac{L}{2} + n \Delta x, \\ k_n &= -\frac{K}{2} + n\Delta k \end{align}$

The values of $K$ and $\Delta k$ will be chosen shortly. The discretized integrals become

: $\begin{align}\psi_m &\approx \sum_{n=0}^{N-1} \frac{\Delta k}{2\pi} \; e^{ik_n x_m}\; \Psi_{k_n}, \\ \left(U_\mathcal{K} \psi\right)_m &\approx \sum_{n=0}^{N-1} \frac{\Delta k}{2\pi} \; e^{-\frac{i\tau}{4}k_n^2} \;e^{ik_n x_m}\; \Psi_{k_n}. \end{align}$

Let us now choose the k-space discretization parameters to be

: $\Delta k = \frac{2\pi}{N\Delta x} = \frac{2\pi}{L} , \quad K = N\Delta k = \frac{2\pi}{\Delta x} \;\;\;\Rightarrow \;\;\; k_n = -\frac{\pi}{\Delta x} + \frac{2\pi n}{N\Delta x}.$

With this choice, we can show with a bit of algebra that the integral for $\psi_m$ reduces to an IDFT:

: $\begin{align}\psi_m &= (-1)^m \frac{1}{N} \sum_{n=0}^{N-1} \left(\frac{1}{\Delta x} e^{iN\pi/2}\;e^{-in\pi}\; \Psi_{k_n}\right) e^{\frac{2\pi i mn}{N}}\\ &= (-1)^m\; \mathrm{IDFT}\left\{\frac{1}{\Delta x} e^{iN\pi/2}\;(-1)^n \;\Psi_{k_n}\right\}_m \\ \Rightarrow \Psi_{k_n} &= \Delta x\, e^{-iN\pi/2}\; (-1)^n \; \;\mathrm{DFT}\Big\{(-1)^m\; \psi_m \Big\}_n \end{align}$

Likewise,

: $(U_\mathcal{K}\; \psi)_m = (-1)^m \; \mathrm{IDFT}\left\{\frac{1}{\Delta x} e^{iN\pi/2}\;(-1)^n e^{-\frac{i\tau}{4}k_n^2}\; \Psi_{k_n}\right\}_m$

Putting these results together, we get

$\begin{align}(U_\mathcal{K}\; \psi)_m = (-1)^m \;\mathrm{IDFT}\Bigg\{ e^{-\frac{i\tau}{4}k_n^2}\; \mathrm{DFT}\Big\{(-1)^{p}\, \psi_{p}\Big\}_n \Bigg\}_m \\ \mathrm{where}\;\; k_n = - \frac{\pi N}{L} + \frac{2\pi n}{L} \end{align}$

Hence, the $U_{\mathcal{K}}$ kinetic stepping operator can be implemented by taking a DFT, multiplying the resulting vector elements by $\exp(-i\tau k_n^2/4)$ phase factors, and taking an IDFT. The runtime of the stepping process is $O(N\log(N))$ . The $m$ , $n$ , and $p$ indices all run over the range $[0, 1, \cdots, N-1]$ , consistent with the standard definition of the DFT and IDFT.