iFFT from FFT

[pyroesp]

Hi guys,

I made an FFT a while back for a project on a TI ARM controller.
I've heard you can make an equalizer using an FFT & iFFT by transforming the samples, applying a "filter" at the output of the FFT and feeding the modified data back into an iFFT.

I've also found that you can make an iFFT, by using the FFT and 2 simple "tricks":

• iFFT(X) = FFT(X*)*/N
• iFFT(X) = swap(FFT(swap(X)))/N

http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Expressing_the_inverse_DFT_in_terms_of_the_DFT

I've tried both tricks and none worked.
I'm attaching my code::blocks project (zip) so you guys can test it.

I'm using the swap trick in the attached zip file by doing :

Code: [Select]

        // swap(x) -> if  x = a - jb then swap(x) = b - ja
        h = data_complex[i].im;
        data_complex[i].im = data_complex[i].re * (h < 0 ? -1 : 1);
        data_complex[i].re = h * (h < 0 ? -1 : 1);
I'm currently using an 8-point FFT with a pre-built sine wave (1kHz, amp 1 + 2kHz, amp 0.5, phase -45°).
The output of my FFT is correct, but not when I'm trying to inverse it and I have no idea why it's not working.

[ejeffrey]

Your swap looks messed up. I can see from the comment that you are trying to flip the signs around in a weird way.  You shouldn't have to do that.  Even if you do, your code for doing it is wrong: for instance, the resulting real component will always be positive.

All your swap should do is:

Code: [Select]

tmp = x.re
x.re = x.im
x.im = tmp

If you want to implement the code in your comments, then you want:

Code: [Select]

tmp = x.re
x.re = -x.im
x.im = -tmp


[pyroesp]

Hmm, I thought if I had "a - jb" it should have been "b - ja" but no, it should be "-b + ja".
My bad, kind of sleepy at 5 am  .

Still doesn't work though, test it out if you can.
Btw, I did this:

Code: [Select]

tmp = x.re
x.re = x.im
x.im = tmp

[ejeffrey]

It also looks like your FFT isn't working properly either.  Here is the output of your FFT program for the input data and forward transform:

Code: [Select]

input:
        X(0) = 0.353553 + 0.000000j
        X(1) = 0.353553 + 0.000000j
        X(2) = 0.646447 + 0.000000j
        X(3) = -1.353553 + 0.000000j
        X(4) = 0.353553 + 0.000000j
        X(5) = -1.060660 + 0.000000j
        X(6) = 1.060660 + 0.000000j
        X(7) = -0.353553 + 0.000000j

output:
        X(0) = -0.000000 + 0.000000j
        X(1) = -0.000000 + 4.000000j
        X(2) = 1.414214 + -1.414213j
        X(3) = 0.000000 + 0.000000j
        X(4) = -0.000000 + 0.000000j
        X(5) = 0.000000 + 0.000000j
        X(6) = 1.414214 + 1.414213j
        X(7) = 0.000000 + -4.000000j

However, when I put that input data into matlab, I get:

Code: [Select]

fft(1j*x).' =

        0 - 0.0000i
  -0.1213 + 1.7071i
   1.0000 - 1.0000i
   0.7071 - 1.7071i
        0 + 4.8284i
  -0.7071 - 1.7071i
  -1.0000 - 1.0000i
   0.1213 + 1.7071i


[pyroesp]

The input samples are 8 samples from a sine wave at 1kHz (ampl = 1) and another one at 2kHz (ampl = 0.5), so the output of the FFT (with sampling frequency of 8kHz) should be a 1 amplitude sine at 1kHz and 0.5 amplitude sine at 2kHz.
And that's what you see in my FFT.

X(0) = 0 => No DC component (0Hz)
X(1) = 0 + j4 = amplitude 4, phase 0 @ 1kHz
X(2) = 1.41 - j1.41 = amplitude 2, phase -45° @ 2kHz
X(3) = 0 => No 3kHz component

You still have to divide the amplitude by N/2, but you'll see it's correct.

I have no idea why your matlab is making those outputs, but that's not right.

[ejeffrey]

Ah.  I was just copying the data that you label as "input", but that label as wrong.  That is not the input data, it is the data after the bit reversed indexing at the beginning of the FFT.  In that case I don't really know what the inputs or outputs mean.  The output of your program is still disagrees with matlab, but only by phases, and I am not really sure whether you are printing the true FFT value or an internal representation.

[pyroesp]

My bad , it's true that I re-arrange the input samples when calling the dataToComplex function.
The "original" input samples should be in the x buffer.

I changed the printfs in the code, everything should be clear now, see attached zip.
- input samples
- complex & bit reversed input
- FFT output (-> complex & amplitude + phase)
- iFFT output (-> complex & amplitude + phase)

