#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <malloc.h>

#define  maxIndex  10000L

/************************************************************************
  fft(int n, double xRe[], double xIm[], double yRe[], double yIm[])
 ------------------------------------------------------------------------
  Input/output:
      int n          transformation length.
      double xRe[]   real part of input sequence.
      double xIm[]   imaginary part of input sequence.
      double yRe[]   real part of output sequence.
      double yIm[]   imaginary part of output sequence.
 ------------------------------------------------------------------------
  Function:
      The procedure performs a fast discrete Fourier transform (FFT) of
      a complex sequence, x, of an arbitrary length, n. The output, y,
      is also a complex sequence of length n.

      y[k] = sum(x[m]*exp(-i*2*pi*k*m/n), m=0..(n-1)), k=0,...,(n-1)

      The largest prime factor of n must be less than or equal to the
      constant maxPrimeFactor defined below.
 ------------------------------------------------------------------------
  Implementation notes:
      The general idea is to factor the length of the DFT, n, into
      factors(2,3,4,5) that are efficiently handled by the routines.

      A number of short DFT's are implemented with a minimum of
      arithmetical operations and using (almost) straight line code
      resulting in very fast execution when the factors of n belong
      to this set. Especially the butterfly-unit of radix-3 and 5 is optimized.

      Any other factors, that are not in the set of {2,3,4,5} will result in error!
 ------------------------------------------------------------------------
  The following procedures are used :
      factorize       :  factor the transformation length.
      transTableSetup :  setup table with sofar-, actual-
	                     remainRadix- and Table_twiddle.
      permute         :  permutation allows in-place calculations.
      twiddleTransf   :  twiddle multiplications and DFT's for one stage.
      fft_4           :  length 4 DFT, a la Nussbaumer.
      fft_5           :  length 5 DFT, a la Nussbaumer.
*************************************************************************/

#define  maxPrimeFactor        5
#define  maxFactorCount        20

/* cos(2*pi/3)-1 not cos(2*pi/3) in order to utilize formal-seat operation */
static double  c3_1 = -1.5000000000000E+00;  /*  c3_1 = cos(2*pi/3)-1;          */
static double  c3_2 =  8.6602540378444E-01;  /*  c3_2 = sin(2*pi/3);            */
                                          
static double  u5   =  1.2566370614359E+00;  /*  u5   = 2*pi/5;                 */

/* (cos(u5)+cos(2*u5))/2-1 not (cos(u5)+cos(2*u5))/2 in order to utilize formal-seat operation*/
static double  c5_1 = -1.2500000000000E+00;  /*  c5_1 = (cos(u5)+cos(2*u5))/2-1;*/
static double  c5_2 =  5.5901699437495E-01;  /*  c5_2 = (cos(u5)-cos(2*u5))/2;  */
static double  c5_3 = -9.5105651629515E-01;  /*  c5_3 = -sin(u5);               */
static double  c5_4 = -1.5388417685876E+00;  /*  c5_4 = -(sin(u5)+sin(2*u5));   */
static double  c5_5 =  3.6327126400268E-01;  /*  c5_5 = (sin(u5)-sin(2*u5));    */

static double   pi;
static int      groupOffset,dataOffset,blockOffset,adr;
static int      groupNo,dataNo,blockNo,twNo;
static double   omega, tw_re,tw_im;
static double  cosw, sinw;
static double   twiddleRe[maxPrimeFactor], twiddleIm[maxPrimeFactor],
                zRe[maxPrimeFactor], zIm[maxPrimeFactor];

 
FILE *fp;
/*factor the transformation length
  Note : this part is just for simulate, it's not cotained in the hardware, 
  all the factors are parameters acquired from outside
 */
void factorize(int n, int *nFact, int fact[])
{
    int i,j;
    int nRadix;
	int radices[5];
    int factors[maxFactorCount];

    nRadix    =  4;  
    radices[1]=  2;
    radices[2]=  3;
    radices[3]=  4;
    radices[4]=  5;

	//factorizing process the sequence is falling order
    if (n==1)
    {
        j=1;
        factors[1]=1;
    }
    else j=0;
    i=nRadix;
    while ((n>1) && (i>0))
    {
      if ((n % radices[i]) == 0)
      {
        n=n / radices[i];
        j=j+1;
        factors[j]=radices[i];
      }
      else  i=i-1;
    }
    
	//change the sequence into the rising order
    for (i=1; i<=j; i++)         
    {
      fact[i] = factors[j-i+1];
    }
    *nFact=j;
}   /* factorize */

