Using fourier coefficients to reconstruct data in matlab

Question

When doing a discrete fourier transform on some data using matlab's fft function, its output is a set of fourier coefficients but I was wondering how do I go about converting these into an and bn so I can reconstruct the signal using sines and cosines.

E.g. I'd have a for loop that continually adds up a_i cos(ix) + b_isin(ix), where i = 1:N

Answer 1

For real signals, it works a little differently than the complex case as there is a symmetry in the DFT that can be exploited.

The DFT and inverse DFT are conventionally defined as:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i n k \frac{2\pi}{N} } $$

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i n k \frac{2\pi}{N} } $$

Now express your bin value in cartesian coordinates (due to the IMO poor convention of using x for the signal and X for the DFT bins, this is a different x):

$$ X[k] = x_k + i y_k $$

$$ X[N-k] = x_k - i y_k $$

Note:

$$ e^{i n (N - k) \frac{2\pi}{N} } = e^{-i n k \frac{2\pi}{N} } $$

Evaluate the kth $(k>0)$ term:

$$ \begin{aligned} T_k &= X[k] e^{i n k \frac{2\pi}{N} } + X[N-k] e^{-i n k \frac{2\pi}{N} } \\ &= ( x_k + i y_k )( \cos( n k \frac{2\pi}{N} ) + i \sin( n k \frac{2\pi}{N} ) ) \\ &+ ( x_k - i y_k )( \cos( n k \frac{2\pi}{N} ) - i \sin( n k \frac{2\pi}{N} ) ) \\ &= 2 x_k \cos( n k \frac{2\pi}{N} ) - 2 y_k \sin( n k \frac{2\pi}{N} ) ) \end{aligned} $$

Put that together, with a shift of one for MATLAB indexing fiasco, and you got what you want.

For the continuous reconstruction, then it is typical to rescale the domain from $ 0 \to N $ to $ 0 \to 2\pi $.

$$ t = \frac{n}{N} 2\pi $$

$$ n = t \frac{N}{2\pi} $$

The $k$th term then becomes:

$$ T_k = 2 x_k \cos( k t ) - 2 y_k \sin( k t ) $$

So for a real valued signal, the DFT real and imaginary are half the Fourier Series coefficients. But you still have to divide them by N if you used the conventional unnormalized DFT.

Now, when we get to the Nyquist frequency when N is even, for real valued signals, to keep the results real:

$$ T_{N/2} = x_{N/2} \cos( \frac{N}{2} t ) $$

For odd N, the last $k$ will be (N-1)/2 and nothing special needs to be done.

Now when you add them up, you still need to divide by N, so in conclusion:

$$ a_k = \frac{2}{N} x_k $$

$$ b_k = -\frac{2}{N} y_k $$

Finally the zeroth bin, $k=0$, aka the DC bin, MATLAB index of one (stupid), for a real signal is the offset you start your series with.

$$ a_0 = \frac{1}{N} x_0 $$ $$ b_0 = 0 $$

For N even:

$$ a_{N/2} = \frac{1}{N} x_{N/2} $$ $$ b_{N/2} = 0 $$

I know this isn't the clearest or cleanest explanation, but you can figure it out.

Ced