PS: The phases are kind of messed up, not really sure why though. For X(1) and X(2) they're correct, but it's calculating a phase of -90° at X(0)...

[ejeffrey]

You are skipping the bit reversal step in the second fft.

[pyroesp]

I thought you didn't need to do that...

I'm going to try that.

EDIT:
It works ! Thanks for the help ejeffrey.

Program output:

input samples:
        x(0) = 0.353553
        x(1) = 0.353553
        x(2) = 0.646447
        x(3) = 1.060660
        x(4) = 0.353553
        x(5) = -1.060660
        x(6) = -1.353553
        x(7) = -0.353553

complex & bit reversed input:
        x(0) = 0.353553 + 0.000000 j
        x(1) = 0.353553 + 0.000000 j
        x(2) = 0.646447 + 0.000000 j
        x(3) = -1.353553 + 0.000000 j
        x(4) = 0.353553 + 0.000000 j
        x(5) = -1.060660 + 0.000000 j
        x(6) = 1.060660 + 0.000000 j
        x(7) = -0.353553 + 0.000000 j

FFT output:
        X(0) = -0.000000 + 0.000000 j = 0.000000 amplitude, 0.000000 phase
        X(1) = -0.000000 + 4.000000 j = 1.000000 amplitude, -0.000002 phase
        X(2) = 1.414214 + -1.414213 j = 0.500000 amplitude, -45.000008 phase
        X(3) = 0.000000 + 0.000000 j = 0.000000 amplitude, 0.000000 phase
        X(4) = -0.000000 + 0.000000 j = 0.000000 amplitude, 0.000000 phase
        X(5) = -0.000000 + 0.000000 j = 0.000000 amplitude, 0.000000 phase
        X(6) = 1.414214 + 1.414213 j = 0.500000 amplitude, 45.000004 phase
        X(7) = 0.000000 + -4.000000 j = 1.000000 amplitude, -0.000002 phase

bit reversed FFT output:
        X(0) = -0.000000 + 0.000000 j
        X(1) = -0.000000 + 0.000000 j
        X(2) = 1.414214 + -1.414213 j
        X(3) = 1.414214 + 1.414213 j
        X(4) = -0.000000 + 4.000000 j
        X(5) = -0.000000 + 0.000000 j
        X(6) = 0.000000 + 0.000000 j
        X(7) = 0.000000 + -4.000000 j

swapped bit reversed FFT output:
        X(0) = 0.000000 + -0.000000 j
        X(1) = 0.000000 + -0.000000 j
        X(2) = -1.414213 + 1.414214 j
        X(3) = 1.414213 + 1.414214 j
        X(4) = 4.000000 + -0.000000 j
        X(5) = 0.000000 + -0.000000 j
        X(6) = 0.000000 + 0.000000 j
        X(7) = -4.000000 + 0.000000 j

swapped & normalized iFFT output:
        X(0) = 0.353553 + 0.000000 j
        X(1) = 0.353553 + -0.000000 j
        X(2) = 0.646447 + 0.000000 j
        X(3) = 1.060660 + 0.000000 j
        X(4) = 0.353553 + 0.000000 j
        X(5) = -1.060660 + 0.000000 j
        X(6) = -1.353553 + 0.000000 j
        X(7) = -0.353553 + -0.000000 j

I'm still having issues with the phase, but it's not like I'm going to use it anyways ...

EDIT 2:

I'm not checking if my imaginary part is 0 when calculating the phase, so I could have a divide by 0 error.
I just added an if statement to the part of the code that calculates the phase:

Code: [Select]

if (imaginary <= -0.000001 || imaginary >= 0.000001) // pseudo code
    // calculate phase
else
    // phase = 0
I edited the quote. No more phase errors .

[ejeffrey]

You can sidestep divide by zero and also correctly get the full -pi to +pi phase by using atan2(x.im, x.re).

[pyroesp]

Oh, I didn't know about that. I'm going to try that.
Thanks.

EDIT:

Yeah, I'm still getting a 180° phase for a complex value of  "-0.0000 + 0.0000j" (FFT output, X(0)) and other stuff like that... :/

[ejeffrey]

The phase of 0+0j is undefined, so you can't expect it to make sense there.  The C library has conventions for positive and negative zeros that respect the limiting behavior as you move towards zero, which avoids NaNs but still isn't physically meaningful.  However, atan2 should handle correctly if either the real or imaginary part is zero but not both.  More importantly, it will handle properly the case where the real part or both parts are negative, which the normal arctan will not.

[pyroesp]

It's seems to work pretty well with atanf(). It's not perfect, but I don't see any use in calculating the phase so it doesn't really bother me that it isn't perfect.
As long as the amplitudes are correct, I'm happy with it ^^.