图

离散傅里叶变换

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

\begin{align} g_p &= \frac{1}{\sqrt{N}}\sum_{q=0}^{N-1} \expRound{-\frac{2\pi\I}{N} p q } f_q \\ f_q &= \frac{1}{\sqrt{N}}\sum_{p=0}^{N-1} \expRound{\frac{2\pi\I}{N} p q} g_p \end{align}

   一个更常见的名词是快速傅里叶变换(Fast Fourier Transform)(FFT), 其定义与离散傅里叶变换一样,只是优化了算法使程序运行更快2

离散傅里叶变换与傅里叶变换

预备知识 傅里叶变换(指数函数)

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

\begin{align} g(k) &= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} f(x)\E^{-\I kx} \dd{x}\\ f(x) &= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} g(k)\E^{\I kx} \dd{k} \end{align}
如果 $f(x)$ 和 $g(k)$ 分别只在区间 $[-L_x/2, L_x, 2]$ 和 $[-L_k/2, L_k/2]$ 内不为零, 积分就可以只在这两个区间内进行. 我们再给两个区间划出 $N$ 个等间距的格点 $\dots, x_{-1}, x_0, x_1,\dots$ 和 $\dots, k_{-1}, k_0, k_1,\dots$ ($x_0 = k_0 = 0$), 并规定相邻格点的间距为
\begin{equation} \Delta x = L_x/N \qquad \Delta k = L_k/N \end{equation}
注意 $x$ 和 $k$ 的首尾格点分别相距 $L_x(N-1)/N$ 和 $L_k(N-1)/N$. 若 $N$ 是奇数, 我们令中间的格点为 $x_0$ 和 $k_0$, 如果 $N$ 是偶数, 我们令中间靠右的格点为 $x_0$ 和 $k_0$.

   现在我们用求和近似上面的积分得3

\begin{align} g(k_p) &= \frac{1}{\sqrt{2\pi}} \sum_q f(x_q) \E^{-\I k_p x_q} \Delta x\\ f(x_q) &= \frac{1}{\sqrt{2\pi}} \sum_p g(k_p) \E^{\I k_p x_q} \Delta k \end{align}
式中 $k_p x_q = (\Delta x \Delta k)pq$, 对比 DFT 中的指数项, 得
\begin{equation} \Delta x\Delta k = \frac{2\pi}{N} \end{equation}
但我们会发现这里的 $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 和 式 6 式 7 , 就只是相差两个常数因子而已了.

   式 8 是 DFT 一个重要的性质. 结合式 5

\begin{equation} N\Delta x \Delta k = L_x \Delta k = L_k \Delta x = \frac{L_x L_k}{N} = 2\pi \end{equation}
离散傅里叶变换只有两个自由度, 只要决定 $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 前后各将数列的左半右半交换一次.

变换矩阵

预备知识 酋矩阵

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

\begin{equation} \bvec g = \mat F \bvec f \qquad \bvec f = \mat F^{-1} \bvec g \end{equation}
其中
\begin{equation} F_{pq} = \frac{1}{\sqrt{N}} \expRound{-\frac{2\pi\I}{N} p q} \end{equation}
下面证明, $\mat F$ 是一个酋矩阵, 所以逆变换矩阵就是 $\mat F$ 的厄米共轭.
\begin{equation} \mat F^{-1} = \mat F\Her \end{equation}
不难验证式 2 的变换矩阵的确等于 $\mat F\Her$.

   根据酋矩阵的定义,我们需要证明

\begin{equation} \sum_{p=0}^{N-1} F^*_{pq_1} F_{pq_2} = 0 \quad (q_1 \ne q_2) \end{equation}
\begin{equation} \sum_{p=0}^{N-1} F^*_{pq_1}F_{pq_2} = \frac1N \sum_{p=0}^{N-1} \exp[\frac{2\pi\I}{N} (q_2-q_1) p] \end{equation}
注意到求和的每一项在复平面上都对应模长为 $1/N$, 幅角为 $(q_2-q_1)p/N$ 个圆周的矢量, 而 $N$ 条矢量恰好向不同方向均匀分布,所以相加为 $0$.证毕.

采样定理(Sampling Theorem)

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

   这里给出傅里叶变换的采样定理, 采样定理的详细表述是: 如果 $g(k)$ 不超出 $[-L_k/2, L_k/2]$(即 $f(x)$ 是有限带宽的), 那么我们只需要用 $\Delta x = 2\pi/L_k$ ( DFT 恰好满足这个条件)来对 $f(x)$ 采样就可以用以下插值公式精确还原出 $f(x)$4

\begin{equation} f(x) = \sum_{n = -\infty}^{\infty} f(x_n)\sinc[\pi(x - x_n)/\Delta x] \end{equation}
其中 $\sinc x = \sin x/x$ (且定义 $\sinc 0 = 1$).

   现在我们将插值公式代入傅里叶变换得

\begin{equation} g(k) = \frac{1}{\sqrt{2\pi}} \sum_{n = -\infty}^{\infty} f(x_n) \int_{-\infty}^{+\infty}\sinc[\pi(x - x_n)/\Delta x] \E^{-\I k x} \dd{x} \end{equation}
结果是
\begin{equation} g(k) = \leftgroup{ \frac{1}{\sqrt{2\pi}} &\sum_{n = -\infty}^{\infty} f(x_n) \E^{-\I kx_n} \Delta x \qquad &(\abs{k} \leqslant L_k/2)\\ &0 & (\abs{k} > L_k/2) } \end{equation}
当 $k = k_p$ 时, 这正是式 6 , 所以这里用求和代替积分是没有误差的. 不仅如此, 该式也可以作为 $g(k)$ 的插值公式.

   有时候我们并不知道 $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)$ 都是有限带宽的? 严格来说这种情况是不存在的5, 但许多时候我们可以近似认为函数在区间外等于零, 其中最常见的就是高斯波包的傅里叶变换.

