Fourier Transform 2 (DTFT, DFT)

离散时间傅里叶变换 (DTFT)

傅里叶级数和傅里叶变换面向的都是连续函数 \(f(x)\)，但在用计算机处理之前，我们必须将其离散化，因此涉及到了采样操作。

冲激串及其傅里叶变换

上一篇文章 Fourier Transform 1 提到，Dirac \(\delta\) 函数能帮助我们采样一个函数值。进一步地，如果要以 \(h\) 为间隔等距采样一系列函数值，就可以用如下冲激串函数： \[ s(x)=\sum_{n=-\infty}^{\infty}\delta(x-nh)\tag{1-1} \] 下图显示了一个冲激函数与冲激串函数：

下面我们求冲激串的傅里叶变换。注意冲激串是周期函数，不是绝对可积的，因此需要用到上一篇文章最后提到的技巧来计算。首先将 \(s(x)\) 展开为傅里叶级数： \[ s(x)=\sum_{n=-\infty}^{\infty}c_ne^{j2\pi nx/h}\tag{1-2} \] 其中： \[ c_n=\frac{1}{h}\int_{-h/2}^{h/2}s(x)e^{-j2\pi nx/h}\mathrm dx=\frac{1}{h}\int_{-h/2}^{h/2}\sum_{m=-\infty}^{\infty}\delta(x-mh)e^{-j2\pi nx/h}\mathrm dx=\frac{1}{h}e^0=\frac{1}{h}\tag{1-3} \]

于是 \(s(x)\) 的傅里叶级数为： \[ s(x)=\frac{1}{h}\sum_{n=-\infty}^{\infty}e^{j2\pi nx/h}\tag{1-4} \] 根据傅里叶变换的线性性，我们只需要求每一项的傅里叶变换： \[ S(u)=\mathscr F\{s(x)\}=\frac{1}{h}\sum_{n=-\infty}^{\infty}\mathscr F\{e^{j2\pi nx/h}\}\tag{1-5} \] 在上一篇文章中我们计算过 Dirac \(\delta\) 函数的傅里叶变换 \(\mathscr F\{\delta(x-x_0)\}=e^{-jux_0}\)，那么根据傅里叶变换的对偶性，有： \[ \mathscr F\{e^{-jxu_0}\}=\delta(-u-u_0)\xrightarrow{a=-u_0}\mathscr F\{e^{jxa}\}=\delta(-u+a)=\delta(u-a)\tag{1-6} \] 因此： \[ S(u)=\frac{1}{h}\sum_{n=-\infty}^{\infty}\mathscr F\{e^{j2\pi nx/h}\}=\frac{1}{h}\sum_{n=-\infty}^{\infty}\delta\left(u-\frac{2\pi n}{h}\right)\tag{1-7}\label{fourier-s} \]

这就是 \(s(x)\) 的傅里叶变换结果。

采样定理

假设每间隔 \(h\) 对 \(f(x)\) 做一次采样，得到序列 \(\{f_n=f(nh)\}\)，其中 \(n\in\{\ldots,-2,-1,0,1,2,\ldots\}\)，如下图所示：

根据上一节的内容，我们可以用冲激串函数将这个离散的序列写作连续的函数的形式： \[ \{f_n=f(nh)\}\longleftrightarrow\tilde f(x)=f(x)s(x)=\sum_{n=-\infty}^{\infty}f(x)\delta(x-nh)\tag{1-8} \] 换句话说，\(\tilde f(x)\) 就是采样后的函数。于是，一个自然的问题就是：采样操作对函数的傅里叶变换有什么影响呢？也就是说，采样后的函数的傅里叶变换与原函数的傅里叶变换有什么关系呢？

根据卷积定理，由于 \(\tilde f(x)\) 是 \(f(x)\) 与 \(s(x)\) 的乘积，因此其傅里叶变换就是 \(F(u)\) 与 \(S(u)\) 的卷积，其中 \(S(u)\) 已经在上一节中计算出来了，即 \(\eqref{fourier-s}\) 式。于是： \[ \begin{align} \tilde F(u)&=F(u)\ast S(u)\\ &=\int_{-\infty}^{\infty}F(\tau)\frac{1}{h}\sum_{n=-\infty}^{\infty}\delta\left(u-\tau-\frac{2\pi n}{h}\right)\mathrm d\tau\\ &=\frac{1}{h}\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}F(\tau)\delta\left(u-\tau-\frac{2\pi n}{h}\right)\mathrm d\tau\\ &=\frac{1}{h}\sum_{n=-\infty}^{\infty}F\left(u-\frac{2\pi n}{h}\right) \end{align}\tag{1-9}\label{DTFT1} \] 由此可以看出，\(\tilde F(u)\) 是一个以 \(\frac{2\pi}{h}\) 为周期的连续周期函数，它是 \(F(u)\) 的多份平移后的副本的叠加。

