C Radix-2 FFT Source
Algorithms in C - fft.c
<source lang="c">
- include <math.h>
- 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>