galois_godelfft.cpp

几个FFT算法

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <complex>
using namespace std;
//*************************************************************************
// *
// * 函数名称：
// *   FFT()
// *
// * 参数:
// *   complex<double> * TD- 指向时域数组的指针
// *   complex<double> * FD- 指向频域数组的指针
// *   r－2的幂数，即迭代次数
// *
// * 返回值:
// *   无。
// *
// * 说明:
// *   该函数用来实现快速付立叶变换。
// *
// ************************************************************************/
void FFT(	complex<double> *TD,  complex<double> *FD, 
			complex<double> *X1,  complex<double> *X2,
			complex<double> *Omega,	  int r)
{
	// 付立叶变换点数
	long count;

	// 循环变量
	int i,j,k;

	// 中间变量
	int bfsize,p;

	// 角度
	double angle;

	complex<double> *X;

	// 计算付立叶变换点数
	count = 1 << r;

	// 分配运算所需存储器
	//Omega  = new complex<double>[count / 2];
	//X1 = new complex<double>[count];
	//X2 = new complex<double>[count];

	// 计算加权系数
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * 3.1415926 * 2 / count;
		Omega[i] = complex<double>(cos(angle), sin(angle));
	}

	// 将时域点写入X1
	memcpy(X1, TD, sizeof(complex<double>) * count);

	// 采用蝶形算法进行快速付立叶变换
	for(k = 0; k < r; k++)
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (r-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * Omega[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}

	// 重新排序
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
		FD[j]=X1[p];
	}

	// 释放内存
	//delete Omega;
	//delete X1;
	//delete X2;
}