Numerical Methods

## Discrete Fourier Transform

Fourier analysis is a form of interpolation that uses periodic functions to interpolate between discrete data points. Consider a set of $n$ data points $(x_j,y_j)$ where the points are equally spaced in $x$: $x_j = j\Delta x$. We seek an interpolation periodic function $I(x)$ with a period $n\Delta x$ that passes through all of the points. This function can be expressed as a Fourier series which can either be written in terms of sines and cosines,

$$$\label{eq:sines} I(x) = \sum\limits_{k=0}^{\infty}\left[a_k\cos\left(\frac{2\pi kx}{n\Delta x}\right)+b_k\sin\left(\frac{2\pi kx}{n\Delta x}\right)\right],$$$

or in terms of complex exponentials,

$$$I(x) = \sum\limits_{k=-\infty}^{\infty}Y_k\exp \left(-\frac{i2\pi kx}{n\Delta x}\right).$$$

The equivalence of these two expressions can be shown using Euler's formula: $e^{ix}=\cos x+i\sin x$. The coefficients can be complex and are related by,

$$$a_k = Y_{k}+Y_{-k}\qquad \text{and}\qquad b_k = Y_{k}-Y_{-k}.$$$

If $I(x)$ is a real function, then $a_k$ and $b_k$ are real and $Y_k = Y_{-k}^*$ so that $a_k = 2\mathcal{Re}(Y_{k})$ and $b_k = 2\mathcal{Im}(Y_{k})$.

Although it is possible to construct many periodic functions that go through all the data points $y_n$ with wavelengths smaller than $\Delta x$, it makes sense to restrict the sum over $k$ to the fewest necessary to go through all the points. There are two common choices: $k$ can be restricted to the range $k=0,1,2,\cdots,n-1$ or $k=-n/2,\cdots,-2,-1,0,1,2,\cdots,n/2$. The first choice appears more often in the literature but the second choice produces the smoothest interpolation function $I(x)$. To keep the discussion general, we restrict the sum over $k$ to the range $k=q,q+1,q+2,\cdots,q+n-1$ and the starting value $q$ can be specified later.

The points $(x_j,y_j)$ are substituted into expression for the Fourier series,

$$$\label{eq:iDFT} y_j = \sum\limits_{k=q}^{q+n-1}Y_k\exp \left(-\frac{i2\pi kx_j}{n\Delta x}\right)=\sum\limits_{k=q}^{q+n-1}Y_ke^{-i2\pi kj/n}.$$$

To determine the Fourier coefficients $Y_k$, multiply by $e^{i2\pi k'j/n}$ and sum over $j$.

$$$\sum\limits_{j=0}^{n-1}y_je^{i2\pi k'j/n} = \sum\limits_{j=0}^{n-1}\sum\limits_{k=q}^{q+n-1}Y_ke^{i2\pi (k'-k)j/n}.$$$

In the double sum on the right side, consider a specific value of $k$ and sum over $j$. This sum is zero unless $k=k'$ so the double sum evaluates to $nY_k$. This yields an expression for the Fourier coefficients,

$$$\label{eq:DFT} Y_k = \frac{1}{n} \sum\limits_{j=0}^{n-1}y_je^{i2\pi kj/n}.$$$

Equation (\ref{eq:DFT}) is called the Discrete Fourier Transform (DFT) of the data series $y_j$ and Eq. (\ref{eq:iDFT}) is known as the inverse discrete Fourier transform.

The form below accepts a sequence of complex numbers $y_j$ and calculates the discrete Fourier transform. It is also possible to input the values of $Y_k$ on the right and calculate the inverse discrete Fourier transform. The smoothest fit is obtained for $q^*=\text{Int}(-n/2+1)$ where $\text{Int} (x)$ rounds down to the nearest integer. The $q^*$ button calculates this value.

 Real space Reciprocal space $y_j$ $Y_k$ $j$ $k$ $q=$ $\text{Re}[y_j]$  $\text{Im}[y_j]$ $\text{Re}[Y_k]$  $\text{Im} [Y_k]$ 0 0 1 0 2 0 3 0 2 0 1 0 0 0

Below is the code that calculates a discrete Fourier transform.

If the data points $y_j$ are real, then real coefficients $a_k$ and $b_k$ can be found for the Fourier series in terms of sines and cosines, Eq. (\ref{eq:sines}).