特别地，假设 \(f(x)\) 是一个有限带宽的函数，即其傅里叶变换 \(F(u)\) 仅在一个有限范围 \([-u_\max,u_\max]\) 内不为零。那么：

• 欠采样：周期 \(\frac{2\pi}{h}<2u_\max\)，平移后的 \(F(u)\) 之间就会发生重叠，导致相加后在重叠区间内的 \(F(u)\) 丢失了，这种现象称作混叠 (aliasing)；
• 过采样：周期 \(\frac{2\pi}{h}>2u_\max\)，那么各个平移后的 \(F(u)\) 不仅没有重叠，还隔开了一段距离，这段距离上 \(\tilde F(u)=0\)；
• 临界采样：周期 \(\frac{2\pi}{h}=2u_\max\)，各个平移后的 \(F(u)\) 刚好不会发生重叠。

下图展示了一个例子，第一行是一个有限带宽函数的傅里叶变换结果，第二、三、四行分别展示了过采样、临界采样和欠采样后，函数的傅里叶变换结果：

因此，要想重建出原本的有限带宽函数，我们要避免欠采样，也即采样率需要满足： \[ f_\text{sampling}=\frac{1}{h}\geq2\cdot\frac{u_\max}{2\pi}=2 f_\max\tag{1-10} \] 其中 \(f_\max\) 表示信号的最大频率，这就是采样定理。当然，如果原信号本身不是有限带宽的，即不存在最大频率 \(f_\max\)，自然就不可能在采样后完美重建。好在现实中大多数信号可以近似为有限带宽的，因此只要采样率充分高，就能将重建误差控制在可接受范围内。

另外，从上面的图可以看出，当 \(\tilde f(x)\) 满足采样定理时，只需要拿出 \(\tilde F(u)\) 的一个周期并放大 \(h\) 倍，就能还原 \(F(u)\)： \[ F(u)=H(u)\tilde F(u),\quad\text{where}\;H(u)=\begin{cases}h&-u_\max\leq u\leq u_\max\\0&\text{otherwise}\end{cases}\tag{1-11}\label{DTFT2} \] 这里 \(H(u)\) 相当于是一个理想低通滤波器。

离散时间傅里叶变换

由于 \(\tilde F(u)\) 是 \(f(x)\) 先采样（即在时域上离散化）、再做傅立叶变换得到的，因此我们称之为离散时间傅里叶变换。在上一节中，我们建立起了 \(\tilde F(u)\) 与 \(F(u)\) 之间的关系，即 \(\eqref{DTFT1}\) 和 \(\eqref{DTFT2}\) 式，发现了 \(\tilde F(u)\) 其实就是 \(F(u)\) 的多份平移副本的叠加，并自然引出了采样定理。这一节中，我们欲显式地写出 \(\tilde F(u)\) 与 \(f_n\) 之间的变换关系。

对 \(\tilde f(x)\) 做傅里叶变换，得： \[ \begin{align} \tilde F(u)&=\int_{-\infty}^{\infty}\tilde f(x)e^{-jux}\mathrm dx=\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}f(x)\delta(x-nh)e^{-jux}\mathrm dx\\ &=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}f(x)\delta(x-nh)e^{-jux}\mathrm dx=\sum_{n=-\infty}^{\infty}f(nh)e^{-jnhu}=\sum_{n=-\infty}^{\infty}f_ne^{-jnhu} \end{align}\tag{1-12}\label{DTFT} \]

这就是离散时间傅里叶正变换。

值得注意的是，由于 \(\tilde F(u)\) 是周期函数，因此不能做傅立叶逆变换，否则结果发散。注意到，如果将 \(\eqref{DTFT}\) 式改写一下： \[ \tilde F(u)=\sum_{n=-\infty}^{\infty}f_{-n}e^{jnhu}\tag{1-13}\label{DTFT3} \] 对比傅里叶级数表达式： \[ f(x)=\sum_{n=-\infty}^{\infty}c_ne^{jn\omega x},\quad\quad c_n=\frac{1}{2l}\int_{-l}^{l}f(x)e^{-jn\omega x}\mathrm dx\tag{1-14}\label{FS} \] 可以发现二者有着相同的形式，也就是说，\(\eqref{DTFT3}\) 式其实就是函数 \(\tilde F(u)\) 的傅里叶级数展开。于是有： \[ f_{-n}=\frac{h}{2\pi}\int_{-\pi/h}^{\pi/h}\tilde F(u)e^{-jnhu}\mathrm du\tag{1-15} \] 也即： \[ f_n=\frac{h}{2\pi}\int_{-\pi/h}^{\pi/h}\tilde F(u)e^{jnhu}\mathrm du\tag{1-16}\label{IDTFT} \] 这就是离散时间傅里叶逆变换。