归一化

   我们知道傅里叶变换可以保持函数的模长不变, 或者说保持函数的归一化. 假设采样定理的条件成立, 我们来看如果通过离散点计算函数的模长. 以 $f(x)$ 为例, 用插值公式来计算模长的平方, 得

\begin{equation}\ali{ \int_{-\infty}^{\infty} &f^*(x) f(x) \dd{x} = \\ &\sum_{m = -\infty}^{\infty}\sum_{n = -\infty}^{\infty} f(x_m)^*f(x_n) \int_{-\infty}^{\infty} \sinc[\pi(x - x_m)/\Delta x]\sinc[\pi(x - x_n)/\Delta x]\dd{x} }\end{equation}
果然, 这个积分的结果是
\begin{equation} \sum_{i = -\infty}^{\infty} \abs{f_i}^2\Delta x \end{equation}
即函数的模方等于列矢量的模方乘以步长, 这对 $g(k)$ 也同样成立. 既然 DFT 是精确的, 我们预期变换矩阵 $\mat F$ 能保证变换前后列矢量的模长相等, 而这恰好是酋矩阵的性质.

插值与添零

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

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

   然而这么做与 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)$. 现在两个函数的地位完全平等, 都可以用式 15 式 17 中的任意一个来精确插值. 而注意补零插值法就相当于式 17 . 所以这时再对 $g(k_j)$ 补任意多的零, 就可以得到 $f(x)$ 任意密度的插值. 在误差范围内, 现在得到的 $f(x)$ 就与 sinc 插值得到的相等了.

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

任意区间的 DFT

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

   我们先定义两个新函数9

\begin{equation} f_1(x) = f(x + x_0)\E^{-\I k_0 x} \qquad g_1(k) = g(k + k_0)\E^{\I (k+k_0) x_0} \end{equation}
其中 $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$. 根据傅里叶变换的平移性质10,$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 仍然可以用式 6 式 7 表示, 只是 $x$ 和 $k$ 的 $N$ 个格点分别落在 $[x_a, x_b]$ 和 $[k_a, k_b]$ 内11. 但在程序中, 最高效的做法还是先将 $\bvec f$ 先变为 $\bvec f_1$, 做 FFT 得到 $\bvec g_1$ 再变为 $\bvec g$.

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

\begin{equation} f(x) = \sum_{n = -\infty}^{\infty} f(x_n)\sinc[\pi(x - x_n)/\Delta x] \E^{\I k_0 (x-x_n)} \end{equation}

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

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

例1 高斯波包

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

\begin{equation} f(x) = \E^{-x^2}\E^{100 \I x} \end{equation}
为了便于验证数值结果, 解析的傅里叶变换结果是
\begin{equation} g(k) = \frac{1}{\sqrt2}\E^{-(k-100)^2/4} \end{equation}
我们分别取 $x, k$ 的区间为 $[-5, 5]$ 和 $[90, 110]$ (区间外的函数值不超过 $1.4\times 10^{-11}$). 所以我们采用式 21 的傅里叶变换, 令 $x_0 = 0$, $k_0 = 100$, 则 $L_x \geqslant 10$, $L_k \geqslant 20$. 为了满足 $L_x L_k = 2\pi N$ (见式 9 ), $N$ 取最小值 32. 所以不妨令 $L_x = 10, L_k = 2\pi N/L_x = 20.106...$. 所以 $\Delta x = L_x/N = 0.313...$, $\Delta k = L_k/N = 2\pi/L_x = 0.628...$. 现在我们就可以生成 $x, k$ 的格点并且做 FFT 了.


1. 工程上的定义常常是正变换没有 $1/\sqrt{N}$ 因子,逆变换的 $1/\sqrt{N}$ 因子变为 $1/N$. 这样的好处是节省运算量.本书中使用的定义好处是变换为幺正变换,有保持归一化的特点.
2. 参考 Numerical Recipes 3ed
3. 事实上, 学习了下文的采样定理会发现这里用求和代替积分是没有任何误差的.
4. 其实我并不确定是否需要 $x_0 = 0$ 或者 $x_i$ 的位置有什么其他要求, 现在姑且认为 $x_0 = 0$ 好了.
5. 我们可以用反证法做一个不太严谨的证明: 如果这种情况存在, 那么 $f(x)$ 就可以用有限个 $\sinc(x)$ 函数的线性组合来表示, 然而从定义上来看 $\sinc(x)$ 函数在 $x$ 的区间外不可能恒为零.
6. 注意实现的时候只能对有限项求和.
7. 不难证明如果新的长度是 $N$ 的整数倍, 反变换后在老格点处仍能得到同样的 $f(x_i)$.
8. 如果 FFT 中只有部分格点不为零, 那么 FFT 理论上可以变得更快, 称为 pruned FFT, 然而 FFTW 的相关页面中介绍, 只要非零格点的个数大于 1%, pruned FFT 都不会有显著的性能提升.
9. 式 20 也可以改为 $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}$
10. 容易证明, DFT 也精确满足这个性质.
11. 同样, 如果 $N$ 是奇数, $x_0, k_0$ 是中间格点, 否则就是中间靠右格点.
12. 对离散的 $f_1(x)$ 使用之前的插值公式, 再由 $f_1(x)$ 求 $f(x)$ 就可以推导出式 21

致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利