利用FFT成功实现拓扑识别（十一）--对FFT算法的理解4

FFT算法

上篇我们用c程序了DFT的算法，从文中我们可以看出，对于N点的DFT，

每计算一个X(k)的值，需要N次的复数乘法运算和（N-1）次的复数加法运算，所以，完成整个DFT运算，复杂度为 ： 当N值较大时，随N增大而急速增大，运算量是相当可观的，同等硬件条件下，运算量大就意味着耗时会增加，实时性会降低。如何解决这个问题昵？

库利－图基FFT算法

库利－图基快速傅里叶变换算法（Cooley-Tukey算法）可能有人不太清楚，可是提起他的另一个通俗名称—蝶形算法，在FFT领域可谓是大名鼎鼎，Cooley-Tukey是最常见的快速傅里叶变换算法。这一方法以分治法为策略，将序列长为N的DFT分区为两个长为N/2的子序列的DFT，并依次递归调用。在FFT算法中，针对输入分组方式的不同，可以使用时域抽取（Decimation-in-time），或是频域抽取（Decimation-in-frequency），两者的区别在于输入和输出的排序不同。下面以最常用的时域抽取DIT FFT算法，对整个流程进行剖析：

算法原理

前面我们已经知道，DFT的每一个输出，都是N个复数乘积的和，即：

按照惯例，
被为旋转因子（twiddle factor）（为什么？请结合前文思考下），用

表示，我们从DFT的计算开始推导：

先按照奇数偶数项进行分组

我们用数学方式表示就是：

从奇数相中，利用指数分解法，对

分解后，为：

将进行等效变换，表示为：的形式，最终结果如下

看出来了吗？原来N点DFT运算，已经可以用两个N/2点的DFT表达出来了，omg！，看个简单的，假设N=8，

对N/2的两个4点DFT递归进行分解，

继续对对N/4的两个2点DFT进行分解

看起来像什么？是不是如一只

蝴蝶翩翩起舞，这就是蝶形算法的来源，至此，递归结束，完整的8点DIT FFT结构图如下：

代码分析

搞清楚了FFT算法的原理，再看代码就简单了，在此要说明的是，原理是统一的，代码如何实现，方法却是多种多样的，此处为贴合上述原理，采用了递归的方法，也有许多算法实现中，将递归算法改写为非递归的形式，总体思想是一样的。fft函数的实现如下：

cplx t; //定义复数临时变量
double a, b, c, d, theta;
printf(“----into fft function para step=%i,n=%i-\n”,step,n);
if (step < n)
{
// 递归调用，计算偶数和奇数组的fft
fft(out, buf, n, step * 2);
fft(out + step, buf + step, n, step * 2);
for (i = 0; i < n; i += 2 * step)
{
//蝶形运算
theta = -PI * i / n;
t.real = cos(theta) * out[i + step].real - sin(theta) * out[i + step].image;
t.image = sin(theta) * out[i + step].real + cos(theta) * out[i + step].image;
// 系数为1
buf[i / 2].real = out[i].real + t.real;
buf[i / 2].image = out[i].image + t.image;
// 系数为-1
buf[(i + n) / 2].real = out[i].real - t.real;
buf[(i + n) / 2].image = out[i].image - t.image;
}
}

您看明白了吗？欢迎交流，共同提高。

效率分析

利用FFT替换DFT算法之后，任何一个N为2的整数幂（即
）
的DFT，都可以通过M次分解，最后成为2点的DFT来计算。经过M次分解计算，每级由N/2个蝶形运算组成。而完成这样一个蝶形计算只需一次乘法和两次复数加法。因此，完成N点的时间抽选FFT计算的总运算量为：

复数乘法次数：

复数加法次数：所以以复数乘法的计算次数来比较DFT与FFT的效率为：以N=1024为例，效率提升：204.8倍