fft.cpp

简单的cdma基础仿真

/*************************************************************************
 * COPYRIGHT 2007     Author：YangTao  
 * The College of Computer Science 
 * Beijing University of Technology
 * Emails：[e-mail]
 ************************************************************************/

#include <math.h>
#include <complex>

using namespace std;
const double PI=3.1415926535;


void fft(complex<double> * TimeDomain, complex<double> * FreqDomain, int ex,int sign)
{
	long Nfft;        //进行傅立叶变换的数据个数
	int	i,j,k;     //定义循环变量
	int	Block;     //划分子块的大小
	int OffSet;    //计算偏移量
	double Angle;  //相位角度
	complex<double> *Wn,*Temp1,*Temp2,*Temp; //定义数据指针

	Nfft = 1 << ex; // 计算付立叶变换点数
	
	// 分配运算所需存储器
	Wn  = new complex<double>[Nfft / 2];
	Temp1 = new complex<double>[Nfft];
	Temp2 = new complex<double>[Nfft];
	
	for(i = 0; i < Nfft / 2; i++) //计算加权系数Wn
	{
		Angle = -i * PI * 2 / Nfft;  //计算相应的角度
		Wn[i] = complex<double> (cos(Angle), sin(Angle)); //计算对应的复数
	}

	for(j = 0; j < Nfft; j++)            //重新排列傅立叶变换结果的次序
	{
		OffSet = 0;                   //初始化原来的数据处理后的新位置
		for(i = 0; i < ex; i++)       //对位置序号的二进制进行倒序反转
		{
			if (j&(1<<i))             //提取二进制的位信息
			{
				OffSet+=1<<(ex-i-1);  //反转二进制位序
			}
		}
	   if(sign==1)	 Temp1[j]=TimeDomain[OffSet]; 
	   else {
		   Temp1[j]=FreqDomain[OffSet];
	   
	      // -imag (Temp1[j])->imag (Temp1[j]);
	    
		     Temp1[j]= complex<double> (real(Temp1[j]), -imag(Temp1[j]))  ;                      //按照新位序将结果重排
	         }     
	   
     	}

	// 采用时间抽取的蝶形算法进行快速傅立叶变换
	for(k = 0; k < ex; k++)  //迭代次数
	{
		for(j = 0; j < 1 << (ex-k-1); j++)  //每次迭代分块
		{
			Block = 1 << (k+1);  //块内数据量
			for(i = 0; i < Block / 2; i++)  //块内进行蝶形计算
			{
				OffSet = j * Block;  //计算偏移量
				Temp2[i + OffSet] = Temp1[i + OffSet] + Temp1[i + OffSet + Block / 2] * Wn[i * (1<<(ex-k-1))]; // A=a+b*Wn
				Temp2[i + OffSet + Block / 2] = Temp1[i + OffSet] - Temp1[i + OffSet + Block / 2] * Wn[i * (1<<(ex-k-1))]; //B=a-b*Wn
			}
		}
		Temp = Temp1;    Temp1 = Temp2;    Temp2 = Temp;    //交换操作和结果数据的指针
	}
	

	for(j = 0; j < Nfft; j++)  //将结果赋给频域数组
	{
		if(sign==1)  FreqDomain[j]=Temp1[j];
		else {
		//	real(TimeDomain)=real(Temp1[j])/N;
        //    imag(TimeDomain)=-imag(Temp1[j])/N;
			
           TimeDomain[j]=complex<double> (real(Temp1[j])/Nfft, -imag(Temp1[j])/Nfft);
		}
		
	}



	//释放内存
	delete Wn;
	delete Temp1;
	delete Temp2;
}