fft_gai.c

标准C语言的高效的fft算法

#include	<stdio.h>
#include	<math.h>
#include	<time.h>
#include	<malloc.h>

void Changeorder_1(float *R, int N);					//倒叙
void PowerM(int N);										//求幂			
void CaculateWn(int N);									//求旋转因子(半长FFT)
void CaculateWn1(int N);								//求旋转因子(对偶组合)			
void Displaydata( float *R, float *I, int N,int flag);	//显示输入输出
void FFT_r(float *R, float *I, int N,int M);			//实序列FFT 

int i=0,N=0,M=0,CT=0;									//循环变量i，序列长,幂，循环次数
clock_t before,after;									//初始时间，结束时间

#define MAX 2048										
const double PI=3.14159265;								

float R2[MAX];		//测试循环用实部
float I2[MAX];		//测试循环用虚部
float R[MAX];		//输入/输出序列的实部
float I[MAX];		//输入/输出序列的虚部
float WR[MAX];		//旋因子实部
float WI[MAX];		//旋因子虚部
float WR1[MAX];		//半长序列用旋因子实部
float WI1[MAX];
float R1[MAX];		//R的偶序列作为半长序列，记为x1
float I1[MAX];		//I的偶序列作为半长序列，记为x2
float xr[MAX];		//x1 FFT后的实部
float xi[MAX];		//x1 FFT后的虚部
float yr[MAX];		//
float yi1[MAX];		//
float or[MAX];		//半长序列实部
float oi[MAX];		//半长序列虚部
int   a[MAX]={0,1};	//记录倒序数下标



//////////////////////////////////////////////////////////////////////////MAIN
void main()
{
	
	int c=0;	//循环测试变量
	
	/////////////////////////////////////////////////////////////////读取信号
	FILE *fp  = NULL;
	char filename[20];
	float in1=0;
	float in2=0;
	
loop:printf("input file name:");
	 scanf("%s",filename);
	 fp=fopen(filename,"r");
	 if (fp==NULL)
	 {
		 printf("file is not exist! please input again:\n");
		 goto loop;
	 }
	 else
	 {
		 fscanf(fp,"%d\n",&N);
		 for (i=0;i<N;i++)
		 {
			 fscanf(fp,"%f\n",&in1);
			 R[i]=in1;
			 I[i]=0;
			 in1=0;
		 }
	 }
	 fclose(fp);	
	 
	 Displaydata(R,I,N,0);
	 printf("input cycle times:");
	 scanf("%d",&CT);
	 
	 ////////////////////////////////////////////////////////////////////////FFT运算
	 
	 before=clock();		//初始时间
	 PowerM(N);			    //求幂
	 Changeorder_1(R,N);	//求倒序数
	 
	 CaculateWn1(N);		//初始旋转因子表1
	 CaculateWn(N);		    //初始旋转因子表2
	 
	 for (c=CT;c>0;c--)
	 {
		 for (i=0;i<N;i++)	//保存并调整倒叙后的输入以循环
		 {	 
			 *(R2+i)=*(R+*(a+i));
			 *(I2+i)=*(I+*(a+i));		 
		 }
		 		 
		 FFT_r(R2,I2,N,M);	//求解FFT		 
	 }
	 
	 after=clock();		//计时结束
	 
	 ////////////////////////////////////////////////////////////////////输出结果
	 
	 Displaydata(R2,I2,N,1);	
	 printf("cycle:%d   time cost:%f\n",CT,(float)(after-before)/CLOCKS_PER_SEC);
	 
}


//////////////////////////////////////////////////////////////////////////

void PowerM(int N)
{
	M=0;
	while (N>1)
	{
		M++;
		N/=2;		
	}	
}

void Changeorder_1(float *R, int N)
{
	int i,j;
	
    for(i=2;i<N;i<<=1)
    {
        for(j=0;j<i;j++)
        {
			*(a+j)=*(a+j)<<1;
			*(a+j+i)=*(a+j)+1;			
        }
	}	
}


void CaculateWn(int N)
{	
	int n = N/4;
	for(i=0;i<n;i++){
		WI[i] =-(float)sin(PI/n*i);
		WR[i] = (float)cos(PI/n*i);
	}		
}


void CaculateWn1(int N)
{		
	int n = N/2;
	for(i=0;i<n;i++){
		WI1[i] =-(float)sin(PI/n*i);
		WR1[i] = (float)cos(PI/n*i);
	}
}


