fft.c

基4的FFT算法程序

#include<stdio.h>
#include<math.h>
#define PI 3.141592653
void fft(x,y,l,f)
double *x,*y,f;
int l;
{
    int i,i0,i1,i2,i3,j,l1,n,ns,n1,n2,n3;
	int arg,arg1,arg2,arg3,*m;
	double s1,c1,s2,c2,s3,c3,sc1,sc2,sc3;
	double x0,y0,x1,y1,x2,y2,x3,y3,t;
	n=1;
	while(l!=0){n*=4;l--;}      /*n=4^l*/
	n1=ns=n/4;n2=2*n1;n3=3*n1;
	sc1=2.0*PI/(double)n;
	sc2=2.0*sc1;
	sc3=3.0*sc1;
	                /*work area for bit operation*/
	if((m=(int*)calioc(n,2))==NULL)
	{
	  printf("\n Out of memory error !!!\n");exit(0);
	}
	while(ns>=1)                /*main loop*/
	{
	    for(l1=0;l1<n;l1+=(4*ns))
		{
		   arg=m[l1]/4;arg1=arg+n1;
		   arg2=arg+n2;
		   arg3=arg+n3;
		   c1=cos(sc1*(double)arg);s1=sin(f*sc1*(double)arg);
		   c2=cos(sc2*(double)arg);s2=sin(f*sc2*(double)arg);
		   c3=cos(sc3*(double)arg);s3=sin(f*sc3*(double)arg);
		   for(i0=l1;i0<l1+ns;++i0)
		   {
		      i1=i0+ns;i2=i1+ns;i3=i2+ns;
			  x0=x[i0];y0=y[i0];
			  x1=x[i1]*c1-y[i1]*s1;
			  y1=y[i1]*c1+x[i1]*s1;
			  x2=x[i1]*c2-y[i1]*s2;
			  y2=y[i1]*c2+x[i1]*s2;
			  x3=x[i1]*c3-y[i1]*s3;
			  y3=y[i1]*c3+x[i1]*s3;
			  x[i0]=x0+x1+x2+x3;
			  y[i0]=y0+y1+y2+y3;
			  x[i2]=x0-x1+x2-x3;
			  y[i2]=y0-y1+y2-y3;
			  if(f<0.0)
			  {
			     x[i1]=x0+y1-x2-y3;
				 y[i1]=y0-x1-y2+x3;
				 x[i3]=x0-y1-x2+y3;
				 y[i3]=y0+x1-y2-x3;
			   }
			   else
			   { 
			     x[i1]=x0-y1-x2+y3;
				 y[i1]=y0+x1-y2-x3;
				 x[i3]=x0+y1-x2-y3;
				 y[i3]=y0-x1-y2+x3;
				}
				m[i0]=arg;m[i1]=arg1;
				m[i2]=arg2;m[i3]=arg3;
			}
		}
		ns/=4;
	}
	if(f<0.0)           /*judge fft or ifft */
	{
	   for(i=0;i<n;++i)
	   { 
	     x[i]/=(double)n;
		 y[i]/=(double)n;
	    }
	}
	for(i=0;i<n;++i)     /*bit reverse operation*/
	{
	   if((j=m[i])>i)
	   {
	     t=x[i];x[i]=x[j];x[j]=t;
		 t=y[i];y[i]=y[j];y[j]=t;
		}
	}
	free(m);
}