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 ai cos(ix) + bisin(ix), where i = 1:N
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
