FFT的C语言算法实现
程序如下:
  /************FFT***********/
  #include  <stdio.h>
  #include  <math.h>
  #include  <stdlib.h>
  #define  N  1000
  typedef  struct
  {
  double  real;
  double  img;
  }complex;
  void  fft(); /*快速傅里叶变换*/
  void  ifft(); /*快速傅里叶逆变换*/
  void  initW();
  void  change();
  void  add(complex  ,complex  ,complex  *);  /*复数加法*/ 
  void  mul(complex  ,complex  ,complex  *);  /*复数乘法*/ 
  void  sub(complex  ,complex  ,complex  *);  /*复数减法*/ 
  void  divi(complex  ,complex  ,complex  *);/*复数除法*/ 
  void  output(); /*输出结果*/
  complex  x[N],  *W;/*输出序列的值*/
  int  size_x=0;/*输入序列的长度,只限2的N次方*/
  double  PI;
  int  main()
  {
  int  i,method;
  system("cls");
  PI=atan(1)*4;
  printf("Please  input  the  size  of  x:\n");
  /*输入序列的长度*/
  scanf("%d",&size_x);
  printf("Please  input  the  data  in  x[N]:(such as:5 6)\n");
  /*输入序列对应的值*/
  for(i=0;i<size_x;i++)
  scanf("%lf %lf",&x[i].real,&x[i].img);
  initW(); 
  /*选择FFT或逆FFT运算*/
  printf("Use  FFT(0)  or  IFFT(1)?\n"); 
  scanf("%d",&method); 
  if(method==0) 
  fft(); 
  else 
  ifft(); 
  output(); 
  return  0; 
  } 
   
  /*进行基-2 FFT运算*/
  void  fft() 
  { 
  int  i=0,j=0,k=0,l=0; 
  complex  up,down,product; 
  change(); 
  for(i=0;i<  log(size_x)/log(2)  ;i++)
  {     
  l=1<<i; 
  for(j=0;j<size_x;j+=  2*l  )
  {                         
  for(k=0;k<l;k++)
  {               
  mul(x[j+k+l],W[size_x*k/2/l],&product); 
  add(x[j+k],product,&up); 
  sub(x[j+k],product,&down); 
  x[j+k]=up; 
  x[j+k+l]=down; 
  } 
  } 
  } 
  } 
   
  void  ifft() 
  { 
  int  i=0,j=0,k=0,l=size_x; 
  complex  up,down; 
  for(i=0;i<  (int)(  log(size_x)/log(2)  );i++) /*蝶形运算*/ 
  {   
  l/=2; 
  for(j=0;j<size_x;j+=  2*l  ) 
  {                       
  for(k=0;k<l;k++)
printf怎么实现的
  {             
  add(x[j+k],x[j+k+l],&up); 
  up.real/=2;up.img/=2; 
  sub(x[j+k],x[j+k+l],&down); 
  al/=2;down.img/=2; 
  divi(down,W[size_x*k/2/l],&down); 
  x[j+k]=up; 
  x[j+k+l]=down; 
  } 
  } 
  } 
  change(); 
  } 
   
  void  initW() 
  { 
  int  i; 
  W=(complex  *)malloc(sizeof(complex)  *  size_x); 
  for(i=0;i<size_x;i++) 
  { 
  W[i].real=cos(2*PI/size_x*i); 
  W[i].img=-1*sin(2*PI/size_x*i); 
  } 
  } 
   
  void  change() 
  { 
  complex  temp; 
  unsigned  short  i=0,j=0,k=0; 
  double  t; 
  for(i=0;i<size_x;i++) 
  { 
  k=i;j=0; 
  t=(log(size_x)/log(2)); 
  while(  (t--)>0  ) 
  { 
  j=j<<1; 
  j|=(k  &  1); 
  k=k>>1; 

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。