fft1.c

FFT实现

#include<stdio.h>
#include<math.h>
//#include<graphics.h>
#define PI 3.142592653589793
#define SWAP(x,y) {double temp=x; x=y; y=temp;}

#define CONSTANT_VALUE 16385

void FFT(double xxx[],int N,int FFTorIFFT)     /*FFT 1 IFFT -1*/
{
   double temp,tempr,tempi;
   double wtemp,wr,wpr,wpi,wi,theta;
   int n,mmax,m,q,istep,p;
   n=N<<1; q=1;
   for(p=1;p<n;p+=2)
   {
      if(q>p)
      {
         SWAP(xxx[q],xxx[p]);
         SWAP(xxx[q+1],xxx[p+1]);
      }
      m=n>>1;
      while(m>=2&&q>m)
      {
         q-=m;
         m>>=1;
      }
   q+=m;
   }
   mmax=2;
   while(n>mmax)
   {
      istep=2*mmax;
      theta=-2.0*PI/(FFTorIFFT*mmax);
      wtemp=sin((double)(0.5*theta));
      wpr=-2.0*wtemp*wtemp;
      wpi=sin((double)theta);
      for(wr=1.0,wi=0.0,m=1;m<mmax;m+=2)
      {
         for(p=m;p<=n;p+=istep)
         {
            q=p+mmax;
            tempr=wr*xxx[q]-wi*xxx[q+1];
            tempi=wr*xxx[q+1]+wi*xxx[q];
            xxx[q]=xxx[p]-tempr;
            xxx[q+1]=xxx[p+1]-tempi;
            xxx[p]+=tempr;
            xxx[p+1]+=tempi;
         }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
      }
      mmax=istep;
   }
}


int main()
{
   double xxx[CONSTANT_VALUE]={0};/*[257]*/
   int i,j;

//   int graphdriver=VGA;
//   int graphmode=VGAHI;
//   initgraph(&graphdriver,&graphmode,"");
   for(i=0;i<8192/*(CONSTANT_VALUE-1)*/;i+=2)
   {
      xxx[i+1]=0/*sin(0.1*i)/(0.1*i)*/;       /*xxx[1,3,5,7...] is real part */
      xxx[i+2]=0.0;     /*xxx[2,4,6,8...] is image part         */
   }
   for(i=10;i<22;i+=2)
   {
      xxx[i+1]=10;
   }
   for(i=0;i<32;i++)
      printf("%f\n",xxx[i+1]);
   printf("\n");
   getchar();
   xxx[0]=0;

   FFT(xxx,(CONSTANT_VALUE-1)/2,1);
   for(i=1;i<=CONSTANT_VALUE-1;i+=2)
   {
      //printf("%f\n",xxx[i]);
      //printf("%f\n",xxx[i+1]);
      j=(i+1)/2;
      xxx[j]=sqrt(pow(xxx[i],2)+pow(xxx[i+1],2));
      printf("%f\n",xxx[j]);
      getchar();
   }
   getchar();
//   cleardevice();
//   for(i=1;i<=(CONSTANT_VALUE-1)/4;i++)
//   {
//      putpixel(i,-(xxx[i])*30+300,1);
//   }
//   getch();
   /*closegragh();*/
}