快速傅里叶变换FFT(一维)

原理不讲，只有代码：

#define PI 3.14159

/*************************************************************
功能：       一维快速傅里叶变换
参数：       complex<double> *TD： 指向时域数组的指针
                 complex<double> *FD：  指向频域数组的指针
			     int r：    2的幂数，即迭代次数
返回值：   无
***************************************************************/

void  Ctry::FFT(complex<double> *TD, complex<double> *FD, int r)
{
	LONG count;         //傅里叶变换点数
	int    i, j, k;              //循环变量
	int    bfsize, p;         //中间变量
	double   angle;        //角度
	complex<double> *w, *x1, *x2, *x;
	//计算傅里叶变换点数
	count = 1 << r;

	//分配运算所需要的存储器
	w = new complex<double>[count / 2];
	x1 = new complex<double>[count];
	x2 = new complex<double>[count];

	//计算加权系数
	for (i = 0; i < count / 2; i++)
	{
		angle = -i*PI * 2 / count;
		w[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])*w[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 w;
	delete x1;
	delete x2;
}