/****************************************************************************
  After N is factored the parameters that control the stages are generated.
  For each stage we have:
    sofar   : the product of the radices so far.
    actual  : the radix handled in this stage.
    remain  : the product of the remaining radices.

	for example : when the N= 3*3*5*5= 225
	then    actual,  sofar(the number of groups),              remain (No. of members in each group)
	          3        1 (only one group in the first stage)    3*5*5  (but 75 butterfly-units in each group with radix=3)   
			  3        3 (3 groups in the second stage)         5*5  (but 25 butterfly-units in each group with radix=3)
			  5        3*3 (9 groups in the third stage)        5   (but 5 butterfly-units in each group with radix=5)
			  5        3*3*5 (45  groups in the fouth stage)    1  (but only one butterfly-unit in each group with radix=5)

    Table_twiddleRe : real part of the twiddle factors
    Table_twiddleIm : image part of the twiddle factors
 ****************************************************************************/

/* this table also input from outside before caculate
 */
void transTableSetup(int sofar[], int actual[], int remain[],
                     int *nFact,
                     int *nPoints)
{
    int i;

    factorize(*nPoints, nFact, actual);
    if (actual[1] > maxPrimeFactor)
    {
        printf("\nPrime factor of FFT length too large : %6d",actual[1]);
        printf("\nPlease modify the value of maxPrimeFactor in mixfft.c");
        exit(1);
    }
    remain[0]=*nPoints;
    sofar[1]=1;
    remain[1]=*nPoints / actual[1];
    for (i=2; i<=*nFact; i++)
    {
        sofar[i]=sofar[i-1]*actual[i-1];
        remain[i]=remain[i-1] / actual[i];
    }
 
}   /* transTableSetup */

/****************************************************************************
  The sequence y is the permuted input sequence x so that the following
  transformations can be performed in-place, and the final result is the
  normal order. This process is sensitive to all the factors.
  Example:
  x[](normal order) -> permute(x[]) (arbitrary order)-> fft() -> y[](normal order)
 ****************************************************************************/

void permute(int nPoint, int nFact,
             int fact[], int remain[],
             double xRe[], double xIm[],
             double yRe[], double yIm[])

{
    int i,j,k;
    int count[maxFactorCount]; 


    for (i=1; i<=nFact; i++) count[i]=0;
    k=0;

    for (i=0; i<=nPoint-2; i++)
    {
        yRe[i] = xRe[k];
        yIm[i] = xIm[k];
        j=1;
        k=k+remain[1];

        count[1] = count[1]+1;
        while (count[j] >= fact[j])
        {
            count[j]=0;
            k=k-remain[j-1]+remain[j+1];
            j=j+1;
            count[j]=count[j]+1;
        }
    }
    yRe[nPoint-1]=xRe[nPoint-1];
    yIm[nPoint-1]=xIm[nPoint-1];
}   /* permute */


/****************************************************************************
  Twiddle factor multiplications and transformations are performed on a
  group of data. The number of multiplications with 1 are reduced by skipping
  the twiddle multiplication of the first stage and of the first group of the
  following stages.
 ***************************************************************************/
