前言
最近学数值分析这门课的时候碰到了 DFT 与 FFT,然而我们学校 PPT 上给出的 FFT 看起来实在是令人费解,教材更是直接摆烂不写了。
于是我想在网上找找教程看看,又发现网上的教程跟我们教材的思路又不完全相同,其看起来似乎是一个更加易懂的改进,所以我想在此简单记录一下。仅解析原理,不提供具体算法实现。
解析
我们考虑 \(f\) 在等距点集 \(\left\{x_j=\frac{2\pi j}N|j=0,1,...,N-1\right\}\) 上用正交函数族 \(1,e^{ix},...,e^{i(N-1)x}\) 得到的最佳平方逼近: \[ S(x)=\sum_{k=0}^{N-1}c_ke^{ikx} \] 其中 \[ c_k=\frac1N \sum_{j=0}^{N-1}f(x_j)e^{-ik\frac{2\pi j}N}\;(k=0,1,...,N-1) \] 被称为 DFT,要想得到所有 \(c_k\) 共需要计算 \(N^2\) 次乘法和加法,FFT 给出了该计算的快速实现。
思路
方便起见,记 \(f_j=f(x_j),\:\omega_N=e^{-i\frac{2\pi}N}\),那么 \(c_k=\frac1N \sum_{j=0}^{N-1}f_j\omega_N^{kj}\)
注意到,\(\omega_N^{kj}\) 事实上只有 \(1,\omega_N,...,\omega_N^{N-1}\) 这 \(N\) 个值,这实际上是复平面上单位圆的 \(N\) 等分点
因而不难发现,若 \(N\) 为偶数,由对称性就有 \(\omega_N^k+\omega_N^{\frac N2+k}=0\)
一个自然的想法是合并公因子,所以我们可以将 \(c_k\) 改写为: \[ c_k=\frac1N\left(\sum_{j=0}^{\frac N2-1}f_j\omega_N^{kj}+\sum_{j=0}^{\frac N2-1}f_{\frac N2+j}\omega_N^{k(\frac N2+j)}\right)=\frac1N\left(\sum_{j=0}^{\frac N2-1}f_j\omega_N^{kj}+\sum_{j=0}^{\frac N2-1}(-1)^kf_{\frac N2+j}\omega_N^{kj}\right)=\frac1N\sum_{j=0}^{\frac N2-1}(f_j+(-1)^kf_{\frac N2+j})\omega_N^{kj} \] 具体考虑 \(k\) 为奇数或者偶数的情况,有: \[ c_{2k}=\frac1N\sum_{j=0}^{\frac N2-1}(f_j+f_{\frac N2+j})\omega_{\frac N2}^{kj}\\ c_{2k+1}=\frac1N\sum_{j=0}^{\frac N2-1}(f_j-f_{\frac N2+j})\omega_N^{j}\omega_{\frac N2}^{kj} \] 从而如果 \(\frac N2\) 仍为偶数,我们就可以对 \(c_{2k}\) 与 \(c_{2k+1}\) 重复上述过程
特别地,若 \(N=2^p\),则该过程可以一直重复下去
改进
一个简单的改进是直接考察 \(c_k\) 中的奇数项和 \(S_k\) 与偶数项和 \(T_k\): \[ c_k=\frac1N \sum_{j=0}^{N-1}f_j\omega_N^{kj}=S_k+T_k\\ S_k=\frac1N \sum_{j=0}^{\frac N2-1}f_{2j+1}\omega_N^k\omega_{\frac N2}^{kj}\\ T_k=\frac1N \sum_{j=0}^{\frac N2-1}f_{2j}\omega_{\frac N2}^{kj} \] 那么由于 \(\frac N2\) 仍为偶数,可以对 \(S_k\) 与 \(T_k\) 反复作相同操作,而这个改进的关键在于: \[ c_{k+\frac N2}=\frac1N \sum_{j=0}^{N-1}f_j\omega_N^{(k+\frac N2)j}=\frac1N \sum_{j=0}^{N-1}(-1)^jf_j\omega_N^{kj}=T_k-S_k \] 由此我们事实上只需要计算 \(c_0,c_1,...,c_{\frac N2-1}\),对于 \(S_k\) 与 \(T_k\) 亦然
参考文献
[1] 李庆扬, 王能超, 易大义. 数值分析第 5 版 [M]. 清华大学出版社, 2009. P87
[2] 快速理解 FFT 算法. 2022. https://zhuanlan.zhihu.com/p/407885496