原理不讲,只有代码:
#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;
}