fft_example.c

// Replaces data[0..2*nn-1] by its discrete Fourier transform,                                                   
// if isign is input as 1; or replaces data[0..2*nn-1] by nn times its inverse                                   
// discrete Fourier transform, if isign is input as −1.                                                          
// data is a complex array of length nn or, equivalently, a real array                                           
// of length 2*nn.                                                                                               
// nn MUST be an integer power of 2 (this is not checked for!).                                                  
void dfft(float data[], unsigned int nn, int isign) {
  initialize_sincos_tables();
  unsigned int n, mmax, m, j, istep, i;
  double wtemp, wr, wpr, wpi, wi, theta;
  float tempr, tempi;
  n = nn << 1;
  j = 1;

  // This is the bit-reversal section of the routine.                                                            
  for (i = 1; i < n; i += 2) {
    if (j > i) {
      SWAP(data[j-1], data[i-1]);
      // Exchange the two complex numbers.                                                                       
      SWAP(data[j], data[i]);
    }
    m = nn;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  // Here  begins  the  Danielson-Lanczos  section  of  the  routine.                                            
  mmax = 2;
  uint theta_index = 1;
  while (n > mmax) {
    // Outer loop executed log2(nn) times.                                                                       
    istep = mmax << 1;
    // Initialize the trigonometric recurrence.                                                                  
    //The angle considered is theta = isign * (6.28318530717959 / mmax)                                          
    wtemp = sin_table[theta_index + 1]; // abs(sin(theta / 2))                                                   
    wpr = -2.0 * wtemp * wtemp;
    wpi = isign * sin_table[theta_index]; // sin(theta)                                                          
    wr = 1.0;
    wi = 0.0;
    for (m = 1; m < mmax; m += 2) { // Here are the two nested inner loops.                                      
      for (i = m; i <= n; i += istep) {
        j = i + mmax;
        // This  is  the  Danielson-Lanczos  formula:                                                            
        tempr = (float)(wr * data[j - 1] - wi * data[j]);
        tempi = (float)(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;
      }
      wtemp = wr;
      wr += wr * wpr - wi * wpi;
      //Trigonometric  recurrence.                                                                               
      wi += wi * wpr + wtemp * wpi;
    }
    mmax = istep;
    theta_index++;
  }
}