I'm attempting to reconstruct a signal from the DFT of the signal. I tried to do it by extracting the individual sinusoids and adding them up, but the answer I get is incorrect.

s = [3 1 4 1 5 9 2 6]#Some random signal(first few digits of pi in this case)
t = 0:(1/7):1;#Sampling at each of these equally spaced time intervals
f = fft(s);
a = abs(f);#Magnitude of each DFT component
an = arg(f);#Phase of each DFT component

recs = (a(1)/8)*cos(0*pi*t + arg(1)) + (a(2)/8)*cos(2*pi*t + arg(2)) + (a(3)/8)*cos(4*pi*t + arg(3)) + (a(4)/8)*cos(6*pi*t + arg(4)) + (a(1)/8)*cos(8*pi*t + arg(5));
#The above equation is derived(hopefully correctly) from the theory given in the question description

This is a program I wrote in GNU Octave(should work on Matlab too). And the basis for using that particular expression for recs is based on some info I got from http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf.

I don't understand if my method for getting the individual components was wrong or is it some basic step I'm missing here.

  • $\begingroup$ i haven't looked at anything in particular, but this appears to me, on the surface, about the difference between a phase vocoder and, what we sometimes call, "sinusoidal modeling". $\endgroup$ – robert bristow-johnson Mar 26 at 0:01

There are a couple of mistakes in your calculations:

  1. The t array is incorrectly specified. It should be 0:1/8:7/8 as the time positions at which the signal is assumed to be sampled for the default FFT calculation performed by fft() are 0, 1, ..., (N-1) where N = 8 in your case, the signal length. Since in your formula you are using t as the ratio n/N (where n is the array of sample times) you get the aforementioned expression to use for t.
  2. You wrote a(1) instead of a(5) at the last term of your recs expression.
  3. You are missing 2* multiplying each term except the first (k=0) and last (k=N/2) ones, where k is the index of the frequency axis of the Fourier transform F (recall that this 2* comes from leveraging the symmetry of the Fourier transform associated to a real signal).

Therefore, your reconstruction expression should read:

recs =  (1*a(1)/8)*cos(0*pi*t + an(1)) + ...
        (2*a(2)/8)*cos(2*pi*t + an(2)) + ...
        (2*a(3)/8)*cos(4*pi*t + an(3)) + ...
        (2*a(4)/8)*cos(6*pi*t + an(4)) + ...
        (1*a(5)/8)*cos(8*pi*t + an(5));

which gives the original signal s as verified by computing:

recs - s

ans =

   1.0e-14 *

         0   -0.1332         0   -0.0444   -0.2665    0.3553    0.1332         0
  • $\begingroup$ Thanks for the answer. Yes, those were my mistakes. I didn't know about when we've to use 2/N and 1/N till I understood the derivation. I did mess up the t array badly, so thanks for that. As for the angle, GNU Octave uses arg and not angle for calculating the angle. So you might want to edit that out of your answer. Otherwise great, thanks $\endgroup$ – Pradyoth Shandilya Mar 26 at 1:16
  • $\begingroup$ @PradyothShandilya Edited the answer as suggested. Thanks for pointing that out about Octave. Also, if you think this answered your question, it would be great if you could mark it as solved or upvote the answer. Thanks! $\endgroup$ – mastropi Mar 26 at 1:31
  • $\begingroup$ Of course. Was just waiting for the edit $\endgroup$ – Pradyoth Shandilya Mar 26 at 1:36
  • $\begingroup$ @PradyothShandilya Thanks! :-) $\endgroup$ – mastropi Mar 26 at 16:21

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.