Friday, October 16, 2009

Radix 2 floating point FFT Implementation

I searched for fft code over the internet but couldn't find a simple code for the FFT. So I am posting the code. I think I have explained everything in the code code is what follows:

/*
The program follows is for the FFT of a one dimensional input
fft algorithm used is radix 2
so works only for input lengths powers of 2
input is complex and output obviously is complex
changes can be made to use the function for your requirement.
I have tried to comment the code anyway there shouldnt be any
problem its such a small program.
*/
//Header Files
#include <stdio.h>
#include<stdlib.h>
#include<math.h>

//Complex number structure
typedef struct complex
{
float real;            //Real part of the number
float imaginary;    //Imaginary part of the number
} complex;

//Function to calculate the FFT Method used is Radix 2
complex* fft(complex *in,int N){
//Local variables
complex *X,*Xg,*Xh;
complex *xe,*xo;
int i=0;
//Memory allocation
xe = (complex*)malloc(sizeof(complex)*(N/2));
xo = (complex*)malloc(sizeof(complex)*(N/2));
X = (complex *)malloc(sizeof(complex)*N);
//every recursion has to stop
//so here is the stopping condition of the recursion
if(N==2){
X[0].real        = in[0].real + in[1].real;
X[0].imaginary    = in[0].imaginary + in[1].imaginary;
X[1].real        = in[0].real - in[1].real;
X[1].imaginary    = in[0].imaginary - in[1].imaginary;
return X;
}
//Breaking the input for recursion
for (i=0;i<N;i++){
if(i%2==0){//even indexes
xe[i/2].real        = in[i].real;
xe[i/2].imaginary    = in[i].imaginary;
}
else{//odd indexes
xo[(i-1)/2].real        = in[i].real;
xo[(i-1)/2].imaginary    = in[i].imaginary;
}
}
//recursive calls
Xg = fft(xe,N/2);
Xh = fft(xo,N/2);
//Combining the two halfs here is where magic happens
for (i=0;i<N;i++){
X[i].real=Xg[i%(N/2)].real + Xh[i%(N/2)].real*cos(2*3.14*i/N)+Xh[i%(N/2)].imaginary*sin(2*3.14*i/N);
X[i].imaginary=Xg[i%(N/2)].imaginary+Xh[i%(N/2)].imaginary*cos(2*3.14*i/N)-Xh[i%(N/2)].real*sin(2*3.14*i/N);
}
//deallocating the waste memory
free(Xg);
free(Xh);
return X;//returning the result
}
//the main function
int main(int argc,char *argv[])
{
int N,i;
complex *x;//input
complex *X;//variable to store FFT of the input
//Note : N has to be a power of 2 otherwise the function doesnt work
//asking for the no of points in the FFT
printf("\nenter N must be a power of 2:");
scanf("%d",&N);
//allocating memory for the input
x = (complex*)malloc(N*sizeof(complex));
//taking the input
for(i=0;i<N;i++){
printf("\nenter %dth element real then imaginary",i);
scanf("%f",&(x[i].real));  
scanf("%f",&(x[i].imaginary));
}
//calling the fft function
X = fft(x,N);
//displaying the result
printf("\nFFT of the input is as follows");
for(i=0;i<N;i++){
printf("\nXr[%d] = %.2f + j %.2f",i,X[i].real,X[i].imaginary);
}
return 0;
}
//end of the program

1 comment:

  1. Thanks for sharing this info here can you let me know how can I use it Kindly share this with all details

    Khurram Shahzad
    Laptop is Much Batter

    ReplyDelete