time_wn.cpp

用C语言实现的对基—2FFT算法的改进

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

#define MAX 2048
#define PI 3.14159

//const int MAX=2048;
//const double PI=3.1415926;


void Changeorder(float *R, float *I, int N);			
void PowerM(int N);												
void CaculateWn(int N);												
void FFT(float *R, float *I, int N,int M);				
void FFT_r(float *R, float *I, int N,int M);		
void Displaydata( float *R, float *I, int N,int flag);	


int i=0,N=0,M=0,CT=0;			
clock_t before,after;


float R[MAX];	
float I[MAX];		
float WR[MAX];	
float WI[MAX];	
float WR1[MAX];	
float WI1[MAX];
float R1[MAX];    
float I1[MAX];	
float xr[MAX];	
float xi[MAX];	//
float yr[MAX];	//
float yi1[MAX];	//
float or[MAX];	//
float oi[MAX];	//



//////////////////////////////////////////////////////////////////////////MAIN
void main()
{
	FILE *fp = NULL;
	char filename[20];
	float in1=0;
//    float FR[MAX];
//    float FI[MAX];
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);
//			FR[i]=in1;
//			FI[i]=0;
           R[i]=in1;
			I[i]=0;
		}
	}
	fclose(fp);	
	
//	for(i=0;i<N;i++)
//	{
//	   R[i]=FR[i]*1000000;
 //   }
    
	
	
//	Displaydata(FR,FI,N,0);
	Displaydata(R,I,N,0);

	printf("input cycle times:");
	scanf("%d",&CT);
	before=clock();
	PowerM(N);	
	Changeorder(R,I,N);	
	CaculateWn(N/2);
	FFT_r(R,I,N,M);
	
//	for(i=0;i<N;i++)
//	{
//	   FR[i]=R[i]/1000000;
//	   FI[i]=I[i]/1000000;
//    }   	
	
	after=clock();
//	Displaydata(FR,FI,N,1);
    Displaydata(R,I,N,0);
	printf("cycle:%d   time cost:%f\n",CT,(float)(after-before)/CLOCKS_PER_SEC);
	goto loop;

/*	free(R);
	free(I);
	free(WR);
	free(WI);
	free(WR1);
	free(WI1);
	free(R1);
	free(I1);
	free(xr);
	free(xi);
	free(yr);
	free(yi1);
	free(or);
	free(oi);	
*/
}

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

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

void Changeorder(float *R, float *I, int N)
{
	int LH=0;		//最高权值
	int	N1=0;
	int	j=0;		//当前倒叙数的10进制数值
	int k=0;
	float tr,ti;
	LH=N/2;
	j=LH;
	N1=N-2;
	for (i=1;i<N1+1;i++)
	{
		if (i<j)
		{
			tr=R[i];
			R[i]=R[j];
			R[j]=tr;

	//		ti=I[i];
	//		I[i]=I[j];
	//		I[j]=ti;
		}
		k=LH;
	
		while(j>=k)
		{
			j=j-k;
			k=(int)(k/2);
		}
		j=j+k;
	}
 /*
	int a[MAX];
//    int *a=(int*)malloc(sizeof(int)*MAX);

	int i,j;
	a[0]=0;
    a[1]=1;
	
    for(i=2;i<N;i<<=1)
    {
        for(j=0;j<i;j++)
        {
            a[j]=2*a[j];
            a[j+i]=a[j]+1;			
        }
	}
	for (j=0;j<N;j++)
	{
		I[j]=R[a[j]];
	
	}
	for (j=0;j<N;j++)
	{
		R[j]=I[j];
	}
//	free(a);
*/
}

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

void FFT_r(float *R, float *I,int N, int M)
{
	int c=0;
	M--;
	N>>=1;


	////////////////////////////////////////////////////////////////////////先倒叙再奇偶抽=先倒叙再前后分
	for (i=0;i<=N;i++)
	{
		WR1[i]=cos(PI*i/N);
		WI1[i]=-sin(PI*i/N);
	//    WR1[i]=WR[i];
	//	WI1[i]=WI[i];

		R1[i]=R[i];
		I1[i]=R[i+N];
	}	

	
//	before=clock();
	for (c=CT;c>0;c--)
	{
		for (i=0;i<N;i++)///保存输入以循环
		{
			R[i]=R1[i];
			I[i]=I1[i];
		}
		FFT(R,I,N,M);//调用求FFT
						
		for (i=0;i<N;i++)
		{
			or[i]=R[i];
			oi[i]=I[i];
		}
		or[N]=or[0];
		oi[N]=oi[0];

		for (i=0;i<N;i++)
		{
			xr[i]=(or[i]+or[N-i])/2;
			xi[i]=(oi[i]-oi[N-i])/2;
			yr[i]=(oi[i]+oi[N-i])/2;
			yi1[i]=-(or[i]-or[N-i])/2;
		}
		
		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];	
		}
		register int TN=N<<1;
		for (i=N+1;i<TN;i++)	//全长后半=前半的共额
		{
			R[i]=R[TN-i];
			I[i]=-I[TN-i];
		}
	}	
//	after=clock();
}
/*
void FFT(float R[], float I[],int N, int M)
{
	register int j=0;
	register int b=0;
	register int k=0;
	int G=1;
	int H=0;
	int C=0;
	float cr=0;
	float ci=0;

	for(b=0; b<M; b++)
	{
		 H=G;
		 G*=2;

		 for(j=0; j<N; j+=G)
		 {
			 for(k=0; k<H; k++)
			 {
				register int long o,p,q;
				o=j+k;
				p=k*N/G;
				q=o+H;
				//乘开

				if (p==0)
				{
					cr=R[q];
					ci=I[q];
				} 
				else
				{cr=WR[p]*R[q]-WI[p]*I[q];
			    ci=WR[p]*I[q]+WI[p]*R[q];
				}
			    
				R[q]=R[o]-cr;
				I[q]=I[o]-ci;

    			R[o] =R[o]+cr;
				I[o]=I[o]+ci;
			 }
		 }
	}
	
}
*/


