CODEWITHSUNDEEP
c#fft
pasq
pasq

Reputation: 3

Calculate values of the spectrum with FFT

I have to calculate the spectrum values of an audio. I used aForge's FFT in Sources/Math/FourierTransform.cs and I used an example of sampling with 16 samples as used in this video to check the results with excel (I tested the results in a spreadsheet like in the video).

FFT:

public enum Direction
{
    Forward = 1,
    Backward = -1
};

private const int minLength = 2;
private const int maxLength = 16384;
private const int minBits = 1;
private const int maxBits = 14;
private static int[][] reversedBits = new int[maxBits][];
private static Complex[,][] complexRotation = new Complex[maxBits, 2][];

static void Main(string[] args)
{
    var Data = new Complex[16];
    Data[0] = new Complex(0, 0);
    Data[1] = new Complex((float)0.998027, 0);
    Data[2] = new Complex((float)0.125333, 0);
    Data[3] = new Complex((float)-0.98229, 0);
    Data[4] = new Complex((float)-0.24869, 0);
    Data[5] = new Complex((float)0.951057, 0);
    Data[6] = new Complex((float)0.368125, 0);
    Data[7] = new Complex((float)-0.90483, 0);
    Data[8] = new Complex((float)-0.48175, 0);
    Data[9] = new Complex((float)0.844328, 0);
    Data[10] = new Complex((float)0.587785, 0);
    Data[11] = new Complex((float)-0.77051, 0);
    Data[12] = new Complex((float)-0.68455, 0);
    Data[13] = new Complex((float)0.684547, 0);
    Data[14] = new Complex((float)0.770513, 0);
    Data[15] = new Complex((float)-0.58779, 0);

    FFT(Data, Direction.Forward);

    for (int a = 0; a <= Data.Length - 1; a++)
    {
        Console.WriteLine(Data[a].Re.ToString());
    }

    Console.ReadLine();
}

public static void FFT(Complex[] data, Direction direction)
{
    int n = data.Length;
    int m = Tools.Log2(n);

    // reorder data first
    ReorderData(data);

    // compute FFT
    int tn = 1, tm;

    for (int k = 1; k <= m; k++)
    {
        Complex[] rotation = GetComplexRotation(k, direction);

        tm = tn;
        tn <<= 1;

        for (int i = 0; i < tm; i++)
        {
            Complex t = rotation[i];

            for (int even = i; even < n; even += tn)
            {
                int odd = even + tm;
                Complex ce = data[even];
                Complex co = data[odd];

                double tr = co.Re * t.Re - co.Im * t.Im;
                double ti = co.Re * t.Im + co.Im * t.Re;

                data[even].Re += tr;
                data[even].Im += ti;

                data[odd].Re = ce.Re - tr;
                data[odd].Im = ce.Im - ti;
            }
        }
    }

    if (direction == Direction.Forward)
    {
        for (int i = 0; i < n; i++)
        {
            data[i].Re /= (double)n;
            data[i].Im /= (double)n;
        }
    }
}

private static int[] GetReversedBits(int numberOfBits)
{
    if ((numberOfBits < minBits) || (numberOfBits > maxBits))
        throw new ArgumentOutOfRangeException();

    // check if the array is already calculated
    if (reversedBits[numberOfBits - 1] == null)
    {
        int n = Tools.Pow2(numberOfBits);
        int[] rBits = new int[n];

        // calculate the array
        for (int i = 0; i < n; i++)
        {
            int oldBits = i;
            int newBits = 0;

            for (int j = 0; j < numberOfBits; j++)
            {
                newBits = (newBits << 1) | (oldBits & 1);
                oldBits = (oldBits >> 1);
            }
            rBits[i] = newBits;
        }
        reversedBits[numberOfBits - 1] = rBits;
    }
    return reversedBits[numberOfBits - 1];
}

