数字图像处理——频域滤波

Posted on 2019-04-29 Edited on 2025-09-18 In 数字图像处理

空间域的算法很多，也容易理解，但是有时候却很难设计。而有时候一些问题，转化到频域上，就变得迎刃而解了。

从空间域到频域，不得不提的是傅里叶变换(一般信号处理上，从时域到频域)。这里首先说明下什么是变换，也就是将原有的函数，转换成另外的一些更容易分析的函数之和。变换给了我们不同看事物的角度。比较有名的变换有傅里叶变换，拉普拉斯变换，小波变换等等。

傅里叶变换

我们这里主要用到的是傅里叶变换。傅里叶变换的原始定义如下：
$$
f(t) \rightarrow F(\omega)=\int f(t) e^{-j \omega t} d t
$$
这里的$t, \omega$可以是标量，也可以是向量。

如果值是离散的，那么有离散傅里叶变化（DFT，Discrete Fourier Transformation）。离散傅里叶变换和连续傅里叶变化形式上差距还是挺大的：
$$
F(u)=\sum_ {x=0}^{N-1} f(x) \cdot e^{-i \frac{2 \pi}{M} x u}
$$
从矩阵视角：
$$
\begin{array}{c}{\vec{v}=A \vec{u} } \ {\vec{u}=[u(0), u[1], \cdots, u(N-1)]^{T}, \vec{v}=[v(0), v[1], \cdots, v(N-1)]^{T} } \ {A=[a(I, m)]_ {N \times N}, \text { where } a(l, m)=\frac{1}{\sqrt{N} } e^{-j \frac{2 \pi}{N}(l-1)(m-1)} }\end{array}
$$
为什么长这个样子，离散傅里叶变换背后的数学推倒是比较复杂的，幸运的是我们也不需要了解太多，因此在这里就不多说了。这里$A$矩阵是酋矩阵（unitary matrix）：

$A^{-1}=A^{* T}=A^{H}.$

有了这个，IDFT（离散傅里叶逆变换）就很容易计算了。因为：
$$
\begin{aligned} u &=A^{-1} v \ u &=A^{H} v \ \Rightarrow u[n] &=\sum_ {k=0}^{N-1} v[k] a^{}(k, n) \end{aligned}
$$
这里，$v(k)$是加权系数，而$a^{}(k, n)$为基函数。

上面提到的是一维的DFT，对于二维的DFT会更麻烦一点。

原始的图像我们用U表示，如下：
$$
U={u[m, n], 0 \leq m, n \leq N-1}
$$
转换系数：
$$
v[k, l]=\sum_ {m=0}^{N-1} \sum_ {n=0}^{N-1} u[m, n] a_ {k l}[m, n]
$$
因为A是酋矩阵：
$$
u[m, n]=\sum_ {k=0}^{N-1} \sum_ {l=0}^{N-1} v[k, l] a_ {k /}^{}[m, n]
$$
基函数$a_ {kl}[m,n]$以及$a_ {kl}^[m,n]$满足：

特殊域正交： $$\sum_ {m=0}^{N-1} \sum_ {n=0}^{N-1} a_ {k l}[m, n] a_ {k^{\prime} / l}^{*}[m, n]=\delta\left[k-k^{\prime}, I-l^{\prime}\right]$$
空间域的完备性： $$\sum_ {k=0}^{N-1} \sum_ {l=0}^{N-1} a_ {k l}[m, n] a_ {k l}^{*}\left[m^{\prime}, n^{\prime}\right]=\delta\left[m-m^{\prime}, n-n^{\prime}\right]$$

可以看到这里任何一个系数估计都是全局计算，需要$N^2$个乘法以及$N^2-1$个加法。而所有系数计算需要$N^4$的复杂度。这个复杂度是难以接受的。因此我们需要快速傅里叶变换（FFT），对于快速傅里叶变换的具体算法这里就不多介绍了，但是它已经被嵌入到各种各样的程序包中了，可以随时调用。

常见的傅里叶变换

为了介绍离散傅里叶变换的性质，我们先认识下面一些函数：

对于傅里叶变换：$f(x,y) \rightarrow F(x,y)$，那么一些常见的函数的傅里叶变换有（$\delta(x, y )$为冲激函数）：

