接上一篇：杨超：fft海面模拟(一)

视频封面

00:15视频封面

00:13

本篇说FFT演算法。

FFT：Fast Fourier Transformation（快速傅里叶变换），是计算DFT的快速方法，而且能在gpu上实现。

一，递归形式的FFT演算法及复杂度

对于如下标准DFT：

X(k)=sum_{n=0}^{N-1}{x(n)e^{-ifrac{2pi k}{N}n}},kin{{0,1,...,N-1}}

注：为了书写方便，通常令 W_{N}^{k}=e^{-ifrac{2pi k}{N}} 。

可以看作是N个输入和N个输出的电器元件（N point DFT calculator），如图：

如果直接按DFT定义式暴力计算，每一个输出都需要计算N次乘法，故N个输出共需乘法N*N次，即演算法复杂度为O(N*N)，是比较高的。

快速傅里叶变换则是使用分治思想对DFT进行计算，可有效降低演算法复杂度。

注：FFT只用于计算N为2的幂的DFT。

考虑如何用两个N/2 point DFT calculator去构造出一个N point DFT calculator（N为2的幂）。

如果将序号为偶数的输入给到第一个N/2 point DFT calculator，序号为奇数的输入给到第二个N/2 point DFT calcuator，如下图所示：

则有：

G(k)=sum_{n=0}^{frac{N}{2}-1}{g(n)e^{-ifrac{2pi k}{frac{N}{2}}n}}=sum_{n=0}^{frac{N}{2}-1}{x(2n)e^{-ifrac{2pi k}{frac{N}{2}}n}},kin{{0,1,...,frac{N}{2}-1}}

H(k)=sum_{n=0}^{frac{N}{2}-1}{h(n)e^{-ifrac{2pi k}{frac{N}{2}}n}}=sum_{n=0}^{frac{N}{2}-1}{x(2n+1)e^{-ifrac{2pi k}{frac{N}{2}}n}},kin{{0,1,...,frac{N}{2}-1}}

如何用G(k)和H(k)得到X(k)呢？

结论是：

X(k)=left{egin{matrix} G(k)+W_{N}^kH(k) &,kin{{0,1,...,frac{N}{2}-1}} \ G(k-frac{N}{2})+W_{N}^kH(k-frac{N}{2}) & ,kin{{frac{N}{2},frac{N}{2}+1,...,N-1}} end{matrix} ight.

推导过程如下：

当 kin{{0,1,...,frac{N}{2}-1}} 时G(k)和H(k)均有定义，有：

X(k)=sum_{n=0}^{N-1}{x(n)e^{-ifrac{2pi kn}{N}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi k(2n)}{N}}}+sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi k(2n+1)}{N}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi kn}{N/2}}}+e^{-ifrac{2pi k}{N}}sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi kn}{N/2}}}

=G(k)+W_{N}^{k}H(k)

当 kin{{frac{N}{2},frac{N}{2}+1,...,N-1}} 时，令K=k-N/2，则 Kin{{0,1,...,frac{N}{2}-1}} ,有：

X(k)=sum_{n=0}^{N-1}{x(n)e^{-ifrac{2pi kn}{N}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi k(2n)}{N}}}+sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi k(2n+1)}{N}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi (K+frac{N}{2})(2n)}{N}}}+sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi (K+frac{N}{2})(2n+1)}{N}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi Kn}{N/2}}}-e^{-ifrac{2pi K}{N}}sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi Kn}{N/2}}}

=sum_{n=0}^{N/2-1}{x(2n)e^{-ifrac{2pi Kn}{N/2}}}+e^{-ifrac{2pi k}{N}}sum_{n=0}^{N/2-1}{x(2n+1)e^{-ifrac{2pi Kn}{N/2}}}

=G(K)+W_{N}^{k}H(K)

=G(k-frac{N}{2})+W_{N}^{k}H(k-frac{N}{2})

根据上面X与H、G的关系，可补全电路图如下：

至此完成了用两个N/2 point DFT calculator构造N point DFT calculator。

以上就是递归形式的FFT演算法。但递归形式一般效率不佳，尤其是不适合在gpu上实现，所以经典的FFT演算法并不是采用这种形式，而是采用展平的形式，所谓蝶形网路（见下一节）。

无论是递归形式还是蝶形网路，演算法复杂度的量级是一样的。下面计算演算法复杂度：

设上面N point DFT calculator的乘法次数为C(N)，则两个N/2 point DFT calculator的乘法次数均为C(N/2)，又由G和H计算X另需N次乘法，故有：

C(N)=2*C(N/2 )+N

此递推方程求解如下：

故演算法复杂度为 O(N*log_2N) ，是不是快了不少。

二，蝶形网路（butterfly diagram）

用上一节的递归电路计算4 point DFT，完整展开如下图所示：

简化得：

此即4 point FFT的蝶形网路。

