fft.c

C语言实现快速傅里叶算法

 #include "conio.h"
 #include "math.h"
 #include "stdio.h"

 #define N 64
 #define PI 3.1415926
 #define w0 (0.125*PI)
 #define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x;
 #define   Cequal(a,b) a.x=b.x;a.y=b.y;
 #define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y;
 #define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y;
 #define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n);
 struct comp
 {
   float x;
   float y;
 };
 void main()
 {
   int i,j,nu2,nm1,n,m,le,le1,k,ip,z;
   int flag,f,n1;
   struct comp a[N],t,t1,w,d;
   float a_ipx,m1;
   printf("\nThis program is about  FFT by DIF way.  ");
   printf("\nplease enter  N : ");
   scanf("%d",&n1);
   n=n1;
   m1=log(n1)/log(2);
   m=log(n1)/log(2);
   if (m!=m1)  n=pow(2,m+1);
   for(i=0;i<n;i++) {a[i].x=a[i].y=0.0;}
   printf("\n");
   for(i=0;i<n1;i++)
    {
     printf("\nplease enter  data(%d)_[Re]: ",i);
     scanf("%f",&a[i].x);
     printf("\nplease enter  data(%d)_[Im]: ",i);
     scanf("%f",&a[i].y);
    }
   for(z=0;z<=1;z++)
   {
    flag=-1;
    for (m=(log(n)/log(2));m>=1;m--)
      {
      le=pow(2,m);
      flag++;
      le1=le/2;
      for( j=0;j<le1;j++)
         { 
          for (i=j;i<=(n-1);i+=le)
             {  
               ip=i+le1;
               Cequal(t,a[i]);
               Cequal(t1,a[ip]);
               f=(int) (i*pow(2,flag))%n;
               Wn(w,f);  
               Cadd(a[i],t,t1);
               Csub(a[ip],t,t1);
               a_ipx=a[ip].x;
               if (z==1)
                 {  
                  w.y*=-1;
                 }
              a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
              a[ip].y=a_ipx*w.y+a[ip].y*w.x;
             }
          }
      }
    nu2=n/2;
    nm1=n-2;
    j=0;i=0;  
   while(i<=nm1)
      {  
       if (i<j)
         { 
          Cequal(d,a[j]);
          Cequal(a[j],a[i]);
          Cequal(a[i],d);
         }
       k=nu2;
       while(k<=j) 
         {
          j=j-k;k=k/2;
         }
        j=j+k;
        i=i+1;
       }
    if(z==0)
       {
        printf("\n the num` fft are:\n\n");
       }
    else
        printf("\n use ifft have the num:\n\n"  ) ;
    for(i=0;i<n;i++)
       if(z==0)
          { 
           printf("    %7.3f",a[i].x);
           if (a[i].y>=0)
           printf("  + %7.3f j \n",a[i].y);
           else
           printf("  - %7.3f j \n",fabs(a[i].y));
           a[i].y= -a[i].y;
           }
       else
          {    
           printf("   %7.3f",a[i].x/n);
           a[i].y=-a[i].y/n;
           if (a[i].y>=0)
           printf("  + %7.3f j \n",a[i].y);
           else
           printf("  - %7.3f j \n",fabs(a[i].y));
          }
    }
  printf("\n");
  getch();
 }