void FFT_r(float *R, float *I,int N, int M)
{	
	register int N1;
	int arg,cntr,pnt0,pnt1;
	int a,b,k;
	float pr,pi;//,tmp;
	
	register float temp1;
	register float temp2;
	
	int p;
	
    register float a1,b1;
	
	////////////分层计算初始化 
	a=8;
	b=4;
	
	N>>=1;
	M--;
	
	p=N/4;   
	
	
	
	///////////////////////////////////////////////////////////后倒叙再奇偶抽=先倒叙再前后分
    for (i=0;i<N;i++)
	{
		R1[i]=R[i];
		I1[i]=R[i+N];
	}
	
	////////////////////////////////////////////////////////////第一级蝶型运算 
	for(i=0;i<N;i=i+2)
	{		
		temp1=*(R1+i+1);
		*(R1+i+1)=*(R1+i)-temp1;
		*(R1+i)=*(R1+i)+temp1;
		
		temp2=*(I1+i+1);
		*(I1+i+1)=*(I1+i)-temp2;
		*(I1+i)=*(I1+i)+temp2;
	}
	
	/////////////////////////////////////////////////////////第二级蝶型运算 

    for(i=0;i<N;i=i+4)
    {	
		a1=WR[p];
		b1=WI[p];
        temp1=*(R1+i+2);
		*(R1+i+2)=*(R1+i)-temp1;
		*(R1+i)=*(R1+i)+temp1;
		
		temp2=*(I1+i+2);
		*(I1+i+2)=*(I1+i)-temp2;
		*(I1+i)=*(I1+i)+temp2;
		
        pr=*(R1+i+3)*a1-*(I1+i+3)*b1;
		pi=*(R1+i+3)*b1+*(I1+i+3)*a1;		
		
        *(R1+i+3)=*(R1+i+1)-pr;
        *(R1+i+1)=*(R1+i+1)+pr;
        *(I1+i+3)=*(I1+i+1)-pi;
        *(I1+i+1)=*(I1+i+1)+pi;	
    }
	
	
	///////////////////////////////////////////////////////////////三-M级正常运算 
	for(cntr=3;cntr<=M;++cntr){
		//pnt0=nfft/a;
		pnt0=N>>cntr;
		pnt1=0;
		for(k=0;k<=b-1;++k){
					
			i=k;
			while(i<N){  //以旋转因子为循环（与教材上实际一致）
				arg=i+b;    //b实际上是2^(cntr-1)的另一算法，pnt1是旋转因子指数
				//////////////////////////////////////////////////////////////////////////
				if(R1[arg]==0&&I1[arg]==0)
				{  if(R1[i]==0&&I1[i]==0)
				{i=i+a;continue;}
                else
                { R1[arg]=R1[i];I1[arg]=I1[i];i+=a;}
				}
				else
				{
					if (R1[i]==0 && I1[i]==0)
					{
						if(k==0){
							pr=R1[arg];
							pi=I1[arg];
						}
						
						else
						{
							register float wa=WR[pnt1];
							register float wb=WI[pnt1];
							pr=wa*R1[arg]-wb*I1[arg];
							pi=wa*I1[arg]+wb*R1[arg];
							
						}
						R1[arg]=-pr;
						I1[arg]=-pi;
						R1[i]=+pr;
						I1[i]=+pi;
						i=i+a;						
					}
					else
					{
						if(k==0){
							pr=R1[arg];
							pi=I1[arg];
						}						
						else
						{							
							register float wa=WR[pnt1];
							register float wb=WI[pnt1];
							pr=wa*R1[arg]-wb*I1[arg];
							pi=wa*I1[arg]+wb*R1[arg];							
						}
						R1[arg]=R1[i]-pr;
						I1[arg]=I1[i]-pi;
						R1[i]=R1[i]+pr;
						I1[i]=I1[i]+pi;
						i=i+a;						
					}				
				}   
			}
			//////////////////////////////////////////////////////////////////////			
			pnt1=pnt1+pnt0;
		}	
		a = a<<1;
		b = b<<1;
	}
	
	//////////////////////////////////////////////////////////////////////////
	
	
	for (i=0;i<N;i++)
	{
		or[i]=R1[i];
		oi[i]=I1[i];		
	}
	or[N]=or[0];
	oi[N]=oi[0];
	
	for (i=0;i<N;i++)
	{
		xr[i]=(or[i]+or[N-i])/2.0f;
		xi[i]=(oi[i]-oi[N-i])/2.0f;
		yr[i]=(oi[i]+oi[N-i])/2.0f;
		yi1[i]=-(or[i]-or[N-i])/2.0f;
	}
		
	xr[N]=xr[0];
	xi[N]=xi[0];
	yr[N]=yr[0];
	yi1[N]=yi1[0];
	
	for (i=0;i<=N;i++)		//全长序列前半（乘开）
	{
		R[i]=xr[i]+yr[i]*WR1[i]-yi1[i]*WI1[i];
		I[i]=xi[i]+yi1[i]*WR1[i]+yr[i]*WI1[i];		
	}
	N1=N<<1;
	for (i=N+1;i<N1;i++)	//全长后半=前半的共额
	{
		R[i]=R[N1-i];
		I[i]=-I[N1-i];
	}
}	




void Displaydata( float *R, float *I,int N,int flag)
{
	
	if (flag==0)
	{
		printf("\n ---------------input data is-------------------\n");
	}
	else
	{
		printf("\n ---------------the answer is-------------------\n");
	}	
	printf("the length of signal:%d\n",N);
	for ( i=0;i<N;i++)
	{
		printf("x[%d]=(%f,%f)\n",i,R[i],I[i]);		
	}
	
}