Inverse FFT in C#

Question

I am writing an application for procedural audiofiles, I have to analyze my new file, get its frequency spectrum and change it in its calculated.

I want to do this with the Fast Fourier Transform (FFT). This is my recursive C# FFT:

void ft(float n, ref Complex[] f)
{
    if (n > 1)
    {
        Complex[] g = new Complex[(int) n / 2];
        Complex[] u = new Complex[(int) n / 2];

        for (int i = 0; i < n / 2; i++)
        {
            g[i] = f[i * 2];
            u[i] = f[i * 2 + 1];
        }

        ft(n / 2, ref g);
        ft(n / 2, ref u);

        for (int i = 0; i < n / 2; i++)
        {
            float a = i;
            a = -2.0f * Mathf.PI * a / n;
            float cos = Mathf.Cos(a);
            float sin = Mathf.Sin(a);
            Complex c1 = new Complex(cos, sin);
            c1 = Complex.Multiply(u[i], c1);
            f[i] = Complex.Add(g[i], c1);

            f[i + (int) n / 2] = Complex.Subtract(g[i], c1);
        }
    }
}

The inspiring example was

I then compared my results with those from wolframalpha for the same input 0.6,0.7,0.8,0.9 but the results aren't be the same. My results are twice as big than Wolfram's and the imaginary part are the -2 times of Wolfram's.

Also, wiki indicates that the inverse of FFT can be computed with

But I compare inputs and outputs and they are different.

Has anyone an idea what's wrong?

Answer 1

Different implementations often use different definitions of the Discrete Fourier Transform (DFT), with correspondingly different results. The correspondence between implementations is usually fairly trivial (such as a scaling factor).

More specifically, your implementation is based on the following definition of the DFT:

On the other hand, Wolfram alpha by default uses a definition, which after adjusting to 0-based indexing looks like:

Correspondingly, it is possible to transform the result of your implementation to match Wolfram alpha's with:

void toWolframAlphaDefinition(ref Complex[] f)
{
  float scaling = (float)(1.0/Math.Sqrt(f.Length));
  for (int i = 0; i < f.Length; i++)
  {
    f[i] = scaling * Complex.Conjugate(f[i]);
  }
}

Now as far as computing the inverse DFT using the forward transform, a direct implementation of the formula

you provided would be:

void inverseFt(ref Complex[] f)
{
  for (int i = 0; i < f.Length; i++)
  {
    f[i] = Complex.Conjugate(f[i]);
  }
  ft(f.Length, ref f);
  float scaling = (float)(1.0 / f.Length);
  for (int i = 0; i < f.Length; i++)
  {
    f[i] = scaling * Complex.Conjugate(f[i]);
  }
}

Calling ft on the original sequence 0.6, 0.7, 0.8, 0.9 should thus get you the transformed sequence 3, -0.2+0.2j, -0.2, -0.2-0.2j.

Further calling inverseFt on this transform sequence should then bring you back to your original sequence 0.6, 0.7, 0.8, 0.9 (within some reasonable floating point error), as shown in this live demo.