离散傅里叶变换

                     

贡献者: 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. 与傅里叶变换的关系

预备知识 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 变换矩阵

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

   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

                     

© 小时科技 保留一切权利