Fourier 方法(01)：从连续到离散

考虑函数\(g(x)\)​的 Fourier 变换以及逆变换 \[ G(k)=\int_{-\infty}^{+\infty}g(x)e^{-2\pi ikx}dx \\ g(x)=\int_{-\infty}^{+\infty}G(k)e^{2\pi ikx}dk \] 数值计算中, 很多时候连续函数也是用离散的形式表达的, 积分不是针对一个解析式进行符号计算, 而是输入一个有限规模的样本来近似替代积分, 例如 \[ \{g(x_0+\frac{n}{f_s})|n\in\mathbb{N}\ ,n<N\} \] 因此我们实际计算的并不是\(g(x)\)本身的变换, 而是如下分段函数的变换 \[ \bar{g}(x)=\begin{aligned}\begin{cases} 0\qquad & x\in(-\infty,x_0) \\ g(x_0)\qquad &\displaystyle x\in[x_0,x_0+\frac{1}{f_s}) \\\displaystyle g(x_0+\frac{1}{f_s})\qquad &\displaystyle x\in[x_0+\frac{1}{f_s}, x_0+\frac{2}{f_s}) \\\vdots\\\displaystyle g(x_0+\frac{N-1}{f_s})\qquad &\displaystyle x\in[x_0+\frac{N-1}{f_s}, x_0+\frac{N}{f_s}] \\0\qquad &\displaystyle x\in(x_0+\frac{N}{f_s},+\infty) \end{cases}\end{aligned} \] 相应的频谱为 \[ \bar{G}(k)=\frac{1-\exp{(\displaystyle-\frac{2\pi ik}{f_s})}}{2\pi ik}e^{-2\pi ikx_0} \sum_{n=0}^{N-1} g(x_0+\frac{n}{f_s})e^{-\frac{2\pi ikn}{f_s}} \] 如果我们将\(g(x_0+\displaystyle\frac{n}{f_s})\)简记为\(g_n\), 同时\(k\)也不取连续变量, 令 \[ k=\frac{mf_s}{N}\ ,\ m=0,1,\cdots,N-1 \\ G_m:=\frac{2\pi ike^{2\pi ikx_0} \bar{G}(k)}{1-\exp{(\displaystyle-\frac{2\pi ik}{f_s})}}\Bigg|_{k=\frac{mf_s}{N}} \] 也就给出 \[ G_m=\sum_{n=0}^{N-1}g_ne^{-\frac{2\pi imn}{N}} \] 这就是所谓的 时域离散 Fourier 变换(DTFT), 它的反变换如何? 可以猜测为 \[ \mathcal{F}^{-1}[G]_n=\sum_{m=0}^{N-1}G_me^{\frac{2\pi imn}{N}}=\sum_{m,n'}g_{n'}\exp{(\frac{2\pi im(n-n')}{N})} \] 幸好时域样本\(g_{n'}\)不依赖于\(m\), 可以按照等比序列的方法对\(m\)求和, 即 \[ \mathcal{F}^{-1}[G]_n=\sum_{n'=0}^{N-1}g_{n'}\frac{1-e^{iN\phi}}{1-e^{i\phi}}\ ,\ \phi=\frac{2\pi(n-n')}{N} \] 显然, \(N\phi\)总是导向\(2\pi\)的整数倍, 因此分子总是\(0\), 只要分数不构成\(\displaystyle\frac{0}{0}\)型未定式, 就给出零结果, 又\(|n-n'|\le N-1\), 因此当且仅当\(n'=n\)时构成未定式, \(n'\ne n\)时都给出零结果.

现在只需要关注\(n'=n\)项, 退回计算等比序列和之前, 立即得到 \[ \mathcal{F}^{-1}[G]_n=\sum_{m=0}^{N-1}g_n=Ng_n \] 我们最终得到完整的变换-反变换式: \[ \begin{aligned} G_m&=\sum_{n=0}^{N-1}g_n\exp{(-\frac{2\pi imn}{N})} \\ g_n&=\frac{1}{N}\sum_{m=0}^{N-1}G_m\exp{(\frac{2\pi imn}{N})} \end{aligned} \] \(N\)因子依然可以随意分配, 也就是说可以改写为 \[ \begin{aligned} G_m&=\frac{1}{N}\sum_{n=0}^{N-1}g_n\exp{(-\frac{2\pi imn}{N})} \\ g_n&=\sum_{m=0}^{N-1}G_m\exp{(\frac{2\pi imn}{N})} \end{aligned} \] 这样写的好处是, 数学上使用反变换更多, 或者说将信号 Fourier 分解为频域上的谱, 频谱前面最好不带有系数, 但是前一种写法在程序中也不少见.

像这样计算一次完整的变换或反变换, 需要对\(m\)和\(n\)二重循环, 给出\(O(N^2)\)的时间复杂度. 下面将要介绍的 快速 Fourier 变换(FFT) 是一种可以简化至\(O(N\log N)\)​的算法.