离散傅里叶变换

                     

贡献者： addis

　　 离散傅里叶变换（Discrete Fourier Transform）（DFT）是一个复数域的线性变换。对两组有序数列 $f_0, f_1, \dots, f_{N-1}$ 和 $g_0,g_2,\dots, g_{N-1}$，正变换和逆变换分别为¹

一个更常见的名词是快速傅里叶变换（Fast Fourier Transform）（FFT），其数学定义与离散傅里叶变换一样，只是优化了算法使程序运行更快²。

1. 与傅里叶变换的关系

　　 在详细分析 DFT 的性质之前，我们先看看它与解析的傅里叶变换如何对应。函数 $f(x)$ 和 $g(k)$ 间的傅里叶变换为

如果 $f(x)$ 和 $g(k)$ 分别只在区间 $[-L_x/2, L_x/2]$ 和 $[-L_k/2, L_k/2]$ 内不为零，积分就可以只在这两个区间内进行。如果函数不为零的区间 $[a,b]$ 并不关于原点对称，那么使用平移公式即可（子节 7 ）。

格点

　　 为了做离散傅里叶变换，我们给两个区间划出 $N$ 个等间距的格点（图 1 ）$\dots, x_{-1}, x_0, x_1,\dots$ 和 $\dots, k_{-1}, k_0, k_1,\dots$ 令 $x_0 = k_0 = 0$，并规定相邻格点的间距为

注意第一个和最后一个格点未必在区间的两端，首尾格点分别相距 $L_x(N-1)/N$ 和 $L_k(N-1)/N$。若 $N$ 是奇数，我们令中间的格点为 $x_0$ 和 $k_0$，如果 $N$ 是偶数，我们令中间偏右的格点为 $x_0$ 和 $k_0$，即

求和代替积分

　　 现在我们用求和近似式 3 式 4 的积分得

式中 $k_p x_q = (\Delta x \Delta k)pq$，对比 DFT（式 1 式 2 ）中的指数项得 但我们会发现这里的 $p, q$ 可以是负整数，而 DFT 中的 $p, q$ 都是非负整数。但稍加计算就会发现当 $p$ 或 $q$ 是负值时，把它们加上 $N$，指数项并不改变。例如 $k_{p+N} x_q = k_p x_q + 2\pi$，并不影响指数项。所以，我们只要将所有小于零的格点编号加上 $N$ 并重新排列即可。例如 $N = 4$ 时，$x$ 格点为 $x_{-2}, x_{-1}, x_0, x_1$，令 $x_2 = x_{-2}, x_3 = x_{-1}$，这四个格点的名字就变为 $x_0, x_1, x_2, x_3$，数值不变。现在再来对比 DFT 和式 7 式 8 ，就只是相差两个常数了：

　　 式 9 是 DFT 一个重要的性质。结合式 5 得

离散傅里叶变换只有两个自由度，只要决定 $N, L_x, \Delta x, L_k, \Delta k$ 中的任意两个，就可以完全决定变换公式。注意 $L_x$ 与 $\Delta k$ 成反比，$L_k$ 与 $\Delta x$ 成反比，$L_xL_k$ 与 $N$ 成正比。

　　 总结起来，要用 DFT 数值近似函数 $f(x)$ 的傅里叶变换，就先用上述方法生成的格点将该函数等间距采样，然后左半和右半调换（负脚标加上 $N$）得到 $f_i$，代入 DFT 公式（在程序中使用 FFT）得到 $g_i$ 再左半和右半调换得到 $g(k)$ 的离散点。反傅里叶变换同理。

　　 下文如无特殊说明，默认 DFT 前后各将数列的左半右半交换一次。

2. 与傅里叶级数的关系

　　 要用 DFT 计算傅里叶级数 的系数（式 2 ），只需要把式 3 乘以 $\sqrt{2\pi}/L_x$ 即可，所以

另外可以发现 DFT 中 $k_n = 2\pi n/L_x$ 的离散值恰好就是傅里叶级数级数中需要的。

3. DFT 变换矩阵

　　 DFT 是一个线性变换，我们来看变换矩阵的性质。把变换和逆变换的系数矩阵用 $ \boldsymbol{\mathbf{F}} $ 和 $ \boldsymbol{\mathbf{F}} ^{-1}$ 来表示，令列矢量 $ \boldsymbol{\mathbf{f}} = (f_0\ f_1\dots f_{N-1}) ^{\mathrm{T}} $，$ \boldsymbol{\mathbf{g}} = (g_0\ g_1\dots g_{N-1}) ^{\mathrm{T}} $，变换和逆变换分别记为

