Click here to Skip to main content
12,250,675 members (42,848 online)
Click here to Skip to main content
Add your own
alternative version

Stats

1.3M views
23K downloads
188 bookmarked
Posted

How to implement the FFT algorithm

, 2 Feb 2005
Rate this:
Please Sign up or sign in to vote.
An article on how to implement the FFT algorithm in C, C++ or C#.

Sample Image

Introduction

Basically, this article describes one way to implement the 1D version of the FFT algorithm for an array of complex samples. The intention of this article is to show an efficient and fast FFT algorithm that can easily be modified according to the needs of the user. I've studied the FFT algorithm when I was developing a software to make frequency analysis on a sample of captured sound.

Background

First of all, 95% of this code wasn't written by me. This is practically the code that is described in the book Numerical Recipes In C of 1982!!! Yes, more than 20 years ago!!! But still, in my opinion, very, very good. But this code is slightly different from the original one. When I was studying the algorithm, I had noticed a pattern that could be exploited, and based on that, I've managed to improve the algorithm with a small change in the code, and the O() (The Big O) (unit measure to the complexity of the algorithm) is reduced in N computations. (After implementing the improved method successfully, I made some research in the web and realized I discovered nothing new. I've just noticed something that someone had already seen :-P)

I will not get "deep in theory", so I strongly advise the reading of chapter 12 if you want to understand "The Why". Other forms of the FFT like the 2D or the 3D FFT can be found on the book too.

The FFT

The Fast Fourier Transform is an optimized computational algorithm to implement the Discreet Fourier Transform to an array of 2^N samples. It allows to determine the frequency of a discreet signal, represent the signal in the frequency domain, convolution, etc... This algorithm has a complexity of O(N*log2(N)). Actually, the complexity of the algorithm is a little higher because the data needs to be prepared by an operation called bit-reversal. This bit-reversal section is presented in the Numerical Recipes In C as a O(2N) complexity. With a small change I've made to the code presented here, it makes it in O(N). This represents something like an 8% improvement of performance.

Example of a signal in the frequency domain.

The FFT is calculated in two parts. The first one transforms the original data array into a bit-reverse order array by applying the bit-reversal method. This makes the mathematical calculations of the second part "much more easy". The second part processes the FFT in N*log2(N) operations (application of the Danielson-Lanzcos algorithm).

Let's start with an array of complex data. This array could be, for example, in this case, an array of floats in witch the data[even_index] is the real part and the data[odd_index] is the complex part. (This can be adapted to an array of real data, just by filling the complex values with 0s, or use the real array FFT implemented on the book.) The size of the array must be in an N^2 order (2, 4, 8, 16, 32, 64, etc...). In case the sample doesn't match that size, just put it in an array with the next 2^N size and fill the remaining spaces with 0s.

Just a small and not very significant consideration: the original code uses data arrays considering the beginning of the information is in index 1 -> data[1], and data[0] is ignored. My code modifies that to start from 0 -> data[0].

First, we define the FFT function:

//data -> float array that represent the array of complex samples
//number_of_complex_samples -> number of samples (N^2 order number) 
//isign -> 1 to calculate FFT and -1 to calculate Reverse FFT
float FFT (float data[], unsigned long number_of_complex_samples, int isign)
{
    //variables for trigonometric recurrences
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

The Bit-Reversal Method

First, the original array must be transformed in order to perform the Danielson-Lanzcos algorithm. For example, the complex[index] will swap places with the complex[bit-reverse-order-index]. If the index (in binary) is 0b00001, the bit-reverse-order-index will be 0b10000. In figure 1, you can see what happens to the data array after the transformation.

The implementation of this method according to Numerical Recipes In C goes like this:

//the complex array is real+complex so the array 
//as a size n = 2* number of complex samples
// real part is the data[index] and the complex part is the data[index+1]
n=number_of_complex_samples * 2;

//binary inversion (note that 
//the indexes start from 1 witch means that the
//real part of the complex is on the odd-indexes
//and the complex part is on the even-indexes
j=1;
for (i=1;i<n;i+=2) {
    if (j > i) {
        //swap the real part
        SWAP(data[j],data[i]);
        //swap the complex part
        SWAP(data[j+1],data[i+1]);
    }
    m=n/2;
    while (m >= 2 && j > m) {
        j -= m;
        m = m/2;
    }
    j += m;
}

The SWAP goes outside the function and it's something like this:

#defineSWAP(a,b)tempr=(a);(a)=(b);(b)=tempr
//tempr is a variable from our FFT function

Figure 1 (8-length data array)

If you pay attention at figure 1, you can see a pattern. Let's see: divide the array in half with a mirror. If you look at the reflecting side of the mirror, you will see exactly the same thing of what's in the other side. This mirrored effect allows you to do the bit-reversal method in the first half of the array and use it almost directly in the second half. But you must be careful with one thing. You can only apply this effect if the change happens in the first half only. This means that if the change is between an index of the first half and an index from the second, this is not valid, otherwise you would be making the swap and then undoing it again (do this on a 16 length data array and you will understand what I'm saying). So the code becomes something like this:

//the complex array is real+complex so the array 
//as a size n = 2* number of complex samples
// real part is the data[index] and 
//the complex part is the data[index+1]
n=number_of_complex_samples * 2;

//binary inversion (note that the indexes 
//start from 0 witch means that the
//real part of the complex is on the even-indexes 
//and the complex part is on the odd-indexes
j=0;
for (i=0;i<n/2;i+=2) {
    if (j > i) {
        //swap the real part
        SWAP(data[j],data[i]);
        //swap the complex part
        SWAP(data[j+1],data[i+1]);
        // checks if the changes occurs in the first half
        // and use the mirrored effect on the second half
        if((j/2)<(n/4)){
            //swap the real part
            SWAP(data[(n-(i+2))],data[(n-(j+2))]);
            //swap the complex part
            SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
        }
    }
    m=n/2;
    while (m >= 2 && j >= m) {
        j -= m;
        m = m/2;
    }
    j += m;
}

The Danielson-Lanzcos

The second half of the code goes exactly like it is described in the book. This applies the N*log2(N) trigonometric recurrences to the data. The code I will show here is my version (data starts at index 0):

//Danielson-Lanzcos routine 
mmax=2;
//external loop
while (n > mmax)
{
    istep = mmax<<  1;
    theta=sinal*(2*pi/mmax);
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    //internal loops
    for (m=1;m<mmax;m+=2) {
        for (i= m;i<=n;i+=istep) {
            j=i+mmax;
            tempr=wr*data[j-1]-wi*data[j];
            tempi=wr*data[j]+wi*data[j-1];
            data[j-1]=data[i-1]-tempr;
            data[j]=data[i]-tempi;
            data[i-1] += tempr;
            data[i] += tempi;
        }
        wr=(wtemp=wr)*wpr-wi*wpi+wr;
        wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
}

How to use the FFT

Let's see now how to use the FFT. Imagine that you want to collect a sample of a signal or a function, and you want to know the fundamental frequency of it. This sample may come from any source: a function that you've inserted in the code, a piece of captured sound, etc...

Let's say that the signal is a real array signal (just like the sound capture buffer), how do I use the FFT??? First of all, you need to choose the FFT variant that you will use. There is a specific variant for real arrays, but in this case, I will use this. It's not the most efficient, but it's easy to use.

Next concern is the amount of data you're going to send to the FFT and the sample rate. The sample rate must be a 2^N number, but you don't need to send an array of 2^N samples to be processed (read the NR for different implementations). You can just send 50 samples, for example, and fill the remaining array with 0s. But remember, the more data you send for calculation, the more precise is the FFT.

After the real array has been passed to a complex array with the complex part equal to 0, you compute the FFT.

And now for the results

After the FFT is calculated, you can use the complex array that resulted from the FFT to extract the conclusions.

If you are interested in knowing the fundamental frequency of the signal, find the absolute maximum of the array, and the frequency will be given by the index of the array. The absolute value of a complex number is the square root of the square of the real plus the square of the complex.

If the absolute maximum occurs in the complex number given by indexes [102][103] (real, complex), then your fundamental frequency will be (102/2)=61Hz. You have to divide it in half because the array is twice longer (remember: real, complex), so the result is not the index position (102), but half (61).

If your intention is to draw the Fourier signal, it goes the same way. Frequency 30 is given by the absolute value of complex [30][31], etc. etc. ...remember. The second half of the computed FFT array must be ignored due to the Nyquist redundancy (the minimum sample rate must twice the highest frequency of the signal). It's only a mirror of the first half. If you want to measure frequencies up to 6000, you will need the next 2^N number next to 2*6000.

See the example were I apply the FFT to a Sine signal. It's on the OnPaint function of the CChildView class. The FFT is implemented on the CFourier class. The example is a stupid example and has a stupid structure, but I think it's easy to understand. Change the parameters, play with it, try different things, and see the results. This way, you will be able to take your own conclusions.

Precautions

If you intend to use this FFT implementation, read the NR license.

There are lots of documentation on this matter. It's not an easy thing to understand, but I think it's a very interesting subject. Wink | ;)

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

A list of licenses authors might use can be found here

Share

About the Author

João Martins
Software Developer
Portugal Portugal
I'm that guy that falls for the impossible...

You may also be interested in...

Comments and Discussions

 
QuestionAm I the only one who noticed.... Pin
DrScruff28-Feb-16 17:41
memberDrScruff28-Feb-16 17:41 
QuestionQuestion sqrt(real*real+img*img) Pin
Member 1172954929-May-15 16:59
memberMember 1172954929-May-15 16:59 
QuestionRe: Question sqrt(real*real+img*img) Pin
g.saldanha9-Jun-15 2:21
professionalg.saldanha9-Jun-15 2:21 
Questionthanks Pin
masoudRajaei7-May-15 7:36
membermasoudRajaei7-May-15 7:36 
GeneralMy vote 5 Pin
mhleecode9-Feb-15 5:35
membermhleecode9-Feb-15 5:35 
GeneralMy vote of 5 Pin
Adam Roderick J16-May-14 0:04
memberAdam Roderick J16-May-14 0:04 
GeneralMy vote of 5 Pin
Hessam Jalali13-May-14 22:08
memberHessam Jalali13-May-14 22:08 
Questioncant open your code Pin
Member 1012084918-Dec-13 5:03
memberMember 1012084918-Dec-13 5:03 
AnswerRe: cant open your code Pin
Eric Hdz20-Aug-15 14:07
memberEric Hdz20-Aug-15 14:07 
BugErrata Pin
live3v1l19-Feb-13 13:57
memberlive3v1l19-Feb-13 13:57 
GeneralMy vote of 5 Pin
xmyy12097-Jan-13 21:33
memberxmyy12097-Jan-13 21:33 
QuestionHow to Use This FFT to Detect Pitch...Pls Help Me... Pin
ngelmusic14-Oct-12 22:47
memberngelmusic14-Oct-12 22:47 
AnswerRe: How to Use This FFT to Detect Pitch...Pls Help Me... Pin
João Martins15-Oct-12 22:58
memberJoão Martins15-Oct-12 22:58 
GeneralRe: How to Use This FFT to Detect Pitch...Pls Help Me... Pin
Mrugrajsinh24-Oct-12 4:43
memberMrugrajsinh24-Oct-12 4:43 
Questioni cant compile your project, there are lot of errors Pin
Aydar Omurbekov2-Jun-12 20:42
memberAydar Omurbekov2-Jun-12 20:42 
GeneralMy vote of 1 Pin
Mehran Taheri2-Apr-12 4:48
memberMehran Taheri2-Apr-12 4:48 
GeneralRe: My vote of 1 Pin
Aydar Omurbekov2-Jun-12 20:48
memberAydar Omurbekov2-Jun-12 20:48 
QuestionKeeps crashing Pin
yellobello31-Aug-11 21:52
memberyellobello31-Aug-11 21:52 
QuestionSign error in Daneilson-Lanczos Summation Pin
RunawayStapler12-Jul-11 8:01
memberRunawayStapler12-Jul-11 8:01 
GeneralReal & Imaginary appear to be switched (in my VB6 version) Pin
JohnAugHarris17-Mar-11 12:59
memberJohnAugHarris17-Mar-11 12:59 
GeneralCannot compile please help Pin
AhmetCakir15-Jun-10 0:01
memberAhmetCakir15-Jun-10 0:01 
GeneralRe: Cannot compile please help Pin
Radu Simionescu5-Dec-10 12:38
memberRadu Simionescu5-Dec-10 12:38 
QuestionResolution Pin
zl3d20-Sep-09 11:02
memberzl3d20-Sep-09 11:02 
GeneralThis doesn't work when I mix two signals together Pin
cbranje12-May-09 10:09
membercbranje12-May-09 10:09 
GeneralRe: This doesn't work when I mix two signals together Pin
Guenter Marx13-Jan-10 23:23
memberGuenter Marx13-Jan-10 23:23 
QuestionHow can I solve this? Pin
catalin7712-Apr-09 23:24
membercatalin7712-Apr-09 23:24 
AnswerRe: How can I solve this? Pin
hnchang23-Dec-10 19:24
memberhnchang23-Dec-10 19:24 
GeneralOk it is not accurate either. But I think I have a fix. Pin
quincy45120-Mar-09 5:08
memberquincy45120-Mar-09 5:08 
GeneralWell it works but it is way too slow to use. Pin
quincy45112-Mar-09 9:58
memberquincy45112-Mar-09 9:58 
GeneralRe: Well it works but it is way too slow to use. Pin
João Martins13-Mar-09 0:24
memberJoão Martins13-Mar-09 0:24 
Generalinterupting FFT output data. Pin
quincy4512-Mar-09 15:56
memberquincy4512-Mar-09 15:56 
Generalaudio data sample how to a load data[] so I get the correct fundamental frequency. Pin
quincy45127-Feb-09 7:31
memberquincy45127-Feb-09 7:31 
GeneralRe: audio data sample how to a load data[] so I get the correct fundamental frequency. Pin
quincy4512-Mar-09 15:51
memberquincy4512-Mar-09 15:51 
GeneralRe: audio data sample how to a load data[] so I get the correct fundamental frequency. Pin
ngelmusic14-Oct-12 22:50
memberngelmusic14-Oct-12 22:50 
GeneralThis article is a ripoff of code from Numerical Recipes in C Pin
pacman12813-Oct-08 9:19
memberpacman12813-Oct-08 9:19 
GeneralRe: This article is a ripoff of code from Numerical Recipes in C Pin
Ahimoth16-Oct-08 12:41
memberAhimoth16-Oct-08 12:41 
GeneralRe: This article is a ripoff of code from Numerical Recipes in C Pin
João Martins24-Nov-08 8:51
memberJoão Martins24-Nov-08 8:51 
Generalorder n Pin
appel45675-Oct-08 10:22
memberappel45675-Oct-08 10:22 
GeneralRe: order n Pin
João Martins9-Oct-08 3:32
memberJoão Martins9-Oct-08 3:32 
GeneralRe: order n Pin
lievenvv2-Nov-09 6:27
memberlievenvv2-Nov-09 6:27 
QuestionNot Very accurate. Please Help Pin
coerciveutopian3-Apr-08 9:11
membercoerciveutopian3-Apr-08 9:11 
AnswerRe: Not Very accurate. Please Help Pin
Guenter Marx13-Jan-10 22:37
memberGuenter Marx13-Jan-10 22:37 
GeneralLP,HP,BP filters sample code Pin
kelvinic28-Jan-08 15:49
memberkelvinic28-Jan-08 15:49 
QuestionCan not compile... Pin
Oner Dikdere (TR)1-Dec-07 11:55
memberOner Dikdere (TR)1-Dec-07 11:55 
GeneralRe: Can not compile... Pin
Oner Dikdere (TR)7-Dec-07 22:31
memberOner Dikdere (TR)7-Dec-07 22:31 
GeneralI Wouldn't Even Think About Using This Code Pin
mirage59930-Jul-07 9:18
membermirage59930-Jul-07 9:18 
GeneralProblem in finding Fundamental frequency Pin
smzhaq20-Jun-07 19:40
membersmzhaq20-Jun-07 19:40 
GeneralSample Rate Pin
smzhaq15-Jun-07 5:00
membersmzhaq15-Jun-07 5:00 
QuestionPlease help me Pin
mrk100009-Jun-07 14:46
membermrk100009-Jun-07 14:46 
GeneralEqualizer Pin
mrk100009-Jun-07 13:30
membermrk100009-Jun-07 13:30 

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.

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.160426.1 | Last Updated 2 Feb 2005
Article Copyright 2005 by João Martins
Everything else Copyright © CodeProject, 1999-2016
Layout: fixed | fluid