贡献者: addis
离散傅里叶变换(Discrete Fourier Transform)(DFT)是一个复数域的线性变换。对两组有序数列 $f_0, f_1, \dots, f_{N-1}$ 和 $g_0,g_2,\dots, g_{N-1}$,正变换和逆变换分别为1
\begin{equation}
g_p = \frac{1}{\sqrt{N}}\sum_{q=0}^{N-1} \exp\left(-\frac{2\pi \mathrm{i} }{N} p q\right) f_q~,
\end{equation}
\begin{equation}
f_q = \frac{1}{\sqrt{N}}\sum_{p=0}^{N-1} \exp\left(\frac{2\pi \mathrm{i} }{N} p q\right) g_p~.
\end{equation}
一个更常见的名词是
快速傅里叶变换(Fast Fourier Transform)(FFT),其数学定义与离散傅里叶变换一样,只是优化了算法使程序运行更快
2。
1. 与傅里叶变换的关系
在详细分析 DFT 的性质之前,我们先看看它与解析的傅里叶变换如何对应。函数 $f(x)$ 和 $g(k)$ 间的傅里叶变换为
\begin{equation}
g(k) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} f(x) \mathrm{e} ^{- \mathrm{i} kx} \,\mathrm{d}{x} ~,
\end{equation}
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} g(k) \mathrm{e} ^{ \mathrm{i} kx} \,\mathrm{d}{k} ~.
\end{equation}
如果 $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$ 格点,$\Delta 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{equation}
\Delta x = L_x/N ~,\qquad \Delta k = L_k/N~.
\end{equation}
注意第一个和最后一个格点未必在区间的两端,首尾格点分别相距 $L_x(N-1)/N$ 和 $L_k(N-1)/N$。若 $N$ 是奇数,我们令中间的格点为 $x_0$ 和 $k_0$,如果 $N$ 是偶数,我们令中间偏右的格点为 $x_0$ 和 $k_0$,即
\begin{equation}
\left\{\begin{aligned}
&x_n = n \Delta x\\
&k_n = n \Delta k~,
\end{aligned}\right. \qquad n = \left\{\begin{aligned}
&-N/2 \dots N/2 - 1 &\quad &(\text{偶数 } N)\\
&-(N-1)/2 \dots (N-1)/2 &\quad &(\text{奇数 } N)~.
\end{aligned}\right. \end{equation}
求和代替积分
现在我们用求和近似式 3 式 4 的积分得
\begin{equation}
g(k_p) \approx \frac{1}{\sqrt{2\pi}} \sum_q f(x_q) \mathrm{e} ^{- \mathrm{i} k_p x_q} \Delta x~,
\end{equation}
\begin{equation}
f(x_q) \approx \frac{1}{\sqrt{2\pi}} \sum_p g(k_p) \mathrm{e} ^{ \mathrm{i} k_p x_q} \Delta k~.
\end{equation}
式中 $k_p x_q = (\Delta x \Delta k)pq$,对比 DFT(
式 1 式 2 )中的指数项得
\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 和
式 7 式 8 ,就只是相差两个常数了:
\begin{equation}
g(k_i) = \sqrt{\frac{N}{2\pi}} \Delta x \cdot g_i, \qquad
f(x_i) = \sqrt{\frac{N}{2\pi}} \Delta k \cdot f_i~.
\end{equation}
式 9 是 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 前后各将数列的左半右半交换一次。
2. 与傅里叶级数的关系
要用 DFT 计算傅里叶级数 的系数(式 2 ),只需要把式 3 乘以 $\sqrt{2\pi}/L_x$ 即可,所以
\begin{equation}
c_n = \frac{\sqrt{2\pi}}{L_x} g(k_i) = \frac{1}{\sqrt{N}} g_i~,
\end{equation}
另外可以发现 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}} $,变换和逆变换分别记为
\begin{equation}
\boldsymbol{\mathbf{g}} = \boldsymbol{\mathbf{F}} \boldsymbol{\mathbf{f}} ~,\qquad
\boldsymbol{\mathbf{f}} = \boldsymbol{\mathbf{F}} ^{-1} \boldsymbol{\mathbf{g}} ~.
\end{equation}
其中
\begin{equation}
F_{pq} = \frac{1}{\sqrt{N}} \exp\left(-\frac{2\pi \mathrm{i} }{N} p q\right) ~.
\end{equation}
下面证明,$ \boldsymbol{\mathbf{F}} $ 是一个
酉矩阵,所以逆变换矩阵就是 $ \boldsymbol{\mathbf{F}} $ 的厄米共轭。
\begin{equation}
\boldsymbol{\mathbf{F}} ^{-1} = \boldsymbol{\mathbf{F}} ^\dagger ~,
\end{equation}
不难验证
式 2 的变换矩阵的确等于 $ \boldsymbol{\mathbf{F}} ^\dagger $。
根据酉矩阵的定义,我们需要证明
\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 \left[\frac{2\pi \mathrm{i} }{N} (q_2-q_1) p \right] ~.
\end{equation}
注意到求和的每一项在复平面上都对应模长为 $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)$ 是有限带宽的),那么我们只需要用步长 $\Delta x \leq 2\pi/L_k$ 来对 $f(x)$ 采样就可以用以下插值公式精确还原出 $f(x)$:
\begin{equation}
f(x) = \sum_{n = -\infty}^{\infty} f(x_n) \operatorname{sinc} [\pi(x - x_n)/\Delta x]~,
\end{equation}
(见
sinc 函数)和
\begin{equation}
g(k) = \left\{\begin{aligned}
&\frac{1}{\sqrt{2\pi}} \sum_{n = -\infty}^{\infty} f(x_n) \mathrm{e} ^{- \mathrm{i} kx_n} \Delta x \quad &( \left\lvert k \right\rvert \leqslant L_k/2)\\
&0 & ( \left\lvert k \right\rvert > L_k/2)
\end{aligned}\right. ~.
\end{equation}
注意把插值公式式 18 代入傅里叶变换式 3 就可直接得到式 19 :
\begin{equation}
g(k) = \frac{1}{\sqrt{2\pi}} \sum_{n = -\infty}^{\infty} f(x_n) \int_{-\infty}^{+\infty} \operatorname{sinc} [\pi(x - x_n)/\Delta x] \mathrm{e} ^{- \mathrm{i} k x} \,\mathrm{d}{x} ~.
\end{equation}
采样定理给出的 $\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)$ 都是有限带宽的?严格来说这种情况是不存在的3,但许多时候我们可以近似认为函数在区间外等于零,其中典型的例子就是高斯分布的傅里叶变换(例 1 )。
5. 归一化
我们知道傅里叶变换可以保持函数的模长不变,或者说保持函数的归一化(式 16 )。
假设采样定理的条件成立,我们来看如果通过离散点计算函数的模长。以 $f(x)$ 为例,用插值公式来计算模长的平方,得
\begin{equation} \begin{aligned}
\int_{-\infty}^{\infty} &f^*(x) f(x) \,\mathrm{d}{x} = \\
&\sum_{m = -\infty}^{\infty}\sum_{n = -\infty}^{\infty} f(x_m)^*f(x_n) \int_{-\infty}^{\infty} \operatorname{sinc} \left[\frac{\pi(x - x_m)}{\Delta x} \right] \operatorname{sinc} \left[\frac{\pi(x - x_n)}{\Delta x} \right] \,\mathrm{d}{x} ~.
\end{aligned} \end{equation}
果然,这个积分的结果是
\begin{equation}
\sum_{i = -\infty}^{\infty} \left\lvert f_i \right\rvert ^2\Delta x~.
\end{equation}
即函数的模方等于列矢量的模方乘以步长,这对 $g(k)$ 也同样成立。既然 DFT 是精确的,我们预期变换矩阵 $ \boldsymbol{\mathbf{F}} $ 能保证变换前后列矢量的模长相等,而这恰好是酉矩阵的性质。
6. 插值与添零
上面我们看到,如果 $f(x)$ 是有限带宽的,就可以使用式 18 进行插值,称为 sinc 插值4。对 $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 变换即可5。
然而这么做与 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,那添零法将会更快6。
7. 任意区间的 DFT
上面提到 DFT 或 FFT 的另一个限制就是只能在以零点为中心的两个区间之间变换 $f(x)$ 和 $g(k)$。但事实上我们只需要把数据稍作处理就可以在两个任意的区间 $[x_a, x_b]$ 和 $[k_a, k_b]$ 进行变换。
我们先定义两个新函数7
\begin{equation}
f_1(x) = f(x + x_0) \mathrm{e} ^{- \mathrm{i} k_0 x}, \qquad
g_1(k) = g(k + k_0) \mathrm{e} ^{ \mathrm{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$。根据傅里叶变换的平移性质(
式 14 )
8,$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]$ 内9。但在程序中,最高效的做法还是先将 $ \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)$ 采样就可以用以下插值公式精确还原出10 $f(x)$。
\begin{equation}
f(x) = \sum_{n = -\infty}^{\infty} f(x_n) \operatorname{sinc} \left[\frac{\pi(x - x_n)}{\Delta x} \right] \mathrm{e} ^{ \mathrm{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) = \mathrm{e} ^{-x^2} \mathrm{e} ^{100 \mathrm{i} x}~.
\end{equation}
为了便于验证数值结果,解析的傅里叶变换结果是(参考
例 1 )
\begin{equation}
g(k) = \frac{1}{\sqrt2} \mathrm{e} ^{-(k-100)^2/4}~,
\end{equation}
我们分别取 $x, k$ 的区间为 $[-5, 5]$ 和 $[90, 110]$(区间外的函数值不超过 $1.4 \times 10^{-11} $)。所以我们采用
式 24 的傅里叶变换,令 $x_0 = 0$,$k_0 = 100$,则 $L_x \geqslant 10$,$L_k \geqslant 20$。为了满足 $L_x L_k = 2\pi N$(见
式 11 ),$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. ^ 我们可以用反证法做一个不太严谨的证明:如果这种情况存在,那么 $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 。