其中 下面证明，$ \boldsymbol{\mathbf{F}} $ 是一个酉矩阵，所以逆变换矩阵就是 $ \boldsymbol{\mathbf{F}} $ 的厄米共轭。 不难验证式 2 的变换矩阵的确等于 $ \boldsymbol{\mathbf{F}} ^\dagger $。

　　 根据酉矩阵的定义，我们需要证明

而 注意到求和的每一项在复平面上都对应模长为 $1/N$，辐角为 $(q_2-q_1)p/N$ 个圆周的矢量， 而 $N$ 条矢量恰好向不同方向均匀分布，所以相加为 $0$。证毕。

4. 采样定理

　　 从以上的分析中，DFT 只是一种近似，且有种种限制，例如我们只能取关于原点对称的区间，又例如变换只有两个自由度。我们不禁想定义更广义的离散傅里叶变换，使式 7 和式 8 中 $x_i$ 和 $k_j$ 可以分别取任意区间的任意多个等间距格点。但事实上，这样的定义并不比 DFT 更精确，而且不能用 FFT 算法提高运算速度。

　　 注意把插值公式式 18 代入傅里叶变换式 3 就可直接得到式 19 ：

　　 采样定理给出的 $\Delta x$ 往往比一般情况下用求和代替积分所需的要小，例如对于 $ \sin\left(ax\right) $ 或 $ \cos\left(ax\right) $，每周期只需要 2 个采样即可！ 所以采样定理一个更简单表述是：

　　 若 $f(t)$ 中不含有 $B$ 赫兹以上的频率，那么只需用至少 $2B$ 的频率对它采样就能将其完全将确定。

　　 有时候我们并不知道 $f(x)$ 的带宽，如何决定采样点的最大步长 $\Delta x$ 呢？一个简单的方法是先选一个较小的 $\Delta x$，使 $f(x_i)$ 散点看起来较为连续。这时解析傅里叶变换就可以用 DFT 近似。得到 $g(k_j)$ 散点后，再判断合理的 $L_k$ 并计算相应的 $\Delta x$ 即可。

　　 由于傅里叶变化和反变换是对称的，所以以上的分析中将 $f(x)$ 和 $g(k)$ 调换同样成立。一个有趣的问题是，是否存在一些情况使 $f(x)$ 和 $g(x)$ 都是有限带宽的？严格来说这种情况是不存在的³，但许多时候我们可以近似认为函数在区间外等于零，其中典型的例子就是高斯分布的傅里叶变换（例 1 ）。

5. 归一化

　　 我们知道傅里叶变换可以保持函数的模长不变，或者说保持函数的归一化（式 16 ）。

　　 假设采样定理的条件成立，我们来看如果通过离散点计算函数的模长。以 $f(x)$ 为例，用插值公式来计算模长的平方，得

果然，这个积分的结果是 即函数的模方等于列矢量的模方乘以步长，这对 $g(k)$ 也同样成立。既然 DFT 是精确的，我们预期变换矩阵 $ \boldsymbol{\mathbf{F}} $ 能保证变换前后列矢量的模长相等，而这恰好是酉矩阵的性质。

6. 插值与添零

　　 上面我们看到，如果 $f(x)$ 是有限带宽的，就可以使用式 18 进行插值，称为 sinc 插值⁴。对 $g(k)$ 的插值也同理。然而在实践中，由于直接用 sinc 插值计算量太大，我们往往用下面介绍的添零法代替。

　　 假设我们要对 $N$ 个散点 $f(x_i)$ 插值，可以先做 DFT 得到 $N$ 个 $g(k_j)$，由式 11 可知，要缩小 $\Delta x$ 而保持 $L_x$ 不变，就要保持 $\Delta k$ 不变且增加 $L_k$ 和 $N$。所以我们可以在 $g(k_j)$ 数列的两边添加相同数量的 0，然后再做反 DFT 变换即可⁵。

　　 然而这么做与 sinc 插值得到的结果并不完全相同。sinc 插值得到的 $f(x)$ 往往会超出原来的区间（因为 sinc 函数只以 $\Delta x/x$ 的速度衰减）。所以我们首先需要在误差能接受的范围内取一个更大 $L_x$（即更小的 $\Delta k$）。由于 sinc 插值函数的带宽严格小于 $L_k = 2\pi/\Delta x$，所以 $L_k$ 不需要改变。现在我们可以认为两个函数都是有限带宽的了，要得到 “包含所有信息” 的 DFT，就要先保持 $\Delta x$ 不变，对 $f(x_i)$ 在新区间补零，做 DFT 得到 $g(k_j)$。现在两个函数的地位完全平等，都可以用式 18 和式 19 中的任意一个来精确插值。而注意补零插值法就相当于式 19 。所以这时再对 $g(k_j)$ 补任意多的零，就可以得到 $f(x)$ 任意密度的插值。在误差范围内，现在得到的 $f(x)$ 就与 sinc 插值得到的相等了。

