// FFT1D.cpp: implementation of the CFFT1D class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FFT1D.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

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

CFFT1D::CFFT1D(unsigned int size)
{
	if(!IsPowerOfTwo(size))
		size = (unsigned int)pow(2, 1 + (int)(log(size)/log(2.0)));

	pInput = new CComplex[size];
	pOutput = new CComplex[size];
	nDimension = size;
}

CFFT1D::~CFFT1D()
{
	if(pInput != NULL)
	{
		delete[] pInput;
		pInput = NULL;
	}

	if(pOutput != NULL)
	{
		delete[] pOutput;
		pOutput = NULL;
	}
}

CComplex CFFT1D::GetOutputAt(unsigned int index)
{
	if(index >=0 && index < nDimension)
		return pOutput[index];
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return NULL;
	}
}

bool CFFT1D::SetOutputAt(unsigned int index, double real/* = 0 */, double imag/* = 0 */)
{
	if(index >=0 && index < nDimension)
	{
		pOutput[index].SetReal(real);
		pOutput[index].SetImag(imag);
		return true;
	}
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return false;
	}
}

CComplex CFFT1D::GetInputAt(unsigned int index)
{
	if(index >=0 && index < nDimension)
		return pInput[index];
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return NULL;
	}
}

bool CFFT1D::SetInputAt(unsigned int index, double real/* = 0 */, double imag/* = 0 */ )
{
	if(index >=0 && index < nDimension)
	{
		pInput[index].SetReal(real);
		pInput[index].SetImag(imag);
		return true;
	}
	else
	{
		#ifdef USE_MESSAGE_BOX_FOR_FAILURE		
			AfxMessageBox("Index out of range!");
		#endif
		return false;
	}
}

void CFFT1D::ForwardTransform()
{
	BitReverseCopy();
	unsigned int numBits = NumberOfBitsNeeded(nDimension);
	int m = 1;
	for(unsigned int s = 1; s <= numBits; s++)
	{
		m *= 2;
		CComplex wm(cos(-2.0 * PI / m), sin(-2.0 * PI / m));
		for(unsigned int k = 0; k < nDimension; k += m)
		{
			CComplex w(1);
			for(int j = 0; j < m/2; j++)
			{
				CComplex t = w * pOutput[k + j + m/2];
				CComplex u = pOutput[k + j];
				pOutput[k + j] = u + t;
				pOutput[k + j + m/2] = u - t;
				w *= wm;
			}
		}
	}
}

void CFFT1D::ReverseTransform()
{
	CComplex* temp = new CComplex[nDimension];
	for(unsigned int i = 0; i < nDimension; i++)
		temp[i] = pOutput[i];

	unsigned int numBits = NumberOfBitsNeeded(nDimension);
	for(i = 0; i < nDimension; i++)
		pOutput[ReverseBits(i, numBits)] = temp[i];

	delete[] temp;

	int m = 1;
	for(unsigned int s = 1; s <= numBits; s++)
	{
		m *= 2;
		CComplex wm(cos(2.0 * PI / m), sin(2.0 * PI / m));
		for(unsigned int k = 0; k < nDimension; k += m)
		{
			CComplex w(1);
			for(int j = 0; j < m/2; j++)
			{
				CComplex t = w * pOutput[k + j + m/2];
				CComplex u = pOutput[k + j];
				pOutput[k + j] = u + t;
				pOutput[k + j + m/2] = u - t;
				w *= wm;
			}
		}
	}
}
////////////////////////////////////
// private helper functions
///////////////////////////////////
void CFFT1D::BitReverseCopy()
{
	unsigned int numBits = NumberOfBitsNeeded(nDimension);
	for(unsigned int i = 0; i < nDimension; i++)
		pOutput[ReverseBits(i, numBits)] = pInput[i] / (double)nDimension;
}

bool CFFT1D::IsPowerOfTwo (unsigned int x)
{
    if ( x < 2 )
        return false;

    if ( x & (x-1) )        // Thanks to 'byang' for this cute trick!
        return false;

    return true;
}

unsigned int CFFT1D::NumberOfBitsNeeded ( unsigned int PowerOfTwo )
{
    unsigned int i;

    if ( PowerOfTwo < 2 )
    {
		//cerr << ">>> Error: argument " << PowerOfTwo << " to NumberOfBitsNeeded is too small." << endl;
        exit(1);
    }

    for ( i=0; ; i++ )
    {
        if ( PowerOfTwo & (1 << i) )
            return i;
    }
}



unsigned int CFFT1D::ReverseBits ( unsigned int index, unsigned int NumBits )
{
    unsigned int i, rev;

    for ( i = rev = 0; i < NumBits; i++ )
    {
        rev = (rev << 1) | (index & 1);
        index >>= 1;
    }

    return rev;
}