Complex Fourier Decomposition

Question

I'm currently working on a python project for which I need to use the complex Fourier decomposition of a complex function. Note that this function is numerical as I don't have the analytical expression, only a list of values. I used the fft from scipy to get the Fourier coefficients and noticed that the modulus of the coefficients were not decreasing. So here is my question: to best approximate the function with a finite number of harmonies, how should I select the harmonies ? After looking up the internet (and my past knowledge) I would have choosen the first harmonies (and if that, should I select the harmonies $[0,1,-1,2,-2,...]$ or $[0,-1,1,-2,2,...]$), but now I wonder wether I should better take the harmonies with the highest amplitudes (it should be the ones with most impact, right?).

Answer 1

The mean squared error between a partial reconstruction $(z_n)_n$ and the true reconstruction $(y_n)_n$ of the samples is given by $$MSE=\frac{1}{N}\sum_{n=0}^{N-1}|y_n-z_n|^2=\frac{1}{N}\sum_{n=0}^{N-1}\bigg|\frac{1}{N}\sum_{k=-(N-1)/2}^{(N-1)/2}\hat{y}_ke^{2\pi i \frac{kn}{N}}-\frac{1}{N}\sum_{j \in J}\hat{y}_je^{2\pi i \frac{jn}{N}}\bigg|^2$$ The set $J$ contains the selected frequencies among $\{-(N-1)/2,...,-1,0,1,...,(N-1)/2\}$ that we want for our partial reconstruction. If we constrain $\textrm{len}(J)<N$ then we have to choose the DFT terms $\hat{y}_j$ that minimize the MSE. Since every $\hat{y}_j$ included in $z_n$cancels out with a component of $y_n$, if we set the length of $J$ we want the MSE to be made of the $N-J$ components that account for the smallest MSE value possible. For example if $\textrm{len}(J)=N-2$ then $$\begin{aligned}MSE&=\frac{1}{N^3}\sum_{n=0}^{N-1}|\hat{y}_{k_1}e^{2\pi i \frac{k_1n}{N}}+\hat{y}_{k_2}e^{2\pi i \frac{k_2n}{N}}|^2=\\ &=\frac{1}{N^3}\sum_{n=0}^{N-1}(|\hat{y}_{k_1}|^2+\hat{y}_{k_1}\hat{y}_{k_2}^*e^{2\pi i\frac{n(k_1-k_2)}{N}}+\hat{y}_{k_1}^*\hat{y}_{k_2}e^{-2\pi i\frac{n(k_1-k_2)}{N}}+|\hat{y}_{k_2}|^2)=\\ &=\frac{1}{N^3}\sum_{n=0}^{N-1}(|\hat{y}_{k_1}|^2+|\hat{y}_{k_2}|^2)=\\ &=\frac{|\hat{y}_{k_1}|^2+|\hat{y}_{k_2}|^2}{N^2}\end{aligned}$$ This extends to any $\textrm{len}(J)<N$. You have to write an algorithm that sorts the energies of the DFT components, and zeros them out in the inverse DFT in order for a given $\textrm{len}(J)$.