Click here to Skip to main content
15,866,398 members
Articles / Programming Languages / C#

Quick FFT

Rate me:
Please Sign up or sign in to vote.
5.00/5 (19 votes)
13 Jul 2013CPOL4 min read 70.1K   6.6K   46   16
An FFT algorithm that runs a bit faster than the standard implementation.

Basic Radix-2-FFT algorithm recursive

Derived from the main clause of Fourier:

 

Image 1

 

the Radix-2-FFT algorithm for N = 2^j samples is based on the formula:

 

Image 2

 

with M = N/2.

With the substitution u[m] = x[2m] and v[m] = x[2m + 1], we get:

 

Image 3

 

In this formula the two sub-DFT’s

 

Image 4

and

Image 5

are visible and that means the basic Fourier transformation can be split in two sub transformations and these two each can be split in two sub-sub transformations again and so on till we have only two samples left per transformation.

If we start with 8 samples for the first splitting this looks like that

Image 6

And if we continue to split the flow finally looks like that:

Image 7

On the first glance it looks quite complicate. The input data is mixed up and there are many arrows in all directions. But if we just look at the first graphic containing one step, we get a simple system. In the formula above we see [2m]-elements are even and [2m+1]-elements are odd elements. That’s why the odd elements are separated from the even ones here. If we look at the first formula and given that the output of one sub transformation overrides the input variables of the transformation, we have to perform the sub transformations and just have to calculate the complex sums after that (given that x[k] is the output of the sub transformation and X[k] the output of this one):

Image 8

Image 9

Image 10

Image 11

Image 12

Image 13

Image 14

Image 15

And with

Image 16 Image 17

Image 18 Image 19

We can replace

Image 20

Image 21

Image 22

Image 23

And that’s what the first graphic shows.

With a recursive approach this can be put into a quite small algorithm that is a bit easier to be understood than a not recursive one. An important detail to recognize here is possibly that we don’t need the so called bit-reversal here. We just have to separate the odd from the even samples in each recursion as can be seen in the graphic. Quite often that’s not mentioned in descriptions of the Radix-2 algorithm.

Of course the bit-reversal can be done here as well. It would have to be done in the beginning of the algorithm…There are many different ways to implement this transformation.

The recursive procedure would look like that:

C#
public void CalcSubFFT(TKomplex[] a, int n, int lo)
{
    int i, m;
    TKomplex w;
    TKomplex v;
    TKomplex h;
    if (n > 1)
    {
        Shuffle(a, n, lo);  // separate odd from even elements
        m = n / 2;
        CalcSubFFT(a, m, lo);  //recursive call for u(k)
        CalcSubFFT(a, m, lo + m); //recursive call for v(k)
        for (i = lo; i < lo + m; i++)
        {
            w.real = Math.Cos(2.0 * Math.PI * (double)(i - lo) / (double)(n));
            w.imag = Math.Sin(2.0 * Math.PI * (double)(i - lo) / (double)(n));
            h = kprod(a[i + m], w);
            v = a[i];
            a[i] = ksum(a[i], h);
            a[i + m] = kdiff(v, h);
        }
    }
}

With the separation of the odd and even samples:

C#
public void Shuffle(TKomplex[] a, int n, int lo)
{
    if (n > 2)
    {
        int i, m = n / 2;
        TKomplex[] b = new TKomplex[m];
        for (i = 0; i < m; i++)
        b[i] = a[i * 2 + lo + 1];
        for (i = 0; i < m; i++)
        a[i + lo] = a[i * 2 + lo];
        for (i = 0; i < m; i++)
        a[i + lo + m] = b[i];
    }
}

At the end we have to put that all in one function and divide every element by 2N and element 0 has to be divided by 4N (see main clause of Fourier).

C#
public void CalcFFT()
{
    int i;
    CalcSubFFT(y, N, 0);
    for (i = 0; i < N; i++)
    {
        y[i].imag = y[i].imag / (double)N * 2.0;
        y[i].real = y[i].real / (double)N * 2.0;
    }
    y[0].imag = y[0].imag / 2.0;
    y[0].real = y[0].real / 2.0;
}