小结

在开启下一节之前，我们先对这一节的内容做一个小结。在从傅里叶变换推导到离散时间傅里叶变换的过程中，我们对函数的时域表示进行了采样，发现它的频域表示从非周期函数变成了周期函数。那么，基于傅里叶变换在形式上的对称性，容易想到，如果我们对函数的频域表示进行采样，那么它的时域表示应该也会从非周期函数变成周期函数——这不就是从傅里叶变换推导回了傅里叶级数吗？事实上，对比傅里叶级数 \(\eqref{FS}\) 和离散时间傅立叶变换 \(\eqref{DTFT},\eqref{IDTFT}\)，可以发现二者在形式上的确是对称的。

值得注意的是，时域和频域在周期性/非周期性和离散性/连续性上有着固定的对应关系——时域上周期（非周期）对应频域上离散（连续）；时域上离散（连续）对应频域上周期（非周期）。对于傅立叶级数、傅立叶变换和离散时间傅立叶变换，分别有：

• 傅立叶级数：时域周期连续、频域离散非周期
• 傅立叶变换：时域非周期连续、频域连续非周期
• 离散时间傅立叶变换：时域非周期离散、频域连续周期

于是，剩下的最后一种组合就呼之欲出了，即时域和频域上均是离散周期函数，这就是离散傅立叶变换。进一步地，我们还能立刻想到两种推导离散傅立叶变换的路径——对离散时间傅立叶变换的频域表示采样、或对傅立叶级数的时域表示采样，如下图所示：

离散傅立叶变换 (DFT)

尽管离散时间傅立叶变换在时域上做了离散化，但是在频域上是连续的，计算机依旧无法处理。同理，傅立叶级数频域是离散的，但时域是连续的，也无法处理。因此，我们希望构建时域、频域均为离散序列之间的变换关系，即离散傅立叶变换。

从离散时间傅立叶变换推导

考虑一个长为 \(M\) 的序列 \(\{f_0,f_1,\ldots,f_{M-1}\}\)，可以等价地将其视作在 \(n>M-1\) 和 \(n<0\) 时取值为零的无限长序列： \[ f_n=\begin{cases}f_n,&0\leq n<M\\0,&\text{otherwise}\end{cases}\tag{2-1} \] 对其进行离散时间傅立叶变换得： \[ F(u)=\sum_{n=-\infty}^{\infty}f_ne^{-jnhu}=\sum_{n=0}^{M-1}f_ne^{-jnhu}\tag{2-2} \] 由前文可知这是一个连续的、以 \(\frac{2\pi}{h}\) 为周期的函数。考虑在每个周期等间隔采样 \(M\) 个点，即间隔为 \(h'=\frac{2\pi}{hM}\)，得到序列 \(\{F_m=F(mh')\}\)，那么代入上式得： \[ F_m=F(mh')=\sum_{n=0}^{M-1}f_ne^{-jmnhh'}=\sum_{n=0}^{M-1}f_ne^{-j2\pi mn/M}\tag{2-3} \] 这就是离散傅立叶正变换。

