/********************************************************

			fft()

**********************************************************

Subroutine does radix two in-place decimation in time fft
Program given on page 367 of Rabiner and Gold,
Digital Signal Processing
*/

#include <math.h>
#include "complex.h"
#include "pi.h"

fft(a,m)

COMPLEX a[];	/* array of complex input values */
int m;		/* dimension of array is n=2**m  */

{
	COMPLEX u,w,t;
	int i,j,k,n,nv2,nm1;
	int le,l,le1,ip;

	n = 1;
	for(i=0; i<m; ++i) n *= 2;
	nv2=n/2;
	nm1=n-1;

/* Apply bit reversal algorithm to prepare for in-place fft */

j=1;
for(i=1; i<=nm1; ++i)
	{
	 if(i < j)
		{
		t = a[j-1];
		a[j-1] = a[i-1];
		a[i-1] = t;
		}
	 k = nv2;
	 while(k < j)
		{
		 j -= k;
		 k /= 2;
		}
	 j += k;
	}

/* Now do stage by stage radix two fft in place */

le = 1;
for(l=1; l<=m; ++l)
	{
	 le *= 2;
	 le1=le/2;

	 u.real = 1.0;
	 u.imag = 0.0;

	 w.real = cos(PI/le1);
	 w.imag = sin(PI/le1);

	 for(j=1; j<=le1; j++)
		{
		 for(i=j; i<=n; i += le)
			{
			 ip=i+le1;
			 t = cmpy(a[ip-1],u);
			 a[ip-1] = csub(a[i-1],t);
			 a[i-1] = cadd(a[i-1],t);
			}
	 	 u = cmpy(u,w);
		}
	}
return;
}
