// Function FFT_func(COMPLEX *X, COMPLEX *W)
// This is a radix 2 decimation-in-time FFT algorithm
// Computations are done in-place, so the
// the output overwrites the input in the array X[].                  
// This algorithm computes an N-point DFT, where N 
// is defined in FFT_header.
// A global declaration of the N/2 twiddle factors, W[]
// must be done by the calling function.
/*1 */ #include "FFT_header.h"
/*2 */ 
/*3 */ void FFT_func(COMPLEX *X, COMPLEX *W)
/*4 */ {
/*5 */   
/*6 */   COMPLEX temp;          // temporary storage complex variable
/*7 */   short i,j,k;           // loop indices
/*8 */   short i_lower;         // Index of lower point in butterfly
/*9 */   short step;
/*10*/   short stage;           // FFT stage
/*11*/   short DFTpts;          // Number of points in sub DFT in current stage
/*12*/                          // and offset to next DFT in stage
/*13*/   short numBF;           // Number of butterflys in one DFT in current
/*14*/                          // stage. Also is offset to lower
/*15*/                          // point in butterfly in current stage
/*16*/ 
/*17*/   // Bit-reversing algorithm
/*18*/   // The index j is the bit reversed value of i.
/*19*/   // k is used to propagate carryover bits to the right.
/*20*/   // Since 0 -> 0 and N-1 -> N-1 under bit-reversal, these two 
/*21*/   // reversals are skipped.
/*22*/ 
/*23*/   j = 0;
/*24*/ 
/*25*/   for(i=1; i<(N-1); i++)
/*26*/   {
/*27*/     k = N2;               // k is 1 in msb, 0 elsewhere (i.e. k=1<<(M-1))
/*28*/ 
/*29*/     while(k<=j)           // Propagate carry to the right if bit is 1
/*30*/     {
/*31*/       j = j - k;          // Bit tested is 1, so clear it.
/*32*/       k = k>>1;           // Carryover binary 1 one bit to right.
/*33*/     }
/*34*/       
/*35*/ 	j = j+k;              // current bit tested is 0, add 1 to that bit
/*36*/ 
/*37*/     // Swap samples if current index is less than bit reversed index.
/*38*/ 
/*39*/ 	if(i<j)
/*40*/     {
/*41*/       temp.real = X[j].real;    // Hold value at bit reversed location.
/*42*/       temp.imag = X[j].imag;
/*43*/       X[j].real = X[i].real;    // Insert value at current location
/*44*/       X[j].imag = X[i].imag;    // at bit reversed location.
/*45*/       X[i].real = temp.real;    // Insert value previous stored at bit reversed
/*46*/       X[i].imag = temp.imag;    // location at current location
/*47*/     }
/*48*/   }
/*49*/ 
/*50*/   // Do M stages of butterflys
/*51*/ 
/*52*/   step = N2;                      // step = N/2, N/4, N/8, ... 1
/*53*/ 
/*54*/   for(stage=1; stage <= M; stage++)
/*55*/   {
/*56*/     DFTpts = 1 << stage;          // DFTpts = 2^stage = points is sub DFT
/*57*/     numBF = DFTpts/2;             // Number of butterflies in sub-DFT
/*58*/     k = 0;                        // initial twiddle factor index
/*59*/     
/*60*/     // Do butterflys for current stage
/*61*/ 
/*62*/     for(j=0; j<numBF; j++)        // Do the numBF butterflys per sub DFT
/*63*/     { 
/*64*/       // Compute butterflys that use same twiddle factor, W[k]
/*65*/       for(i=j; i<N; i += DFTpts)
/*66*/       {
/*67*/         i_lower = i + numBF;      // Index of lower point in butterfly
/*68*/         temp.real = X[i_lower].real * W[k].real + X[i_lower].imag * W[k].imag;
/*69*/         temp.imag = X[i_lower].imag * W[k].real - X[i_lower].real * W[k].imag;
/*70*/ 
/*71*/         X[i_lower].real = X[i].real - temp.real;
/*72*/         X[i_lower].imag = X[i].imag - temp.imag;
/*73*/ 
/*74*/         X[i].real  = X[i].real + temp.real;
/*75*/         X[i].imag  = X[i].imag + temp.imag;
/*76*/ 	  }
/*77*/ 	  
/*78*/ 	  k += step;                  // increment twiddle index
/*79*/ 	}
/*80*/ 	
/*81*/ 	step = step/2;                // calculate step for next stage
/*82*/   }
/*83*/   return;
/*84*/}