fft.c

C++下的乘法 阶乘 fft运算

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#define PI 3.1415926

float x[1024],y[1024],am[1024];
float c[1024][1024],s[1024][1024];

/*旋转因子表*/
  void MakeTable(int l)
  {
	  int i,r,N;
	  int m1,m2;

	  N = (int)(1<<l);

	  for(i=0;i<l;i++)
	  {
          m1=1<<i;
          m2=2*m1;
          for(r=0;r<m1;r++)
		  {
	    	  c[r][m2] = cos(2*PI*r/m2);
              s[r][m2] = sin(2*PI*r/m2);
		  }
	  }
  }
 
 /*位反转*/
  int BitReverse(int src, int size)
  {
	  int tmp = src;
	  int des = 0;
	  int i;
	  for(i=size-1; i>= 0; i--)
	  {
		  des = ((tmp&0x1)<<i)|des;
		  tmp = tmp >> 1;
	  }
	  return des;
  }

 /*倒序重排*/
  void re_bit(int N, int l)
  {
     int j, M;
     float a;
	 M = N/2;
     for(j=1;j<M;j++)
     {
		 a=x[j];
	     x[j]=x[BitReverse(j,l)];
	     x[BitReverse(j,l)]=a;
         a=y[j];
	     y[j]=y[BitReverse(j,l)];
	     y[BitReverse(j,l)]=a;		
	 }
  }

/*蝶形运算*/
  void butterfly(int l, int flag)
  {
      int i,j,r,m1,m2,m3,m4;
      int k1,k2,k3;
      int N;
      float u,v,s_PI;
	  int Sign = (flag==1)?(1):(-1);
	  
      N=1<<l;
      for(i=0;i<l;i++)
	  {
          m1=1<<i;
          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*PI*r/m2)+y[k2]*sin(2*s_PI*r/m2);//不查表
                  //v=y[k2]*cos(2*PI*r/m2)-x[k2]*sin(2*s_PI*r/m2);

				  u=x[k2]*c[r][m2]+y[k2]*s[r][m2]*Sign;//查表,加快运算速度
                  v=y[k2]*c[r][m2]-x[k2]*s[r][m2]*Sign;//Sign在fft变换时为+1,在ifft时为-1,如果内存足够大,可以再做一张逆变换表格
				  x[k2]=x[k1]-u;
                  y[k2]=y[k1]-v;
                  x[k1]=x[k1]+u;
                  y[k1]=y[k1]+v;
			  }
		  }
	  }
	  if(flag!=1)
		  for(k3=0; k3<N; k3++)
		  {
			  x[k3] /= N;
			  y[k3] /= N;
		  }
  }

  /*fft 变换*/
  void fft(int l)
  {
      int N;
	  int flag = 1;
      N=1<<l;

	  re_bit(N,l);
      butterfly(l,flag);
 }

  /*fft反变换*/
  void ifft(int l)
  {
	  int N;
	  int flag = 0;
      N=1<<l;

	  re_bit(N,l);
	  butterfly(l,flag);	  
  }

  int main()
  {
	  int i,l,N;	  
      printf("input l(l<=10)=");
      scanf("%d",&l);

	  N = 1<<l;

	  //输入数据，当然也可以用别的方法输入，这几个只是测试用
	  x[0]=1;y[0]=5;
	  x[1]=2;y[1]=6;
	  x[2]=3;y[2]=7;
	  x[3]=4;y[3]=8;
      
	  MakeTable(l);

	  fft(l);

      printf("\nfft Result:\n");
      for(i=0;i<N;i++)
		  printf("%6.4f%s%6.4fi\n",x[i],(y[i]>=0)?"+":"-",(y[i]>=0)?y[i]:-y[i]);
      printf("\n");

	  ifft(l);

	  printf("\nifft Result:\n");
      for(i=0;i<N;i++)
		  printf("%6.4f%s%6.4fi\n",x[i],(y[i]>=0)?"+":"-",(y[i]>=0)?y[i]:-y[i]);
      printf("\n");
  }