1024fft.cpp

一个基-2对时间抽取的1024点的快速傅立叶C语言算法,位倒序输入

 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 #include<process.h>
 #include<complex.h>
 #define p 3.1416
 float x[1024],y[1024],am[1024];
 void main()
 {
 int l,N,i;
 float a;
 void put_in(int N);
 void re_bit(int N);
 void butterfly(int l);
 printf("input l(1<=10)=");
 scanf("%d",&l);
 N=int(pow(2,1));
 put_in(N);
 re_bit(N);
 butterfly(1);
 for(i=0;i<N;i++)
   printf("%f+j%f  ",x[i],y[i]);
   printf("\n\n  ");
   getchar();
   for(i=0;i<N;i++)
   {
   am[i]=sqrt(x[i]*x[i]+y[i]*y[i]);
   printf("%f  ",am[i]);
   }
   printf("\n\n  ");
   getchar();
   }
   void butterfly(int l)
   {int i,j,r,m1,m2,m3,m4,N,k1,k2;
   float u,v;
   N=int(pow(2,1));
   for(i=1;i<l+1;i++)
   {m1=int(pow(2,i-1));
   m2=2*m1;
   m3=N/m2;
   for(j=0;j<m3;j++)
   {m4=j*m2;
   for(r=0;r<m1;r++)
   {
   k1=m4+r;
   k2=k1+m1;
   u=x[k2]*cos(2*p*r/m2)+y[k2]*sin(2*p*r/m2);
   v=y[k2]*cos(2*p*r/m2)-x[k2]*sin(2*p*r/m2);
   x[k2]=x[k1]-u;
   y[k2]=y[k1]-v;
   x[k1]=x[k1]+u;
   y[k1]=y[k1]+v;
       }
      }
     }
    }
   void re_bit(int N)
   {
    int i,j,M,s;
    float a;
    M=N/2;i=0;
    for(j=1;j<N;j++)
    {
     s=M;
     while(i>=s)
     {i=i-s;s=s/2;}
     i=i+s;
     if(j<i)
      {
      a=x[j];x[j]=x[i];x[i]=a;
      a=y[j];y[i]=y[i];y[i]=a;
	 }
	}
      }
      void put_in(int N)
      {
      int i;
      for(i=0;i<N;i++)
      {x[i]=cos(p*i*i/N);
      y[i]=sin(p*i*i/N);
	}
      }