傅氏变换的快速算法.cpp

几个c语言的程序

	   ImageTemp=-CosF*bReal-SinF*cImage;	   
	   DataOdd[j]=ImageTemp+dImage;
	   DataOdd[k]=ImageTemp-dImage;
	   //加权的角度要变化
      temp=CosA*CosF-SinA*SinF;//新的Cos余弦值
	  SinF=CosF*SinA+SinF*CosA;//新的正弦值
	  CosF=temp;
	}

	//这样就得出了其实部和虚部,再将其返回到X数组中,

for(int i=0;i<M+1;i++)
	{	
		X[2*i]=DataEven[i];//存贮的是实部
        X[2*i+1]=DataOdd[i];//存贮的是虚部		
	}	
}

//一种为实数列的逆变换,利用卷积来完成

void Fourea_ForRealNum_IFFT(double *X,long Length)
{
	long i,kk,ImageJ,AnImageJ,ImageAJ,ImageBJ;	
	long LengthM=Length-2;
    long M=LengthM>>1;//通过移位最快
	long MM=LengthM>>2;//作为角度循环的终值
	long MMM=LengthM>>3;//计算角度时用到的循环终值
	double *DataRealEven,*DataImageEven,*CtanValue;//设四个数组,一次计算偶数和奇数	
	double EvenNum=(X[0]+X[LengthM])/LengthM,temp,temp1,sum;//偶数时用到的常数值
	double OddNum=(X[0]-X[LengthM])/LengthM;//奇数时用的常数值
	double CosA=cos(M_PI/LengthM),SinA=sin(M_PI/LengthM);//定值为PI除去2N的正弦和余弦值
	double CosB=CosA*CosA-SinA*SinA,SinB=2.0*CosA*SinA;//同时也为一个正弦和余弦值
	double CosF=1.0,SinF=0.0;
	DataRealEven=new double[M];
	DataImageEven=new double[M];//数组的长度是原来的一半
    CtanValue=new double[MM];//用来存贮卷积时用到的虚部值
	

	
	//要将数组中的数据取共轭
	for( i=1;i<M;i++)
	{
		DataRealEven[i]=X[i*2]/M;		
		DataImageEven[i]=-X[i*2+1]/M;
	}
	DataRealEven[0]=0.0;
	DataImageEven[0]=0.0;
	//先计算,先反序,计算方法与傅立叶变换大致一样
  // Fourea(DataRealEven,DataImageEven,M);
	Fourea(DataRealEven,DataImageEven,M);//这里是用此程序来试一下.
 
 
//算余切值以减少函数调用,同时也为了减少运算
  for(kk=0;kk<MM;kk++)
  {   
	  //算余切值
	  CtanValue[kk]=(CosA*CosF-SinA*SinF)/(SinA*CosF+CosA*SinF);//余切值的表达式
	  //CtanValue[MM-1-kk]=1.0/CtanValue[kk];//根据对称关系得出
	  //对CosF,SinF进行重新赋值
	  temp=CosF;
	  CosF=CosF*CosB-SinF*SinB;
	  SinF=temp*SinB+SinF*CosB;
  }
 
  
 //直接赋值以及算奇数部分来赋值以减少内存占用
  for(i=0;i<M;i++)
  {
	  sum=0.0;
	  X[2*i]=DataRealEven[i]+EvenNum;//偶数部分	
	 for( kk=0;kk<MM;kk++)
	 {//卷积的算法还是;较慢的
		ImageJ=(kk+i+1)%M;
		AnImageJ=(M-kk+i)%M;
		//ImageAJ=(MM-kk+i)%M;
		//ImageBJ=(MM+kk+i+1)%M;
		//temp1=(DataImageEven[ImageJ]-DataImageEven[AnImageJ])*CtanValue[kk];
		//temp=(DataImageEven[ImageAJ]-DataImageEven[ImageBJ])/CtanValue[kk];
        //sum-=temp1+temp;	
		sum-=(DataImageEven[ImageJ]-DataImageEven[AnImageJ])*CtanValue[kk];		
	  } 
	  
	  X[2*i+1]=OddNum+sum/M;//奇数部分	  
  }
}
//下面采用逆变换程序来做

