I was given a task to obtain a Fourier Series approximation from the DFT (Discrete Fourier Transorm), more exactly fft function from MATLAB. I know that a Fourier Series has the form $\frac{a_0}{2}+\sum_{k=1}^{k=\infty}(a_k\cos(\dfrac{\pi kx}{L})+b_k\sin(\dfrac{\pi kx}{L}))$ if the period is $2L$. Also, the Discrete Fourier Transform for a given set of values $\{x_1,x_2,\dotsc x_N\}$ gives a set of complex numbers ${X_1, X_2, \dotsc X_N}$. While doing this, we assume the function has period of length $N-1$. From the disscusion I had, I've drawn the conclusion that the real and imaginary part of $X_k$ give us $a_k$ and $b_k$. In the end, this was proved wrong as the approximation I get from the Fourier series obtained does not resemble my initial function.
Therefore, is there anyone who could clarify the relation between the values from the Discrete Fourier Transform and the coefficients from a Fourier Series approximation?
Thank you!
Best Answer
A $1$-periodic function, discretized (sampled) at step $h=1/n$, is a collection of $n$ values $y_0,\dots,y_{n-1}$. The functions $f_k = \exp(2\pi i k/n)$, $k=0,\dots,n-1$, form an orthogonal basis among such functions, so there is an expansion $f=\sum c_k f_k$ where $$ c_k = \frac{1}{n} \sum_{\ell} y_\ell \exp(-2\pi i k\ell/n) \tag{1} $$ Except for $1/n$, this computation amounts to the multiplication of vector $y$ by the Fourier matrix $W$ with the entries $\exp(-2\pi i kj)$ (if we index the rows and columns starting with 0). The Discrete Fourier transform is this multiplication (algorithmically,
fftarranges it in a clever way, reusing previous intermediate results to speed it up.)However, one should not jump to the conclusion that $$ y = \sum_{k=0}^{n-1} c_k \exp(2\pi i k t) \tag{2} $$ is the right way to interpolate the sampled values. Although the function on the right does interpolate $y_0,\dots, y_{n-1}$, it has unnatural oscillations in between. One should consider the aliasing of frequencies and fold them appropriately.
For example, take $n=5$, so that (2) says $$ c_0+ c_1e^{2\pi i t} + c_2e^{4\pi i t} + c_3e^{6\pi i t} + c_4e^{8\pi i t} \tag{3} $$ The aliasing of frequencies is the fact that $6 \pi t \equiv -4 \pi t \bmod 2\pi$ when $t$ is a grid point (one of $k/5$ points). Similarly, $8\pi t \equiv -2\pi t \bmod 2\pi $. So, we can get smoother curve by folding the frequencies, replacing (3) with $$ c_0+ c_1e^{2\pi i t} + c_2e^{4\pi i t} + c_3e^{ - 4\pi i t} + c_4e^{-2\pi i t} \tag{4} $$ The equation (4) is the correct interpolant, and it can be turned into sine-cosine Fourier series (rather, its partial sum) using Euler's formula.
Note that for real functions, the vector of coefficients is conjugate-symmetric: $c_{n-k}=\overline {c_k}$. So one can take the real part of the terms with $c_1, c_2$ and double them; this accounts for $c_3, c_4$. Separate treatment is required for the constant term $c_0$ (which has no counterpart) and, when $n$ is even, for the Nyquist frequency $k=n/2$, which also has no counterpart.