void scin_table (int sofar_mul_Radix)
{
	/* 72 different probabilities*/
     switch(sofar_mul_Radix) {
         case     3 : cosw= -0.500000; break;
         case    12 : cosw= 0.866025; break;
         case    60 : cosw= 0.994522; break;
         case   300 : cosw= 0.999781; break;
         case  1500 : cosw= 0.999991; break;
         case     9 : cosw= 0.766044; break;
         case    36 : cosw= 0.984808; break;
         case   180 : cosw= 0.999391; break;
         case   900 : cosw= 0.999976; break;
         case    27 : cosw= 0.973045; break;
         case   108 : cosw= 0.998308; break;
         case   540 : cosw= 0.999932; break;
         case  2700 : cosw= 0.999997; break;
         case    81 : cosw= 0.996993; break;
         case   324 : cosw= 0.999812; break;
         case  1620 : cosw= 0.999992; break;
         case   243 : cosw= 0.999666; break;
         case   972 : cosw= 0.999979; break;
         case   729 : cosw= 0.999963; break;
         case  2916 : cosw= 0.999998; break;
         case     2 : cosw= -1.000000; break;
         case     6 : cosw= 0.500000; break;
         case    24 : cosw= 0.965926; break;
         case   120 : cosw= 0.998630; break;
         case   600 : cosw= 0.999945; break;
         case  3000 : cosw= 0.999998; break;
         case    18 : cosw= 0.939693; break;
         case    72 : cosw= 0.996195; break;
         case   360 : cosw= 0.999848; break;
         case  1800 : cosw= 0.999994; break;
         case    54 : cosw= 0.993238; break;
         case   216 : cosw= 0.999577; break;
         case  1080 : cosw= 0.999983; break;
         case   162 : cosw= 0.999248; break;
         case   648 : cosw= 0.999953; break;
         case   486 : cosw= 0.999916; break;
         case  1944 : cosw= 0.999995; break;
         case    48 : cosw= 0.991445; break;
         case   240 : cosw= 0.999657; break;
         case  1200 : cosw= 0.999986; break;
         case   144 : cosw= 0.999048; break;
         case   720 : cosw= 0.999962; break;
         case   432 : cosw= 0.999894; break;
         case  2160 : cosw= 0.999996; break;
         case  1296 : cosw= 0.999988; break;
         case    96 : cosw= 0.997859; break;
         case   480 : cosw= 0.999914; break;
         case  2400 : cosw= 0.999997; break;
         case   288 : cosw= 0.999762; break;
         case  1440 : cosw= 0.999990; break;
         case   864 : cosw= 0.999974; break;
         case  2592 : cosw= 0.999997; break;
         case   192 : cosw= 0.999465; break;
         case   960 : cosw= 0.999979; break;
         case   576 : cosw= 0.999941; break;
         case  2880 : cosw= 0.999998; break;
         case  1728 : cosw= 0.999993; break;
         case   384 : cosw= 0.999866; break;
         case  1920 : cosw= 0.999995; break;
         case  1152 : cosw= 0.999985; break;
         case   768 : cosw= 0.999967; break;
         case  2304 : cosw= 0.999996; break;
         case  1536 : cosw= 0.999992; break;
         case  3072 : cosw= 0.999998; break;
         case     8 : cosw= 0.707107; break;
         case    32 : cosw= 0.980785; break;
         case   128 : cosw= 0.998795; break;
         case     4 : cosw= 0.000000; break;
         case    16 : cosw= 0.923880; break;
         case    64 : cosw= 0.995185; break;
         case   256 : cosw= 0.999699; break;
         case   512 : cosw= 0.999925; break;
         case  1024 : cosw= 0.999981; break;	   
         default : cosw= 0; break;
	 }
	 /* 72 different probabilities*/
     switch(sofar_mul_Radix) {
         case     3 : sinw= -0.866025; break;
         case    12 : sinw= -0.500000; break;
         case    60 : sinw= -0.104528; break;
         case   300 : sinw= -0.020942; break;
         case  1500 : sinw= -0.004189; break;
         case     9 : sinw= -0.642788; break;
         case    36 : sinw= -0.173648; break;
         case   180 : sinw= -0.034899; break;
         case   900 : sinw= -0.006981; break;
         case    27 : sinw= -0.230616; break;
         case   108 : sinw= -0.058145; break;
         case   540 : sinw= -0.011635; break;
         case  2700 : sinw= -0.002327; break;
         case    81 : sinw= -0.077492; break;
         case   324 : sinw= -0.019391; break;
         case  1620 : sinw= -0.003878; break;
         case   243 : sinw= -0.025854; break;
         case   972 : sinw= -0.006464; break;
         case   729 : sinw= -0.008619; break;
         case  2916 : sinw= -0.002155; break;
         case     2 : sinw= -0.000000; break;
         case     6 : sinw= -0.866025; break;
         case    24 : sinw= -0.258819; break;
         case   120 : sinw= -0.052336; break;
         case   600 : sinw= -0.010472; break;
         case  3000 : sinw= -0.002094; break;
         case    18 : sinw= -0.342020; break;
         case    72 : sinw= -0.087156; break;
         case   360 : sinw= -0.017452; break;
         case  1800 : sinw= -0.003491; break;
         case    54 : sinw= -0.116093; break;
         case   216 : sinw= -0.029085; break;
         case  1080 : sinw= -0.005818; break;
         case   162 : sinw= -0.038775; break;
         case   648 : sinw= -0.009696; break;
         case   486 : sinw= -0.012928; break;
         case  1944 : sinw= -0.003232; break;
         case    48 : sinw= -0.130526; break;
         case   240 : sinw= -0.026177; break;
         case  1200 : sinw= -0.005236; break;
         case   144 : sinw= -0.043619; break;
         case   720 : sinw= -0.008727; break;
         case   432 : sinw= -0.014544; break;
         case  2160 : sinw= -0.002909; break;
         case  1296 : sinw= -0.004848; break;
         case    96 : sinw= -0.065403; break;