离散傅里叶变换

贡献者： addis

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

\begin{matrix} (1) & g_{p} = \frac{1}{\sqrt{N}} \sum_{q = 0}^{N - 1} \exp (- \frac{2 π i}{N} p q) f_{q}, \end{matrix}

\begin{matrix} (2) & f_{q} = \frac{1}{\sqrt{N}} \sum_{p = 0}^{N - 1} \exp (\frac{2 π i}{N} p q) g_{p} . \end{matrix}

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

1. 与傅里叶变换的关系

预备知识 1　傅里叶变换（指数函数）

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

\begin{matrix} (3) & g (k) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{+ \infty} f (x) e^{- i k x} d x, \end{matrix}

\begin{matrix} (4) & f (x) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{+ \infty} g (k) e^{i k x} d k . \end{matrix}

如果

f (x)

和

g (k)

分别只在区间

[- L_{x} / 2, L_{x} / 2]

和

[- L_{k} / 2, L_{k} / 2]

内不为零，积分就可以只在这两个区间内进行。如果函数不为零的区间

[a, b]

并不关于原点对称，那么使用平移公式即可（子节 7 ）。

格点

图 1：

N = 4

和

N = 5

时的

x

格点，

Δ x = L_{x} / N

。

k

轴格点同理。

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

\begin{matrix} (5) & Δ x = L_{x} / N, Δ k = L_{k} / N . \end{matrix}

注意第一个和最后一个格点未必在区间的两端，首尾格点分别相距

L_{x} (N - 1) / N

和

L_{k} (N - 1) / N

。若

N

是奇数，我们令中间的格点为

x_{0}

和

k_{0}

，如果

N

是偶数，我们令中间偏右的格点为

x_{0}

和

k_{0}

，即

