本文详细介绍了数字图像处理中的频率域滤波，包括傅里叶变换原理、一维与二维傅里叶变换，以及在Python中实现傅里叶级数。通过低通与高通滤波的实验，展示了如何在频率域中分析和处理图像，以改变图像的频率特征，进而平滑图像或突出边缘细节。
数字图像处理——第4章 频率域滤波
文章目录
频率域
上一章学的灰度变换和空间滤波，主要目的是减少噪声和平滑图像，同时也是在空间域进行的图像增强操作。而这章主要是频率域滤波，滤波的概念在上一章节有涉及到，所以先了解下频率域。频率域是指从函数的频率角度出发分析函数，和频率域相对的是时间域。简单说就是如果从时间域分析信号时，时间是横坐标，振幅是纵坐标。而在频率域分析的时候则是频率是横坐标，振幅是纵坐标。那么下一个问题就是为什么要在频率域中进行图像处理？ 1). 可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务，在频率域中变得非常普通；2). 滤波在频率域更为直观，它可以解释空间域滤波的某些性质； 3).可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为空间域滤波器的指导。
频率域滤波是对图像进行傅里叶变换，将图像由图像空间转换到频域空间，然后在频率域中对图像的频谱作分析处理，以改变图像的频率特征。
1.傅里叶级数原理
傅里叶是18世纪法国的一位伟大的数学家。他指出任何周期函数都可以表示为不同频率的正弦和或者余弦和的形式，每个正弦或者余弦乘以不同的系数就是傅里叶级数。无论函数有多复杂，只要它是周期性的，并且满足一定的数学条件，就一定可以用这样的正弦和或者余弦和的形式来表示。翻开考完研没扔的高数，回看到傅里叶级数，公式如下：
f
(
x
)
=
a
0
2
+
∑
n
=
1
∞
(
a
n
cos
n
π
l
x
+
b
n
sin
n
π
l
x
)
f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \frac{n \pi}{l} x+b_{n} \sin \frac{n \pi}{l} x\right)
f(x)=2a0+n=1∑∞(ancoslnπx+bnsinlnπx)
其中：
a n = 1 l ∫ − l l f ( x ) cos n π l x d x , n = 0 , 1 , 2 , ⋯ , 其 中 a 0 为 n = 0 时 a_{n}=\frac{1}{l} \int_{-l}^{l} f(x) \cos \frac{n \pi}{l} x \mathrm{~d} x, n=0,1,2, \cdots,其中a_{0}为n=0时 an=l1∫−llf(x)coslnπx dx,n=0,1,2,⋯,其中a0为n=0时
b n = 1 l ∫ − l l f ( x ) sin n π l x d x , n = 0 , 1 , 2 , ⋯ , b_{n}=\frac{1}{l} \int_{-l}^{l} f(x) \sin \frac{n \pi}{l} x \mathrm{~d} x, n=0,1,2, \cdots, bn=l1∫−llf(x)sinlnπx dx,n=0,1,2,⋯,
除此之外，若展开式中只有正弦函数，则称其为正弦函数；相反，若展开式中只有余弦函数，则称其为余弦级数。
1.1.一维傅里叶变换
傅里叶变换：
F
(
μ
)
=
∫
−
∞
+
∞
f
(
t
)
e
−
j
2
π
μ
t
d
t
F(\mu)=\int_{-\infty}^{+\infty} f(t) e^{-j 2 \pi \mu t} d t
F(μ)=∫−∞+∞f(t)e−j2πμtdt
傅里叶反变换：
f
(
t
)
=
∫
−
∞
+
∞
F
(
μ
)
e
j
2
π
μ
t
d
μ
f(t)=\int_{-\infty}^{+\infty} F(\mu) e^{j 2 \pi \mu t} d \mu
f(t)=∫−∞+∞F(μ)ej2πμtdμ
其
中
，
2
π
μ
代
表
频
率
，
t
代
表
时
间
。
其中，2 \pi \mu代表频率，t代表时间。
其中，2πμ代表频率，t代表时间。傅里叶变换认为一个周期函数(信号)包含多个频率分量，任意函数（信号）f(t)可通过多个周期函数（基函数）相加而合成。从物理角度理解傅里叶变换是以一组特殊的函数（三角函数）为正交基，对原函数进行线性变换，物理意义便是原函数在各组基函数的投影。
1.2.二维傅里叶变换
同样连续的傅里叶变换为：
F
(
μ
,
v
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
f
(
t
,
z
)
e
−
j
2
π
(
μ
t
+
v
z
)
d
t
d
z
F(\mu, v)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(t, z) e^{-j 2 \pi(\mu t+v z)} d t d z
F(μ,v)=∫−∞+∞∫−∞+∞f(t,z)e−j2π(μt+vz)dtdz
它的反变换：
f
(
t
,
z
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
F
(
μ
,
v
)
e
j
2
π
(
μ
t
+
v
z
)
d
μ
d
v
f(t, z)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} F(\mu, v) e^{j 2 \pi(\mu t+v z)} d \mu d v
f(t,z)=∫−∞+∞∫−∞+∞F(μ,v)ej2π(μt+vz)dμdv
傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数傅里叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，通常用一个二维矩阵表示空间上各点，记为z=f(x,y)。又因空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就必须由梯度来表示，这样我们才能通过观察图像得知物体在三维空间中的对应关系。
傅里叶频谱图上我们看到的明暗不一的亮点，其意义是指图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅里叶变换后的频谱图，也叫功率图，我们就可以直观地看出图像的能量分布：如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小）；反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的、边界分明且边界两边像素差异较大的。
2.python×傅里叶级数
理论过于抽象，还是实实在在的用代码折腾一下比较好理解。首先我们需要一个原函数：
f
(
x
)
=
{
exp
(
x
)
,
−
π
≤
x
<
0
1
,
0
≤
x
<
π
f(x)=\left\{\begin{array}{lr} \exp (x), & -\pi \leq x<0 \\ 1, & 0 \leq x<\pi \end{array}\right.
f(x)={exp(x),1,−π≤x<00≤x<π
函数图像如下：
可以用numpy和plt画出，代码如下：
import os
import cv2
from pylab import *
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(-3,pi)# 定义x轴
cond = [True if (i>-pi and i<0) else False for i in x] #使用列表解析的方法
y=1*(x>0)+np.exp(x)*cond # y的函数
plt.plot(x,y)
plt.show()
接着算出该函数的a0，an，bn，最后自定义函数画出傅里叶级数变换，图像如下：
由上图看出，这图像就有正弦余弦内味了，周期函数都可以表示为不同频率的正弦和或者余弦和的形式。代码如下，建议做这种小实验使用jupyter，方便快捷
import os
import cv2
from pylab import *
import matplotlib.pyplot as plt
import numpy as np
x = mgrid[-10:10:0.2]
n = arange(1, 1000)
def fourier_transform():
a0 = (1-exp(-pi))/pi+1
s = a0 / 2
for i in range(1, 100, 1):
s0 = ((1-(-1)**i*exp(-pi))/(pi*(1+i**2))*cos(i*x)+1/pi*( (-i*(1-(-1)**i*exp(-pi)))/(1+i**2) + (1-(-1)**i)/i ) * sin(i*x))
s = s + s0
plot(x,s,'orange',linewidth=0.6)
title('fourier_transform')
show()
fourier_transform()
2.1.傅里叶变换后的频谱图
使用傅里叶变换获得的频谱图。图像的主要成分是低频信息，它形成了图像的基本灰度等级，对图像结构的决定作用较小；高频信息形成了图像的边缘和细节，是在中频信息上对图像内容的进一步强化。
如上图所示，中间的亮点代表着图像的低频，其余周围则表示图像的高频信息。若不将低频信息移至中间，则低频信息过于分散，频谱未中心化如下图所示：
代码如下所示：
import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读取图片
img = cv2.imread(' ', 0)
# 进行float32形式转换
float32_img = np.float32(img)
# 使用cv2.dft进行傅里叶变化
dft_img = cv2.dft(float32_img, flags=cv2.DFT_COMPLEX_OUTPUT)
# 使用np.fft.shiftfft()将变化后的图像的低频转移到中心位置
dft_img_ce = np.fft.fftshift(dft_img)
# 使用cv2.magnitude将实部和虚部转换为实部，乘以20是为了使得结果更大
# img_dft = 20 * np.log(cv2.magnitude(dft_img_ce[:, :, 0], dft_img_ce[:, :, 1]))
# 频谱未中心化
img_dft = 20 * np.log(cv2.magnitude(dft_img[:, :, 0], dft_img[:, :, 1]))
# 进行画图操作
plt.subplot(121) # 我是第一行第二列的第一个
plt.imshow(img, cmap='gray')
plt.subplot(122)
plt.imshow(img_dft, cmap='gray')
plt.show()
3.频率域滤波
3.1.低频与高频
低频对应图像内变换缓慢的灰度分两。例如在一幅大草原的图像中，低频对应着广袤的颜色趋于一致的草原。总结就是低频是图像中平坦的，灰度变化不大的点，图像中的大部分区域。
高频对应图像内变化越来越快的灰度分量，是由灰度的尖锐过度造成的，例如，还是大草原那幅图像，若图像中有一只狮子，则狮子的边缘信息就是高频。总结就是高频是图像中灰度变化剧烈的点，一般是图像轮廓或者是噪声。
根据图像的高频与低频的特征，我们可以设计相应的高通与低通滤波器。通过低频的滤波器叫做低通滤波器，通过高频的滤波器叫做高通滤波器。高通滤波可以检测图像中尖锐、变化明显的地方；低通滤波可以让图像变得光滑，滤除图像中的噪声。
3.2.频率域滤波步骤
为了方便下述低通高通滤波器的实验，必须先明确频率域滤波的基本步骤。
- 首先将图片读取成灰度图且转化为np.float32()
- 使用cv2.dft函数进行傅里叶变化
- 使用np.fft.fftshift函数将低频转移到图像中心
- 定义掩模：生成的掩模中间为1周围为0
- 将掩模与傅里叶变化后图像相乘，保留中间部分
- 使用np.fft.ifftshift将低频移动到原来的位置
- 使用cv2.idft进行傅里叶的反变化
- 使用cv2.magnitude转化为空间域内
- 使用plt进行绘图展示
由上述叙述可知，滤波的关键取决于滤波函数，也称为滤波器，因为它在滤波中抑制或除去了频谱中的某些分量，而保留其他的一些频率不受影响，从而达到滤波的目的。而下述的 低通滤波器和高通滤波器的不同点，便在于定义不同的滤波，即步骤中的掩模，不同的滤波将产生不同的效果。
3.3.低通滤波器
理想的低通滤波器模板为：
H
(
u
,
v
)
=
{
1
,
D
(
u
,
v
)
≤
D
0
0
,
D
(
u
,
v
)
>
D
0
H(u, v)=\left\{\begin{array}{ll} 1, & D(u, v) \leq D_{0} \\ 0, & D(u, v)>D_{0} \end{array}\right.
H(u,v)={1,0,D(u,v)≤D0D(u,v)>D0
其中，D0表示通带半径，D(u，v)是到频谱中心的距离（欧式距离），计算公式如下：
D
(
u
,
v
)
=
(
u
−
M
/
2
)
2
+
(
v
−
N
/
2
)
2
\mathrm{D}(u, v)=\sqrt{(u-M / 2)^{2}+(v-N / 2)^{2}}
D(u,v)=(u−M/2)2+(v−N/2)2
M和N表示频谱图像的大小，（M/2，N/2）即为频谱中心。低通滤波器代码如下所示：
# 低通滤波器
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 第一步读入图片
img = cv2.imread(' ', 0)# filepath
# 使用cv2.dft进行傅里叶变化
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
# 使用np.fft.fftshift将低频转移到图像中心
dft_center = np.fft.fftshift(dft)
# 定义掩模：生成的掩模中间为1周围为0
crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
mask = np.zeros((img.shape[0], img.shape[1], 2), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1
# 将掩模与傅里叶变化后图像相乘，保留中间部分
mask_img = dft_center * mask
# 使用np.fft.ifftshift(将低频移动到原来的位置
img_idf = np.fft.ifftshift(mask_img)
# 使用cv2.idft进行傅里叶的反变化
img_idf = cv2.idft(img_idf)
# 使用cv2.magnitude转化为空间域内
img_idf = cv2.magnitude(img_idf[:, :, 0], img_idf[:, :, 1])
# 进行绘图操作
plt.subplot(121)
plt.imshow(img, cmap='gray')
plt.subplot(122)
plt.imshow(img_idf, cmap='gray')
plt.show()
由实验可以看出，图片经过低通滤波器，使低频通过而使高频衰减 ，降低了图像变换的幅度。平滑且模糊掉了尖锐和边缘部分。保留了图像的整体概貌（低频部分），去掉了图像的细节（高频部分）。
3.4.高通滤波器
和低通滤波器相反，理想的高通滤波器为：
H
(
u
,
v
)
=
{
0
,
D
(
u
,
v
)
≤
D
0
)
1
,
D
(
u
,
v
)
>
D
0
)
H(u, v)=\left\{\begin{array}{ll} 0, & \left.\mathrm{D}(\mathrm{u}, \mathrm{v}) \leq D_{0}\right) \\ 1, & \left.\mathrm{D}(\mathrm{u}, \mathrm{v})>D_{0}\right) \end{array}\right.
H(u,v)={0,1,D(u,v)≤D0)D(u,v)>D0)
其中，D0表示通带半径，D(u，v)是到频谱中心的距离（欧式距离），计算公式如下：
D
(
u
,
v
)
=
(
u
−
M
/
2
)
2
+
(
v
−
N
/
2
)
2
\mathrm{D}(u, v)=\sqrt{(u-M / 2)^{2}+(v-N / 2)^{2}}
D(u,v)=(u−M/2)2+(v−N/2)2
常见的巴特沃斯高通滤波器(BHPF)公式如下：
H
(
u
,
v
)
=
1
1
+
[
D
0
/
D
(
u
,
v
)
]
2
n
H(u, v)=\frac{1}{1+\left[D_{0} / D(u, v)\right]^{2 n}}
H(u,v)=1+[D0/D(u,v)]2n1
还有高斯高通滤波器（GHPF)：
H
(
u
,
v
)
=
1
−
e
−
D
2
(
u
,
v
)
/
2
σ
2
,
σ
=
D
0
H(u, v)=1-e^{-D^{2}(u, v) / 2 \sigma^{2}}, \sigma=D_{0}
H(u,v)=1−e−D2(u,v)/2σ2,σ=D0
# 实现高通滤波器
import cv2
import numpy as np
from matplotlib import pyplot as plt
#读取图像
img = cv2.imread(' ', 0) # filepath
#傅里叶变换
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)
fshift = np.fft.fftshift(dft)
#设置高通滤波器
rows, cols = img.shape
crow,ccol = int(rows/2), int(cols/2) #中心位置
mask = np.ones((rows, cols, 2), np.uint8) # 低通滤波器这是np.zeros
mask[crow-30:crow+30, ccol-30:ccol+30] = 0 # 低通滤波器这就是1
#掩膜图像和频谱图像乘积
f = fshift * mask
#傅里叶逆变换
ishift = np.fft.ifftshift(f)
iimg = cv2.idft(ishift)
res = cv2.magnitude(iimg[:,:,0], iimg[:,:,1])
#显示原始图像和高通滤波处理图像
plt.subplot(121), plt.imshow(img, 'gray')
plt.subplot(122), plt.imshow(res, 'gray')
plt.show()
实现效果如下图所示：
由上述实验可以看出：高通滤波器使高频通过而使低频衰减，被高通滤波的图像比原始图像少灰度级的平滑过渡而突出边缘等细节部分，保留了图像的细节，模糊了概貌。也就是边缘信息更加突出，例如轮廓。
3.5.低通与高通滤波器实验总结
由上述实验可看出，低通滤波器的滤波定义为：
crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
mask = np.zeros((img.shape[0], img.shape[1], 2), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1
高通滤波器的滤波为：
crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
mask = np.ones((img.shape[0], img.shape[1], 2), np.uint8) # 低通滤波器这是np.zeros
mask[crow-30:crow+30, ccol-30:ccol+30] = 0 # 低通滤波器这就是1
相同点在于定义的掩模大小和图片大小相同，不同点在于低通滤波器中间为1，那么经过傅里叶变换的图像和掩模相乘之后，便只保留中间部分；相反，高通滤波器中间为0，经过傅里叶变换的图像和掩模相乘之后，便只保留边缘部分。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
