fft.c

可以配置点数的fft算法源代码

#include "math.h"

int n,k,l,il;                      //n输入的点数;k满足n=2^k（2的k次方）;若l＝0，表示要求本函数计算傅立叶变换；若l＝1，表示要求本函数计算逆傅立叶变换。。
                                   //若il＝0，表示不要求本函数计算傅立叶变换或逆傅立叶变换的模与幅角；若il＝1，表示要求本函数计算傅立叶变换或逆傅立叶变换的模与幅角。
double pr[256],pi[256],fr[256],fi[256];        //pr长度为n,存放n个采样输入的实部，返回时存放离散傅立叶变换的模
                                   //pi长度为n,存放n个采样输入的虚部，返回时存放离散傅立叶变换的幅角
                                     //fr长度为n,返回傅立叶变换的实部
                                   //fi长度为n,返回傅立叶变换的虚部



void kbfft(double *pr,double *pi,int n,int k,double *fr,double *fi,int l,int il)   //傅立叶变换,令l=0

{
      int it,m,is,i,j,nv,t;
      double p,q,s,vr,vi,poddr,poddi;
      for(it=0;it<=n-1;it++)
      {
         m=it;
         is=0;
         for(i=0;i<=k-1;i++)
         {
            j=m/2;
            is=2*is+(m-2*j);
            m=j;
         }
         fr[it]=pr[is];
		 fi[it]=pi[is];
      }
      pr[0]=1.0;
      pi[0]=0.0;
      p=6.283185306/(1.0*n);
      pr[1]=cos(p);
      pi[1]=-sin(p);
      if(l!=0) 
      pi[1]=-pi[1];
      for(i=2;i<n-1;i++)
      {
           p=pr[i-1]*pr[1];
           q=pi[i-1]*pi[1];
           s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
           pr[i]=p-q;
           pi[i]=s-p-q;
      }
      for(it=0;it<=n-2;it=it+2)
      {
           vr=fr[it];
	       vi=fi[it];
           fr[it]=vr+fr[it+1];
           fi[it]=vi+fi[it+1];
           fr[it+1]=vr-fr[it+1];
           fi[it+1]=vi-fi[it+1];
      } 
      m=n/2;
      nv=2;
      for(t=k-2;t>=0;t--)
      {
         m=m/2;
         nv=2*nv;
         for(it=0;it<=(m-1)*nv;it=it+nv)      
         for(j=0;j<=(nv/2)-1;j++)
         {
                p=pr[m*j]*fr[it+j+nv/2];
                q=pi[m*j]*fi[it+j+nv/2];
                s=pr[m*j]+pi[m*j];
                s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
                poddr=p-q;
                poddi=s-p-q;
                fr[it+j+nv/2]=fr[it+j]-poddr;
                fi[it+j+nv/2]=fi[it+j]-poddi;
                fr[it+j]=fr[it+j]+poddr;
                fi[it+j]=fi[it+j]+poddi;
         }
      }
      if(l!=0)
      for(i=0;i<=n-1;i++)
      {
            fr[i]=fr[i]/(1.0*n);
            fi[i]=fi[i]/(1.0*n);
      }
      if(il!=0)
      for(i=0;i<=n-1;i++)
      {
             pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
             if(fabs(fr[i])<0.000001*fabs(fi[i]))
             {
                 if((fi[i]*fr[i])>0) 
                     pi[i]=90.0;
                 else  
                     pi[i]=-90.0;
             }
             else
                 pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
 	      
         } 
      return;    
}
void main()
{ 
	int i;
	for(i=0;i<256;i++) 
	{
		pr[i]=cos(i*(6.283185306/256)); 
		pi[i]=0;
	}
	
	 kbfft(pr,pi,256,8,fr,fi,0,1);
	 while(1);
}