void Fourea_IFFT(double * DataReal,double * DataImage,int Length)
	{//我们要对原来的数据先进行共轭,而后再进行运算,其实共轭只针对虚
	//部进行,所以就直接对此进行运算,我们直接在原数据上进行.
	int i;
	//下面对虚部取共轭
	for(i=0;i<Length;i++)
		DataImage[i]=-DataImage[i];
	//然后调用基4算法
	Fourea_Base4(DataReal,DataImage,Length);
	//求时域点的共轭,并要除以一个系数
	for(i=0;i<Length;i++)
	{
		DataReal[i]=DataReal[i]/Length;
		DataImage[i]=-DataImage[i]/Length;
	}
   }//此逆变换


//直接用原数据来做实数列的逆变换,程序较易编制
void Fourea_ForRealNumBer1_IFFT(double *X,long Length)
{   
	int i;
	long LengthM=Length-2;
	double *DataReal,*DataImage;
	DataReal=new double[LengthM];
	DataImage=new double[LengthM];
	//将数据按规律赋值
	for( i=0;i<LengthM;i++)
	{//这里是将其先共轭
		if(i>LengthM/2)
		{
			DataReal[i]=X[(LengthM-i)*2];
			DataImage[i]=-X[(LengthM-i)*2+1];
			//DataImage[i]=X[(LengthM-i)*2+1];
		}
		else
		{
			DataReal[i]=X[i*2];
			DataImage[i]=X[i*2+1];
			//DataImage[i]=-X[i*2+1];
		}
	}
	//先度一下
	
	//
 /* DataReal[0]=0.0;
  DataImage[0]=0.0;
  DataReal[LengthM/2]=0.0;
  DataImage[LengthM/2]=0.0;
  double EvenNum=(X[0]+X[LengthM])/LengthM;
  double OddNum=(X[0]-X[LengthM])/LengthM;*/
  //下面调用基4的程序一试
  Fourea_IFFT(DataReal,DataImage,LengthM);
  for(i=0;i<LengthM;i++)
	  X[i]=DataReal[i];

	//Fourea(DataReal,DataImage,LengthM,true);
	//再将Datareal的值赋予X数组
	/* for(int i=0;i<LengthM;i++)
	 {
		//X[2*i]=EvenNum+DataReal[i*2]/LengthM;
		//X[2*i+1]=OddNum+DataReal[i*2+1]/LengthM;
		X[i]=DataReal[i]/LengthM;
	 }*/
}

//下面是一个最快的

