Discrete Fourier Transformation (DFT)
In the following, we consider a periodic function that has been discretely sampled on an evenly distributed set of \(N\) sample points contained within a single period \(L\). Such points will have the coordinates:
$$ x_{k}=k\times\frac{L}{N}\qquad k=0,\cdots,N-1 $$
The trigonometric polynomial approximation evaluated at the sample points should reproduce the function values:
$$ f(x)=\sum_{j=0}^{N-1}F_{j}\exp\left(\mathrm{i}\frac{2\pi j}{L}x\right)\rightarrow f_{k}=\sum_{j=0}^{N-1}F_{j}\exp\left(\mathrm{i}\frac{2\pi j}{L}x_{k}\right) $$
This gives rise to the following relation between the function values at the sample points and the Fourier coefficients:
$$\require{color}$$ $$ \fcolorbox{red}{#f2d5d3}{ $f_{k}=\sum_{j=0}^{N-1}F_{j}\exp\left(\mathrm{i}\frac{2\pi jk}{N}\right)$ } $$
Which defines the Discrete Fourier Transformation (DFT)!
Our main task now it to numerically calculate Fourier coefficients \(F_{j}\) for a given set of discrete function values.
Before we proceed we introduce a shorthand notation that will make the expressions more transparent.
We introduce the shorthand notation for the complex \(N-\)th root of unity:
$$W_{N}=\exp\left[\frac{2\pi \mathrm{i}}{N}\right]$$
We furthermore define a matrix \(\mathbf M\) with elements:
$$ M_{kj} = W_{N}^{kj} = \exp\left[\frac{2\pi \mathrm{i} jk}{N}\right]. $$
Using this notation the discrete Fourier Transform can be written as:
$$ f_{k}=\sum_{j=0}^{N-1}M_{kj}F_{j} $$ which is nothing else than a matrix-vector product.
Thus the discrete Fourier transformation is a linear matrix transformation of a given vector or function values written compactly as:
$$\mathbf {f} = \mathbf {M}\mathbf {F}$$
The matrix transformation transforms the vector of Fourier coefficients \(\mathbf{F} \) into the vector of function values \( \mathbf{f} \) at the sample points!
In order to find the values of the Fourier coefficients we now have two distinct possibilities:
- Inverting the matrix \( \mathbf{M} \) such that the Fourier coefficients can be determined as \( \mathbf {F} = \mathbf {M}^{-1}\mathbf {f} \).
- Calculate Fourier coefficients by some discrete approximation for the integral $$ c_{k}=\frac{1}{L}\int_{0}^{L}\mathrm{d}x\ f(x)\exp\left(-\mathrm{i}\frac{2\pi k}{L}x\right) $$ We first follow the route 1 and prove by explicit calculation that the matrix \( \mathbf M \) is up to a constant of \( N \) unitary.
For this purpose we calculate a general matrix element as follows:
$$\mathbf{M}^{\dagger}\mathbf{M}=?$$
$$\left(\mathbf{M}^{\dagger}\mathbf{M}\right)_ {ij}=\sum_{k=0}^{N-1}\left(\mathbf{M}^{\dagger}\right)_ {ik}\left(\mathbf{M}\right)_ {kj}=\sum_{k=0}^{N-1}\left(W_{N}^{ * }\right)^{ik}W_{N}^{kj} $$
By exploiting properties of the complex roots of unity we can easily prove that:
$$(W_{N}^{ * })^{ik}=\left(W_{N}^{ik}\right)^{ * }=W_{N}^{-ik}$$
This can be utilized to rewrite the matrix element as:
$$\left(\mathbf{M}^{\dagger}\mathbf{M}\right)_ {ij}=\sum_{k=0}^{N-1}W_{N}^{k(j-i)}$$
We now consider two cases:
Case 1: \( i = j \):
$$\left(\mathbf{M}^{\dagger}\mathbf{M}\right)_ {ii}=\sum_{k=0}^{N-1}1=N$$
Case 2: \( i \ne j\):
$$\left(\mathbf{M}^{\dagger}\mathbf{M}\right)_ {ij}=\sum_{k=0}^{N-1}\left[W_{N}^{(j-i)}\right]^{k}=\frac{1-\left(W_{N}^{N}\right)^{j-i}}{1-W_{N}^{(j-i)}}=0$$
Here, we have utilized the formula for the partial sum of a geometric series as well as the following identity:
$$W_{N}^{N}=1$$
This shows that the matrix $\mathbf M$ has an inverse determined by the relation:
$$\mathbf{M}^{\dagger}\mathbf{M}=N\mathbf{I}\rightarrow \mathbf{M}^{-1}=\frac{1}{N}\mathbf{M}^{\dagger}$$
This can now be utilized in order to perform the discrete Fourier transformation and calculate the Fourier coefficients:
$$\mathbf{f}=\mathbf{M}\mathbf{F}\rightarrow\mathbf{F}=\mathbf{M}^{-1}\mathbf{f}=\frac{1}{N}\mathbf{M}^{\dagger}\mathbf{f}$$
Explicitly in terms of the input vector components the Fourier coefficients read:
$$\fcolorbox{red}{#f2d5d3} {$F_{k}=\frac{1}{N}\sum_{j=0}^{N-1}f_{j}\exp\left(-\mathrm{i}\frac{2\pi jk}{N}\right)$}$$
Together with the inverse relation:
$$ \fcolorbox{green}{#d3f2d9}{$ f_{k}=\sum_{j=0}^{N-1}F_{j}\exp\left(\mathrm{i}\frac{2\pi jk}{N}\right)$}$$
These expressions fully specify the discrete Fourier transformation and its inverse!
What if we would have followed the route 2 and calculate the Fourier coefficients by approximating the integral:
$$F_{k}=\frac{1}{L}\int_{0}^{L}\mathrm{d}x\ f(x)\exp\left(-\mathrm{i}\frac{2\pi k}{L}x\right)$$
using some of the numerical integration procedures (e. g. trapezoidal or Simpson's rule)?
Most numerical integration schemes can be recast in the following general from by introducing the integration weights for each sample point:
$$\int_{0}^{L}f(x)\mathrm{d}x=\frac{L}{N}\sum_{j=0}^{N-1}w_{j}f(x_{j})$$
If we now apply this to the Fourier coefficients we get:
$$ F_{k} = \frac{1}{L}\int_{-\frac{L}{2}}^{+\frac{L}{2}}\mathrm{d}x\ f(x)\exp\left(-\mathrm{i}\frac{2\pi k}{L}x\right) \rightarrow F_{k}\approx\frac{1}{N}\sum_{j=0}^{N-1}w_{j}f_{j}\exp\left(-\mathrm{i}\frac{2\pi jk}{N}\right). $$
Thus, different numerical integration schemes would generally lead to different Fourier coefficients.
If we choose all weights to be equal to one, which corresponds to the usual trapezoidal rule we recovere the Fourier coefficients that we have determined along the route 1 exploiting the unitarity of the \(\mathbf M\)-matrix!
The advantage of this choice in contrast to all other possible integration scheme is that in this way the values of the input function at the sample points are exactly recovered if the inverse transformation is performed. In order to show this let us calculate the value of the trigonometric polynomial at the point \(x_{i}\):
$$f(x_{i})=\sum_{k=0}^{N-1}F_{k}W_{N}^{ik}$$
By inserting the Fourier coefficient obtained by a general numerical integration scheme:
$$F_{k}=\frac{1}{N}\sum_{j=0}^{N-1}w_{j}f_{j}W_{N}^{-jk}$$
we get:
$$f(x_{i})=\frac{1}{N}\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}w_{j}f_{j}W_{N}^{(i-j)k}$$
which results in the following value of the trigonometric polynomial at the sampling point \(x_i\):
$$f(x_{i})=\frac{1}{N}\sum_{j=0}^{N-1}w_{j}f_{j}\sum_{k=0}^{N-1}W_{N}^{(i-j)k}=\sum_{j=0}^{N-2}w_{j}f_{j}\delta_{ij} = w_{i} f_{i}$$
We see that only in the case where \(w_{i}=1\) the trigonometric interpolation will be exact at the sampling points.
$$f(x_{i})=f_{i}$$
Since we anyway do not know anything about the function between the sampling points it is reasonable to require from the interpolation scheme that it reproduces the exact values at the sampling points!