fft.cpp

FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。按照蝶形运算的构成不同可分为基2、基4、基8以

#include<stdio.h>
#include<string.h>
#include<math.h>
#define PI 3.1415926 
#define MAX 100

void main()
{int i=0;
 int *m=NULL;
 int *N=NULL;
  

 float x[MAX],y[MAX];                        /* y第一步帮助二进制翻转，然后第二步保存虚部*/
 void zero_add(float *,int *,int * );                      /*补零*/
 void FFT(float *,float *,int *,int *);
 void flip(float *,int *,int *,float *);/*翻转*/
m=new int;                               /*m级数*/
N=new int;                               /*点的个数*/

  
  printf("输入数据:\n");
 scanf("%f",&x[i]);
  while(getchar()!='\n')
  {i++;
   scanf("%f",&x[i]);
   if(i>MAX)
	   printf("超过最大个数点");
     }
   *N=i+1;
  zero_add( x, N,m);/*补零*/
flip(x,N,m,y);           /*翻转 重新排序*/ 
for(i=0;i<*N;i++)
{printf("%f\n",x[i]);
}

 FFT(x,y,m,N);
for(i=0;i<*N;i++)
{printf("X[%d]=(%f)+(%f)i\n",i,x[i],y[i]);
}


}  
 
void zero_add(float *x,int *N,int *m) /*补零*/
{int i,t;
 *m=(int)(log(*N)/log(2));
	t=(int)(pow(2,*m));
if(*N!=t)
  {(*m)++;
  }

  t=(int)(pow(2,*m));
  for(i=*N;i<t;i++)
  {x[i]=0;
  }
  *N=t;
}
void flip(float*x,int *N ,int *m,float *y)  /*排序*/
{int i,j;
int k;
float t;
for(i=0;i<*N;i++)
{j=i;
 k=*m;
y[i]=0;
	while(j!=0)
	{ k--;
		y[i]=y[i]+pow(2,k)*(j%2);
		j=j/2;
	}
}
  for(i=0;i<*N;i++)
  {  if(y[i]!=0)
  {  t=x[i];
      k=(int)(y[i]);
     x[i]=x[k];
     x[k]=t;
     y[i]=0;
     y[k]=0;
  }
  }
}
void FFT(float *x,float *y,int *m,int *N)
{int i,r,h,k,j;
float xp,xq,xp1,xq1,w;
 for(i=0;i<*m;i++)     /*实现m级运算*/
 {h=(int)(pow(2,i));   /*计算步长*/
	 for(r=0;r<pow(2,i);r++)/*r相同的分一组计算*/
	 { w=2*PI*r/pow(2,i+1); /*旋转因子的值*/
		 for(j=0;j<(*N/pow(2,i+1));j++)/*r相同的值的组的计算*/
	 {k=(int)( r+j*pow(2,i+1));          /*每一组成员的下标*/
		 
		 xp=x[k]+cos(w)*x[k+h]+sin(w)*y[k+h];/*蝶形计算*/
		 xp1=y[k]+cos(w)*y[k+h]-sin(w)*x[k+h];
		 xq=x[k]-cos(w)*x[k+h]-sin(w)*y[k+h];
		 xq1=y[k]-cos(w)*y[k+h]+sin(w)*x[k+h];
		 
		 x[k]=xp;
		 x[k+h]=xq;
		 y[k]=xp1;
		 y[k+h]=xq1;
			
			
	 }
	 }
 }
}