#include <stdio.h>
#include <strings.h>
#include <math.h>

FILE *fopen(), *fdata;

void main()
{
#define pi 3.14159265358979323846
#define dB(x) ( ( (fabs(x)) > (0.0) ) ? (20.*log10(fabs(x))) : (-1000.) )
#define Q 8
float ar[Q], ai[Q];    /* array of points */
float freq_test, mag;
int i;

fdata = fopen("data", "w");

/* generate original signal */
freq_test = (float)Q/(float)4;
for (i=0; i<Q; i++)
 { ar[i] = cos(freq_test*2.*pi*((float)i/(float)Q));
   ai[i] = sin(freq_test*2.*pi*((float)i/(float)Q)); }

fft(Q, ar, ai);

for(i=0; i<Q; i++)
 { mag = sqrt(ar[i]*ar[i]+ai[i]*ai[i]);
fprintf(fdata, "%6.4lf %20.16f %8.2f\n",                                                     (float)i/(float)Q, mag, dB(mag) ); }

} /* main */

fft(int N, float *ar, float *ai)
/*
in-place complex fft

After Cooley, Lewis, and Welch; from Rabiner & Gold (1975)

 program adapted from FORTRAN 
    by K. Steiglitz  (ken@princeton.edu)
       Computer Science Dept. 
       Princeton University 08544          */
{

#define verbose 2     /* 2 = print input and output */

#define pi 3.14159265358979323846

int i, j, k, L;            /* indexes */
int M, TEMP, LE, LE1, ip;  /* M = log N */
int NV2, NM1;
float t;               /* temp */
float Ur, Ui, Wr, Wi, Tr, Ti;
float Ur_old;

NV2 = N/2;
NM1 = N-1;
TEMP = N; /* get M = log N */
M = 0;
while(TEMP>1){
 if(TEMP%2 != 0){
  printf("Error: N not a power of 2, stopping\n");
  exit(0); }
 TEMP /=2;  
 M ++; }
printf("FFT with N = %d, M = log N = %d\n", N, M); 

if(verbose >= 2)
{ printf("input:\n");
  printf("real: ");
  for (i=0; i<N; i++)
   printf("%6.4lf ", ar[i]);
  printf("\nimag: ");
  for (i=0; i<N; i++)
   printf("%6.4lf ", ai[i]);
  printf("\n\n"); }

/* shuffle */
j = 1;
for (i=1; i<=NM1; i++)
 { if(i<j)              /* swap a[i] and a[j] */
    { t = ar[j-1];     
      ar[j-1] = ar[i-1];
      ar[i-1] = t;
      t = ai[j-1];
      ai[j-1] = ai[i-1];
      ai[i-1] = t; }

   k = NV2;             /* bit-reversed counter */
   while(k<j)
    { j -= k;
      k /= 2; }

   j += k; }

if(verbose >= 2)
{ printf("shuffled input:\n");
  printf("real: ");
  for (i=0; i<N; i++)
   printf("%6.4lf ", ar[i]);
  printf("\nimag: "); 
  for (i=0; i<N; i++)
   printf("%6.4lf ", ai[i]);
  printf("\n\n"); }

LE = 1.;
for (L=1; L<=M; L++)                 /* stage L */
 { LE1 = LE;                         /* LE1 = LE/2 */ 
   LE  *= 2;                         /* LE  = 2^L */
   Ur = 1.0;
   Ui = 0.; 
   Wr = cos(M_PI/(float)LE1);
   Wi = -sin(M_PI/(float)LE1); /* Cooley, Lewis, and Welch have "+" here */
   for (j=1; j<=LE1; j++)

    { for (i=j; i<=N; i+=LE) /* butterfly */
       { ip = i+LE1;
         Tr = ar[ip-1]*Ur-ai[ip-1]*Ui;
         Ti = ar[ip-1]*Ui+ai[ip-1]*Ur;
         ar[ip-1] = ar[i-1] - Tr;
         ai[ip-1] = ai[i-1] - Ti;
         ar[i-1]  = ar[i-1] + Tr;
         ai[i-1]  = ai[i-1] + Ti; }  /* end of butterfly */

      Ur_old = Ur;
      Ur = Ur_old*Wr-Ui*Wi;
      Ui = Ur_old*Wi+Ui*Wr; }         /* end of j loop */

 } /* end of stage L */

if (verbose>=2)
 { printf("output:\n");
   printf("real: ");
   for(i=0; i<N; i++)
    printf("%6.4lf ", ar[i]);
   printf("\nimag: ");
   for(i=0; i<N; i++)
    printf("%6.4lf ", ai[i]);
   printf("\n"); }

} /* fft */