void FFT(float *R, float *I,int N, int M)
{
	register int GN=N/2;    //每级多少组蝶型GN=GroupNum
	register int BN=1;		//每组多少个蝶型BM=ButterflyNum
//	register int l;			//l=0....M-1级
//	register int j;			//j=0....GN组
//	register int k;			//每组蝶型k=0...BN个
//   register  int l,j,k;

//     int l,j,k;
    int l;
	for (l=0;l<M;l++)
	{ int j;
	   for (j=0;j<GN;j++)

		{   int k;
			for (k=0;k<BN;k++)

			{
				register int p1=BN*j*2+k;  //前半位置p1=position1
				register int p2=p1+BN;		//p2=position2
				register int p=k*GN;				//旋转因子系数WnP
				register float cr,ci;
/////////////////////////////////////////////////////  判断方式1	
/*	
				if(R[p1]==0&&I[p1]==0    &&     R[p2]==0&&I[p2]==0)         //pi=p2=0
                    break;
                else if(R[p2]==0&&I[p2]==0  )
                    // && R[p1]!=0&&I[p1]!=0)     //p1!=0,p2=0
                {
                    R[p2]=R[p1];I[p2]=I[p1];
                    
                    //break;
                }
               else if(R[p1]==0&&I[p1]==0   )
                    //&&R[p2]!=0&&I[p2]!=0)     //p1=0,p2!=0
                {
                    if (p==0)        //wr[p]=1,wi[p]=0
				    {
					   cr=R[p2];
					   ci=I[p2];
				    }
		/*		    else if(p==N/2)	//wr[p]=-1,wi[p]=0
				    {	
					cr=-R[p2];
					ci=-I[p2];		
				    }
		          	else if(p==N/4)    //wr[p]=0,wi[p]=-1
				    {
				    cr=I[p2];
					ci=-R[p2];
			         }
			     	else if(p==3*N/4)   //wr[p]=0,wi[p]=1
				    {
					cr=-I[p2];
					ci=R[p2];
				    }
				    else
				    {
					   cr=WR[p]*R[p2]-WI[p]*I[p2];
					   ci=WR[p]*I[p2]+WI[p]*R[p2];
				    }

                  //  cr=WR[p]*R[p2]-WI[p]*I[p2];
				  //  ci=WR[p]*I[p2]+WI[p]*R[p2];
				    
				    R[p2]=-cr;
				    I[p2]=-ci;
				    
				    R[p1]=cr;
				    I[p1]=ci;
				  //  break;
                }
                else
               {
                    
				
              //  cr=WR[p]*R[p2]-WI[p]*I[p2];
              //  ci=WR[p]*I[p2]+WI[p]*R[p2];
				

				    if (p==0)        //wr[p]=1,wi[p]=0
				    {
					   cr=R[p2];
					   ci=I[p2];
				    }
			/*	    else if(p==N/2)	//wr[p]=-1,wi[p]=0
				    {	
					cr=-R[p2];
					ci=-I[p2];		
				    }
		          	else if(p==N/4)    //wr[p]=0,wi[p]=-1
				    {
				    cr=I[p2];
					ci=-R[p2];
			         }
			     	else if(p==3*N/4)   //wr[p]=0,wi[p]=1
				    {
					cr=-I[p2];
					ci=R[p2];
				    }
				    else
				    {
					   cr=WR[p]*R[p2]-WI[p]*I[p2];
					   ci=WR[p]*I[p2]+WI[p]*R[p2];
				    }
				
			     	R[p2]=R[p1]-cr;
			     	I[p2]=I[p1]-ci;

    		      	R[p1] =R[p1]+cr;
		      		I[p1]=I[p1]+ci;
                }
                
*/
/////////////////////////////////判断法2

            if(R[p2]==0&&I[p2]==0)
                if(R[p1]==0&&I[p1]==0)break;
                else
                {
                    R[p2]=R[p1];I[p2]=I[p1];
                }
            else
                if(R[p1]==0&&I[p1]==0)
                {
                    if (p==0)        //wr[p]=1,wi[p]=0
				    {
					cr=R[p2];
					ci=I[p2];
				    }
				    else
				    {
					cr=WR[p]*R[p2]-WI[p]*I[p2];
					ci=WR[p]*I[p2]+WI[p]*R[p2];
				
                    R[p2]=-cr;
				    I[p2]=-ci;
				
				    R[p1]=cr;
				    I[p1]=ci;

                    }
				
                }
                else
                {
                  //  cr=WR[p]*R[p2]-WI[p]*I[p2];
                  //  ci=WR[p]*I[p2]+WI[p]*R[p2];
				
		      		if (p==0)        //wr[p]=1,wi[p]=0
		      		{
		      			cr=R[p2];
		      			ci=I[p2];
		      		}
	       			else
	       			{
					cr=WR[p]*R[p2]-WI[p]*I[p2];
					ci=WR[p]*I[p2]+WI[p]*R[p2];
	       			}
		  		
	       			R[p2]=R[p1]-cr;
	       			I[p2]=I[p1]-ci;

           			R[p1] =R[p1]+cr;
		      		I[p1]=I[p1]+ci;
                }

			}
		}
		GN>>=1;
		BN<<=1;
	}
}


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]);		
	}
	
}