DSPRelated.com
Forums

IFFT is not working

Started by Il Totem April 13, 2009
I am struggling to make IFFT work, but I have had a few problems. I have
read over the Internet that it is relatively simple to compute an IFFT of a
frequency domain: conjugate the frequency domain's samples, than compute a
direct FFT on them, and therefore take the conjugate of the output.
I have tried with a 4-points FFT with a simple sinusoid time domain signal
(only real parts set: 1, 0, -1, 0) and this method worked well. I therefore
tried an 8-points FFT, where the real time domain is:
ReX(i) = Sin(2Pi * i / 8)  i = 0..7
But the IFFT produces a completely wrong output. Here it is the code:

Dim x(7) As FourierTransforms.ComplexSample
For i As Int16 = 0 To 7
  x(i).RealPart = Math.Sin(2 * Math.PI * i / 8)
Next

Dim c() As FourierTransforms.ComplexSample =
FourierTransforms.ComplexFFTComplexOutput(x, False)
Dim d(c.Length - 1) As FourierTransforms.ComplexSample

For i As Int16 = 0 To c.Length - 1
  d(i).RealPart = c(i).ImaginaryPart
  d(i).ImaginaryPart = c(i).RealPart
Next

Dim k() As FourierTransforms.ComplexSample =
FourierTransforms.ComplexFFTComplexOutput(d, False)
Dim x2(k.Length - 1) As Double

For i As Int16 = 0 To k.Length - 1
  x2(i) = k(i).ImaginaryPart / k.Length
Next

Moreover, there is a thing that I haven't understood yet. On the pages I
have read there is written that conjugation consists in swapping real and
imaginary part of a complex sample. But the complex conjugate of a+ib is
a-ib and not b+ia.


