傅氏变换的快速算法.cpp

几个c语言的程序

// 这是使用应用程序向导生成的 VC++ 
// 应用程序项目的主项目文件。

// 这是使用应用程序向导生成的 VC++ 
// 应用程序项目的主项目文件。

// 这是使用应用程序向导生成的 VC++ 
// 应用程序项目的主项目文件。

#include "stdafx.h"
#define _USE_MATH_DEFINES //使用系统定义的PI等常量
#using <mscorlib.dll>
#include <tchar.h>
#include "math.h"
#include <complex>
#include "iostream.h"


void ReOrder_For_Fourea(double *DataReal,double *DataImage ,int ClassNum);
int * Reorder(int Length);
void ReOrder_For_FoureaAnother(double *DataReal,double *DataImage,long Length);
using namespace std;


void Fourea(double * DataReal,double *DataImage,long Length=32,bool IsIFFT=false)//以指针操作
{
	
	int NN=1,ButterFlyNum=0;
	int i,j=0,K;
	int q,kk,x;
    int Fswitch=1;
    int LPointNum=1,LDistance,NextNum;//每一次计算的点数及点距
	double WReal,WImage,WCReal=1.0,WCImage=0.0;//计算时用到的系数值
	double tReal,tImage,temp;//计算时用的中间存贮数
//先检验长度是否是基2的
	do{
		NN*=2;   
	if(NN>Length)
		break;
    ButterFlyNum++;//蝶形的级数
	 }while(ButterFlyNum<32);

	
    Length=NN>>1;//获得以保证长度是Length为2的N次方

	
/*	for(int i=1;i<Length;i++)
	{ //先按位反转,注意到这里的技巧,是按位模
		 K=Length/2;//
		 q=i;
		 j=0;
         for(kk=ButterFlyNum;kk>=1;kk--)
		 {
			 x=q%2;//取位			
			 q=q/2;//为下一次取位做准备
			 j +=x*K;//取反后的值
			 K/=2;//下一次取位时要用的值
		  }	
	 //按位取反后,调换,但要注意到如果I>J时,又将换回去,所以规定I<J
	   if(i<j)
		  {
           double tempReal=DataReal[j];
		   double tempImage=DataImage[j];
		   DataReal[j]=DataReal[i];
		   DataImage[j]=DataImage[i];
		   DataReal[i]=tempReal;
		   DataImage[i]=tempImage;		
		}
	   }*/
	ReOrder_For_Fourea(DataReal,DataImage,ButterFlyNum);

	//下面进入傅立叶变换
    if(IsIFFT==true)
		Fswitch=-1;

	//逐级进行变换
	for(kk=1;kk<=ButterFlyNum;kk++)
	{
      LDistance=LPointNum;//间隔点数,上一次的LPointNum
      LPointNum*=2;//=pow(2,kk) 去掉平方函数的调用	 
      WCReal=1.0;//初始系数
	  WCImage=0.0;
	  WReal=cos(M_PI/LDistance);//本次计算的加权实部
	  WImage=Fswitch*sin(M_PI/LDistance);//虚部	
     for(j=0;j<LDistance;j++)
	 {		 
		 i=j;		 
		 while(i<Length)
		 {			
			 NextNum=i+LDistance;			
			 tReal=DataReal[NextNum]*WCReal+DataImage[NextNum]*WCImage;
			 tImage=-DataReal[NextNum]*WCImage+DataImage[NextNum]*WCReal;
			 DataReal[NextNum]=DataReal[i]-tReal;
			 DataImage[NextNum]=DataImage[i]-tImage;//注意到赋值先后顺序,
			 DataReal[i]+=tReal;
			 DataImage[i]+=tImage;
			 i+=LPointNum;//遍历整个数组中的同一个J下的数,即用同一个系数乘
		 }
		 temp=WCReal*WReal-WCImage*WImage;//获得下一次的加权系数实部		
		 WCImage=WCReal*WImage+WCImage*WReal;//虚部
		 WCReal=temp;
	 }
   }

	//如果是逆傅立叶变换,则将运算结果除Length;

	if(IsIFFT==true)
	{
		for(i=0;i<Length;i++)
		{
			DataReal[i]/=Length;
			DataImage[i]/=Length;
		}
	}
//输出傅立叶系数结果,由于是利用指针,故可在主函数中直接显示	
	
/*for(i=0;i<Length;i++)
	{  
		//cout<<"经傅立叶变换后,其实部的系数为: "<<endl;
		//cout<<DataReal[i]<<endl;
        cout<<"经傅立叶变换后,其虚部的系数为: "<<endl;
		cout<<DataImage[i]<<endl;
	}*/	
}



