0

Hi guy's I'm working on simple signals and I want to calculate the Fourier transform of a signal, get the magnitude and phases, then reconstruct the original signal from that.

I'm basing my code on this thread.

Code:

>> n=0:99;
>> N=length(n);
>> x = sin((2*pi/N).*n).*cos((pi/N).*n);
>> F = fft(x);
>> mag =  sqrt(real(F).^2 + imag(F).^2);
>> phase = atan2(imag(F),real(F));
>> re = mag .* cos(phase);
>> im = mag .* sin(phase);
>> F_i = re + 1i*im;
>> x_i = ifft(F_i);
>> figure;stem(x);figure;stem(x_i);

I completely get different graphs.

EDIT: I'm actually doing this to test what would happen to a signal if the phase changes. So with this I will need the phase angle to construct the signal again.

I'm still new to both Fourier + Matlab so I apologize if I'm making some random stupid mistake. I'll appreciate it if you guys can point me in the right direction. Thank you.

Community
  • 1
  • 1
kir
  • 581
  • 1
  • 6
  • 22
  • The imaginary part of `F_i` is mainly roundoff error, if you use `F_i=re` the figure is reproduced correctly. I don't know enough about Fourier transforms to tell you why this is so prone to numerical instability though. – David Nov 13 '13 at 00:15
  • This is true. I'm actually doing this for an experiment on where I change the phase angle and reconstruct the signal and see what's the change. Is there any other way I can do this then? – kir Nov 13 '13 at 00:22

2 Answers2

2

The problem is caused by round-off errors in phase, so don't use them when calulating the sine and consine of the phase angles. Instead, use trig identities cos(atan(A))=(1+A^2)^(-1/2), and sin(atan(A))=A*(1+A^2)^(-1/2), and so

re = mag .* real(F)./sqrt(real(F).^2+imag(F).^2);
im = mag .* imag(F)./sqrt(real(F).^2+imag(F).^2);

EDIT: I think that if you want to change the phase angle by S, this will do the trick:

re = mag .* (real(F)*cos(S)-imag(F)*sin(S))./sqrt(real(F).^2+imag(F).^2);
im = mag .* (real(F)*sin(S)+imag(F)*cos(S))./sqrt(real(F).^2+imag(F).^2);

EDIT2: You will still sometimes get bad results with non-zero imaginary part (e.g. if S=pi), and you will need to plot either stem(real(x_i)) or stem(1:length(x_i),x_i) as Luis suggested.

David
  • 8,449
  • 1
  • 22
  • 32
  • Thank you this properly reproduces the signal. If I want to modify the im to lets say adding a constant "10" do I need to do any calculations or would `im = im + 10` suffice? – kir Nov 13 '13 at 00:36
  • Yeah if you want to change `im` just do it. I have edited my answer to include what I think is the formulae for `re` and `im` with a phase shift of `S` radians. – David Nov 13 '13 at 00:54
  • @David Suppose if I want to change the entire phase-matrix to something say phase-new and then want to apply ifft() using mag and phase-new, then how should I proceed? – Ayush Nov 16 '14 at 16:19
  • 1
    @Jerky Might be best to ask a new question, but I don't think my method will work, as it's just a constant phase change. Just replace your calculated `phase` with a different matrix and it ought to work. – David Nov 16 '14 at 21:40
2

There is no such numerical instability in the FFT. The problem is that roundoff errors give a very small imaginary part in the recovered signal x_i. That's normal. The real part of x_i correctly reproduces x, whereas the imaginary part of x_i is very small. You can check with stem(real(x_i)) and stem(imag(x_i))

Now stem(x_i) with complex x_i plots the imaginary part versus the real part. On the othe hand, stem(1:length(x_i),x_i) (or of course stem(real(x_i))) plots the real part of x_i. Note that this is consistent with plot's behaviour (see also the answer to this question)

So just plot the real imaginary parts of x_i separately with stem. You'll see that the real part reproduces x as expected, and the imaginary part is negligible (of the order of eps).

Community
  • 1
  • 1
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • Yes this is true. I wonder why my solution fixes the problem? – David Nov 13 '13 at 00:41
  • @David Somehow it doesn't cause roundoff errors, hence no imaginary part in `stem` and no problem in plotting. The weird part is why `stem` does that with a complex argument. It would make more sense if it ignored the imaginary part, as `plot(x,y)` does – Luis Mendo Nov 13 '13 at 00:42
  • @David In fact, `stem(x,y)` _does_ behave like `plot(x,y)`: it ignores imaginary part of `y`. – Luis Mendo Nov 13 '13 at 00:46