fft.c

fft快速变换

#include"stdio.h"
#include"alloc.h"
#include"stdlib.h"
#include"math.h"
#include"graphics.h"
void reSort(double array[],int length);
void main(){
   int N,                           /*采样点数目*/
       F,                           /*正弦信号频率*/
       layers,                      /*fft运算的层数*/
       length,tag=0,                /*采样点数目、是否追加0标记*/
       i,j,c,LE,LEI,IP,I,t,         /*int型中间变量*/
       init=DETECT,mode=0,          /*屏幕画图的参数*/
       step=15;                     /*画频谱的步长*/
   char ch;
   double *real,*image,*mod,        /*实部、虚部、模*/
         ModMax=0,                  /*最大模*/
         T,                         /*采样周期*/
         PI=3.141592653589793,      /*常量*/
         d,treal,timage,wreal,
         wimage,ureal,uimage;       /*doublu型中间变量*/
   initgraph(&init,&mode,0);  
   printf("Append zeros?(Y/N)\n");
   ch=getchar();
   printf("Please input the number of nodes:\n");
   scanf("%d",&N);
   length=N;
   printf("Please input frequency of the signal:\n");
   scanf("%d",&F);
   printf("please input T:\n");
   scanf("%lf",&T);
   if(ch=='y'||ch=='Y'){ 
	   length=N*2;
       tag=1;
   }
   real=(double*)malloc(length*sizeof(double));/*动态创建数组*/
   image=(double*)malloc(length*sizeof(double));
   mod=(double*)malloc(length*sizeof(double));
   d=T*F*2*PI;
   for(i=0;i<N;i++){
	    real[i]=sin(i*d);
	}
	if(tag)
		for(j=N;j<2*N;j++)
			real[j]=0;  
    for(c=0;c<length;c++){
	   image[c]=0;
	}
	reSort(real,length);
/****************************************************/
/*      蝶型算法                                     */
/****************************************************/
   layers=(int)(log(length)/log(2));
   for(i=1;i<=layers;i++){
       LE=(int)pow(2,i);
       LEI=LE/2;
       ureal=1,uimage=0;
       wreal=cos(PI/LEI),wimage=-sin(PI/LEI);
	for(j=1;j<=LEI;j++){
          for(I=j;I<=length;I=I+LE){
	   IP=I+LEI;			
	   treal=real[IP-1]*ureal-image[IP-1]*uimage;
           timage=real[IP-1]*uimage+image[IP-1]*ureal;
	   real[IP-1]=real[I-1]-treal;
	   image[IP-1]=image[I-1]-timage;
           real[I-1]=real[I-1]+treal;
           image[I-1]=image[I-1]+timage;
	   }
           treal=ureal;timage=uimage;
	   ureal=treal*wreal-timage*wimage;
     	   uimage=treal*wimage+timage*wreal;
	   }
   }
   printf("The anwser is:\n");
 /*  for(t=0;t<length;t++){*/
 /*  printf("%f+j(%f)\n",real[t],image[t]);}*/
/*****************************************************************/
/*                模归一化                               */
/*****************************************************************/
   	for(c=0;c<length;c++){
         mod[c]=sqrt(real[c]*real[c]+image[c]*image[c]);
	     if(ModMax<mod[c])
			 ModMax=mod[c];
	}  
    printf("the max value of mod is %lf\n\n",ModMax);
    for(c=0;c<length;c++)
	   mod[c]=mod[c]/ModMax;
     printf("\t\t\tthe graphic of signal's frequency");
     moveto(50,200);
     lineto(50,380);
     moveto(30,360);
     lineto(600,360);
     setcolor(12);
     printf("\n\n\n\n\n\n\n\n\n\n\n\n");  /*在坐标轴上标k值*/
     if(length==32)printf("\t1          7              15             23             31");
     else if(length==64)printf("       1     7       15      23      31      39      47      55     63");
     if(length==64)step=8;            /*当length为64时把画线步长设小为8*/
     for(c=0;c<length;c++){
       moveto(50+c*step,360);       /*画频谱线*/
       lineto(50+c*step,360-mod[c]*130); 
       
     }
     getch();
     closegraph();
}
/*****************************************************************/
/*      排序子函数reSort                                    */
/*****************************************************************/
void reSort(double *array,int length){
int NV2=length/2;
int NM1=length-1;
int j=1,i=0;
double t;
int k;
while(i<=NM1){
   i++;
  if(i<j)
  {  
     i--;j--;
     t=array[j];
	 array[j]=array[i];
	 array[i]=t;
	 i++;j++;
  }
      k=NV2;
	while(k<j){
	   j=j-k;
	   k=k/2;
	}
    j=j+k;if(j==length) return;/*当j=length时跳出，防止死循环*/
}
}