void Fourea_for_realnumber_IFFT(double *X,long Length)
{
	long LengthM=Length-2;
	long i,j;
	long M=LengthM>>1;
	double Angle=M_PI/M;//一个角度
	double OddNum=(X[0]-X[LengthM])/LengthM;
	double EvenNum=(X[0]+X[LengthM])/2.0;
	double *DataReal,*DataImage;
	DataReal=new double[M];
	DataImage=new double[M];

	//赋值
	for(i=1;i<M;i++)
	{
		DataReal[i]=(X[2*i]*cos(Angle*i)-X[2*i+1]*sin(Angle*i));
		DataImage[i]=-(X[2*i]*sin(Angle*i)+X[2*i+1]*cos(Angle*i));
	}
	DataReal[0]=0.0;
	DataImage[0]=0.0;
   Fourea_Base4(DataReal,DataImage,M); 

}
//下面是将获取实序列变换后的数据
//下面是获得实序列的变换后的数据.
	void GetDataOfReal(double* Xreal,int Length, double* DataReal, 
		double * DataImage)
	{//还是以实部,虚部分开存贮,当然可以用复数来存贮.
		int i;
		int LengthM=Length-2;
		int LengthMM=LengthM/2;
		//对输出数组进行分配内存赋值
		Xreal[LengthM]=0.0;
		Xreal[LengthM+1]=0.0;
		
		//先调用上面的程序获得数据
		Fourea_For_RealNumber(Xreal,Length);
		//下面是将原始数据进行重新赋值
		for(i=0;i<LengthM;i++)
		{
			if(i>LengthMM)
			{//这里后面是将其共轭
				DataReal[i]=Xreal[(LengthM-i)*2];
				DataImage[i]=-Xreal[(LengthM-i)*2+1];
				
			}
			else
			{
				DataReal[i]=Xreal[i*2];
				DataImage[i]=Xreal[i*2+1];				
			}
		}
		//重新赋值完毕
	}//程序完毕

   //下面是二维图像的傅氏变换程序
	void FFT_2D(double * PData, int Length , double* DataReal,  double *DataImage)
	{//在这里我们的图像数据为实序列,所以我们将先调用此实序列的
	//来得到第一次数据,而后用此数据来进行第二次变换.
		int x,y,TempN,Ktemp;//这是用来循环的控制变量.
       
		int nTransWidth=8,nTransHeight=8;


		//下面进行第一次傅立叶变换,是X方向,即行方向的.
		double* Xreal=new double [nTransWidth+2];//这里因为是对实数据进行的
		double* DataRealTemp1=new double[nTransWidth];
		double* DataImageTemp2=new double[nTransWidth];//以上用来暂存数据用的.
		double* DataRealTemp3=new double[nTransWidth];
		double* DataImageTemp4=new double[nTransWidth];

		TempN=0;
		Ktemp=0;
		for(y=0;y<nTransHeight;y++)
		{
          //先赋值给Xreal数组
			for(x=0;x<nTransWidth;x++)
			{
				Xreal[x]=PData[y*nTransWidth+x];
			
			}
			//多赋两个值是为了便于计算
			Xreal[nTransWidth]=0.0;
			Xreal[nTransWidth+1]=0.0;
			//下面调用一维傅氏变换来进行处理.
		   //GetDataOfReal(Xreal,10, DataRealTemp1, DataImageTemp2);
			Fourea_For_RealNumber(Xreal,10);
			//下面将得来的数据重新赋给总的DataReal,DataImage
			//Fourea_Base4(DataRealTemp3,DataImageTemp4,nTransWidth);

			
			for(x=0;x<nTransWidth;x++)
			{
				//DataReal[y*nTransWidth+x]=DataRealTemp[x];
				//DataImage[y*nTransWidth+x]=DataImageTemp[x];
				if(x>nTransWidth/2)
			    {//这里后面是将其共轭
				DataReal[y*nTransWidth+x]=Xreal[(nTransWidth-x)*2];				
				DataImage[y*nTransWidth+x]=-Xreal[(nTransWidth-x)*2+1];				
			    }
			    else
			    {
				DataReal[y*nTransWidth+x]=Xreal[x*2];				
				DataImage[y*nTransWidth+x]=Xreal[x*2+1];
				
			    }
			}

		}//
		//DataReal[],DataImage[]目前存贮了经过行变换的结果,为了直接利用一维FFT,需要将
		//此数据进行转置,再一次利用一维FFT来进行变换
		//下面将转置与变换放到一起做,只是要明白其中的原理就可以了,我们直接取每列的数        //据来进行变换,然后再将此变换后的数据放回,但要注意的是,这时的变换是复
       //数变换了.
		//先设两个数组来存贮变换的数据
		double *TempReal=new double [nTransHeight];
		double *TempImage=new double [nTransHeight];
		for(x=0;x<nTransWidth;x++)
		{//需要调用的次数
			//先将原始要变换的数据赋予上面所设的两个数组
			for(y=0;y<nTransHeight;y++)
			{
				TempReal[y]=DataReal[x+y*nTransWidth];//注意此处的数据是如何赋予的,根                                                      //据转置本身的定义来的.
				TempImage[y]=DataImage[x+y*nTransWidth];
			}
			//得到数据后就将此数据进行变换
			Fourea_Base4(TempReal,TempImage,nTransHeight);
			//变换后的数据重新赋予给原数组
			for(y=0;y<nTransHeight;y++)
			{
				DataReal[x+y*nTransWidth]=TempReal[y];
				DataImage[x+y*nTransWidth]=TempImage[y];
				
			}

		}
	}//完毕*/
	//

	//下面是图像二维傅氏逆变换.
	void IFFT_2D(double* DataReal,double * DataImage, double *PData)
	{//PData是最后用于接受数据,即转变回来的数据,只接受变换回来的实部
		//为了和上面的相一致,先做列方向上的一维傅氏逆变换
		int i,j;//循环变量
		int nTransWidth=8,nTransHeight=8;
		//下
		//PData=new double[TotalLength];//这里的数据用来存

		double * TempReal=new double[nTransHeight];
		double *TempImage=new double[nTransHeight];
		//下面是进行列逆变换.我们这里不能直接利用逆变换
		for(i=0;i<nTransWidth;i++)
		{//调用次数
			//先将原始数据赋予给临时数组
			for(j=0;j<nTransHeight;j++)
			{
				TempReal[j]=DataReal[i+j*nTransWidth];
				TempImage[j]=DataImage[i+j*nTransWidth];
			}
			//然后进行逆变换
			Fourea_IFFT(TempReal,TempImage,nTransHeight);
			//逆变换完之后将数据写入原始数据中
			for(j=0;j<nTransHeight;j++)
			{
				DataReal[i+j*nTransWidth]=TempReal[j];
				DataImage[i+j*nTransWidth]=TempImage[j];
			}

		}
		//然后再进行行的逆变换
		double * tempReal=new double [nTransWidth];
		double* tempImage=new double[nTransWidth];
		for(i=0;i<nTransHeight;i++)
		{//先赋予数据给临时数组
			for(j=0;j<nTransWidth;j++)
			{
				tempReal[j]=DataReal[i*nTransWidth+j];
				tempImage[j]=DataImage[i*nTransWidth+j];
			}
			//然后进行逆变换
			Fourea_IFFT(tempReal,tempImage,nTransWidth);
			//最后将变换的数据的实部赋予PData
			for(j=0;j<nTransWidth;j++)
				PData[i*nTransWidth+j]=tempReal[j];
		}
	}





