Numerical Methods
|
|
Numerical Methods | |
|
|
Discrete Fourier TransformFourier 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, \[ \begin{equation} \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], \end{equation} \]or in terms of complex exponentials, \[ \begin{equation} I(x) = \sum\limits_{k=-\infty}^{\infty}Y_k\exp \left(-\frac{i2\pi kx}{n\Delta x}\right). \end{equation} \]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, \[ \begin{equation} a_k = Y_{k}+Y_{-k}\qquad \text{and}\qquad b_k = Y_{k}-Y_{-k}. \end{equation} \]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, \[ \begin{equation} \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}. \end{equation} \]To determine the Fourier coefficients $Y_k$, multiply by $e^{i2\pi k'j/n}$ and sum over $j$. \[ \begin{equation} \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}. \end{equation} \]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, \[ \begin{equation} \label{eq:DFT} Y_k = \frac{1}{n} \sum\limits_{j=0}^{n-1}y_je^{i2\pi kj/n}. \end{equation} \]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}).
| |||||||||||||||||||||||||||||||||||||||||||||||||