（数字图像处理MATLAB+Python）第六章图像平滑-第二节：频域平滑滤波_fimage程序

频域平滑滤波：是数字图像处理中常用的一种滤波方法，它的主要思想是在频域中对图像进行滤波处理。其步骤如下

• 对原始图像进行傅里叶变换，得到频域图像
• 对频域图像进行滤波操作，可以选择常用的低通滤波器，例如高斯低通滤波器、均值滤波器等
• 对滤波后的频域图像进行傅里叶逆变换，得到平滑后的图像

频域低通滤波的核心关键为设计合适的低通滤波器$H(u,v)$

$$G(u,v)=H(u,v)F(u,v)$$

频域平滑滤波有以下几类

• 理想低通滤波
• 巴特斯沃滤波
• 指数低通滤波
• 梯形低通滤波

一：理想低通滤波

（1）概述

理想低通滤波：是一种在频域中将高频信号滤除，保留低频信号的滤波器。其数学表达式为

$$H(u,v)=\{\begin{array}{ll}1& D(u,v)\le D_0\\ 0& D(u,v)>D_0\end{array}$$

其中，$H(u,v$)表示频域滤波器在$(u,v)$处的值，$D(u,v)=\sqrt{u^2+v^2}$表示频域中点$(u,v)$到中心点的距离，$D_0$为截止频率，即在此频率以下的信号通过，此频率以上的信号被滤除

具体来说，可以通过对图像进行傅里叶变换，将图像转换到频域中。然后，将理想低通滤波器与频域图像做乘积，得到滤波后的频域图像。最后，对滤波后的频域图像进行反傅里叶变换，得到空域中的滤波结果

理想低通滤波器具有截止频率清晰、滤波效果好等优点，但由于其截止频率突变，会导致滤波结果出现“振铃”现象，同时也不适用于实际应用中，因为实际信号中的高频分量很难完全截止

• “振铃”现象：$H(u,v)$在$D_9$处由1突变到0，其对应的冲激响应$h(x,y)$在空域中表现为同心环形式；同心环半径与$D_0$成反比；$D_0$越小，同心环半径越大，模糊越厉害

如下图为“振铃”现象

（2）程序

实现如下效果：截断频率不同的理想低通滤波器

---

matlab实现：

% 读入图像并显示
Image=imread('lena.bmp');
imshow(Image);

% 对图像进行傅里叶变换，并进行频谱搬移
FImage=fftshift(fft2(double(Image)));  

% 获取频域图像的大小
[N M]=size(FImage);

% 初始化一个全零的数组
g=zeros(N,M);

% 计算频谱的中心点坐标
r1=floor(M/2);  r2=floor(N/2);

% 在傅里叶频谱图上绘制4个圆形，并显示傅里叶频谱图像
figure;
imshow(log(abs(FImage)+1),[]),title('傅里叶频谱');
hold on
circle(r1,r2,5)
circle(r1,r2,11)
circle(r1,r2,45)
circle(r1,r2,68)

% 初始化截止频率数组
d0=[5 11 45 68];

% 循环对不同截止频率进行低通滤波并显示
for i=1:4
    for x=1:M
        for y=1:N
            % 计算当前像素点到频谱中心的距离
            d=sqrt((x-r1)^2+(y-r2)^2);
            % 根据距离和截止频率确定滤波器的值
            if d<=d0(i)
                h=1;
            else
                h=0;
            end
            % 用滤波器乘以频域图像
            g(y,x)=h*FImage(y,x);
        end
    end
    % 对滤波后的频域图像进行反傅里叶变换得到空域的图像，并显示
    g= real(ifft2(ifftshift(g)));
    figure,imshow(uint8(g)),title(['理想低通滤波D0=',num2str(d0(i))]);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读入图像并显示
Image = cv2.imread('lena.bmp', cv2.IMREAD_GRAYSCALE)
plt.imshow(Image, cmap='gray')
plt.show()

# 对图像进行傅里叶变换，并进行频谱搬移
FImage = np.fft.fftshift(np.fft.fft2(Image.astype(float)))

# 获取频域图像的大小
N, M = FImage.shape

# 初始化一个全零的数组
g = np.zeros((N, M))

# 计算频谱的中心点坐标
r1 = int(M/2)
r2 = int(N/2)

# 在傅里叶频谱图上绘制4个圆形，并显示傅里叶频谱图像
plt.imshow(np.log(np.abs(FImage)+1), cmap='gray')
plt.title('傅里叶频谱')
plt.hold(True)
circle1 = plt.Circle((r1,r2), 5, color='r', fill=False)
circle2 = plt.Circle((r1,r2), 11, color='r', fill=False)
circle3 = plt.Circle((r1,r2), 45, color='r', fill=False)
circle4 = plt.Circle((r1,r2), 68, color='r', fill=False)
plt.gca().add_artist(circle1)
plt.gca().add_artist(circle2)
plt.gca().add_artist(circle3)
plt.gca().add_artist(circle4)
plt.show()

# 初始化截止频率数组
d0 = [5, 11, 45, 68]

# 循环对不同截止频率进行低通滤波并显示
for i in range(4):
    for x in range(M):
        for y in range(N):
            # 计算当前像素点到频谱中心的距离
            d = np.sqrt((x-r1)**2 + (y-r2)**2)
            # 根据距离和截止频率确定滤波器的值
            if d <= d0[i]:
                h = 1
            else:
                h = 0
            # 用滤波器乘以频域图像
            g[y,x] = h * FImage[y,x]
    # 对滤波后的频域图像进行反傅里叶变换得到空域的图像，并显示
    g = np.real(np.fft.ifft2(np.fft.ifftshift(g)))
    plt.imshow(g.astype(np.uint8), cmap='gray')
    plt.title('理想低通滤波D0=' + str(d0[i]))
    plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

二：巴特沃斯低通滤波

（1）概述

巴特沃斯低通滤波：是一种用于去除图像中高频噪声的数字滤波方法。它是一种优化了幅频特性的滤波器，能够保留图像中更多的低频信息。该滤波器是基于巴特沃斯函数设计的，巴特沃斯函数是一种理想的低通滤波器，其频率响应具有非常陡峭的截止特性，但在实际应用中会出现一些缺陷。因此，巴特沃斯低通滤波器是基于巴特沃斯函数的特点，通过添加一些适当的补偿来消除其缺陷的一种优化滤波器

$$\{\begin{array}{l}H(u,v)=\frac{1}{1+{\left[D(u,v)/D_0\right]}^{2n}}\\ H(u,v)=\frac{1}{1+(\sqrt{2}-1){\left[D(u,v)/D_0\right]}^{2n}}\end{array}$$

当$D(u,v)=D_0$时，$H(u,v)$降为最大值的$\frac{1}{2}$或$\frac{1}{\sqrt{2}}$

与其他低通滤波器相比，巴特沃斯低通滤波器有以下特点：

• 具有更陡峭的截止特性，可以更好地去除高频噪声
• 可以通过调整阶数和截止频率来灵活地控制滤波器的性能
• 在频域中的实现相对简单，通常使用二阶滤波器级联实现

如下图所示，随着阶数$n$的增加，图像振铃效应越来越明显

（2）程序

实现如下效果：阶数不同的巴特沃斯低通滤波器

---

matlab实现：

Image=imread('lena.bmp'); % 读取lena.bmp图像
Image=imnoise(Image,'gaussian'); % 加入高斯噪声
imshow(Image); % 显示噪声图像
FImage=fftshift(fft2(double(Image)));  % 对加噪声图像进行傅里叶变换及频谱搬移
[N M]=size(FImage); % 获取傅里叶变换后频率矩阵的大小
g=zeros(N,M); % 用0矩阵初始化滤波后的频率矩阵
r1=floor(M/2);  r2=floor(N/2); % 计算中心点位置
figure;
imshow(log(abs(FImage)+1),[]),title('傅里叶频谱'); % 显示傅里叶频谱
hold on
circle(r1,r2,30) % 在频谱图上绘制圆形滤波器的半径
d0=30; % 设置圆形滤波器的半径
n=[1 2 3 4]; % 不同的Butterworth滤波器阶数
for i=1:4
    for x=1:M
        for y=1:N
            d=sqrt((x-r1)^2+(y-r2)^2); % 计算像素点到中心点的距离
            h=1/(1+(d/d0)^(2*n(i))); % 计算Butterworth滤波器的值
            g(y,x)=h*FImage(y,x); % 对傅里叶变换后的频率矩阵进行滤波
        end
    end
    g=ifftshift(g); % 频谱搬移
    g=real(ifft2(g)); % 逆傅里叶变换得到滤波后的图像
    figure,imshow(uint8(g)),title(['Butterworth低通滤波n=',num2str(n(i))]); % 显示滤波后的图像
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像并加入高斯噪声
Image = cv2.imread('lena.bmp', 0)
Image = cv2.GaussianBlur(Image, (5, 5), 0)
Image = Image.astype(np.float32)
Image += np.random.normal(0, 20, Image.shape)
Image = np.uint8(np.clip(Image, 0, 255))

# 显示噪声图像
plt.imshow(Image, cmap='gray')
plt.title('加噪声的图像')
plt.show()

# 对加噪声图像进行傅里叶变换并进行频谱搬移
FImage = np.fft.fftshift(np.fft.fft2(Image))

# 获取频域图像的大小
N, M = FImage.shape

# 初始化一个全零的数组
g = np.zeros((N, M))

# 计算频谱的中心点坐标
r1 = M // 2
r2 = N // 2

# 在傅里叶频谱图上绘制一个圆形，并显示傅里叶频谱图像
plt.imshow(np.log(np.abs(FImage) + 1), cmap='gray')
plt.title('傅里叶频谱')
plt.plot(r1, r2, 'wo', markerfacecolor='none', markersize=60)
plt.show()

# 初始化截止频率数组和Butterworth滤波器阶数数组
d0 = 30
n = [1, 2, 3, 4]

# 循环对不同Butterworth滤波器阶数进行低通滤波并显示
for i in range(4):
    for x in range(M):
        for y in range(N):
            # 计算当前像素点到频谱中心的距离
            d = np.sqrt((x - r1) ** 2 + (y - r2) ** 2)
            # 根据距离和截止频率计算Butterworth滤波器的值
            h = 1 / (1 + (d / d0) ** (2 * n[i]))
            # 用滤波器乘以频域图像
            g[y, x] = h * FImage[y, x]
    # 对滤波后的频域图像进行反傅里叶变换得到空域的图像，并显示
    g = np.real(np.fft.ifft2(np.fft.ifftshift(g)))
    plt.imshow(np.uint8(g), cmap='gray')
    plt.title('Butterworth低通滤波 n=' + str(n[i]))
    plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

三：指数低通滤波

（1）概述

指数低通滤波：是一种在频域中对图像进行滤波的方法，主要用于去除高频噪声。该滤波器的核心思想是将频域中高频分量的值进行衰减，从而达到降低图像的高频部分，保留图像的低频部分的效果

$$H(u,v)=e^{-\frac{D^2(u,v)}{2D_0^2}}$$

其中，$D(u,v)$是频域中像素点$(u,v)$到中心点的距离，$D_0$是滤波器的截止频率，也就是滤波器开始起作用的频率。$H(u,v)$是滤波器在频域的值，它与原始图像在频域中的傅里叶变换进行卷积，即可得到滤波后的图像。当$D(u,v)=D_0$时，$H(u,v)$降为最大值的0.607处

如下图所示，随着截止频率$D_0$的增加，图像振铃效应越来越不明显

（2）程序

实现如下效果：截止频率不同的指数低通滤波器

---

matlab实现：

Image=imread('lena.bmp');  %读取lena.bmp图像
Image=imnoise(Image,'gaussian');   %添加高斯噪声
imshow(Image);  %显示加噪声后的图像
FImage=fftshift(fft2(double(Image)));  %对图像进行傅里叶变换，并对频谱进行移位
[N M]=size(FImage);  %获取变换后频谱的大小
g=zeros(N,M);  %初始化滤波后的频谱
r1=floor(M/2);  r2=floor(N/2);  %获取频谱中心位置
figure;
imshow(log(abs(FImage)+1),[]),title('傅里叶频谱');  %显示傅里叶频谱
d0=[20 40];  %设置指数低通滤波器的参数
n=2;
for i=1:2  %遍历两个不同的截止频率
    for x=1:M
        for y=1:N
            d=sqrt((x-r1)^2+(y-r2)^2);  %计算距离
            h=exp(-0.5*(d/d0(i))^n);   %计算滤波器响应
            g(y,x)=h*FImage(y,x);  %对频谱进行滤波
        end
    end
    g=ifftshift(g);  %对滤波后的频谱进行移位
    g=real(ifft2(g));  %对滤波后的频谱进行逆傅里叶变换
figure,imshow(uint8(g)),title(['指数低通滤波D0=',num2str(d0(i))]);  %显示滤波后的图像
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取lena.bmp图像
Image = cv2.imread('lena.bmp', cv2.IMREAD_GRAYSCALE)

# 添加高斯噪声
Image = cv2.GaussianBlur(Image, (3, 3), 0)
Image = Image + np.uint8(np.random.normal(0, 20, size=Image.shape))

# 显示加噪声后的图像
plt.imshow(Image, cmap='gray')
plt.show()

# 对图像进行傅里叶变换，并对频谱进行移位
FImage = np.fft.fftshift(np.fft.fft2(np.double(Image)))

# 获取变换后频谱的大小
N, M = FImage.shape

# 初始化滤波后的频谱
g = np.zeros((N, M))

# 获取频谱中心位置
r1, r2 = M // 2, N // 2

# 显示傅里叶频谱
plt.imshow(np.log(np.abs(FImage) + 1), cmap='gray')
plt.title('傅里叶频谱')
plt.show()

# 设置指数低通滤波器的参数
d0 = [20, 40]
n = 2

# 遍历两个不同的截止频率
for i in range(2):
    for x in range(M):
        for y in range(N):
            # 计算距离
            d = np.sqrt((x - r1)**2 + (y - r2)**2)
            # 计算滤波器响应
            h = np.exp(-0.5 * (d / d0[i])**n)
            # 对频谱进行滤波
            g[y, x] = h * FImage[y, x]

    # 对滤波后的频谱进行移位
    g = np.fft.ifftshift(g)
    # 对滤波后的频谱进行逆傅里叶变换
    g = np.real(np.fft.ifft2(g))

    # 显示滤波后的图像
    plt.imshow(np.uint8(g), cmap='gray')
    plt.title('指数低通滤波D0=' + str(d0[i]))
    plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

四：梯形低通滤波

（1）概述

梯形低通滤波：是一种可以在一定程度上同时平衡频率响应的振铃和带宽的低通滤波器。它的频率响应类似于一个梯形，因此得名。梯形低通滤波器可以用于图像处理中的频域滤波，通过控制梯形上下两条边界的斜率和距离，可以调整其截止频率和带宽。梯形低通滤波器相比于其他低通滤波器，可以提供更平滑的过渡区域和更宽的过渡带宽，从而减少了图像的失真和信息的丢失

$$H(u,v)=\{\begin{array}{cc}1,& D(u,v)\le D_0\\ \frac{D(u,v)-D_1}{D_0-D_1},& D_0<D(u,v)\le D_1\\ 0& D(u,v)>D_1\end{array}$$

其中$D_0$为截止频率，$D_0、D_1$需要满足：$D_0<D_1$

如下图所示

（2）程序

实现如下效果截止频率不同的梯度低通滤波器

---

matlab实现：

Image=imread('lena.bmp');
Image=imnoise(Image,'gaussian');   %加入噪声
imshow(Image);
FImage=fftshift(fft2(double(Image)));  %傅里叶变换及频谱搬移
[N M]=size(FImage);
g=zeros(N,M);
r1=floor(M/2);  r2=floor(N/2);
%figure;
%imshow(log(abs(FImage)+1),[]),title('傅里叶频谱');
d0=[5 30];
d1=[45 70];

for i=1:2
    for x=1:M
        for y=1:N
            d=sqrt((x-r1)^2+(y-r2)^2);
            
            if d>d1  
                h=0;
            else
                if d>d0
                    h=(d-d1)/(d0-d1);
                else
                    h=1;
                end
            end            
            g(y,x)=h*FImage(y,x);
        end
    end
    g=ifftshift(g);
    g=real(ifft2(g));
    figure,imshow(uint8(g)),title(['梯度低通滤波D0=',num2str(d0(i)),',D1=',num2str(d1(i))]);
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

Python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像并添加高斯噪声
Image = cv2.imread('lena.bmp', cv2.IMREAD_GRAYSCALE)
Image = cv2.GaussianBlur(Image, (5, 5), 0)
Image = cv2.addWeighted(Image, 1.5, np.zeros(Image.shape, Image.dtype), 0, -0.5)

# 显示原始图像
plt.imshow(Image, cmap='gray')
plt.title('Original Image')
plt.show()

# 进行傅里叶变换及频谱搬移
FImage = np.fft.fftshift(np.fft.fft2(Image.astype(np.float64)))

# 获取图像尺寸
N, M = FImage.shape

# 初始化梯度低通滤波器参数
g = np.zeros((N, M), dtype=np.complex64)
r1 = int(M / 2)
r2 = int(N / 2)
d0 = [5, 30]
d1 = [45, 70]

# 循环处理两组参数
for i in range(2):
    # 计算滤波器系数
    for x in range(M):
        for y in range(N):
            d = np.sqrt((x - r1)**2 + (y - r2)**2)
            
            if d > d1[i]:
                h = 0
            else:
                if d > d0[i]:
                    h = (d - d1[i]) / (d0[i] - d1[i])
                else:
                    h = 1
                    
            g[y, x] = h * FImage[y, x]
    
    # 进行傅里叶逆变换并显示结果
    g = np.fft.ifftshift(g)
    g = np.real(np.fft.ifft2(g))
    plt.imshow(g, cmap='gray')
    plt.title(f'梯度低通滤波D0={d0[i]},D1={d1[i]}')
    plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51