// 这是此应用程序的入口点
int _tmain(void)
{
    // TODO: 请用您自己的代码替换下面的示例代码。;//用指针来较为方便
	int i;
	double Data1[64]={0.22925599,0.76687497,
		              0.68317699,0.50919088,
					  0.87456000,0.64464098,
					  0.84746796,0.35396299,
					  0.39889199,0.45709398,
					  0.23630899,0.13318199,
					  0.16605200,0.22602700,
					  0.66245896,0.25021198,
					  0.61769694,0.26246497,
					  0.51266795,0.93920696,
					  0.62402797,0.42238200,
					  0.93970597,0.28206798,
					  0.46921799,0.054879200,
					  0.51983094,0.39682698,
					  0.11315700,0.60751694,
					  0.70150697,0.88705498,					  
	                 };
	for(i=32;i<64;i++)
		Data1[i]=i*0.01;
	double Data2[64],Data3[64],Data4[64];
	/*double Data2[10],Data3[8],Data4[8],Data5[8],Data6[8];
	for(i=0;i<8;i++)
	{
		Data2[i]=Data1[i];
		Data5[i]=Data2[i];
		Data6[i]=0.0;
	}
	Data2[8]=0.0;
	Data2[9]=0.0;

	GetDataOfReal(Data2,10,Data3,Data4);
	Fourea_Base4(Data5,Data6,8);
	for(i=0;i<8;i++)
	{
		if(fabs(Data3[i]-Data5[i])<0.001)
			cout<<"0"<<endl;

	}*/



	FFT_2D(Data1,64,Data2,Data3);
	IFFT_2D(Data2,Data3,Data4);
	for(i=0;i<64;i++)
		cout<<Data4[i]<<endl;

	


	
	 

    return 0;
}