// FFT.cpp: implementation of the FFT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FFT.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

FFT::FFT()
{

}
FFT::FFT(int len, double inRe[],double inIm[])
{
	//預設參數初始化
	c3_1 = -1.5000000000000E+00;  /*  c3_1 = cos(2*pi/3)-1;          */
	c3_2 =  8.6602540378444E-01;  /*  c3_2 = sin(2*pi/3);            */
	
	u5   =  1.2566370614359E+00;  /*  u5   = 2*pi/5;                 */
	c5_1 = -1.2500000000000E+00;  /*  c5_1 = (cos(u5)+cos(2*u5))/2-1;*/
	c5_2 =  5.5901699437495E-01;  /*  c5_2 = (cos(u5)-cos(2*u5))/2;  */
	c5_3 = -9.5105651629515E-01;  /*  c5_3 = -sin(u5);               */
	c5_4 = -1.5388417685876E+00;  /*  c5_4 = -(sin(u5)+sin(2*u5));   */
	c5_5 =  3.6327126400268E-01;  /*  c5_5 = (sin(u5)-sin(2*u5));    */
	c8   =  7.0710678118655E-01;  /*  c8 = 1/sqrt(2);    */
	
	radixLen=6;
	radices[1]=  2;
    radices[2]=  3;
    radices[3]=  4;
    radices[4]=  5;
    radices[5]=  8;
    radices[6]= 10;
	
	setLength(len);

	setX(inRe,inIm);
}

FFT::~FFT()
{

}

void FFT::setX(double inRe[], double inIm[])
{
	int i;

	for( i = 0 ; i <length ; i++ )
	{
		xRe[i]=inRe[i];
		xIm[i]=inIm[i];
	}
}

void FFT::getY(double outRe[], double outIm[])
{
	int i;

	for( i = 0 ; i <= maxIndex ; i++ )
	{
		outRe[i]=yRe[i];
		outIm[i]=yIm[i];
	}
}

void FFT::setLength(int len)
{
	if(len>=0)
		length=len;
	else{
		cout<<"FAIL! The FFT length cannot be negative!"<<endl;
		
	}
}
int FFT::getLength(void)
{
	return length;
}

void FFT::setRadices(int rLen, int r[])
{
	int i;
	radixLen=rLen;
	for( i=1 ; i <= radixLen ; i++ )
		radices[i]=r[i];
}


void FFT::fft(void)
{
	int sofarRadix[maxFactorCount], 
		actualRadix[maxFactorCount], 
		remainRadix[maxFactorCount];
  
    int count;
	factorize();
    transTableSetup(sofarRadix, actualRadix, remainRadix);
    permute(remainRadix);
	
    for (count=1; count<=factLen; count++)
		twiddleTransf(sofarRadix[count], actualRadix[count], remainRadix[count]);

}


void FFT::factorize(void)
{
	int i,j,k,ll;
    int factors[maxFactorCount];
	
	ll=length;

    if (ll==1)
    {
        j=1;
        factors[1]=1;
    }
    else j=0;
    i=radixLen;
    while ((ll>1) && (i>0))
    {
		if ((ll % radices[i]) == 0)
		{
			ll=ll / radices[i];
			j=j+1;
			factors[j]=radices[i];
		}
		else  i=i-1;
    }
    if (factors[j] == 2)   /*substitute factors 2*8 with 4*4 */
    {   
		i = j-1;
		while ((i>0) && (factors[i] != 8)) i--;
		if (i>0)
		{
			factors[j] = 4;
			factors[i] = 4;
		}
    }
    if (ll>1)
    {
        for (k=2; k<sqrt(ll)+1; k++)
            while ((ll % k) == 0)
            {
                ll= ll / k;
                j=j+1;
                factors[j]=k;
            }
			if (ll>1)
			{
				j=j+1;
				factors[j]=ll;
			}
    }               
    for (i=1; i<=j; i++)         
    {
		fact[i] = factors[j-i+1];
    }
    factLen=j;
}