//这是一个基四的快速傅立叶变换
void Fourea_Base4(double *DataReal,double *DataImage,long Length)
{
	int ButterflyClass=0;
	int ID,IS;
	long j,I0,I1,I2,I3,NFourth;
	long NTwo=Length<<1;//2*Length;
	long i=1;
	double EAngle,Angle;
	double CosFirst,SinFirst,CosThree,SinThree;
	double Real1,Real2,Image1,Image2,Image3;
	//先求运算级数
	do{
        i*=2;
		ButterflyClass++;
		if(i==Length)
			break;
	  }while(i<Length);

	

  //直接运算而不需要反序
 for(int k=0;k<ButterflyClass-1;k++)
 {
   NTwo=NTwo>>1;//=NTwo/2;这里用移位是为了加速  
   NFourth=NTwo>>2;//NTwo/4;//这里要注意,将NTwo序列分为四段   
   EAngle=2.0*M_PI/NTwo;//公式中的值
   CosFirst=1.0;//旋转因子的初始值
   SinFirst=0.0;
   CosThree=1.0;
   SinThree=0.0;
   for( j=0;j<NFourth;j++)
   { 	
	 IS=j;//初始值
	 ID=NTwo*2;//在同一个旋转因子时要用到的跳过的序列的标号数,这里很重要

    do{
	   I0=IS;
       while(I0<Length)
	   {//在此序列中用相同的旋转因子进行运算
		 I1=I0+NFourth;//这里的I0,I1,I2,I3都是下标,I0为第一个数,I1为I0+N/4,
		 I2=I1+NFourth;//同理,I2为I0+N/2;
		 I3=I2+NFourth;//I3为I0+3N/4;
         Real1=DataReal[I0]-DataReal[I2];
		 DataReal[I0]=DataReal[I0]+DataReal[I2];
		 Real2=DataReal[I1]-DataReal[I3];
		 DataReal[I1]=DataReal[I1]+DataReal[I3];
         Image1=DataImage[I0]-DataImage[I2];
		 DataImage[I0]=DataImage[I0]+DataImage[I2];
		 Image2=DataImage[I1]-DataImage[I3];
		 DataImage[I1]=DataImage[I1]+DataImage[I3];
         Image3=Real1-Image2;
		 Real1=Real1+Image2;
		 Image2=Real2-Image1;
		 Real2=Real2+Image1;
		 DataReal[I2]=Real1*CosFirst-Image2*SinFirst;
		 DataImage[I2]=-Image2*CosFirst-Real1*SinFirst;//以上都是按公式进行计算
		 DataReal[I3]=Image3*CosThree+Real2*SinThree;//所以没有什么技巧可言
		 DataImage[I3]=Real2*CosThree-Image3*SinThree;
		 I0+=ID;//重新赋值为循环 ,这里要仔细体会     
	    }
	 IS=2*ID-NTwo+j;//这里很重要
	 ID=4*ID;//这里的很重要	 
	 }while(IS<Length);
	//当相同的旋转因子遍历整个序列后,才进入下一次的循环,下面是计算新的旋转因子
     Angle=(j+1)*EAngle;//这里的赋值很有技巧
     CosFirst=cos(Angle);
	 SinFirst=sin(Angle);
     CosThree=CosFirst*(CosFirst*CosFirst-3.0*SinFirst*SinFirst);//cos(Angle3);为了减少函数调用
	 SinThree=SinFirst*(4.0*CosFirst*CosFirst-1);//sin(Angle3);算出3倍角度正,余弦值
	 }
    //整个循环要做到4点化为两点为止   
   }
 

//,以下是计算两点的傅立叶变换
IS=0;
ID=4;//这里用4为步长是根据规律性来的
do{
  I0=IS;
  while(I0<Length-1)
   {
	I1=I0+1;
	Real1=DataReal[I0];
	DataReal[I0]=Real1+DataReal[I1];
	DataReal[I1]=Real1-DataReal[I1];
	Real1=DataImage[I0];
	DataImage[I0]=Real1+DataImage[I1];
	DataImage[I1]=Real1-DataImage[I1];
	I0+=ID;
   }
IS=2*(ID-1);//注意到这里的重新赋值
ID=4*ID;
}while(IS<Length);
//反序,直接调用函数
//整序调用了最快的一种算法。
ReOrder_For_Fourea(DataReal,DataImage,ButterflyClass); 	
}




//这个程式用递归的办法来获得逆的序列下标号
int * Reorder(int N)
{//N为序列的标号最大的值,采用递归算法来获得逆序列的标号
 //然后与需要整序的序列进行比较,如果
//相同则不需要交换.
//当然只存偶数部分也可,只是算法以及最后的
	int *Data,*Data1;	
	int DataTwo,IS;
	int NSecond=N>>1;//相当于除2运算	
	Data=new int[N];
	Data1=new int[NSecond];
	
	 if(N==2)
	{
		Data[0]=0;
		Data[1]=1;
	}
	else
	{
		Data1=Reorder(NSecond);//递归算法
		IS=NSecond;//为下面的后半部分运算减少运算量
		for(int i=0;i<NSecond;i++)
		{
			DataTwo=Data1[i]<<1;//移位运算相当于乘2运算			
			Data[i]=DataTwo;
			Data[IS]=DataTwo+1;
            IS++;//此下标的计算是后面的奇数部分
		}
	}
	return Data;  
}
//下面是利用上面的递归算法获得的逆序列来进行整序