　　 如果直接按照 DFT 公式计算添零法，程序的速度未必比 sinc 插值要快，但如果用 FFT，那添零法将会更快⁶。

7. 任意区间的 DFT

　　 上面提到 DFT 或 FFT 的另一个限制就是只能在以零点为中心的两个区间之间变换 $f(x)$ 和 $g(k)$。但事实上我们只需要把数据稍作处理就可以在两个任意的区间 $[x_a, x_b]$ 和 $[k_a, k_b]$ 进行变换。

　　 我们先定义两个新函数⁷

其中 $x_0 = (x_a + x_b)/2$， $k_0 = (k_a + k_b)/2$。另外定义 $L'_x = x_b - x_a$，$L'_k = k_b-k_a$，$\Delta x' = 2\pi/L'_k$，$\Delta k' = 2\pi/L'_x$。根据傅里叶变换的平移性质（式 14 ）⁸，$f_1(x)$ 的傅里叶变换就是 $g_1(k)$。注意 $f_1(x)$ 和 $g_1(k)$ 的 DFT 区间 $[-L_x/2, L_x/2]$ 和 $[-L_k/2, L_k/2]$ 分别关于原点对称。现在我们只需要做 $f_1(x)$ 和 $g_1(k)$ 间的 DFT，就能得到 $f(x)$ 和 $g(k)$ 间的离散变换，不妨称为广义 DFT。

　　 经过简单的推导可以发现广义 DFT 仍然可以用式 7 和式 8 表示，只是 $x$ 和 $k$ 的 $N$ 个格点分别落在 $[x_a, x_b]$ 和 $[k_a, k_b]$ 内⁹。但在程序中，最高效的做法还是先将 $ \boldsymbol{\mathbf{f}} $ 先变为 $ \boldsymbol{\mathbf{f}} _1$，做 FFT 得到 $ \boldsymbol{\mathbf{g}} _1$ 再变为 $ \boldsymbol{\mathbf{g}} $。

　　 由采样定理易证，如果 $g(k)$ 只在 $[k_a, k_b]$ 不为零，那么只需用 $\Delta x' = 2\pi/L'_k$ 对 $f(x)$ 采样就可以用以下插值公式精确还原出¹⁰ $f(x)$。

　　 实际应用中同样也可以用补零法来插值。

　　 广义 DFT 的好处是什么呢？例如有时候非零区间 $[x_a, x_b]$ 或（和）$[k_a, k_b]$ 很窄但却远离原点，则 $L_k$ 或 $L_x$ 仍需取较大的值，这时 DFT 就没有发挥最大效率，因为大部分格点的函数值都是 0。但如果用广义 DFT，就可以大大减少格点数。

---

1. ^ 工程上的定义常常是正变换没有 $1/\sqrt{N}$ 因子，逆变换的 $1/\sqrt{N}$ 因子变为 $1/N$。这样的好处是节省运算量。我们使用的定义好处是变换为幺正变换，有保持归一化的特点。
2. ^ 参考 Numerical Recipes 3ed
3. ^ 我们可以用反证法做一个不太严谨的证明：如果这种情况存在，那么 $f(x)$ 就可以用有限个 $ \operatorname{sinc} (x)$ 函数的线性组合来表示，然而从定义上来看 $ \operatorname{sinc} (x)$ 函数在 $x$ 的区间外不可能恒为零。
4. ^ 注意实现的时候只能对有限项求和。
5. ^ 不难证明如果新的长度是 $N$ 的整数倍，反变换后在老格点处仍能得到同样的 $f(x_i)$。
6. ^ 如果 FFT 中只有部分格点不为零，那么 FFT 理论上可以变得更快，称为 pruned FFT，然而 FFTW 的相关页面中介绍，只要非零格点的个数大于 1%，pruned FFT 都不会有显著的性能提升。
7. ^ 式 23 也可以改为 $f_1(x) = f(x + x_0) \mathrm{e} ^{- \mathrm{i} k_0 (x+x_0)}$ 和 $g_1(k) = g(k + k_0) \mathrm{e} ^{ \mathrm{i} k x_0}$。
8. ^ 容易证明，DFT 也精确满足这个性质。
9. ^ 同样，如果 $N$ 是奇数，$x_0, k_0$ 是中间格点，否则就是中间靠右格点。
10. ^ 对离散的 $f_1(x)$ 使用之前的插值公式，再由 $f_1(x)$ 求 $f(x)$ 就可以推导出式 24 。