fft.txt

关于傅立叶变换的C语言代码, 可用于信息处理

  函数名：_FFT   功能：FFT变换   用法：void _FFT(double *fr, double *fi, int n, int flag)   参数：fr：采样点的实数表，fi：采样点的虚数表，n：采样点个数，flag：flag=0表示求 Fourier变换，flag=1表示求逆Fourier变换   返回：无 


C的FFT:       void _FFT(double *fr, double *fi, int n, int flag) {  int mp,arg,q,cntr,p1,p2;  int i,j,a,b,k,m;  double sign,pr,pi,harm,x,y,t;  double *ca,*sa;    ca=(double *)calloc(n,sizeof(double));  sa=(double *)calloc(n,sizeof(double));  if(ca==NULL||sa==NULL) Error("calloc error in _FFT!");  j=0;  if(flag!=0)  {   sign=1.0;   for(i=0;i<=n-1;i++)   {    fr[i]=fr[i]/n;    fi[i]=fi[i]/n;   }  }  else   sign=-1.0;  for(i=0;i<=n-2;i++)  {   if(i<j)   {    t=fr[i];    fr[i]=fr[j];    fr[j]=t;    t=fi[i];    fi[i]=fi[j];    fi[j]=t;   }   k=n/2;   while(k<=j)   {    j-=k;    k/=2;   }   j+=k;  }  mp=0;  i=n;  while(i!=1)  {   mp+=1;   i/=2;  }  harm=2*M_PI/n;  for(i=0;i<=n-1;i++)  {   sa[i]=sign*sin(harm*i);   ca[i]=cos(harm*i);  }  a=2;  b=1;  for(cntr=1;cntr<=mp;cntr++)  {   p1=n/a;   p2=0;   for(k=0;k<=b-1;k++)   {    i=k;    while(i<n)    {     arg=i+b;     if(k==0)     {      pr=fr[arg];      pi=fi[arg];     }     else     {      pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];      pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];     }     fr[arg]=fr[i]-pr;     fi[arg]=fi[i]-pi;     fr[i]+=pr;     fi[i]+=pi;     i+=a;    }    p2+=p1;   }   a*=2;   b*=2;  }  free(ca);  free(sa); }