Reducing calculation time

As I already found in my article about the quick DFT, we have to calculate many sine and cosine values in these algorithms too. If we want to do a FFT repeatedly in an application, it makes some sense to put all the sine and cosine values into al look up table and just get them from there. That speeds the calculation up quite well.

We create the look up table like

C#
for (i = 0; i < (N / 2); i++)
    we[i].real = Math.Cos(2.0 * Math.PI * (double)(i) / (double)(N));
    we[i].imag = Math.Sin(2.0 * Math.PI * (double)(i) / (double)(N));
}

And replace the lines

C#
w.real = Math.Cos(2.0 * Math.PI * (double)(i - lo) / (double)(n));
w.imag = Math.Sin(2.0 * Math.PI * (double)(i - lo) / (double)(n));

by

C#
w.real = we[(i - lo)*N/n].real;
w.imag = we[(i - lo)*N/n].imag;

This modification of the algorithm reduced the calculation time for 4096 samples from 5 ms to 3 ms on my computer.

At the DFT algorithm we had to perform a modulo division to get the correct index. Here this is not necessary as we anyway have just sine and cosine values for 0 to 2Pi.

Basic Radix-2-FFT algorithm not recursive

The disadvantage of the easy to be understood recursive approach is that not every environment allows recursive function calls and a bit higher time consumption. All together there are more operations to be done with these shuffles and recursions. Therefore I implemented a not recursive algorithm too. That looks a bit more complicate but it’s worthwhile.

Here we need the Bitreversal function:

C#
public void BitInvert(TKomplex[] a, int n)
{
    int i, j, mv = n/2;
    int k, rev = 0;
    TKomplex b;
    for (i = 1; i < n; i++) // run tru all the indexes from 1 to n
    {
        k = i;
        mv = n / 2;
        rev = 0;
        while (k > 0) // invert the actual index
        {
            if ((k % 2) > 0)
            rev = rev + mv;
            k = k / 2;
            mv = mv / 2; 
        }
        {  // switch the actual sample and the bitinverted one
            if (i < rev)
            {
                b = a[rev];
                a[rev] = a[i];
                a[i] = b;
            }
        }
    }
}

For each index “I” runs through all the bits by “k” and inverts them. The modulo division tells me whether the least significant bit is 1 or 0 and has to be put in reversed index or not. Starting at the most significant bit. Dividing k by 2 moves the second least significant bit one position to the right. In this way we get all the bits from least significant trough most significant and have the reversal in “rev”. Now we just have to switch the index “i” and “rev”. But just replace if “I” is smaller the “rev”. Otherwise we change all the indexes back if “I” becomes bigger than “rev”.

After this bit reversal we can calculate the sub FFT by:

C#
public void CalcSubFFT(TKomplex[] a, int n)
{
    int i, k, m;
    TKomplex w;
    TKomplex v;
    TKomplex h;
    k = 1;
    while (k <= n/2)  // we start at stage 0 whit 2 sample per transformation
    {
        m = 0;
        while (m <= (n-2*k))
        {
            for (i = m; i < m + k; i++)
            {
                w.real = Math.Cos( Math.PI * (double)(i-m)/(double)(k));
                w.imag = Math.Sin( Math.PI * (double)(i-m)/(double)(k));
                h = kprod(a[i + k], w);
                v = a[i];
                a[i] = ksum(a[i], h);
                a[i + k] = kdiff(v, h);
            }
            m = m + 2 * k;
        }
        k = k * 2;
    }
}

These two functions have to be put in main call and all the elements have to be divided by 2N and element 0 has to be divided by 4N.

C#
public void CalcFFT()
{
    int i;
    BitInvert(y, N);
    CalcSubFFT(y, N);
    for (i = 0; i < N; i++)
    {
        y[i].imag = y[i].imag / (double)N * 2.0;
        y[i].real = y[i].real / (double)N * 2.0;
    }
    y[0].imag = y[0].imag / 2.0;
    y[0].real = y[0].real / 2.0;
}