//
//
//下面是另一种逆序标号算法的设计,主要是为了减少函数的调用而设计
int *Reorder_A(int Length)
{//Length为此序列的长度
//先求出需要运算的级数
	short ClassNum;
	int *Data;
	Data=new int[Length];
	Data[0]=0;
	Data[1]=1;
    int i,Bknum,kpow=1,j=1;

	//
	for(int i=1;i<32;i++)
	{
		j<<=1;//用移位运算来获得快速的算法
		ClassNum++;
		if(j==Length)
			break;
	}	
	//ClassNum就是运算的级数

	for(i=2;i<=ClassNum;i++)
	{
		kpow<<=1;//通过移位运算来获得pow(2,i-1)的效果
		Bknum=kpow;//通过这种赋值来获得下标的快速运算
		for(j=0;j<kpow;j++)
		{
			Data[j]<<=1;//原位乘2再赋值
			Data[Bknum]=Data[j]+1;//递推归律
			Bknum++;
		}
	}	
 return Data;
}

void ReOrder_For_Fourea(double *DataReal,double *DataImage,int ClassNum)
{
	int *Serial;
	int k;
	double temp;
	long Length=pow(2,ClassNum);//数据长度

	Serial=Reorder_A(Length);//调用前面的函数来获得逆序的标号的数组的指针.
	
	for(int i=0;i<Length;i++)
	{
		if(i<Serial[i])
		{
			k=Serial[i];
			temp=DataReal[i];
			DataReal[i]=DataReal[k];
			DataReal[k]=temp;
			temp=DataImage[i];
			DataImage[i]=DataImage[k];
			DataImage[k]=temp;
		}
	}
}


//下面是另一个反位的子程序
void ReOrder_For_FoureaAnother(double *DataReal,double *DataImage,long Length)
{
	int i,j,x,K,q=1 ,ButterFlyNum;

	//先求出蝶形运算的级数
	for(i=1;i<=32;i++)
	{
		ButterFlyNum=i;
		q*=2;
		if(q==Length)
			break;
	}

	for(i=1;i<Length;i++)
	{
		 K=Length/2;
         q=i;
		 j=0;     
      for(int kk=ButterFlyNum;kk>=1;kk--)
		 {
			 x=q%2;//取位			
			 q=q/2;//为下一次取位做准备
			 j +=x*K;//取反后的值
			 K/=2;//下一次取位时要用的值
		  }	
	 //按位取反后,调换,但要注意到如果I>J时,又将换回去,所以规定I<J
	   if(i<j)
		  {
           double tempReal=DataReal[j];
		   double tempImage=DataImage[j];
		   DataReal[j]=DataReal[i];
		   DataImage[j]=DataImage[i];
		   DataReal[i]=tempReal;
		   DataImage[i]=tempImage;		
		  }
	  }
  }


//下面的一个程序是为了当输入全为实数时的计算方法1

void Fourea_For_RealNumber(double *X,long Length,bool IsIFFT=false)
{
	double *DataOdd, *DataEven;//存放X数下标为奇数和偶数
	long M=Length/2-1;
	long MM=M/2;//下面循环时的终值

	DataOdd=new double[M+1];
	DataEven=new double[M+1];
	//下面的数据为后面的计算
	double aReal,bReal,cImage,dImage,RealTemp,ImageTemp;
	double Angle=M_PI/M;
    double CosA=cos(Angle),SinA=sin(Angle);
	double CosF=1.0,SinF=0.0,temp;	
	int j,k;

	for(int i=0;i<M;i++)
	{
		DataEven[i]=X[i*2];//存放偶数点号,存实部		
		DataOdd[i]=X[i*2+1];//存虚部		
	}	

   //下调用傅立叶变换
	//Fourea(DataEven,DataOdd,M);
	Fourea_Base4(DataEven,DataOdd,M);
   //将最后一个数据设为第一个数据的实部和虚部,以便于编程式
 
	DataEven[M]=DataEven[0];
	DataOdd[M]=DataOdd[0];
	for(j=0;j<=MM;j++)
	{  // 公式中的数据计算
		k=M-j;
		aReal=(DataEven[j]+DataEven[k])/2.0;
		bReal=(DataEven[j]-DataEven[k])/2.0;
		cImage=(DataOdd[j]+DataOdd[k])/2.0;
		dImage=(DataOdd[j]-DataOdd[k])/2.0;
		//先算实部
       RealTemp=CosF*cImage-SinF*bReal;	  
	   DataEven[j]=aReal+RealTemp;
	   DataEven[k]=aReal-RealTemp;
	   //再算虚部