$f(x,y)$	$F(x,y)$
$\delta(x, y)$	$1$
$\delta\left(x \pm x_ {0}, y \pm y_ {0}\right)$	$e^{ \pm j 2 \pi\left(x_ {0} \xi_ {1}+y_ {0} \xi_ {2}\right)}$
$e^{-\pi\left(x^{2}+y^{2}\right)}$	$e^{-\pi\left(\xi_ {1}^{2}+\xi_ {2}^{2}\right)}$
$\operatorname{rect}(x, y)=\operatorname{rect}(x) \operatorname{rect}(y)$	$\operatorname{sinc}\left(\xi_ {1}, \xi_ {2}\right)=\operatorname{sinc}\left(\xi_ {1}\right) \operatorname{sinc}\left(\xi_ {2}\right)$
$\operatorname{tri}(x, y)$	$\sin c^{2}\left(\xi_ {1}, \xi_ {2}\right)$
$\operatorname{comb}(x, y)$	$\operatorname{comb}\left(\xi_ {1}, \xi_ {2}\right)$

Translation

${f(x, y) \rightarrow F(u, v)}$
${f(x, y) e^{\left(j u_ {0} x+j v_ {0} y\right)} \rightarrow F\left(u-u_ {0}, v-v_ {0}\right)}$
${f\left(x-x_ {0}, y-y_ {0}\right) \rightarrow F(u, v) e^{\left(-j u x_ {0}-j y_ {0}\right)} }$

Multiplication and Convolution

${f(x, y) \star g(x, y) \Rightarrow F(u, v) G(u, v)} $
${f(x, y) g(x, y) \Rightarrow F(u, v) \star G(u, v)} $

Scaling

f(a t) \Longleftrightarrow \frac{1}{|a|} F\left(\frac{\omega}{a}\right), \quad a \in \mathbb{R}, a \neq 0

对于傅里叶，必须要提的是高斯函数。它有很多优秀的性质：

在空间和频域上都有一定的跨度。
高斯的傅里叶变换依然是高斯函数
高斯函数乘以高斯函数还是一个高斯函数，只是标准差不同
高斯函数的标准差是一个尺度参数
我们经常需要在不同尺度上观察图片，而高斯金字塔是一个很好的工具

高斯金字塔是著名的SIFT特征必不可少的预处理步骤。不同尺寸的高斯金字塔，模拟了场景离我们眼睛越来越远时候看到的效果，想对高斯金字塔有详细了解的可以参考论文：Gaussian Pyramid.

频域滤波

在频域上滤波，能够多个看图片的视角，也很容易设计滤波器的形状。它的缺点是计算量较大，由于数字化，使得算法也较为复杂，比如上面的傅里叶变换，我就没完全理解。但是有一些问题，在频域下就非常直观。比如下图：

可以看到原图的去噪（很有规律）在空间域上是几乎无法去除的，但是在频域上就可以看出，是多余的一些异常点造成了这样的噪声。频域滤波的一般过程如下：

在空间域上的卷积，对应着频域的乘积，因此频域中使用的是乘积操作。频域上一般有高通滤波，低通滤波以及带通滤波。

下面是一张图片的原始图和频域图：

低通滤波

低通滤波可以用来平滑或者模糊图像。关于低通滤波，这里介绍三种滤波器，理想的低通滤波，巴特沃斯滤波，高斯低通滤波。它们都是各向同性的，它们之间的关系是：

