快速傅立叶变换(FFT)算法(原来这就是蝶形变换)

为了实现FFT的海面模拟，不得不先撸个FFT算法实现。

离散傅立叶变换(DFT)

学习FFT之前，首先要先了解什么是DFT，我们都知道傅立叶变换是将时域转换为频域。但是我们计算机是没办法处理连续的点，因此就有了离散傅立叶变换DFT。

标准DFT公式：

我们令：

W的一些性质

证明：

证明方法同上。

快速傅立叶变换(FFT)

我们将X(k)按照奇偶组合，在，有：

因为是常数，所以X(k)可以更改为：

即：

我们令

，

可得在时，

我们这里发现F，G为X的递归函数。

然后计算，：

根据以上公式：，

我们发现想要求与，只需要要求出与就可以了，这两条公式也就是我们所说的蝶形运算式。

我们令N=4，这时候出现了一幅熟悉的图，4-point FFT蝴蝶变换图：

说实话，我第一次看这个看了半天没看懂，就像瞎比划出来的东西，然而，我们只要举个例子就明白了：

我们首先计算k=0的情况：

其中，为F的递归函数。，为G的递归函数。

所以我们的为：

首先公式：

，

过程如下图所示(其中)：

下面两条公式：

，

图解为：

由于我们求出与也能求出的值，其完成的图为：

到这里，就是我们标准的傅立叶变换了，我们的算法复杂度从变成了。

bitreverse算法

然而，我们的FFT还可以使用非递归的版本，如我们的4-point FFT计算X[0],X[1],X[2],X[3]最后对应的输入为x[0],x[2],x[1],x[3]。

再看8-point FFT计算过程x的位置变化：

0 1 2 3 4 5 6 7

0 2 4 6 1 3 5 7

0 4 2 6 1 5 3 7

bitreverse算法则是从中发现了x[i]在几号位，只要将i转化为(N-point FFT)位二进制，然后将bit反序再转回10进制即可，所得的数就是在第几位。比如在8-point FFT中我们想知道x(4)在几号位，我们只要将2转化为位二进制100，然后反序得001，即10进制的1。所以x(4)最后的位置为1号位(从0开始)。

所以我们的FFT便可以写成迭代的版本了(即上面蝶形变换图从左往右推)。贴上代码：

//其中a为x[i],omg为WNkvoid fft(cp *a, cp *omg){int lim = 0;while((1 << lim) < n) lim++;//logN //二进制反序for(int i = 0; i < n; i++){int t = 0;for(int j = 0; j < lim; j++)//每当该位上的值为1时，t和该位反转后的值进行或运算if((i >> j) & 1) t |= (1 << (lim - j - 1));if(i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次}//FFT计算for(int l = 2; l <= n; l *= 2){int m = l / 2;for(cp *p = a; p != a + n; p += l)for(int i = 0; i < m; i++){//蝶形变换cp t = omg[n / l * i] * p[i + m];p[i + m] = p[i] - t;p[i] += t;}}}

其中FFT计算的第一层循环为下图所示分块计算：

第二层循环为如下分块计算：

我们的蝶形变换就是按下图：

每次迭代到下一次的g'(x)都可由当前的g(x)自计算获得。

到这里我们的FFT就全部结束啦，接下来可以拿他做FFT海面模拟啦，先上个图：