倘若将序列 \(\{F_m=F(mh')\}\) 写作连续函数形式： \[ \{F_m=F(mh')\}\longleftrightarrow \tilde F(u)=F(u)s(u)=\sum_{m=-\infty}^{\infty}F(u)\delta(u-mh') \] 那么代入离散时间傅立叶逆变换 \(\eqref{IDTFT}\) 式： \[ \begin{align} f_n&=\frac{h}{2\pi}\int_{-\pi/h}^{\pi/h}\tilde F(u)e^{jnhu}\mathrm du\\ &=\frac{h}{2\pi}\int_{-\pi/h}^{\pi/h}\sum_{m=-\infty}^{\infty}F(u)\delta(u-mh')e^{jnhu}\mathrm du\\ &=\frac{h}{2\pi}\sum_{m=-\infty}^{\infty}\int_{-\pi/h}^{\pi/h}F(u)\delta(u-mh')e^{jnhu}\mathrm du\\ &=\frac{h}{2\pi}\sum_{-\frac{\pi}{hh'}\leq m\leq\frac{\pi}{hh'}}F(mh')e^{jmnhh'}\\ &=\frac{1}{M}\sum_{-\frac{M}{2}\leq m\leq\frac{M}{2}}F_me^{j2\pi mn/M}\\ &=\frac{1}{M}\sum_{m=0}^{M-1}F_me^{j2\pi mn/M}\\ \end{align}\tag{2-4} \] 这就是离散傅立叶逆变换。

为了与后文的符号统一，我们用 \(F(u)\) 和 \(f(x)\) 来代替 \(F_m\) 和 \(f_n\)，就得到了离散傅里叶变换对： \[ \begin{align} &F(u)=\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M},\quad u=0,\ldots,M-1\tag{2-5}\label{DFT}\\ &f(x)=\frac{1}{M}\sum_{u=0}^{M-1}F(u)e^{j2\pi ux/M},\quad x=0,\ldots,M-1\tag{2-6}\label{IDFT} \end{align} \]

矩阵形式

将 \(F(u)\) 和 \(f(x)\) 都视作 \(m\) 维向量，则式 \(\eqref{DFT},\eqref{IDFT}\) 可以写作矩阵形式。具体而言，记 \(\omega=e^{-j2\pi/M}\) 为 \(M\) 次单位复数根，构建矩阵： \[ \mathcal F=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega^{1}&\omega^{2}&\cdots&\omega^{(M-1)}\\ 1&\omega^{2}&\omega^{4}&\cdots&\omega^{2(M-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega^{(M-1)}&\omega^{2(M-1)}&\cdots&\omega^{(M-1)^2} \end{bmatrix}\tag{2-7} \] 简记作 \(\mathcal F_{nm}=\omega^{nm}\)，其中 \(n,m\) 从 0 开始计数。于是离散傅里叶变换写作： \[ F=\mathcal F f\tag{2-8} \] 对于逆变换，可以验证 \(\mathcal F\mathcal F^H=MI\)，其中 \(\mathcal F^H\) 表示 \(\mathcal F\) 的共轭转置矩阵，因此 \(\mathcal F^{-1}=\frac{1}{M}\mathcal F^H\)，故离散傅里叶逆变换写作： \[ f=\frac{1}{M}\mathcal F^HF\tag{2-9} \] 当然，这也可以直接从式 \(\eqref{IDFT}\) 中看出来。

离散傅里叶变换的性质

记 \(\mathscr F\{f(x)\}=F(u),\,\mathscr F\{g(x)\}=G(u)\)，DFT 具有以下常用性质。

线性性：

\[ \begin{align} &\mathscr F\{af(x)+bg(x)\}=a\mathscr F\{f(x)\}+b\mathscr F\{g(x)\}=F(u)+G(u)\\ &\mathscr F^{-1}\{aF(u)+bG(u)\}=a\mathscr F^{-1}\{F(u)\}+b\mathscr F^{-1}\{G(u)\} \end{align} \]

平移性质：

\[ \begin{align} &\mathscr F\{f(x-x_0)\}=e^{-j2\pi ux_0/M}F(u)\\ &\mathscr F^{-1}\{F(u-u_0)\}=e^{j2\pi u_0x/M}f(x) \end{align} \] 特别地，当 \(x_0=u_0=M/2\) 时，有： \[ \begin{align} &\mathscr F\{f(x-M/2)\}=(-1)^uF(u)\\ &\mathscr F^{-1}(F(u-M/2))=(-1)^xf(x) \end{align} \]

周期性：

\[ \begin{align} &F(u)=F(u+kM)&&k\in\mathbb Z\\ &f(x)=f(x+kM)&&k\in\mathbb Z \end{align} \]

共轭对称性：若 \(f(x)\) 是实函数，则 \(F(u)=\mathscr F\{f(x)\}\) 是共轭对称的（实部偶函数，虚部奇函数），即： \[ F^\ast(u)=F(-u) \] 若 \(f(x)\) 是虚函数，则 \(F(u)=\mathscr F\{f(x)\}\) 是共轭反对称的（实部奇函数，虚部偶函数），即： \[ F^\ast(u)=-F(-u) \]

循环卷积定理：

\[ \begin{align} &\mathscr F\{f(x)\ast g(x)\}=F(u)G(u)\\ &\mathscr F\{f(x)g(x)\}=\frac{1}{M}F(u)\ast G(u) \end{align} \] 值得注意的是，由于离散傅里叶变换具有周期性，因此这里的卷积是循环卷积，即先将 \(f(x)\) 与 \(g(x)\) 进行周期延拓后再做卷积。在下一篇文章中我们将看到，在图像处理中这种卷积模式并不是我们想要的，因此需要做一些额外的操作。

参考资料

1. Rafael C. Gonzalez. Digital Image Processing, Fourth Edition. ↩︎
2. Wikipedia. Discrete-time Fourier transform. https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform ↩︎
3. 彻底搞懂傅里叶变换之实用干货分享(四)-离散傅里叶变换(DFT) - anders的文章 - 知乎 https://zhuanlan.zhihu.com/p/405143684 ↩︎
4. Half_Kettle. 从傅里叶级数（Fourier series）到离散傅里叶变换（Discrete Fourier transform）. https://www.cnblogs.com/yang-ding/p/15925430.html ↩︎
5. 《数字图像处理》图像表征：离散傅里叶变换（DFT）、离散余弦变换（DCT）、主成分分析（PCA） - zhiwei的文章 - 知乎 https://zhuanlan.zhihu.com/p/563668048 ↩︎