在一维信号处理中,我们利用傅里叶变换实现信号从时域到频域的变换。现在在二维图像信号中,我们学习了空间频率的概念,自然可以想到,可以利用二维傅里叶变换实现空间到空间频率的转换。注意关于频谱中心化和可视化在文末有提及到。
1. 1D Fourier Transform
回顾一下一维傅里叶变换: \[
\begin {aligned}
X(j \omega) &= \int_{-\infty}^{\infty} x(t) e^{-j \omega t}dt \\
X(t) &= \cfrac {1} {2 \pi} \int_{-\infty}^{\infty} X(j \omega) e^{j \omega x} d\omega
\end {aligned}
\] 在实际计算中,我们一般没有用上面的模拟角频率 \(\omega\) ,单位为 radians/second ,而是用频率 \(f = \omega / 2\pi\) 表示 (单位为 : cycles/seconds),如下图所示,非常完美的对称性,除了符号的不同。 \[
\begin {aligned}
X(f) &= \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft}dt \\
x(t) &= \int_{-\infty}^{\infty} X(f) e^{j 2 \pi ft} df
\end {aligned}
\]
2. 2D Fourier Transform
\[
\begin {aligned}
F(u, v) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-j2\pi (ux + vy)} dx dy \\
f(x, y) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u, v) e^{j 2 \pi (ux + vy)} du dv
\end {aligned}
\]
这里 \(u, v\) 均为空间频率 (spatial frequency)
利用 \(\delta(x, y) = \delta(x) \delta(y)\) 以及 \(\delta\) 函数的欧拉公式积分表达 证明:
\[
\begin {aligned}
f(x, y) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u, v) e^{j 2 \pi (ux + vy)} du dv \\
& = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \big ( \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(m, n) e^{-j2\pi (um + vn)} dm dn \big )e^{j 2 \pi (ux + vy)} du dv \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(m, n) \big ( \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{j 2 \pi [u(x-m)+v(y-n)]} du dv \big ) dm dn \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(m, n) \big ( \int_{-\infty}^{\infty} e^{j 2 \pi u(x-m)}du \int_{-\infty}^{\infty} e^{j 2\pi v(y-n)} dv \big ) dm dn \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(m, n) \delta(x-m) \delta(y-n)dm dn \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(m, n) \delta(x-m, y-n)dm dn
\end {aligned}
\] 性质:
- \(F(u, v)\) 通常为复数,即 \(F(u, v) = F_R(u, v) + j F_I(u, v)\)
-
\(|F(u, v)| = \sqrt {F_R^2(u, v) + F_I^2(u, v)}\) 是幅值频谱 (magnitude spectrum)
-
\(\arctan(F_I(u, v)/F_R(u, v))\) 是相位频谱 (phase spectrum)
-
功率谱 (power spectrum): \(P(u, v) = |F(u, v)|^2 = F_R^2(u, v) + F_I^2(u, v)\)
-
共轭 (conjugacy ): \(f^{\ast}(x, y) \Leftrightarrow F(-u, -v)\)
\(F(u, v) = F^{\ast}(-u, -v)\) -
线性 (linearity) \[
\alpha f(x, y) + \beta g(x, y) \Leftrightarrow \alpha F(u, v) + \beta G(u, v)
\] -
伸缩特性 (scaling) \[
f(ax, by) \Leftrightarrow \cfrac {1}{|ab|} F(\cfrac {u}{a}, \cfrac {v}{b})
\] -
平移特性 (shift) \[
f(x-a, y-b) \Leftrightarrow e^{j 2\pi (au + bv)} F(u, v) \\
f(x, y) e^{j 2 \pi (u_0 x + v_0 y)} \Leftrightarrow F(u-u_0, v- v_0)
\] -
旋转 (rotation)
二维傅里叶变换和一维相比,最大的不同就是二维坐标 \((x, y)\) 可以旋转
利用极坐标系,可以获得一个优良性质: \[
x = r \cos \theta \qquad y=r \sin \theta \\
u = w \cos \varphi \qquad v = w \sin \varphi
\] 则: \(f(r, \theta + \theta_0) \Leftrightarrow F(w, \varphi + \theta_0)\)
原始二维空间信号绕中心旋转 \(\theta_0\) 角度,则对应傅里叶变换后的空间频率信号绕相同方向旋转 \(\theta_0\) 角度
证明时只需要注意 积分 \(dxdy = r dr d\theta\) -
微分 (differentiation) \[
\cfrac {\partial ^n f(x, y)} {\partial x^n} \Leftrightarrow (ju)^n F(u, v) \\
(jx)^n f(x,y) \Leftrightarrow \cfrac{\partial ^n F(u, v)} {\partial u^n}
\] 由此到处著名的 拉普拉斯算子 (laplacian) \[
\nabla^2 f(x, y) = \cfrac {\partial ^2 f(x, y)} {\partial x^2} + \cfrac {\partial ^2 f(x, y)} {\partial y^2}
\] 则: \(\nabla^2 f(x, y) \Leftrightarrow -(u^2+v^2) F(u, v)\)
2.1 常见的 FT 变换对
-
rectangle centered at origin with sides of length \(X\) and \(Y\)
主要利用可分性 \[
\begin {aligned}
F(u, v) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-j2\pi (ux + vy)} dx dy \\
&= \int_{-X/2}^{X/2} e^{-j2 \pi ux}dx \int_{-Y/2}^{Y/2} e^{-j2 \pi vy}dy \\
&= \bigg [\cfrac {e^{-j 2\pi ux}} {-j 2 \pi u} \bigg ] ^{X/2}_{-X/2} \bigg [\cfrac {e^{-j 2\pi vy}} {-j 2 \pi v} \bigg ] ^{Y/2}_{-Y/2} \\
&= \cfrac {1} {-j 2 \pi u} [e^{-j \pi uX} – e^{j \pi uX} ] \cfrac {1} {-j 2 \pi v} [e^{-j \pi vY} – e^{j \pi uY} ] \\
&= XY \bigg[ \cfrac {sin(\pi Xu)} {\pi Xu} \bigg] \bigg[ \cfrac {sin(\pi Yv)} {\pi Yv} \bigg] \\
&= XY \hbox{sinc}(\pi Xu) \hbox{sinc}(\pi Yv)
\end {aligned}
\] -
Gaussian centered on origin \[
f(x, y) = \cfrac {1} {2 \pi \sigma^2} e^{- (x^2+y^2)/ 2 \sigma^2}
\] 可以计算得到傅里叶变换对: \[
F(u, v) = F(\rho) = e^{-2 \pi^2 \rho^2 \sigma^2}
\] 这里 \(\rho^2 = u^2 + v^2\)
高斯函数的一维和二维傅里叶变换仍然是高斯函数
但我们需要注意尺度的反比关系 (inverse scale relation)。 -
Circular disk unit height and radius centered on origin \[
f(x, y) =
\begin{cases}
1, & |r| \lt a, \\
0, & |r| \ge a.
\end{cases}
\]
\[
F(u ,v) = F(\rho) = a J_1(\pi a \rho)/ \rho
\]
这里 \(J_1(x)\) 是 Bessel function -
\(\delta(x, y)\) 函数 \[
f(x, y) = \delta(x, y) = \delta(x) \delta(y) \\
F(u, v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \delta(x, y) e^{-j2\pi (ux + vy)} dx dy =1
\] 如果: \(f(x, y) = \frac {1}{2} (\delta(x, y-a) + \delta(x , y +a))\)
则: \[
\begin {aligned}
F(u, v) &=\cfrac {1}{2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (\delta(x, y-a) + \delta(x , y +a))e^{-j2\pi (ux + vy)} dx dy \\
&= \cfrac {1}{2} (e^{-j 2\pi av} + e^{j 2\pi av}) = \cos 2\pi av
\end {aligned}
\]
2.3 二维傅里叶变换用于滤波的例子
基本流程如下图所示:
和传统的数字信号处理非常类似。
上图展示了
- 低通滤波,图片模糊 (blur),空间频域只剩下低频成分;
- 高通滤波,图片边缘 (edge),空间频域将低频成分消除。
下面这幅图非常重要,展示了图像处理中常见的 5 种滤波方法: mean filter, gaussian filter, laplacian filter, soble filter, scharr filter 的频域分布:
下面是经验结论:
Image with periodic structure, FT has peaks at spatial frequencies of repeated texture.
应用上面的结论,我们可以用于消除背景的 periodic structure
- Forensic application
-
Lunar image
3. 2D Discrete Fourier Transform
对于大小为 \(M \times N\) 的2D DFT,如下: \[
\begin {aligned}
F[k, l] &= \cfrac {1}{MN} \sum_{m =0}^{M-1} \sum_{n=0}^{N-1} f[m, n] e^{-j 2 \pi ( {k \over M} m + {l \over N}n ) } \\
f[m, n] &= \sum_{k =0}^{M-1} \sum_{l=0}^{N-1} F[k, l] e^{-j 2 \pi ( {m \over M} k + {n \over N}l ) }
\end {aligned}
\] (Note: spatial coordinates = \(k,l\), frequency coordinates: \(u,v\))
3.1 性质
FT 的性质,DFT 基本都有,这里介绍几个重要的性质:
- 周期性,一个 \(M \times N\) 矩阵的 2D DFT 的周期性为 \([M, N]\) \[
\begin {aligned}
F[k+M, l+N] &= \cfrac {1}{MN} \sum_{m =0}^{M-1} \sum_{n=0}^{N-1} f[m, n] e^{-j 2 \pi ( {k+M \over M} m + {l+N \over N}n ) } \\
&= \cfrac {1}{MN} \sum_{m =0}^{M-1} \sum_{n=0}^{N-1} f[m, n] e^{-j 2 \pi ( {k \over M} m + {l \over N}n ) } e^{-j 2 \pi ( {M \over M} m + {N \over N}n ) } \\
&= F[k, l]
\end {aligned}
\] -
可分性 (separability) \[
\begin {aligned}
F[k, l] &= \cfrac {1}{MN} \sum_{m =0}^{M-1} \sum_{n=0}^{N-1} f[m, n] e^{-j 2 \pi ( {km \over M} + {ln \over N} ) } \\
&= \cfrac {1}{N} \sum_{n=0}^{N-1} e^{-j 2\pi {ln \over N}} \Big [\sum_{m =0}^{M-1} \cfrac {1}{M} f[m, n] e^{-j 2 \pi {km \over M}} \Big ] \\
\end {aligned}
\] 定义: \[
w_M[m, k] \triangleq \cfrac {1}{M}e^{-j 2 \pi {mk \over M}}; \quad w_N[n, l] \triangleq \cfrac {1}{N}e^{-j 2 \pi {nl \over N}}
\] 进一步定义: \[
F'[k, n] \triangleq \sum_{m =0}^{M-1} f[m, n] w_M[m, k] \quad (k=0, 1, \ldots, M-1)
\] 则上面的 2D 转换为: \[
F[k, l] = \sum_{n =0}^{N-1} F'[k, n] w_N[n,l], \quad (l=0, 1, \ldots, N-1)
\] 解释上面的步骤:- column transform
考虑表达式 \(F'[k, n]\) ,我们将列索引 \(n\) 是为参数,那么这个表达式其实是在对列向量进行 1D DFT 计算 \[
\begin {bmatrix}
F'[0, n] \\
\vdots \\
F'[M-1, n]
\end {bmatrix}_{M \times 1} =
\begin {bmatrix}
\cdots & \cdots & \cdots \\
\cdots & w_M[m, k] & \cdots \\
\cdots & \cdots & \cdots
\end {bmatrix}_{M \times M}
\begin {bmatrix}
f[0, n] \\
\vdots \\
f[M-1, n]
\end {bmatrix}_{M \times 1}
\] 更简洁些: \(\mathbf {F}'\) 的第 \(n\) 列向量是 \(\mathbf {f}\) 的第 \(n\) 列向量进行 1D DFT 计算的结果,即: \[
\mathbf {F}'_n = \mathbf {W}_M \mathbf {f}_n, \quad (n =0, \cdots, N-1)
\] 将所有 \(N\) 列合并起来: \[
\begin {bmatrix} \mathbf {F}'_0 & \cdots & \mathbf {F}'_{N-1} \end {bmatrix} = \mathbf {W}_M \begin {bmatrix} \mathbf {f}_0 & \cdots & \mathbf {f}_{N-1} \end {bmatrix}
\] 即:
\[
\mathbf {F}' = \mathbf {W}_M \mathbf {f}
\] -
row transform
和上面类似,我们按行划分计算 \[
\begin {bmatrix}
F[k, 0] & \cdots & F[k, N-1]
\end {bmatrix} =
\begin {bmatrix}
F'[k, 0] & \cdots & F'[k, N-1]
\end {bmatrix}
\begin {bmatrix}
\cdots & \cdots & \cdots \\
\cdots & w_N[n, l] & \cdots \\
\cdots & \cdots & \cdots
\end {bmatrix}
\] \(\mathbf {F}\) 的第 \(k\) 行向量是 \(\mathbf {F}'\) 的第 \(k\) 行向量进行 \(N\) 点 DFT 运算得到的
如果我们直接计算 \(M \times N\) 像素点的图片的 DFT,时间复杂度为 \(\mathcal {O}(M^2N^2)\), 现在只需要 \(N\) 次 \(M\) 点DFT,以及 \(M\) 次 \(N\) 点 DFT 就可以, 复杂度为 \(\mathcal {O}(NM^2 + MN^2)\)
- column transform
3.3 频谱中心化 (spectrum centralization)
在 数字信号处理中聊过 DFT的频谱搬移问题 ,这里同样存在这个问题
DC 分量位于 \(F[0, 0]\) (左上角),最高频分量位于 \(F[M/2, N/2]\) (位于中心)。也就是说: 高频分量在中心,低频分量在四个角落。
这并不符合我们的认知习惯。因为我们喜欢把中心原点放在图片的中心,而不是左上角。我们通常会将频谱中心化,将频谱垂直移动 \(M/2\) ,水平移动 \(N/2\) 。
和 1D DFT 类似,利用 2D DFT 的平移特性 (上面有提到过) \[
x[m, n] e^{j 2\pi (mM/2M + nN/2N) } \Leftrightarrow X[k-M/2, l – N/2]
\] 即: \[
x[m, n] (-1)^{m+n} \Leftrightarrow X[k-M/2, l – N/2]
\] 如下面所示: \[
\begin {bmatrix}
x[0,0] & -x[0,1] & x[0,2] & \cdots \\
-x[1,0] &x[1,1] & -x[1,2] & \cdots \\
x[2,0] & -x[2,1] & x[2,2] & \cdots \\
\cdots & \cdots & \cdots & \cdots
\end {bmatrix}
\]
3.4 频谱可视化
前面放了很多张图的例子,是不是觉得频谱图片很好看 ?
但是如果你现在调用 OpenCV 中的 fft 函数,并且将幅频图中心化,可能得到的效果并没有我们显示的那么好。
一个重要原因就是我们并没有考虑 幅频值的范围, 或者说低频对应的值可能是高频对应的值的 1000 倍,如果在非常大的范围内显示幅频图,那么由于低频强度太高,会导致大部分区域是黑色的 (相比之下,值接近 0)
我们通常会对数化处理,即: \[
F'[k ,l] = \log (1 + F[k, l])
\] 此时效果就很棒了。
如下图所示,效果对比非常鲜明: