fft.c

FFT源代码

//========================================
//函数名：_FFT   
//功能：FFT变换    
//用法：void   _FFT(double   *fr,   double   *fi,   int   n,   int   flag)      
//参数：fr：采样点的实数表，fi：采样点的虚数表，n：采样点个数，
//    flag：flag=0表示求Fourier变换，flag=1表示求逆Fourier变换   
//返回：无   
//========================================    
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);   
}