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,

or in terms of complex exponentials,

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,

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,

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

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,

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.

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}$).