C Radix-2 FFT Source

Algorithms in C - fft.c

<source lang="c">

1. include <math.h>

1. define FFT_LENGTH 32

/* Constants */ const double TWO_PI = 2.0 * 3.1415926535902069414;

/* Global Variables */ static double sine_table_real[ FFT_LENGTH ]; static double sine_table_imag[ FFT_LENGTH ]; static int current_sine_table_length = 0;

/* Function declarations */ void fft( double* const values_real, double* const values_imag ); void inverse_fft( double* const values_real, double* const values_imag ); void compute_sine_table();

/* Calculate an in-place fourier transform using radix-2 algorithm */ /* Based on Matlab source code by M.D.Macleod, CUED (13/11/2001) */ void fft( double* const values_real, double* const values_imag ) {

  double temp_real, temp_imag;
  int i1, i2;

  double omega_real, omega_imag;
  double product_real, product_imag;
  int half_length, max_half_len, time_step, sine_table_index;
  
  /* Compute the sine table if necessary */
  if( current_sine_table_length == 0 )
  {
     compute_sine_table();
  }

  /* Now swap values with the element at the bit reversed location */
  /*  e.g. If there are 8 elements then 0-0, 1-4, 2-2, 3-6, 7-7    */
  i2 = 0;
  for( i1 = 0; i1 < FFT_LENGTH; i1++ )
  {
     if( i2 > i1 )
     {
        temp_real = values_real[ i1 ];
        temp_imag = values_imag[ i1 ];

        values_real[ i1 ] = values_real[ i2 ];
        values_imag[ i1 ] = values_imag[ i2 ];

        values_real[ i2 ] = temp_real;
        values_imag[ i2 ] = temp_imag;
     }

     /* Find the next location for i2 */
     half_length = FFT_LENGTH / 2;

     while( half_length >= 1 && i2 >= half_length )
     {
        i2 -= half_length;
        half_length /= 2;
     }
     i2 += half_length;
  }

  /* FIRST STAGE Omega is 0 or 1 : i1=i1+i2 , i2=i1-i2   */
  for( i1 = 0; i1 < FFT_LENGTH; i1 += 2 )
  {
     i2 = i1 + 1;

     temp_real = values_real[ i1 ];
     temp_imag = values_imag[ i1 ];

     values_real[ i1 ] = temp_real + values_real[ i2 ];
     values_imag[ i1 ] = temp_imag + values_imag[ i2 ];

     values_real[ i2 ] = temp_real - values_real[ i2 ];
     values_imag[ i2 ] = temp_imag - values_imag[ i2 ];
  }

  /* SECOND STAGE Omega is 0 or 1 real                   */
  for( i1 = 0; i1 < FFT_LENGTH; i1 += 4 )
  {
     i2 = i1 + 2;

     temp_real = values_real[ i1 ];
     temp_imag = values_imag[ i1 ];

     /* i1 = i1 + i2 */
     values_real[ i1 ] = temp_real + values_real[ i2 ];
     values_imag[ i1 ] = temp_imag + values_imag[ i2 ];

     /* i2 = i1 - i2 */
     values_real[ i2 ] = temp_real - values_real[ i2 ];
     values_imag[ i2 ] = temp_imag - values_imag[ i2 ];
  }

  /* SECOND STAGE Omega is 0 or 1 real minus 1 imaginary */
  for( i1 = 1; i1 < FFT_LENGTH; i1 += 4 )
  {
     i2 = i1 + 2;

     temp_real = values_real[ i1 ];
     temp_imag = values_imag[ i1 ];

     /* i1 = i1 + -i * i2 */
     values_real[ i1 ] = temp_real + values_imag[ i2 ];
     values_imag[ i1 ] = temp_imag - values_real[ i2 ];

     /* i2 = i1 - -i * i2 */
     product_imag = values_real[ i2 ];
     values_real[ i2 ] = temp_real - values_imag[ i2 ];
     values_imag[ i2 ] = temp_imag + product_imag;
  }

  /* REMAINING STAGES Omega must be calcuated          */
  /* time_step = 2 * FFT_LENGTH / ( 2 * max_half_len ) */
  time_step = FFT_LENGTH / 4;
  
  /* We've done the simple Omega values already, so start from 4 */
  for( max_half_len = 4; FFT_LENGTH > max_half_len; max_half_len *= 2)
  {
     /* As max_half_len is doubled, the time_step is halved */
     time_step /= 2;

     /* Step through at the max Half length              */
     /* This is no longer the original data in the array */
     for( half_length = 0; half_length < max_half_len; half_length++ )
     {
        sine_table_index = half_length * time_step;
        omega_real = sine_table_real[ sine_table_index ];
        omega_imag = sine_table_imag[ sine_table_index ];

        /* Compute once for each half_length in the data array */
        for( i1 = half_length; i1 < FFT_LENGTH; i1 += 2*max_half_len )
        {
           i2 = i1 + max_half_len;

           temp_real = values_real[ i1 ];
           temp_imag = values_imag[ i1 ];

           product_real = omega_real * values_real[ i2 ] -
                          omega_imag * values_imag[ i2 ];

           product_imag = omega_real * values_imag[ i2 ] +
                          omega_imag * values_real[ i2 ];

           /* i1 = i1 + omega * i2 */
           values_real[ i1 ] = temp_real + product_real;
           values_imag[ i1 ] = temp_imag + product_imag;

           /* i2 = i1 - omega * i2 */
           values_real[ i2 ] = temp_real - product_real;
           values_imag[ i2 ] = temp_imag - product_imag;
        }
     }
  }

}

/* Inverse Fast Fourier Transform */ void inverse_fft( double* const values_real, double* const values_imag ) {

  int i;

  /* Replace the elements with their complex conjugate */
  for( i = 0; i < FFT_LENGTH; i++ )
  {
     values_imag[ i ] = -values_imag[ i ];
  }

  /* Use the FFT function to compute the Fourier Transform */
  fft( values_real, values_imag );

  /* Replace the elements with their complex conjugate */
  for( i = 0; i < FFT_LENGTH; i++ )
  {
     values_real[ i ] =  values_real[ i ] / (double)FFT_LENGTH;
     values_imag[ i ] = -values_imag[ i ] / (double)FFT_LENGTH;
  }

}

/* Precompute the sine table to avoid calculating it repeatedly */ void compute_sine_table() {

  int i;

  /* Compute the cos and -sin of angles from 0 to 2 pi radians */
  /* -ve sign for Forward FFT ( i.e. exp() = cos() - i*sin() ) */
  for( i = 0; i < FFT_LENGTH; i++ )
  {
     double t = TWO_PI * i / (double)FFT_LENGTH;
     sine_table_real[ i ] =  cos( t );
     sine_table_imag[ i ] = -sin( t );
  }
  
  current_sine_table_length = FFT_LENGTH;

} </source>