理想的低通滤波数学形式如下：
$$
H(u, v)=\left{\begin{array}{ll}{1} & {\text { if } D(u, v) \leq D_ {0} } \ {0} & {\text { if } D(u, v)>D_ {0} }\end{array}\right.
$$
$D(u, v )$代表的是到原点的欧几里得距离，低通滤波如下图：

这里$D_0$是截断频率，如果确定截断频率是非常重要的一件事，它和我们想滤掉多少频分以及图像的分布有关。下面是一个低通滤波的例子：

低通滤波很明显可以模糊图像，但是它会带来振铃效应（ringing effect）:

**高斯低通滤波（Gaussian Low-Pass Filter）**的做法是用$D_0$代替高斯滤波中的$\delta$：
$$
H(u, v)=e^{-\frac{D^{2}(u, v)}{2 D_ {0}^{2} } }
$$

高斯低通滤波的好处在于它在空间域内也容易得到解析形式，而不需要傅里叶变换，而且与理想的低通滤波器相比，在空间域它没有振铃效应。

但是高斯低通滤波的缺点是有时候它不够sharp。因为我们知道，理想低通滤波是直接截断，而高斯滤波则是有比较缓慢的过渡过程。

巴特沃斯低通滤波(Butterworth Low-Pass Filter), 是介于低通滤波与高斯滤波中间的滤波器。它的数学形式如下：
$$
H(u, v)=\frac{1}{1+\left[\frac{D(u, v)}{D_ {0} }\right]^{2 n} }
$$

它的振铃效应和次数n有关。

效果图如下：

高通滤波

有时候我们需要锐化特征的滤波器，高通滤波可以做到这一点。类似于低通滤波，它也有三种类型，理想高通滤波（Ideal High Pass Filter），高斯高通滤波（Gaussian High Pass Filter）以及巴特沃斯高通滤波（Butterworth High Pass）。

理想高通滤波数学形式：
$$
H(u, v)=\left{\begin{array}{ll}{0} & {\text { if } D(u, v) \leq D_ {0} } \ {1} & {\text { if } D(u, v)>D_ {0} }\end{array}\right.
$$
高斯高通滤波数学形式：
$$
H(u, v)=1-e^{-\frac{D^{2}(u, v)}{2 D_ {0}^{2} } }
$$
巴特沃斯高通滤波数学形式：
$$
H(u, v)=\frac{1}{1+\left[\frac{D_ {0} }{D(u, v)}\right]^{2 n} }
$$
下图分别是IHPF，GHPF，BHPF的图像：

滤波效果：

空间域和频域

在空间域的滤波可以转换成设计在频域的滤波，比如拉普拉斯算子。在空间域中，它的推导如下：
$$
\begin{aligned} \nabla^{2} f &=\frac{\partial^{2} f}{\partial x^{2} }+\frac{\partial^{2} f}{\partial y^{2} } \ \frac{\partial^{2} f}{\partial x^{2} } &=f(x+1, y)+f(x-1, y)-2 f(x, y) \ \frac{\partial^{2} f}{\partial y^{2} } &=f(x, y+1)+f(x, y-1)-2 f(x, y) \ \nabla^{2} f &=[f(x+1, y)+f(x-1, y)+f(x, y+1)+f(x, y-1)]-4 f(x, y) \end{aligned}
$$
由此得到下面的滤波器：

而在频域中，我们可以计算：
$$
\begin{aligned} \mathfrak{F}\left[\frac{d^{n} f(x)}{d x^{n} }\right] &=(j u)^{n} F(u) \ \widetilde{\mathfrak{F} }\left[\frac{\partial^{2} f(x, y)}{\partial x^{2} }+\frac{\partial^{2} f(x, y)}{\partial y^{2} }\right] &=(j u)^{2} F(u, v)+(j v)^{2} F(u, v) \ &=-\left(u^{2}+v^{2}\right) F(u, v) \ \mathfrak{F}\left[\nabla^{2} f(x, y)\right] &=-\left(u^{2}+v^{2}\right) F(u, v) \end{aligned}
$$
因此在频域中，滤波器形状为：$$H(u, v)=-\left(u^{2}+v^{2}\right)$$.

在之前，我们介绍到锐化特征时候的做法，就是先通过原图与模糊图相减得到边缘，然后再加到原始图像上。使用拉普拉斯算子也可以这样做，将使用拉普拉斯算子滤波后的图加到原始图像上，可以得到锐化的效果。数学形式如下：
$$
g(x, y)=\left{\begin{array}{ll}{f(x, y)-\nabla^{2} f(x, y)} & {\text { if the center coefficient of the Laplacian mask is negative} } \ {f(x, y)+\nabla^{2} f(x, y)} & {\text { if the center coefficient of the Laplacian mask is positive.} }\end{array}\right.
$$
效果如图：

在频域上我们这样做，也可以实现锐化的效果：

在空间域上，我们做非锐化屏蔽（unsharp masking）是这样的：
$$
f_ {\mathrm{hp} }(x, y)=f(x, y)-f_ {\mathrm{lp} }(x, y)
$$
在频域上，需要下面这样做：
$$
F_ {\mathrm{hp} }(u, v)=F(u, v)-F_ {\mathrm{lp} }(u, v)
$$
因为$F_ {\mathrm{lp} }(u, v)=H_ {\mathrm{lp} }(u, v) F(u, v)$，所以我们得到：
$$
H_ {\mathrm{hp} }(u, v)=1-H_ {\mathrm{lp} }(u, v)
$$
一般来说，空间域易于计算，频域更容易理解和设计算法，但是相对复杂，计算量也较大，DFT是实现频域算法的关键步骤。

其他一些频域滤波

高通量滤波

上一篇文章介绍了，在空间域中，高通量滤波（High-boosting filtering）的做法，首先计算非锐化屏蔽（unsharp masking）的图：
$$
f_ {s}(x, y)=f(x, y)-\overline{f}(x, y)
$$
然后让原图加上计算的图（这里的A可能和之前的k相关）：
$$
\begin{aligned} f_ {\mathrm{hb} }(x, y) &=A f(x, y)-\overline{f}(x, y) \ f_ {\mathrm{hb} }(x, y) &=(A-1) f(x, y)+f(x, y)-\overline{f}(x, y) \ f_ {\mathrm{hb} }(x, y) &=(A-1) f(x, y)+f_ {s}(x, y) \end{aligned}
$$
而在频域上的高通量滤波如下：
$$
H_ {\mathrm{hb} }(u, v)=(A-1)+H_ {\mathrm{hp} }(u, v)
$$

高频强调滤波

除了高通量滤波，这里再介绍一个类似于高通量滤波的滤波器：高频强调滤波（High-Frequency Emphasis Filtering）。它可以增强高频成分。数学形式如下：
$$
H_ {\mathrm{hfe} }(u, v)=a+b H_ {\mathrm{hp} }(u, v)
$$
可以看到的是，当上式中$a = (A-1),b=1$时，它与高通量滤波一致。一般来说这里的$b$会大于1，可以强调高频的成分。

同态滤波

同态滤波（Homomorphic Filtering）想法来源于illumination-reflectance模型：
$$
f(x, y)=i(x, y) r(x, y)
$$
但是在频域上是很难实现的，因为空间域中的乘积意味着频域中的卷积。解决方法是，在转化到频域之前，先利用log函数，然后经过滤波之后回到空间域，再使用log的逆达到想要的效果。流程如下：

同态滤波应用在于减少光照变化并且锐化边缘，它可以减少低频增加高频。因为显示中光照变化缓慢，属于低频，而反射变化比较突然，属于高频信息。同态滤波的函数曲线应该如下图，低频的信息给较低的权重，而高频的信息给较高的权重。

下面是同态滤波效果图：

实现上图的效果的滤波器数学形式如下：
$$
\begin{array}{l}{H(u, v)=\left(\gamma_ {H}-\gamma_ {L}\right)\left[1-e^{-c\left(\frac{D^{2}(u, y)}{D_ {0}^{2} }\right)}\right]+\gamma_ {L} } \ {\text { where } \gamma_ {L}=0.5, \gamma_ {H}=2.0}\end{array}
$$
实现中一些值得注意的问题

FS，FT，DFT

对于不同的信号，空间域中，非周期的连续信号，使用FT，周期连续信号，使用FT，FS（傅里叶级数），周期离散信号，使用DFT。而使用DFT时，我们暗示假设了周期性，因此总是会使用循环卷积。

Zero padding

补零操作细节如下：
$$
f_ {e}(x)=\left{\begin{array}{ll}{f(x)} & {0 \leq x \leq A-1} \ {0} & {A \leq x \leq P}\end{array}\right. g_ {e}(x)=\left{\begin{array}{ll}{g(x)} & {0 \leq x \leq B-1} \ {0} & {B \leq x \leq P}\end{array}\right.
$$
上式中，$P\ge A +B - 1$。

DC component

DFT中直流分量占据了大量的位置，一般来说，我们会把直流分量放在图的中间。如果实现这个？在matlab中使用fftshift，它利用了DFT的性质：
$$
f\left(x-x_ {0}, y-y_ {0}\right) \longleftrightarrow F(u, v) e^{-j 2 \pi\left(u \frac{x_ {0} }{M}+v \frac{y_ {0} }{N}\right)}
$$
当$x_0=\frac{M}{2},y_0=\frac{N}{2}$：
$$
f\left(x-\frac{M}{2}, y-\frac{N}{2}\right) \longleftrightarrow F(u, v)(-1)^{(u+v)}
$$
类似的：
$$
f(x, y) e^{j 2 \pi\left(u \frac{x_ {0} }{M}+v \frac{y_ {0} }{N}\right)} \longleftrightarrow F\left(u-u_ {0}, v-v_ {0}\right)
$$
当$u_0=\frac{M}{2},v_0 = \frac{N}{2}$，
$$
f(x, y)(-1)^{(x+y)} \longleftrightarrow F\left(u-\frac{M}{2}, v-\frac{N}{2}\right)
$$