Reducing calculation time

And if we finally put the sine and cosine values into a look up table:

C#
we = new TKomplex[N / 2];
for (i = 0; i < (N / 2); i++)  // Init look up table for sine and cosine values
{
    we[i].real = Math.Cos(2* Math.PI * (double)(i) / (double)(N));
    we[i].imag = Math.Sin(2* Math.PI * (double)(i) / (double)(N));
}

And replace

C#
w.real = Math.Cos( Math.PI * (double)(i-m)/(double)(k));
w.imag = Math.Sin( Math.PI * (double)(i-m)/(double)(k));

by

C#
w.real = we[((i-m)*N / k/ 2)].real;
w.imag = we[((i-m)*N / k / 2)].imag;

That reduces the calculation time for the transformation by approx. 30 % as well. There is no difference in time consumption visible compared to the recursive algorithm as long as we just run the transformation once. But if we run it 10 times the difference becomes visible:

Comparison of time consumption while 10 transformations whit 4096 samples:

 

Calculation time

Recursive standard

45 ms

Recursive with look up table

30 ms

Not recursive standard

41 ms

Not recursive with look up table

24 ms

For more math visit www.mosismath.com :-)

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Tester / Quality Assurance Annax Switzerland AG
Switzerland Switzerland
Computers are very straight... They always do exactly what we tell them to do... Only, much too often what we tell them to do is not really what we want them to do Smile | :)

Writing Software is one of the most creative tings one can do. I have been doing this for more than ten years now and still having a lot of fun with it. Besides doing software for HMI's on C# for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time Smile | :)

For more detailed descriptions and math visit me on my own page

www.mosismath.com

Comments and Discussions

 
QuestionCannot download the sample Pin
Member 1394275718-May-19 9:11
Member 1394275718-May-19 9:11 
AnswerRe: Cannot download the sample Pin
Southmountain25-Feb-23 13:56
Southmountain25-Feb-23 13:56 
QuestionCan't Download the Sample Pin
gongchaaa30-Jun-16 16:41
gongchaaa30-Jun-16 16:41 
AnswerRe: Can't Download the Sample Pin
Mosi_6220-Aug-16 23:51
professionalMosi_6220-Aug-16 23:51 
Questionresult Pin
raoul123418-Jun-15 5:41
raoul123418-Jun-15 5:41 
AnswerRe: result Pin
Mosi_6230-Sep-15 9:09
professionalMosi_6230-Sep-15 9:09 
QuestionCross spectrum Pin
Ray17023-Sep-13 4:13
Ray17023-Sep-13 4:13 
AnswerRe: Cross spectrum Pin
Mosi_6226-Sep-13 8:47
professionalMosi_6226-Sep-13 8:47 
GeneralRe: Cross spectrum Pin
Ray17026-Sep-13 11:38
Ray17026-Sep-13 11:38 
QuestionSource Code for FFT code Pin
Alex Toh18-Jul-13 20:51
Alex Toh18-Jul-13 20:51 
AnswerRe: Source Code for FFT code Pin
Mosi_6220-Jul-13 5:24
professionalMosi_6220-Jul-13 5:24 
QuestionComparision Pin
Steffen Villadsen13-Jul-13 2:00
Steffen Villadsen13-Jul-13 2:00 
SuggestionFile Upload? Pin
firmwaredsp12-Jul-13 20:49
firmwaredsp12-Jul-13 20:49 
GeneralRe: File Upload? Pin
Mosi_6213-Jul-13 3:00
professionalMosi_6213-Jul-13 3:00 
GeneralRe: File Upload? Pin
firmwaredsp13-Jul-13 13:02
firmwaredsp13-Jul-13 13:02 
GeneralRe: File Upload? Pin
Mosi_6213-Jul-13 22:49
professionalMosi_6213-Jul-13 22:49 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.