/*************************************************************************
*
* 函数名称:
*   FFT()
*
* 参数:
*   complex<double> * TD - 指向时域数组的指针
*   complex<double> * FD - 指向频域数组的指针
*   r -2的幂数,即迭代次数
*
* 返回值:
*   无。
*
* 说明:
*   该函数用来实现快速付立叶变换。
*
************************************************************************/
VOID WINAPI 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;
}
今天看了有关快速傅里叶变换的算法(如上),红色字体部分是我不明白的地方。根据书上原理介绍,我感觉该部分应该是
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2]* W[i * (1<<k)];
X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2] * W[i * (1<<k)];
但是我看了好多参考书,结果里面介绍的却是同红字部分的相同,这是为什么。请高手指点。

解决方案 »

  1.   

    FFT至少是十几年前的事了,忘了
    但LZ有没有意识到实际上X2和X1一直在交换,也就是X1X2可以交换红字部分,如果将X2与X1交换X1[i + p] = X2[i + p] + X2[i + p + bfsize / 2]; 
    X1[i + p + bfsize / 2] = (X2[i + p] - X2[i + p + bfsize / 2]) * W[i * (1 < <k)]; } 
    重新调整公式 得到
    X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2]/(2*W[i * (1 < <k)]);
    X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2]/(2*W[i * (1 < <k)]);
    看来LZ的公式从结构上与标红一致,问题出现在W的定义上
      

  2.   

    LZ没仔细看FFT,FFT快就快在他通过分立实虚部,达到实部只有+-运算(权值+-1),虚步才有权值参与