/*************************************************************************
 **  FILE NAME="fft001.c"                                               **
 **       <<< Fast Fourier Transformation ( radix 2 ) >>>               **
 **     programed by Motoi Fujita   1993/05/24 (translated from FORTRAN)**
 **                        revision 1993/11/09 (IFFT)                   **
 **     calling seqence:                                                **
 **         #include <math.h>                                           **
 **         struct complex x[LENGTH_OF_DATA];                           **
 **         int m;                              LENGTH_OF_DATA=2**m     **
 **                                                                     **
 **           fft001( x, m, 0 )                 make bit tables         **
 **                ...                                                  **
 **           fft001( x, m, 1 )                 perform FFT             **
 **                ...                                                  **
 **           fft001( x, m,-1 )                 perform IFFT            **
 *************************************************************************/

#include <math.h>

#define NMAX 1024

static int inv[NMAX],*q1,*q2;
static struct complex w[NMAX],*r1,*r2;

fft001( x, m, sw )
struct complex x[];
int m, sw;
{
	int n,k,l,p,q,flag;
	int nob,lob,lhb,mm,ks,kf,i,j,jj;
	double pi=3.14159265358979323;
	double pn,t;
	struct complex tocomplex(), add_complex(), sub_complex(),
		       mul_complex(), t1,t2,wf;

	n=1<<m;
	if(sw==-1) {sw=1;flag=-1;}
	else flag=1;
	switch(sw)
	{
case 0:                                /* make bit inversed table */
	inv[0]=0;
	k=1<<(m-1); l=1;
	for( p=1; p<=m; p++,k=k>>1,l=l<<1 )
		for( q=0,q2=inv+q,q1=q2+l; q<l; q++,q1++,q2++ ) *q1=*q2+k;

				       /* make twiddle factor table */
	pn=(double)2*pi/(double)n;
	for( p=0,t=(double)0; p<n; p++,t+=pn )
		*(w+p)=tocomplex( cos(t), -sin(t) );
	break;

case 1:                                /* perform fft */
	nob=1; lob=n; lhb=lob>>1;
	for( mm=0; mm<m; mm++,nob<<=1,lob>>=1,lhb>>=1 )
	{
		for( i=0; i<nob; i++ )
		{
		  ks=i*lob; kf=inv[i<<1];
		  if(flag==-1) kf=(n-kf)%n; wf=*(w+kf);
		  for(j=0,r1=x+ks,r2=r1+lhb;j<lhb;j++,r1++,r2++)
		  {
		    t1=*r1;
		    t2=mul_complex( wf, *r2 );
		    *r1=add_complex( t1, t2 );
		    *r2=sub_complex( t1, t2 );
		  }
		}
	}

				       /* bit reverse operation */
	for( i=0; i<n; i++ )
	  if( i<(j=inv[i]) )
	  {
		  r1=x+i; r2=x+j;
		  t1=*r1; *r1=*r2; *r2=t1;
	  }
	break;

default:
	printf("\n FFT-illegal switch %d",sw);
	break;
	}
}

/** complex operations **/

struct complex tocomplex( x1, x2 )
double x1, x2;
{
	struct complex x3;
	x3.x=x1; x3.y=x2;
	return( x3 );
}

struct complex add_complex( x1, x2 )
struct complex x1, x2;
{
	struct complex x3;
	x3.x=x1.x+x2.x; x3.y=x1.y+x2.y;
	return( x3 );
}

struct complex sub_complex( x1, x2 )
struct complex x1, x2;
{
	struct complex x3;
	x3.x=x1.x-x2.x; x3.y=x1.y-x2.y;
	return( x3 );
}

struct complex mul_complex( x1, x2 )
struct complex x1, x2;
{
	struct complex x3;
	x3.x=x1.x*x2.x-x1.y*x2.y;
	x3.y=x1.x*x2.y+x1.y*x2.x;
	return( x3 );
}






 fft2d1(x,m,ln)
 struct complex x[];
 int m,ln;
 {
	int i,j,k,n;
	struct complex w[1024];
	extern int fft001();
	n=1;
	for(k=0;k<m;n*=2,k++);
	fft001(w,m,0);
	for(j=0;j<n;j++)
	{
		for(i=0;i<n;i++) w[i]=x[i*ln+j];
		fft001(w,m,1);
		for(i=0;i<n;i++) x[i*ln+j]=w[i];
	}
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++) w[j]=x[i*ln+j];
		fft001(w,m,1);
		for(j=0;j<n;j++) x[i*ln+j]=w[j];
	}
	return(0);
 }