0
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.