fft.h

利用最佳一致逼近设计数字滤波器的程序

#include"math.h"
#include "stdlib.h"
#include "malloc.h"
#define PI 2*asin(1)

/*基２时间抽选ＦＦＴ算法Ｔurbo Ｃ源程序*/

void swap(float *s1,float *s2)    /*交换两数据值的函数*/
 {float temp;
  temp=(*s1);
  (*s1)=(*s2);
  (*s2)=temp;
 }

void fft(float *xreal,float *yimag,int numdat,int flag)

     /*FFT和逆FFT函数:xreal和yimag分别为存放输入序列的实部和虚部数组的指针;
       flag=0,正FFT;flag非0,逆FFT;numdat为序列点数*/

{int maxpower,arg,cntr,pnt0,pnt1,i;
 int j,a,b,k;
 float sign,prodreal,prodimag;double harm;
 float *cosary;
 float *sinary;
 cosary=(float *)calloc(numdat+1,sizeof(float));
 if(cosary==NULL)
   exit(1);
 sinary=(float *)calloc(numdat+1,sizeof(float));
 if(sinary==NULL)
   exit(1);

 if(flag!=0){
    sign=1.0;
    
    for(i=0;i<=numdat-1;++i){
       xreal[i]=xreal[i]/numdat;
       yimag[i]=yimag[i]/numdat;
       }
    }
 else{
      sign=-1.0;
    }

    j=0;                              /*二进制反向重排*/
    for(i=0;i<=numdat-2;i++){
      if(i<j){
      swap(&xreal[i],&xreal[j]);
      swap(&yimag[i],&yimag[j]);
      }
      k=numdat/2;
      while(k<=j){
      j=j-k;
      k=k/2;
      }
    j=j+k;
  }
maxpower=0;      /*求N=(2的M次方)中的M:M=maxpower*/
i=numdat;
while(i!=1){
    maxpower=maxpower+1;
    i=i/2;
}
harm=2.*PI/numdat;                  /***原程序为:harm=2*M-PI/numdat***/
for(i=0;i<numdat;i++)
  {sinary[i]=(float)(sign*sin(harm*i));
   cosary[i]=(float)cos(harm*i);
   }
   a=2;
   b=1;
   for(cntr=1;cntr<=maxpower;++cntr){    /*此循环进行M轮蝶型运算*/
      pnt0=numdat/a;
      pnt1=0;
      for(k=0;k<=b-1;++k){        /*b为每一轮蝶型运算中具有相同旋转因子的组数*/
	i=k;
	while(i<numdat){              /*while循环进行蝶型运算*/
	  arg=i+b;
	  if(k==0){
	    prodreal=xreal[arg];
	    prodimag=yimag[arg];
	  }
	  else{
	   prodreal=xreal[arg]*cosary[pnt1]-yimag[arg]*sinary[pnt1];
	   prodimag=xreal[arg]*sinary[pnt1]+yimag[arg]*cosary[pnt1];
	   }
	   xreal[arg]=xreal[i]-prodreal;
	   yimag[arg]=yimag[i]-prodimag;
	   xreal[i]=xreal[i]+prodreal;
	   yimag[i]=yimag[i]+prodimag;
	   i=i+a;
	   }
	   pnt1=pnt1+pnt0;
	   }
	   a=2*a;
	   b=2*b;
	   }
	   free(cosary);
	   free(sinary);
	   }