private static Complex[] GetComplexRotation(int numberOfBits, Direction direction)
{
    int directionIndex = (direction == Direction.Forward) ? 0 : 1;

    // check if the array is already calculated
    if (complexRotation[numberOfBits - 1, directionIndex] == null)
    {
        int n = 1 << (numberOfBits - 1);
        double uR = 1.0;
        double uI = 0.0;
        double angle = System.Math.PI / n * (int)direction;
        double wR = System.Math.Cos(angle);
        double wI = System.Math.Sin(angle);
        double t;
        Complex[] rotation = new Complex[n];

        for (int i = 0; i < n; i++)
        {
            rotation[i] = new Complex(uR, uI);
            t = uR * wI + uI * wR;
            uR = uR * wR - uI * wI;
            uI = t;
        }

        complexRotation[numberOfBits - 1, directionIndex] = rotation;
    }
    return complexRotation[numberOfBits - 1, directionIndex];
}

// Reorder data for FFT using
private static void ReorderData(Complex[] data)
{
    int len = data.Length;

    // check data length
    if ((len < minLength) || (len > maxLength) || (!Tools.IsPowerOf2(len)))
            throw new ArgumentException("Incorrect data length.");

    int[] rBits = GetReversedBits(Tools.Log2(len));

    for (int i = 0; i < len; i++)
    {
        int s = rBits[i];

        if (s > i)
        {
            Complex t = data[i];
            data[i] = data[s];
            data[s] = t;
        }
    }
}

These are the results after the transformation:

Output FFT results:             Excel FFT results:
0,0418315622955561              0,669305
0,0533257974328085              0,716163407
0,137615673627316               0,908647001         
0,114642731070279               1,673453043
0,234673940537634               7,474988602
0,0811255020953362              0,880988382          
0,138088891589122               0,406276784
0,0623766891658306              0,248854492
0,0272978749126196              0,204227
0,0124250144575261              0,248854492
0,053787064184711               0,406276784
0,00783331226557493             0,880988382
0,0884368745610118              7,474988602
0,0155431246384978              1,673453043
0,0301093757152557              0,908647001
0                               0,716163407

The results are not at all similar. Where is it wrong? Is the implementation of complex (Data) wrong or is the FFT method wrong or other?

Thanks in advance!

Upvotes: 0

Views: 245

Answers (1)

Jeff
Jeff

Reputation: 7674

First, the resulting FFT is a complex function in general. You're only displaying the real parts in your code but the thing you're comparing to is displaying the magnitudes, so of course they're going to be different: you're comparing apples to oranges.

When you use magnitudes and compare apples to apples, you should get this:

for (int a = 0; a <= Data.Length - 1; a++)
{
    Console.WriteLine(Data[a].Magnitude.ToString());
}

...

0.0418315622955561
0.0447602132472683
0.0567904388057513
0.104590813761862
0.46718679147454
0.0550617784710375
0.025392294285886
0.0155534081359397
0.0127641875296831
0.0155534081359397
0.025392294285886
0.0550617784710375
0.46718679147454
0.104590813761862
0.0567904388057513
0.0447602132472683

That looks a little better -- it has the same symmetry property as the Excel output and there appear to be peaks in the same locations.

It almost looks like the scale is off. If I divide each element by the corresponding element from the Excel output, I get:

16
16
16
16
16
16
16
16
16
16
16
16
16
16
16
16

So your results are pretty much correct, just off by a scaling factor.

You're dividing everything by n in the last step of your FFT:

if (direction == Direction.Forward)
{
    for (int i = 0; i < n; i++)
    {
        data[i].Re /= (double)n;
        data[i].Im /= (double)n;
    }
}

This is conventionally done for the inverse transform, not the forward transform.

In summary, changing the output from Data[a].Re to Data[a].Magnitude and changing the condition at the end of FFT from if (direction == Direction.Forward) to if (direction == Direction.Backward), I get this output:

0.669304996728897
0.716163411956293
0.908647020892022
1.67345302018979
7.47498866359264
0.880988455536601
0.406276708574176
0.248854530175035
0.20422700047493
0.248854530175035
0.406276708574176
0.880988455536601
7.47498866359264
1.67345302018979
0.908647020892022
0.716163411956293

which matches the Excel output.

Upvotes: 1

Related Questions