FAST FOURIER TRANSFORM USING C PROGRAM IN DSP - Code

Latest

programs

Thursday, November 10, 2011

FAST FOURIER TRANSFORM USING C PROGRAM IN DSP


 FAST FOURIER TRANSFORM
AIM
To find the Fast Fourier Transform for the realtime  samples.
ALGORITHM
Step 1     sample the input (N)  of any desired frequency. Convert it to fixed-point format and scale the input to avoid overflow during manipulation.
Step 2     Declare four buffers namely real input, real exponent, imaginary exponent and imaginary input.
Step 3     Declare three counters for stage, group and butterfly.
Step 4     Implement the Fast Fourier Transform for the input signal.
Step 5     Store the output (Real and Imaginary) in the output buffer.
Step 6     Decrement the counter of butterfly. Repeat from the Step 4 until the counter reaches zero.
Step 7     If the butterfly counter is zero, modify the exponent value.
Step 8     Repeat from the Step 4 until the group counter reaches zero.
Step 9     If the group counter is zero, multiply the butterfly value by two and divide the group value by two.
Step 10   Repeat from the Step 4 until the stage counter reaches zero.
Step 11   Transmit the FFT output through line out port.


PROGRAM: -


#include <math.h>                         
#define PTS 128    //# of points for FFT
#define PI 3.14159265358979
typedef struct {float real,imag;} COMPLEX;
void FFT(COMPLEX *Y, int n);     //FFT prototype
float iobuffer[PTS];                          //as input and output buffer
float x1[PTS],x[PTS];                                   //intermediate buffer
short i;                                   //general purpose index variable                 
short buffercount = 0;           //number of new samples in iobuffer        
short flag = 0;                                    //set to 1 by ISR when iobuffer full  
float y[128];
COMPLEX w[PTS];                        //twiddle constants stored in w
COMPLEX samples[PTS];                           //primary working buffer                                              

main()
{

   float j,sum=0.0 ;
   int n,k,i,a;
 for (i = 0 ; i<PTS ; i++)        // set up twiddle constants in w
  {
   w[i].real = cos(2*PI*i/(PTS*2.0)); //Re component of twiddle constants
   w[i].imag =-sin(2*PI*i/(PTS*2.0));  /*Im component of twiddle constants*/
  }




 
/****************InputSignalX(n)***************************************/  
 
      for(i=0,j=0;i<PTS;i++)
     { x[i] = sin(2*PI*5*i/PTS);      // Signal x(Fs)=sin(2*pi*f*i/Fs);
       samples[i].real=0.0;
       samples[i].imag=0.0;
     }
   


/********************** FFT of R(t) *****************************/

   for (i = 0 ; i < PTS ; i++)    //swap buffers
    {
     samples[i].real=iobuffer[i]; //buffer with new data
    }
   
   for (i = 0 ; i < PTS ; i++)
       samples[i].imag = 0.0;                //imag components = 0

   FFT(samples,PTS);              //call function FFT.c

/******************** PSD *******************************************/
  
   for (i = 0 ; i < PTS ; i++)    //compute magnitude
    {
     x1[i] = sqrt(samples[i].real*samples[i].real
         + samples[i].imag*samples[i].imag);         
      }
 
}                                                        //end of main

void FFT(COMPLEX *Y, int N)      //input sample array, # of points    
{
 COMPLEX temp1,temp2;             //temporary storage variables         
 int i,j,k;                       //loop counter variables              
 int upper_leg, lower_leg;                 //index of upper/lower butterfly leg  
 int leg_diff;                    //difference between upper/lower leg  
 int num_stages = 0;              //number of FFT stages (iterations)
 int index, step;                 //index/step through twiddle constant
 i = 1;                           //log(base2) of N points= # of stages 
 do
  {
   num_stages +=1;
   i = i*2;
  }while (i!=N);
 leg_diff = N/2;                                 //difference between upper&lower legs
 step = (PTS*2)/N;                                        //step between values in twiddle.h   // 512          
 for (i = 0;i < num_stages; i++)  //for N-point FFT                





 {

   index = 0;
   for (j = 0; j < leg_diff; j++)
    {
     for (upper_leg = j; upper_leg < N; upper_leg += (2*leg_diff))
      {  
       lower_leg = upper_leg+leg_diff;
       temp1.real = (Y[upper_leg]).real + (Y[lower_leg]).real;
       temp1.imag = (Y[upper_leg]).imag + (Y[lower_leg]).imag;
       temp2.real = (Y[upper_leg]).real - (Y[lower_leg]).real;
       temp2.imag = (Y[upper_leg]).imag - (Y[lower_leg]).imag;
       (Y[lower_leg]).real = temp2.real*(w[index]).real
                                              -temp2.imag*(w[index]).imag;
       (Y[lower_leg]).imag = temp2.real*(w[index]).imag
                                              +temp2.imag*(w[index]).real;
       (Y[upper_leg]).real = temp1.real;
       (Y[upper_leg]).imag = temp1.imag;
      }
   index += step;
    }
     leg_diff = leg_diff/2;
     step *= 2;
      }
       j = 0;
       for (i = 1; i < (N-1); i++)     //bit reversal for resequencing data
      {
     k = N/2;
     while (k <= j)
      {
       j = j - k;
       k = k/2;
      }
     j = j + k;
     if (i<j)
      {
       temp1.real = (Y[j]).real;
       temp1.imag = (Y[j]).imag;
       (Y[j]).real = (Y[i]).real;
       (Y[j]).imag = (Y[i]).imag;
       (Y[i]).real = temp1.real;
       (Y[i]).imag = temp1.imag;
      }
    }
   return;
}




OUTPUT:
  

                             DFT or FFT spectrum of sinusoidal signal f= 10 Hz

No comments:

Post a Comment