Fourier Transform 3 (2D DFT, Image)

二维离散傅里叶变换

分量形式

设 \(f(x,y)\) 是一个二元函数，其中 \(x\in\{0,\ldots,M-1\},\,y\in\{0,\ldots,N-1\}\). 类似于一维的离散傅里叶变换和逆变换，二维离散傅里叶变换为： \[ F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(ux/M+vy/N)},\quad u=0,\ldots,M-1,\;v=0,\ldots,N-1\label{2D-DFT} \] 逆变换为： \[ f(x,y)=\frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(ux/M+vy/N)},\quad x=0,\ldots,M-1,\;y=0,\ldots,N-1\label{2D-IDFT} \]

矩阵形式

同一维情形类似，我们也可以将二维 DFT 写作矩阵形式。视 \(F(u,v),f(u,v)\) 为 \(M\times N\) 的矩阵，记 \(\omega_M=e^{-j2\pi/M},\,\omega_N=e^{-j2\pi/N}\) 分别为 \(M\) 次单位复数根和 \(N\) 次单位复数根，构建矩阵： \[ \mathcal F_M=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_M^{1}&\omega_M^{2}&\cdots&\omega_M^{(M-1)}\\ 1&\omega_M^{2}&\omega_M^{4}&\cdots&\omega_M^{2(M-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_M^{(M-1)}&\omega_M^{2(M-1)}&\cdots&\omega_M^{(M-1)^2} \end{bmatrix},\quad \mathcal F_N=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega_N^{1}&\omega_N^{2}&\cdots&\omega_N^{(N-1)}\\ 1&\omega_N^{2}&\omega_N^{4}&\cdots&\omega_N^{2(N-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_N^{(N-1)}&\omega_N^{2(N-1)}&\cdots&\omega_N^{(N-1)^2} \end{bmatrix} \] 则二维离散傅里叶变换可以写作： \[ F=\mathcal F_Mf\mathcal F_N^T \] 进一步地，利用拉直算子 \(\text{vec}(\cdot)\) 以及 Kronecker 积的性质 \(\text{vec}(AXB)=(B^T\otimes A)\text{vec}(X)\)，上式也可写作： \[ \text{vec}(F)=(\mathcal F_N\otimes \mathcal F_M)\text{vec}(f) \] 该形式与一维情形具有统一的形式，能方便后续的推导。

性质

二维 DFT 的大部分性质可类比一维 DFT 的性质推广而来。

平移性质：设 \(\mathscr F\{f(x,y)\}=F(u,v)\)，则： \[ \begin{align} &\mathscr F\{f(x-x_0,y-y_0)\}=e^{-j2\pi(ux_0/M+vy_0/N)}F(u,v)\\ &\mathscr F^{-1}\{F(u-u_0,v-v_0)\}=e^{j2\pi (u_0x/M+v_0y/N)}f(x,y) \end{align} \]

旋转性质：使用极坐标：\(x=r\cos\theta,\,y=r\sin\theta,\,u=\omega\cos\varphi,\,v=\omega\sin\varphi\)，可得如下变换对： \[ f(r,\theta+\theta_0)\iff F(\omega,\varphi+\theta_0) \] 即若 \(f(x,y)\) 旋转 \(\theta_0\) 角度，则 \(F(u,v)\) 也旋转相同的角度，反之亦然。

周期性： \[ \begin{align} &F(u,v)=F(u+k_1M,v+k_2N)&&k_1,k_2\in\mathbb Z\\ &f(x,y)=f(x+k_1M,y+k_2N)&&k_1,k_2\in\mathbb Z \end{align} \]

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

频谱和功率谱：由于 DFT 是复函数，因此可以用极坐标表示： \[ \begin{align} &F(u,v)=R(u,v)+jI(u,v)=|F(u,v)|e^{j\phi(u,v)}\\ \text{where}\quad&|F(u,v)|=\sqrt{R^2(u,v)+I^2(u,v)}\\ &\phi(u,v)=\arctan\frac{I(u,v)}{R(u,v)} \end{align} \] 其中 \(|F(u,v)|\) 称为幅度谱，\(\phi(u,v)\) 称为相位谱，二者统称为频谱。另外，功率谱定义为幅度谱的平方：

\[ P(u,v)=|F(u,v)|^2=R^2(u,v)+I^2(u,v) \] 前文说到，对于实函数 \(f(x,y)\)，其傅里叶变换 \(F(u,v)\) 是共轭对称的，于是容易知道其频谱是关于原点偶对称的，而相位谱是关于原点奇对称的： \[ |F(u,v)|=|F(-u,-v)|\quad\quad\phi(u,v)=-\phi(-u,-v) \]

二维循环卷积定理：设 \(\mathscr F\{f(x,y)\}=F(u,v),\,\mathscr F\{g(x,y)\}=G(u,v)\)，则： \[ \begin{align} &\mathscr F\{f(x,y)\ast g(x,y)\}=F(u,v)\cdot G(u,v)\\ &\mathscr F\{f(x,y)\cdot g(x,y)\}=\frac{1}{MN}F(u,v)\ast G(u,v) \end{align} \]

其中卷积依旧是循环卷积。

图像的频域滤波

图像的空间域与频域

前文中我们研究的一元函数可以看作是时域信号，但是对于图像而言，时域信号这个说法显然不合适了。那应该怎么理解呢？

首先，一幅数字图像其实是一个二元函数 \(f=f(x,y)\)，其中自变量 \((x,y)\) 是像素的坐标，函数值是像素的灰度级。因此，将图像在三维空间中绘制出来，就好像一个连绵不断的山脉，类似于信号波形的二维版本。只不过这里自变量表示空间位置而非时间，因此我们称之为空间域。

与时域信号类似，我们可以用无数个二维的不同角频率的正余弦函数去近似一幅图像。其中，图像变化比较剧烈的地方——比如边缘、纹理等细节，需要依靠高频的正余弦函数去拟合；相反，图像大体的构图就对应低频的正余弦函数。因此，我们经常把图像的细节称作高频分量，而整体构图称作低频分量。

对图像做二维离散傅里叶变换就得到了其频域表示，将频域表示的频谱可视化出来就是频谱图。频域滤波指在频谱图上使用一个滤波器对图像做滤波的过程。然而，出于可视化、浮点误差等考虑，下面我们关注对图像做频域滤波时一些值得注意的方面。

循环卷积与交叠误差

根据循环卷积定理，在空间域中做循环卷积等价于在频域中做乘法。然而，我们在图像处理中使用的卷积通常是 padding 0 的卷积，两种卷积方式导致的结果差别就是所谓的交叠误差。

幸运的是交叠误差很容易解决。设 \(f(x,y)\) 和 \(h(x,y)\) 分别是 \(A\times B\) 和 \(C\times D\) 的图像，则选取 \(P\geq A+C-1\) 和 \(Q\geq B+D-1\)，将 \(f\) 和 \(h\) 都往右下填充 0 直到 \(P\times Q\) 大小，那么在频域中处理（对应循环卷积）与直接 padding 0 的卷积结果就一样了。

频谱图的可视化：中心化与对数变换

我们知道 \(F(0,0)\) 表示频率最低的分量，其值其实是图像所有位置的亮度之和： \[ F(0,0)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(0x/M+0y/N)}=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y) \] 因此在频谱图上，这个点的亮度是最大的。又由于周期性，所以事实上四个角的亮度都很大，这不利于我们观察频谱图。因此，我们常常对频谱图做中心化，即将 \(F(0,0)\) 平移到 \(F(M/2,N/2)\) 的位置。根据傅里叶变换的平移性质，这相当于在做傅里叶变换之前先对 \(f(x,y)\) 乘以 \((-1)^{x+y}\).

另外，由于 \(F(0,0)\) 实在是太大了，当我们可视化出频谱图后其他地方的亮度很低，不利于观察，因此我们常常在可视化之前做如下对数变换： \[ F_\text{visualize}(u,v)=\log(1+|F(u,v)|) \]

浮点误差

理论上，傅里叶变换和逆变换是一对可逆操作，但是由于计算机存在浮点误差，实际操作中数值也许会稍有扰动。例如，数字图像原本是实值函数，在傅里叶变换到频域后执行某些滤波操作，然后傅里叶逆变换后理应得到实值函数，但实际中可能存在很小的虚部。因此，我们常常需要在逆变换后增加一个取实部的操作。

频域滤波步骤

根据上文的讨论，在对图像进行频域滤波时有如下步骤：

1. 对于 \(M\times N\) 的图像 \(f(x,y)\)，取 \(P=2M,Q=2N\)；
2. 将原图零填充到 \(P\times Q\) 大小，得到 \(f_p(x,y)\)；
3. 对 \(f_p(x,y)\) 乘以 \((-1)^{x+y}\)，为频谱图的中心化做准备；
4. 计算二维离散傅里叶变换 \(F(u,v)\)；
5. 构造大小为 \(P\times Q\) 的频域滤波器 \(H(u,v)\)；
6. 进行频域滤波：\(G(u,v)=H(u,v)F(u,v)\)；
7. 计算傅里叶逆变换，取实部并将 \((-1)^{x+y}\) 乘回来：\(g_p(x,y)=\left(\text{real}\left[\mathscr F^{-1}\{G(u,v)\}\right]\right)(-1)^{x+y}\)；
8. 取 \(g_p(x,y)\) 左上角 \(M\times N\) 区域 \(g(x,y)\)，得到最终结果。

参考资料

1. Rafael C. Gonzalez. Digital Image Processing, Fourth Edition. ↩︎