\begin{matrix} (6) & {\begin{aligned} x_{n} = n Δ x \\ k_{n} = n Δ k, \end{aligned} n = {\begin{aligned} - N / 2 \dots N / 2 - 1 & (偶数 N) \\ - (N - 1) / 2 \dots (N - 1) / 2 & (奇数 N) . \end{aligned} \end{matrix}

求和代替积分

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

\begin{matrix} (7) & g (k_{p}) \approx \frac{1}{\sqrt{2 π}} \sum_{q} f (x_{q}) e^{- i k_{p} x_{q}} Δ x, \end{matrix}

\begin{matrix} (8) & f (x_{q}) \approx \frac{1}{\sqrt{2 π}} \sum_{p} g (k_{p}) e^{i k_{p} x_{q}} Δ k . \end{matrix}

式中

k_{p} x_{q} = (Δ x Δ k) p q

，对比 DFT（式 1 式 2 ）中的指数项得

\begin{matrix} (9) & Δ x Δ k = \frac{2 π}{N} . \end{matrix}

但我们会发现这里的

p, q

可以是负整数，而 DFT 中的

p, q

都是非负整数。但稍加计算就会发现当

p

或

q

是负值时，把它们加上

N

，指数项并不改变。例如

k_{p + N} x_{q} = k_{p} x_{q} + 2 π

，并不影响指数项。所以，我们只要将所有小于零的格点编号加上

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 ，就只是相差两个常数了：

\begin{matrix} (10) & g (k_{i}) = \sqrt{\frac{N}{2 π}} Δ x \cdot g_{i}, f (x_{i}) = \sqrt{\frac{N}{2 π}} Δ k \cdot f_{i} . \end{matrix}

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

\begin{matrix} (11) & N Δ x Δ k = L_{x} Δ k = L_{k} Δ x = \frac{L_{x} L_{k}}{N} = 2 π . \end{matrix}

离散傅里叶变换只有两个自由度，只要决定

N, L_{x}, Δ x, L_{k}, Δ k

中的任意两个，就可以完全决定变换公式。注意

L_{x}

与

Δ k

成反比，

L_{k}

与

Δ x

成反比，

L_{x} L_{k}

与

N

成正比。

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

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

2. 与傅里叶级数的关系

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

\begin{matrix} (12) & c_{n} = \frac{\sqrt{2 π}}{L_{x}} g (k_{i}) = \frac{1}{\sqrt{N}} g_{i}, \end{matrix}

另外可以发现 DFT 中

k_{n} = 2 π n / L_{x}

的离散值恰好就是傅里叶级数级数中需要的。

3. DFT 变换矩阵

预备知识 2　正交矩阵、酉矩阵

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

\begin{matrix} (13) & g = F f, f = F^{- 1} g . \end{matrix}

其中

\begin{matrix} (14) & F_{p q} = \frac{1}{\sqrt{N}} \exp (- \frac{2 π i}{N} p q) . \end{matrix}

下面证明，

F

是一个酉矩阵，所以逆变换矩阵就是

F

的厄米共轭。

\begin{matrix} (15) & F^{- 1} = F^{†}, \end{matrix}

不难验证式 2 的变换矩阵的确等于

F^{†}

。

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

\begin{matrix} (16) & \sum_{p = 0}^{N - 1} F_{p q_{1}}^{*} F_{p q_{2}} = 0 (q_{1} \neq q_{2}), \end{matrix}

而

\begin{matrix} (17) & \sum_{p = 0}^{N - 1} F_{p q_{1}}^{*} F_{p q_{2}} = \frac{1}{N} \sum_{p = 0}^{N - 1} \exp [\frac{2 π i}{N} (q_{2} - q_{1}) p] . \end{matrix}

注意到求和的每一项在复平面上都对应模长为

1 / N

，辐角为

(q_{2} - q_{1}) p / N

个圆周的矢量，而

N

条矢量恰好向不同方向均匀分布，所以相加为

0

。证毕。

4. 采样定理

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

定理 1　奈奎斯特—香农采样定理（Nyquist–Shannon sampling theorem）

　　如果函数 $f (x)$ 的傅里叶变换 $g (k)$ 不为零的区间不超出 $[- L_{k} / 2, L_{k} / 2]$ （即 $f (x)$ 是有限带宽的），那么我们只需要用步长 $Δ x \leq 2 π / L_{k}$ 来对 $f (x)$ 采样就可以用以下插值公式精确还原出 $f (x)$ ：

\begin{matrix} (18) & f (x) = \sum_{n = - \infty}^{\infty} f (x_{n}) sinc [π (x - x_{n}) / Δ x], \end{matrix}

（见 sinc 函数）和

\begin{matrix} (19) & g (k) = {\begin{aligned} \frac{1}{\sqrt{2 π}} \sum_{n = - \infty}^{\infty} f (x_{n}) e^{- i k x_{n}} Δ x & (| k | ⩽ L_{k} / 2) \\ 0 & (| k | > L_{k} / 2) \end{aligned} . \end{matrix}

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

\begin{matrix} (20) & g (k) = \frac{1}{\sqrt{2 π}} \sum_{n = - \infty}^{\infty} f (x_{n}) \int_{- \infty}^{+ \infty} sinc [π (x - x_{n}) / Δ x] e^{- i k x} d x . \end{matrix}

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

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

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

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

5. 归一化

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

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

\begin{matrix} (21) & \begin{aligned} \int_{- \infty}^{\infty} & f^{*} (x) f (x) d x = \\ \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} f (x_{m})^{*} f (x_{n}) \int_{- \infty}^{\infty} sinc [\frac{π (x - x_{m})}{Δ x}] sinc [\frac{π (x - x_{n})}{Δ x}] d x . \end{aligned} \end{matrix}

果然，这个积分的结果是

\begin{matrix} (22) & \sum_{i = - \infty}^{\infty} {| f_{i} |}^{2} Δ x . \end{matrix}

即函数的模方等于列矢量的模方乘以步长，这对

g (k)

也同样成立。既然 DFT 是精确的，我们预期变换矩阵

F

能保证变换前后列矢量的模长相等，而这恰好是酉矩阵的性质。

6. 插值与添零

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

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

　　然而这么做与 sinc 插值得到的结果并不完全相同。sinc 插值得到的 $f (x)$ 往往会超出原来的区间（因为 sinc 函数只以 $Δ x / x$ 的速度衰减）。所以我们首先需要在误差能接受的范围内取一个更大 $L_{x}$ （即更小的 $Δ k$ ）。由于 sinc 插值函数的带宽严格小于 $L_{k} = 2 π / Δ x$ ，所以 $L_{k}$ 不需要改变。现在我们可以认为两个函数都是有限带宽的了，要得到 “包含所有信息” 的 DFT，就要先保持 $Δ 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}]$ 进行变换。

