#include<stdio.h>
#include<math.h>
#define nFFT 16
#define nFFT3_4 ((nFFT*3)/4)

double wn_FFT[nFFT3_4];  /* for twiddle factor */
short br_FFT[nFFT];     /* for bit reversal   */
double xr[nFFT], xi[nFFT], yr[nFFT], yi[nFFT];

void fft_table(double wn_FFT[], short br_FFT[], int N_FFT);
void fft(double xr[], double xi[], double wn_FFT[], short br_FFT[], int N_FFT);
void swap(double *a, double *b);

/***************************/
/* Genarate tables for FFT */
/***************************/
void fft_table(double wn_FFT[], short br_FFT[], int N_FFT)
{
    int i, n_half, ne, jp;
    double arg;
/* Calculation of twiddle factor */
    arg = 6.283185307f/N_FFT;
    for (i=0; i<((N_FFT*3)>>2); i++) wn_FFT[i] = cos(arg*i);
/* Calculation of bit reversal table */
    n_half = N_FFT>>1;
    br_FFT[0] = 0;
    for (ne=1; ne<N_FFT; ne=ne<<1)
    {
        for (jp=0; jp<ne; jp++) br_FFT[jp+ne] = br_FFT[jp] + n_half;
        n_half = n_half>>1;
    }
}

/**************************/
/* FFT for complex signal */
/**************************/
void fft(double xR[], double xI[], double wn_FFT[], short br_FFT[], int N_FFT)
{
    double xtmpR, xtmpI;
    int   j, jnh, k, jxC, jxS, ne, n_half, n_half2;

    n_half = N_FFT>>1;
    for (ne=1; ne<N_FFT; ne=ne<<1)
    {
        n_half2 = n_half<<1;
        for (k=0; k<N_FFT; k=k+n_half2)
        {
            jxC = 0;
            jxS = N_FFT>>2;
            for (j=k; j<(k+n_half); j++)
            {
                jnh = j + n_half;
/* beginning of butterfly operations */
                xtmpR = xR[j];
                xtmpI = xI[j];
                xR[j] = xtmpR + xR[jnh];
                xI[j] = xtmpI + xI[jnh];
                xtmpR = xtmpR - xR[jnh];
                xtmpI = xtmpI - xI[jnh];
                xR[jnh] = xtmpR*wn_FFT[jxC] - xtmpI*wn_FFT[jxS];
                xI[jnh] = xtmpR*wn_FFT[jxS] + xtmpI*wn_FFT[jxC];
/* end of butterfly operations */
                jxC = jxC + ne;
                jxS = jxS + ne;
            }
        }
        n_half = n_half>>1;
    }
/* Bit reverse */
    for (j=0; j<N_FFT; j++)
        if (j<br_FFT[j])
        {
             swap(&xR[j], &xR[br_FFT[j]]);
             swap(&xI[j], &xI[br_FFT[j]]);
        }
}

/* used in FFT procedure */
void swap(double *a, double *b)
{
    double tmp;
    tmp = *a;
    *a = *b;
    *b = tmp;
}

void main()
{
 short i;

    fft_table(wn_FFT, br_FFT, nFFT);       /* ganarate tables for FFT */
/* genarate rectangular wave */

    for (i=0; i<nFFT; i++)
    {                                   
            xr[i] =sin(6.28*i/3);
    }

    for (i=0; i<nFFT; i++)
    {
        yr[i] = xr[i];
        yi[i] = xi[i];
    }

    fft(yr, yi, wn_FFT, br_FFT, nFFT);     /* FFT */

 for(i=0;i<nFFT;i++)
 {
  printf("%f ",xr[i]);
 }
 printf("\n\n ");
 for(i=0;i<nFFT;i++)
 {
  printf("%f ",yi[i]);
 }
}