/* -----------------------------------------------------------------------------
	GSFramework
	Copyright 2001-2013 Emmanuel Julien. All Rights Reserved.
----------------------------------------------------------------------------- */


		#include <cmath>
		#include "picture/pict.h"


//------------------------------------------------------------------------------
void				FFT(int size, bool inverse, float *inReal, float *inIm, float *outReal, float *outIm)
{
	// Calculate m = log_2(n).
	int m = 0, p = 1;
	for	(; p < size; ++m)
		p *= 2;

	// Bit reversal.
	outReal[size - 1] = inReal[size - 1];
	outIm[size - 1] = inIm[size - 1];

	int j = 0;
	for	(int i = 0; i < size - 1; ++i)
	{
		outReal[i] = inReal[j];
		outIm[i] = inIm[j];

		int k = size / 2;
		while (k <= j)
		{
			j -= k;
			k /= 2;
		}

		j += k;
	}

	// Calculate the FFT.
	float ca = -1.0, sa = 0.0;
	int l1 = 1, l2 = 1;

	for	(int l = 0; l < m; ++l)
	{
		l1 = l2;
		l2 *= 2;

		float u1 = 1.0, u2 = 0.0;

		for(int j = 0; j < l1; j++)
		{
			for(int i = j; i < size; i += l2)
			{
				int i1 = i + l1;

				float t1 = u1 * outReal[i1] - u2 * outIm[i1],
					  t2 = u1 * outIm[i1] + u2 * outReal[i1];

				outReal[i1] = outReal[i] - t1;
				outIm[i1] = outIm[i] - t2;
				outReal[i] += t1;
				outIm[i] += t2;
			}

			double z =  u1 * ca - u2 * sa;
			u2 = u1 * sa + u2 * ca;
			u1 = (float)z;
		}

		sa = (float)sqrt((1.f - ca) / 2.f);
		if	(!inverse)
			sa = -sa;
		ca = (float)sqrt((1.f + ca) / 2.f);
	}

	// Divide through n if it isn't the IDFT.
	if	(!inverse)
		for	(int i = 0; i < size; ++i)
		{
			outReal[i] /= size;
			outIm[i] /= size;
		}
}

void		testFFT()
{
	float inR[8], inI[8], outR[8], outI[8];

	for (int n = 0; n < 8; ++n)
	{
		inR[n] = 5;
		inI[n] = 2;
	}

	FFT(8, false, inR, inI, outR, outI);

	for (int n = 0; n < 8; ++n)
	{
		inR[n] = 0;
		inI[n] = 0;
	}

	FFT(8, true, outR, outI, inR, inI);
}
//------------------------------------------------------------------------------