Fourier 方法(02)：FFT 及其 python 实现

本节介绍基于 DTFT 的算法——快速 Fourier 变换, 以及用于实现这一算法的 python 模块.

快速 Fourier 变换(FFT)

FFT 的基本原理其实很简单, 基本上就是 分而治之 (divide and conquer, D&C) 的方法. 为了方便, 将相因子记为 \[ \phi_N:=\exp{(\frac{2\pi i}{N})} \] 那么规模为\(2N\)的 DTFT 可以写为 \[ G_m=\frac{1}{2N}\sum_{n=0}^{2N-1}\phi_{2N}^{mn}g_n \] 将它按照奇数项和偶数项拆开: \[ \begin{aligned} G_m&=\frac{1}{2N}\sum_{j=0}^{N-1}\phi_{2N}^{2mj}g_{2j}+\frac{1}{2N}\sum_{k=0}^{N-1}\phi_{2N}^{m(2k+1)}g_{2k+1} \\ &=\frac{1}{2N}\sum_{j=0}^{N-1}\phi_{N}^{mj}g_{2j}+\frac{1}{2N}\sum_{k=0}^{N-1}\phi_{N}^{mk}g_{2k+1}\exp{(\frac{i\pi m}{N})} \end{aligned} \] 这就是著名的 Danielson-Lanczos 引理, 类似于著名的归并排序算法, 它实际上给出两个规模为\(N\)​的 DTFT, 可以不断向下递归, 每次二分时域样本集, 直到集合很小以至于极其好算, 这样就把\(O(N^2)\)的问题转化为了\(O(N\log N)\)的问题.

numpy.fft模块简介

著名的科学计算工具numpy有一套成熟的函数模块用于计算 FFT. 当然, 它与常使用的公式略有区别. 上节提到, 数学上喜欢将因子全部放在变换公式中, 让反变换不带有系数, 即 \[ \begin{aligned} G_m&=\frac{1}{N}\sum_{n=0}^{N-1}g_n\exp{(-\frac{2\pi imn}{N})} \\ g_n&=\sum_{m}G_m\exp{(\frac{2\pi imn}{N})} \end{aligned} \] 这样的\(G_m\)正是按照复指数函数分解后的系数, 更具有"频谱"的意义. 如果\(g_n\)就是\(x_0+\displaystyle\frac{n}{f_s}\)处的函数值, 那么\(G_m\)代表着空间频率为\(k=\displaystyle\frac{mf_s}{N}\)的分量的复振幅(这里其实还相差一个相因子\(\exp({2\pi ikx_{center}})\), 其中\(x_{center}\)是空时域采样位置的中位数), 科学计算工具未必采用这种约定, 因此编写代码时反复乘除\(N\)因子并不是值得惊奇的事情.

上述式子中, 我们没有指定\(m\)​的求和范围, 是打算说明空间频率的选取具有一定的任意性.

空间频率的选取, 采样定理

按照复指数函数的周期性, 很容易看出 \[ G_{m\pm N}=G_m \] 因此, 所谓的空间频率为\(\displaystyle0,\frac{f_s}{N},\frac{2f_s}{N}, \cdots,\frac{N-1}{N}f_s\)和\(\displaystyle-f_s,-\frac{N-1}{N}f_s,-\frac{N-2}{N}f_s,\cdots,-\frac{f_s}{N}\)是等价的. 从这个无限序列中, 任选\(N\)个在不同等价类中的空间频率, 都能根据反变换复原出全套的\(g_n\). 如果我们执意从\(0\)频率出发向两侧选取频率, 那么频率绝对值中最大的也不会超出\(\displaystyle\frac{f_s}{2}\)​, 这就是著名的 采样定理.

正/反变换

numpy.fft.fft(a, n=None, axis=-1, norm="backward")
numpy.fft.ifft(a, n=None, axis=-1, norm="backward")

它们是最基本的 FFT 函数, 其中

• a: numpy.ndarray表示变换/反变换的对象,
  • 它可能是一个多维数组, axis指定变换所处的维度, 默认取-1即最后一维.
• n: int指定样本长度, 输入None表示根据a的长度来决定. 如果输入具体值, 长则切片, 短则用0补齐.
• norm: str表示归一化, "forward"/"backward"分别指定对变换/反变换得到的矢量归一化, "ortho"则会将\(N\)因子均分给正/反变换.

如果不打算对后三个 params 多加干涉, 那么以下两点需要重点注意:

1. 如果输入的g正好是空时域采样, 则numpy.fft.fft(g)给出的是\(N\cdot G\), 因为默认的归一化方式是backward, 它会把\(N\)因子全部放到反变换里; 同样地, 如果输入的G是特定频率的复振幅, 则numpy.fft.ifft(G)输出的是\(\displaystyle\frac{g}{N}\); 将两个函数复合显然仍能够回到数组自身.

2. 令G = numpy.fft.fft(g) / len(g), 则G[m]未必代表\(\displaystyle\frac{mf_s}{N}\)对应的复振幅, 这是由于返回的数组采用相当奇葩的方式. 给定正指标m, 如果m <= (N - 1) // 2, 那么G[m]的确具有这种含义; 否则, G[m]实际上对应着\(\displaystyle\frac{(m-N)f_s}{N}\). 频率序列详解请看下一部分.

频率序列

numpy.fft.fftfreq(n, d=1.0)

这是numpy.fft.fft对应的频率序列.

• n仍表示样本长度, 这时必须准确地指定.

• d是采样间隔, 即\(\displaystyle\frac{1}{f_s}\), 因此\(\displaystyle\frac{mf_s}{n}\)​​即是m / (n * d).

频率序列的具体形式见如下代码与注释.

d = 1 / fs
freqs = numpy.fft.fftfreq(N, d)
# N = 2 * k: freqs = [0, 1, ..., k - 1, - k, - k + 1, ..., -1] / (N * d)
# N = 2 * k + 1: freqs = [0, 1, ..., k,  - k, - k + 1, ..., -1] / (N * d)

即大约前一半是零频率以及递增的正频率, 大约后一半是递增的负频率(绝对值递减), 如果N是偶数, 负频率项要比正频率项多一项.

有时候, 将序列按照频率从负到正从小到大顺次重排是更有必要的, 这有赖于所谓的 shift 函数.

shift 函数

numpy.fft.fftshift(x, axes=None)
numpy.fft.ifftshift(x, axes=None)

fftshift是能够重排频谱的函数, ifftshift不是ifft的 shift 函数, 而是fftshift的反函数. 下面以前者为例说明参数:

• x输入一个数组, 最好是numpy.fft.fftshift返回的结果
• axes: int | tuple规定作用轴, 默认None表示每根轴都作用

fftshift将\(\{0, 1, ..., k - 1(,- k),- k + 1, ..., -1\}\)重排为\(\{- k + 1, ..., -1,0, 1, ..., k - 1(,- k)\}\), 而ifftshift正相反.