　　我们先定义两个新函数⁷

\begin{matrix} (23) & f_{1} (x) = f (x + x_{0}) e^{- i k_{0} x}, g_{1} (k) = g (k + k_{0}) e^{i (k + k_{0}) x_{0}} . \end{matrix}

其中

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}

，

Δ x^{'} = 2 π / L_{k}^{'}

，

Δ k^{'} = 2 π / 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}]$ 内⁹。但在程序中，最高效的做法还是先将 $f$ 先变为 $f_{1}$ ，做 FFT 得到 $g_{1}$ 再变为 $g$ 。

　　由采样定理易证，如果 $g (k)$ 只在 $[k_{a}, k_{b}]$ 不为零，那么只需用 $Δ x^{'} = 2 π / L_{k}^{'}$ 对 $f (x)$ 采样就可以用以下插值公式精确还原出¹⁰ $f (x)$ 。

\begin{matrix} (24) & f (x) = \sum_{n = - \infty}^{\infty} f (x_{n}) sinc [\frac{π (x - x_{n})}{Δ x}] e^{i k_{0} (x - x_{n})}, \end{matrix}

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

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

例 1　高斯波包

　　我们现在来用 FFT 数值计算以下高斯波包的傅里叶变换。

\begin{matrix} (25) & f (x) = e^{- x^{2}} e^{100 i x} . \end{matrix}

为了便于验证数值结果，解析的傅里叶变换结果是（参考例 1 ）

\begin{matrix} (26) & g (k) = \frac{1}{\sqrt{2}} e^{- (k - 100)^{2} / 4}, \end{matrix}

我们分别取

x, k

的区间为

[- 5, 5]

和

[90, 110]

（区间外的函数值不超过

1.4 \times 10^{- 11}

）。所以我们采用式 24 的傅里叶变换，令

x_{0} = 0

，

k_{0} = 100

，则

L_{x} ⩾ 10

，

L_{k} ⩾ 20

。为了满足

L_{x} L_{k} = 2 π N

（见式 11 ），

N

取最小值 32。所以不妨令

L_{x} = 10, L_{k} = 2 π N / L_{x} = 20.106 . . .

。所以

Δ x = L_{x} / N = 0.313 . . .

，

Δ k = L_{k} / N = 2 π / L_{x} = 0.628 . . .

。现在我们就可以生成

x, k

的格点并且做 FFT 了。

未完成：比较数值和解析结果

1. ^ 工程上的定义常常是正变换没有 $1 / \sqrt{N}$ 因子，逆变换的 $1 / \sqrt{N}$ 因子变为 $1 / N$ 。这样的好处是节省运算量。我们使用的定义好处是变换为幺正变换，有保持归一化的特点。
2. ^ 参考 Numerical Recipes 3ed
3. ^ 我们可以用反证法做一个不太严谨的证明：如果这种情况存在，那么 $f (x)$ 就可以用有限个 $sinc (x)$ 函数的线性组合来表示，然而从定义上来看 $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}) e^{- i k_{0} (x + x_{0})}$ 和 $g_{1} (k) = g (k + k_{0}) e^{i k x_{0}}$ 。
8. ^ 容易证明，DFT 也精确满足这个性质。
9. ^ 同样，如果 $N$ 是奇数， $x_{0}, k_{0}$ 是中间格点，否则就是中间靠右格点。
10. ^ 对离散的 $f_{1} (x)$ 使用之前的插值公式，再由 $f_{1} (x)$ 求 $f (x)$ 就可以推导出式 24 。

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。