fft2.cpp

使用C语言详细实现了两种不同的快速傅立叶变换算法

/*************************************************************************
 * 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;

/*************************************************************************
 * 函数名称：TD_FFT2()
 * 参数:
 *   complex<double> * TimeDomain	- 指向时域数组的指针
 *   complex<double> * FreqDomain	- 指向频域数组的指针
 *   ex				                - 2的幂数，即迭代次数
 * 返回值:void
 * 说明:该函数实现时间抽取的以2为基的快速付立叶变换。
 ************************************************************************/

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

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

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

	// 采用时间抽取的蝶形算法进行快速傅立叶变换
	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 < N; j++)  //将结果赋给频域数组
	{
		FreqDomain[j]=Temp1[j];
	}

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

/*************************************************************************
 * 函数名称：FD_FFT2()
 * 参数:
 *   complex<double> * TimeDomain	- 指向时域数组的指针
 *   complex<double> * FreqDomain	- 指向频域数组的指针
 *   ex				                - 2的幂数，即迭代次数
 * 返回值:void
 * 说明:该函数实现频率抽取的以2为基的快速付立叶变换。
 ************************************************************************/

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

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

	// 采用频率抽取的蝶形算法进行快速傅立叶变换
	for(k = 0; k < ex; k++)  //迭代次数
	{
		for(j = 0; j < 1 << k; j++)  //每次迭代分块
		{
			Block = 1 << (ex-k);  //块内数据量
			for(i = 0; i < Block / 2; i++)  //块内进行蝶形计算
			{
				OffSet = j * Block;  //计算偏移量
				Temp2[i + OffSet] = Temp1[i + OffSet] + Temp1[i + OffSet + Block / 2];  // A=(a+b)
				Temp2[i + OffSet + Block / 2] = (Temp1[i + OffSet] - Temp1[i + OffSet + Block / 2]) * Wn[i * (1<<k)];  // B=(a-b)Wn
			}
		}
		Temp = Temp1;    Temp1 = Temp2;    Temp2 = Temp;    //交换操作和结果数据的指针
	}
	
	for(j = 0; j < N; j++)            //重新排列傅立叶变换结果的次序
	{
		OffSet = 0;                   //初始化原来的数据处理后的新位置
		for(i = 0; i < ex; i++)       //对位置序号的二进制进行倒序反转
		{
			if (j&(1<<i))             //提取二进制的位信息
			{
				OffSet+=1<<(ex-i-1);  //反转二进制位序
			}
		}
		FreqDomain[j]=Temp1[OffSet];  //按照新位序将结果重排
	}
	
	//释放内存
	delete Wn;
	delete Temp1;
	delete Temp2;
}
void main()
{
    complex<double> *TD=new complex<double>[8];
    complex<double> *FD=new complex<double>[8];
    for(int i=0;i<8;i++)
		TD[i]=complex<double>(8-i, 0);
	printf("The data in Time Domain :\n");
	for(int i=0;i<8;i++)
		printf("%f,%fi\n",TD[i].real(),TD[i].imag());
	TD_FFT2(TD,FD,3);
	printf("The data in Frequency Domain after TD_FFT2 function :\n");
	for(int i=0;i<8;i++)
		printf("%f,%fi\n",FD[i].real(),FD[i].imag());
	FD_FFT2(TD,FD,3);
	printf("The data in Frequency Domain after FD_FFT2 function :\n");
	for(int i=0;i<8;i++)
		printf("%f,%fi\n",FD[i].real(),FD[i].imag());
	getchar();
}