On Apr 13, 9:47&#4294967295;am, "Il Totem" <nicolo1...@yahoo.it> wrote:
> I am struggling to make IFFT work, but I have had a few problems. I have > read over the Internet that it is relatively simple to compute an IFFT of a > frequency domain: conjugate the frequency domain's samples, than compute a > direct FFT on them, and therefore take the conjugate of the output. > I have tried with a 4-points FFT with a simple sinusoid time domain signal > (onlyrealparts set: 1, 0, -1, 0) and this method worked well. I therefore > tried an 8-points FFT, where therealtime domain is: > ReX(i) = Sin(2Pi * i / 8) &#4294967295;i = 0..7 > But the IFFT produces a completely wrong output. Here it is the code: > > Dim x(7) As FourierTransforms.ComplexSample > For i As Int16 = 0 To 7 > &#4294967295; x(i).RealPart = Math.Sin(2 * Math.PI * i / 8) > Next > > Dim c() As FourierTransforms.ComplexSample = > FourierTransforms.ComplexFFTComplexOutput(x, False) > Dim d(c.Length - 1) As FourierTransforms.ComplexSample > > For i As Int16 = 0 To c.Length - 1 > &#4294967295; d(i).RealPart = c(i).ImaginaryPart > &#4294967295; d(i).ImaginaryPart = c(i).RealPart > Next > > Dim k() As FourierTransforms.ComplexSample = > FourierTransforms.ComplexFFTComplexOutput(d, False) > Dim x2(k.Length - 1) As Double > > For i As Int16 = 0 To k.Length - 1 > &#4294967295; x2(i) = k(i).ImaginaryPart / k.Length > Next > > Moreover, there is a thing that I haven't understood yet. On the pages I > have read there is written that conjugation consists inswappingrealandimaginarypart of a complex sample. But the complex conjugate of a+ib is > a-ib and not b+ia.
you can do it either way use, use a conjugate or swap the real and imaginary, or to save computations for a real input you can do neither and get the correct result but backwards!
> >you can do it either way use, use a conjugate or swap the real and >imaginary, or to save computations for a real input you can do >neither and get the correct result but backwards! >
Indeed, nothing seems to work. I have tried every method explained on the Internet, including reverting input to receive the original output, but the results are constantly wrong. In the resulting array there are some correct values, but they are less than an half and are located in apparently random places in the array.
On Apr 14, 5:03&#4294967295;am, "Il Totem" <nicolo1...@yahoo.it> wrote:
> >you can do it either way use, use a conjugate or swap the real and > >imaginary, or &#4294967295;to save computations for a real input you can do > >neither and get the correct result but backwards! > > Indeed, nothing seems to work. I have tried every method explained on the > Internet, including reverting input to receive the original output, but the > results are constantly wrong. In the resulting array there are some correct > values, but they are less than an half and are located in apparently random > places in the array.
Works find in matlab... for an input of 1 2 3 4 5 6 7 8 the conj method
>> fft(conj (fft ([1 2 3 4 5 6 7 8])))/8
ans = 1 2 3 4 5 6 7 8
>>
and the swap real and imag method >> fft ([1 2 3 4 5 6 7 8]) ans = Columns 1 through 5 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 Columns 6 through 8 -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
>> imag(fft(imag(ans) + j*real(ans))/8)
ans = 1 2 3 4 5 6 7 8
>>
>Works find in matlab... > >for an input of 1 2 3 4 5 6 7 8 > >the conj method > >>> fft(conj (fft ([1 2 3 4 5 6 7 8])))/8 > >ans =3D > > 1 2 3 4 5 6 7 8 > >>> > >and the swap real and imag method > > >> fft ([1 2 3 4 5 6 7 8]) > >ans =3D > > Columns 1 through 5 > > 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + >1.6569i -4.0000 > > Columns 6 through 8 > > -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i > >>> imag(fft(imag(ans) + j*real(ans))/8) > >ans =3D > > 1 2 3 4 5 6 7 8 > >>> >
I am sure that with matlabs it works well, but this does not help me. I will show you and exmplae of what I receive as output. Suppose that the time domain is composed by real-only sample, with these values: 0, 7.07, 10, 7.07, 0, -7.07, -10, -7.07 I take the fft of those, I swap the real and imaginary part of freq domain's samples, then I compute a direct fft of the latter and therefore take the imaginary part of the output. It gives me the following values: 0, 9.57, 7.5, -3.96, -3.54, -1.04, -7.5, -1.04
>Works find in matlab... > >for an input of 1 2 3 4 5 6 7 8 > >the conj method > >>> fft(conj (fft ([1 2 3 4 5 6 7 8])))/8 > >ans =3D > > 1 2 3 4 5 6 7 8 > >>> > >and the swap real and imag method > > >> fft ([1 2 3 4 5 6 7 8]) > >ans =3D > > Columns 1 through 5 > > 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + >1.6569i -4.0000 > > Columns 6 through 8 > > -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i > >>> imag(fft(imag(ans) + j*real(ans))/8) > >ans =3D > > 1 2 3 4 5 6 7 8 > >>> >
I am sure that with matlabs it works well, but this does not help me. I will show you and exmplae of what I receive as output. Suppose that the time domain is composed by real-only sample, with these values: 0, 7.07, 10, 7.07, 0, -7.07, -10, -7.07 I take the fft of those, I swap the real and imaginary part of freq domain's samples, then I compute a direct fft of the latter and therefore take the imaginary part of the output. It gives me the following values: 0, 9.57, 7.5, -3.96, -3.54, -1.04, -7.5, -1.04
On Apr 14, 1:16&#4294967295;pm, "Il Totem" <nicolo1...@yahoo.it> wrote:
> >Works find in matlab... > > >for an input of 1 2 3 4 5 6 7 8 > > >the conj method > > >>> fft(conj (fft ([1 2 3 4 5 6 7 8])))/8 > > >ans =3D > > > &#4294967295; &#4294967295; 1 &#4294967295; &#4294967295; 2 &#4294967295; &#4294967295; 3 &#4294967295; &#4294967295; 4 &#4294967295; &#4294967295; 5 &#4294967295; &#4294967295; 6 &#4294967295; &#4294967295; 7 &#4294967295; &#4294967295; 8 > > >and the swap real and imag method > > > >> fft ([1 2 3 4 5 6 7 8]) > > >ans =3D > > > &#4294967295;Columns 1 through 5 > > > &#4294967295;36.0000 &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;-4.0000 + 9.6569i &#4294967295;-4.0000 + 4.0000i &#4294967295;-4.0000 + > >1.6569i &#4294967295;-4.0000 > > > &#4294967295;Columns 6 through 8 > > > &#4294967295;-4.0000 - 1.6569i &#4294967295;-4.0000 - 4.0000i &#4294967295;-4.0000 - 9.6569i > > >>> imag(fft(imag(ans) + j*real(ans))/8) > > >ans =3D > > > &#4294967295; &#4294967295; 1 &#4294967295; &#4294967295; 2 &#4294967295; &#4294967295; 3 &#4294967295; &#4294967295; 4 &#4294967295; &#4294967295; 5 &#4294967295; &#4294967295; 6 &#4294967295; &#4294967295; 7 &#4294967295; &#4294967295; 8 > > I am sure that with matlabs it works well, but this does not help me. I > will show you and exmplae of what I receive as output. Suppose that the > time domain is composed by real-only sample, with these values: > &#4294967295;0, 7.07, 10, 7.07, 0, -7.07, -10, -7.07 > I take the fft of those, I swap the real and imaginary part of freq > domain's samples, then I compute a direct fft of the latter and therefore > take the imaginary part of the output. It gives me the following values: > &#4294967295;0, 9.57, 7.5, -3.96, -3.54, -1.04, -7.5, -1.04- Hide quoted text - > > - Show quoted text -
something must be wrong with your FFT, what do you get for your very first FFT? Should be 0,-40i,0,0,0,0,0,40i
I just did the following:

Set real inputs r[0 to 7] to: 0., 7.07, 10., 7.07, 0., -7.07, -10.,
-7.07, and imaginary inputs i[0 to 7] to zero. Then I did an 8 point
forward FFT.  Then I divided all outputs by N (= 8 in this case).
Then I inversed by exchanging the r[ ] and i[ ] in the forward call:

fft( N , r, i )  // forward fft

Divide results by N ( = 8 in this case)

fft( N, i, r )  //inverse fft via a forward fft

The output result in the array r[ ] is the same as what I started
with:

0., 7.07, 10., 7.07, 0., -7.07, -10., -7.07

Note that there's no real swapping of anything going on.  It's all
done due the way you give the array names to the function that
matters.  After the inverse, the r[ ] array is your real output and
the i[ ] array is your imaginary output.

This technique has been around for at least 20+ years (perhaps much
longer), and it works. If it doesn&#4294967295;t work for you, then check your
code.  And, as noted in steve&#4294967295;s previous post, perhaps your fft is not
working correctly. Note that to do the above, you have to have an fft
routine that accepts and works with both real and imaginary inputs/
outputs.  Some routines for real inputs won&#4294967295;t do that.

Kevin



> >something must be wrong with your FFT, what do you get for your very >first FFT? Should be 0,-40i,0,0,0,0,0,40i >
The results of the first FFT are: 0, -40i, 0, 0, 0, 0, -14.14+34.13i, 14.14+5.87i The FFT code is: Public Shared Function ComplexFFTComplexOutput(ByRef TimeDomain() As ComplexSample, Optional ByVal HalfResult As Boolean = True) As ComplexSample() n = TimeDomain.Length nu = CInt((Math.Log(n) / Math.Log(2))) Dim n2 As Integer = n / 2 Dim nu1 As Integer = nu - 1 Dim xre(n - 1) As Double Dim xim(n - 1) As Double Dim tr, ti, p, arg, c, s As Double For i As Integer = 0 To n - 1 xre(i) = TimeDomain(i).RealPart xim(i) = TimeDomain(i).ImaginaryPart Next Dim k As Integer = 0 For l As Integer = 1 To nu While (k < n) For i As Integer = 1 To n2 p = BitReverse(k >> nu1) arg = 2 * CDbl(Math.PI) * p / n c = CDbl(Math.Cos(arg)) s = CDbl(Math.Sin(arg)) tr = xre(k + n2) * c + xim(k + n2) * s ti = xim(k + n2) * c - xre(k + n2) * s xre(k + n2) = xre(k) - tr xim(k + n2) = xim(k) - ti xre(k) += tr xim(k) += ti k += 1 Next k += n2 End While k = 0 nu1 -= 1 n2 = n2 / 2 Next k = 0 Dim r As Integer While (k < n) r = BitReverse(k) If r > k Then tr = xre(k) ti = xim(k) xre(k) = xre(r) xim(k) = xim(r) xre(r) = tr xim(r) = ti End If k += 1 End While Dim Result() As ComplexSample If HalfResult Then ReDim Result(n / 2 - 1) For i As Int64 = 0 To n / 2 - 1 Result(i).RealPart = xre(i) Result(i).ImaginaryPart = xim(i) Next Else ReDim Result(n - 1) For i As Int64 = 0 To n - 1 Result(i).RealPart = xre(i) Result(i).ImaginaryPart = xim(i) Next End If xre = Nothing xim = Nothing Return Result End Function HalfResult determines wheter to take n/2 samples from the result (if true) or the whole output (if false). Since the intermediate results seems not to be correct (quite disturbing), could anyone look at that code and tell me if there's something wrong, put that HalfResult is always false?
>HalfResult determines wheter to take n/2 samples >from the result (if true)or the whole output (if false).
>Since the intermediate results seems not to >be correct (quite disturbing)
Your intermediate result is definitely not correct if you want all N outputs. But it would seem to be correct for the &#4294967295;half results&#4294967295; case for outputs 0 to N/2. Perhaps the value of HalfResult isn't what you think it is within the function? Find out by printing it out while within the function. Your routine is definitely based on a negative exponent, so your result should be: Im: 0,-40i,0,0,0,0,0,40i Re: all zeros It looks like you're getting only half the results with some leftover intermediate data in the upper part. And note that if you added (-14.14+34.13i) to (14.14+5.87i), you'd get 40i. It's almost as if you&#4294967295;re missing the upper part of the last column of butterflies. Given that your forward results are not correct for all N outputs, it's not surprising that an inverse won't work. You need to get the outputs shown above before doing the inverse. Kevin