另外可以利用公式 W_{N}^{k+frac{N}{2}}=-W_N^k 对蝶形网路权重作如下变形：

得：

这是4 point FFT蝶形网路的另一种形式。

类似的，8 point FFT蝶形网路（第二种形式）为：

对于给定的point数，蝶形网路是固定，可预计算。在fft的gpu实现中，通常将其生成为lut。

上图中还标出了stage。使用蝶形网路计算FFT，是按stage推进的：N个输入经过第一个stage得到N个中间结果，再输入第二个stage...直至得到最终N个输出。N point FFT有 log_2N 个stage。

三，bitreverse演算法

注意到蝶形网路的N个输入的顺序是打乱的，以8 point蝶形网路为例，可以看到：

x(0)在0号位，x(1)在4号位，x(2)在2号位，x(3)在6号位，x(4)在1号位，x(5)在5号位，x(6)在3号位，x(7)在7位号。

对于一般N point的情况，这个顺序是否可以直接算出来呢？

答案是肯定的，有称为bitreverse的演算法，说：

对于N point蝶形网路，求x(k)在几号位，只需将k化为 log_2N 位二进位数，然后将bit反序，再转回十进位，所得结果即为x(k)所在位号。

以8 point蝶形网路为例，我们求x(3)在几号位，将3化为 log_28=3 位二进位数得011，bit反序得110，将110化回十进位得6，所以x(3)在6号位。

下面是完整列表：（图取自：OpenStax CNX）

此演算法看起来很神奇，但其实是比较容易理解的。

作为事后诸葛，我估计它是这么想出来的：

以上就是用于快速计算DFT的FFT演算法。

四，IFFT

回到海面模型，遗憾的是它并非DFT而是IDFT，所以无法套用FFT演算法。

不过没事儿，比较标准DFT和标准IDFT的表达式：

DFT：

X(k)=sum_{n=0}^{N-1}{x(n)e^{-ifrac{2pi kn}{N}}},kin{{0,1,...,N-1}}

IDFT：

x(n)=frac{1}{N}sum_{k=0}^{N-1}{X(k)e^{ifrac{2pi kn}{N}}},nin{{0,1,...,N-1}}

我们发现，两者很像（脸盲的同学甚至会搞混）。

所以，虽然无法直接套用，前者思路仍可运用于后者。

模仿FFT演算法推导过程重来一遍，可以得到IFFT演算法：

用两个N/2 point IDFT calculator去构造一个N point IDFT calculator。将序号为偶数的输入给到第一个N/2 point IDFT calculator，序号为奇数的输入给到第二个N/2 point IDFT calcuator，如下图所示：

则有：

G(n)=frac{1}{N}sum_{k=0}^{frac{N}{2}-1}{g(k)e^{ifrac{2pi kn}{frac{N}{2}}}}=frac{1}{N}sum_{k=0}^{frac{N}{2}-1}{x(2k)e^{ifrac{2pi kn}{frac{N}{2}}}},nin{{0,1,...,frac{N}{2}-1}}

H(n)=frac{1}{N}sum_{k=0}^{frac{N}{2}-1}{h(k)e^{ifrac{2pi kn}{frac{N}{2}}}}=frac{1}{N}sum_{k=0}^{frac{N}{2}-1}{x(2k+1)e^{ifrac{2pi kn}{frac{N}{2}}}},nin{{0,1,...,frac{N}{2}-1}}

如何用G(n)和H(n)得到x(n)呢？

与前面类似方法可推得：

x(n)=left{egin{matrix} G(n)+W_{N}^{-n}H(n) &,nin{{0,1,...,frac{N}{2}-1}} \ G(n-frac{N}{2})+W_{N}^{-n}H(n-frac{N}{2}) & ,nin{{frac{N}{2},frac{N}{2}+1,...,N-1}} end{matrix} ight.

于是补全电路图：

与前面相同，取N=4将电路彻底展开并简化，得到4 point IFFT蝶形网路：

利用公式 W_{N}^{-n-frac{N}{2}}=-W_N^{-n} 可变形为第二种形式：

8 point IFFT蝶形网路（第二种形式）：

另外，IFFT的bitreverse与FFT相同。

最后由于DFT/IDFT是线性的，所以常数因子并不会影响演算法。

故适用于标准IDFT：

x(n)=frac{1}{N}sum_{n=0}^{N-1}{X(k)e^{ifrac{2pi kn}{N}}},nin{{0,1,...,N-1}}

的IFFT演算法，可不加任何修改地应用于未归一化的IDFT：

x(n)=sum_{n=0}^{N-1}{X(k)e^{ifrac{2pi kn}{N}}},nin{{0,1,...,N-1}}

海面的IDFT模型更接近于后者。

--

下一篇说如何用IFFT计算海面IDFT模型。

下一篇：杨超：fft海面模拟（三）


推荐阅读：
相关文章