void FFT::transTableSetup(int sofar[], int actual[], int remain[])
{
	int i;
	

    if (fact[1] > maxPrimeFactor)
    {
        cout<<"\nPrime factor of FFT length is overrun : "<< fact[1]
			<<"\nPlease modify the value of maxPrimeFactor in StdAfx.h"<<endl;
        exit(1);
    }

	for( i = 1 ; i <= factLen ; i++ )
		actual[i] = fact[i];

    remain[0]=length;
    sofar[1]=1;
    remain[1]=length / actual[1];
    for (i=2; i<=factLen; i++)
    {
        sofar[i]=sofar[i-1]*actual[i-1];
        remain[i]=remain[i-1] / actual[i];
    }
}

void FFT::permute(int remain[])
{
	int i,j,k;
    int count[maxFactorCount]; 
	
    for (i=1; i<=factLen; i++) 
		count[i]=0;
    k=0;
    for (i=0; i<=length-2; i++)
    {
        yRe[i] = xRe[k];
        yIm[i] = xIm[k];
        j=1;
        k=k+remain[j];
        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[length-1]=xRe[length-1];
    yIm[length-1]=xIm[length-1];
}

void FFT::twiddleTransf(int sofarRadix, int radix, int remainRadix)
{
	double  cosw, sinw, gem;
    double  t1_re,t1_im, t2_re,t2_im, t3_re,t3_im;
    double  t4_re,t4_im, t5_re,t5_im;
    double  m2_re,m2_im, m3_re,m3_im, m4_re,m4_im;
    double  m1_re,m1_im, m5_re,m5_im;
    double  s1_re,s1_im, s2_re,s2_im, s3_re,s3_im;
    double  s4_re,s4_im, s5_re,s5_im;


    initTrig(radix);
    omega = 2*pi/(double)(sofarRadix*radix);
    cosw =  cos(omega);
    sinw = -sin(omega);
    tw_re = 1.0;
    tw_im = 0;
    dataOffset=0;
    groupOffset=dataOffset;
    adr=groupOffset;
    for (dataNo=0; dataNo<sofarRadix; dataNo++)
    {
        if (sofarRadix>1)
        {
            twiddleRe[0] = 1.0; 
            twiddleIm[0] = 0.0;
            twiddleRe[1] = tw_re;
            twiddleIm[1] = tw_im;
            for (twNo=2; twNo<radix; twNo++)
            {
                twiddleRe[twNo]=tw_re*twiddleRe[twNo-1]
                               - tw_im*twiddleIm[twNo-1];
                twiddleIm[twNo]=tw_im*twiddleRe[twNo-1]
                               + tw_re*twiddleIm[twNo-1];
            }
            gem   = cosw*tw_re - sinw*tw_im;
            tw_im = sinw*tw_re + cosw*tw_im;
            tw_re = gem;                      
        }
        for (groupNo=0; groupNo<remainRadix; groupNo++)
        {
            if ((sofarRadix>1) && (dataNo > 0))
            {
                zRe[0]=yRe[adr];
                zIm[0]=yIm[adr];
                blockNo=1;
                do {
                    adr = adr + sofarRadix;
                    zRe[blockNo]=  twiddleRe[blockNo] * yRe[adr]
                                 - twiddleIm[blockNo] * yIm[adr];
                    zIm[blockNo]=  twiddleRe[blockNo] * yIm[adr]
                                 + twiddleIm[blockNo] * yRe[adr]; 
                    
                    blockNo++;
                } while (blockNo < radix);
            }
            else
                for (blockNo=0; blockNo<radix; blockNo++)
                {
                   zRe[blockNo]=yRe[adr];
                   zIm[blockNo]=yIm[adr];
                   adr=adr+sofarRadix;
                }
            switch(radix) {
              case  2  : gem=zRe[0] + zRe[1];
                         zRe[1]=zRe[0] -  zRe[1]; zRe[0]=gem;
                         gem=zIm[0] + zIm[1];
                         zIm[1]=zIm[0] - zIm[1]; zIm[0]=gem;
                         break;
              case  3  : t1_re=zRe[1] + zRe[2]; t1_im=zIm[1] + zIm[2];
                         zRe[0]=zRe[0] + t1_re; zIm[0]=zIm[0] + t1_im;
                         m1_re=c3_1*t1_re; m1_im=c3_1*t1_im;
                         m2_re=c3_2*(zIm[1] - zIm[2]);