数字图像处理第五次作业——频域滤波器

---

一、基本概念及原理

1. 理想低通滤波器：

理想低通滤波器在以原点为圆心、$D_0$为半径的园内，通过所有的频率，而在圆外截断所有的频率。（圆心的频率最低，为变换的直流分量）。函数如下：
$$H(s,v)=\{\begin{array}{r}1,D(u,v)\le D_0\\ 0,D(u,v)>D_0\end{array}$$

根据理想低通滤波器的函数画出其图像如下：

可以看出，理想低通滤波器的过渡非常急剧。

2. 巴特沃斯低通滤波器：

巴特沃斯低通滤波器的函数表达式如下，其中n称为其阶数：
$$H(u,v)=\frac{1}{1+(D(u,v)/D_0{)}^{2n}}$$

根据函数画出其图像如下：

从图中，我们可以看出它的过渡没有理想低通滤波器那么剧烈。同时阶数越高，过渡越剧烈，越像理想低通滤波器；当阶数越低时，过渡越平缓，越像高斯低通滤波器。

3. 高斯低通滤波器：

高斯低通滤波器的函数表达式如下：
$$H(u,v)=e^{\frac{-D^2(u,v)}{{}_{{2D_{{}_0}}^2}}}$$

根据函数画出其图像如下：

4. 高通滤波器：

高通滤波与低通滤波正好相反，是频域图像的高频部分通过而抑制低频部分。在图像中图像的边缘对应高频分量，因此高通滤波的效果是图像锐化。通常，我们用1减去一种低通滤波器，就可以得到相应的高通滤波器，反之亦然。

同样最简单的高通滤波器是理想高通滤波器。通过设置一个频率阈值，将高于该阈值的频率部分通过，而低于阈值的低频部分设置为0。

对于巴特沃斯滤波器，它的函数表达式除了可以用1减去其低通滤波的函数表达式之外，还可以直接将原低通表达式中的$D_0$和$D(u,v)$分子分母上下颠倒即可。

5. 拉普拉斯高通滤波器：

频域的拉普拉斯算子可以由如下的滤波器实现：

$$H(u,v)=-4{\pi }^2(u^2+v^2)$$

前提是$F(u,v)$的原点在进行图像变换之前已经通过执行运算$f(x,y)(-1{)}^{x+y}$中心化了，使得变换中心$(u,v)=(0,0)$就是频率矩形的中点$(M/2,N/2)$。否则

$$H(u,v)=-4{\pi }^2[(u-\frac{M}{2}{)}^2+(v-\frac{N}{2}{)}^2]$$

6. Unmask高通滤波器：

钝化模板由如下的表达式给出：

$$\begin{gathered} g_{max}(x,y)=f(x,y)-f_{LP}(x,y) \\ {f}_{LP}(x,y)={\Im }^{-1}[H_{LP}(u,v)F(u,v)] \end{gathered}$$

最后图像由下式给出：

$$g(x,y)=f(x,y)+k\cdot g_{max}(x,y)$$

当$k=1$时，为钝化模板；$k>1$时，为高频提升滤波器。

二、实现过程和结果分析

第一题：

频域低通滤波器：设计低通滤波器包括 butterworth and Gaussian (选择合适的半径，计算功率谱比),平滑测试图像test1和2;分析各自优缺点；

1. 巴特沃斯低通滤波：

Matlab代码：

function p=Butterworth(Img_in,D0,n)
[M,N]=size(Img_in);    M=2*M;N=2*N;
u = -M/2:(M/2-1);   v = -N/2:(N/2-1);
[u,v] = meshgrid(u,v);
D = sqrt(u.^2+v.^2);
%设计滤波器：
H = 1./(1+(D./D0).^(2*n)); 
Img_fft= fft2(Img_in,size(H,1),size(H,2));
Img_fft_shift = fftshift(Img_fft);
Img_Butterworth = Img_fft_shift.*H;
Img_out = ifft2(ifftshift(Img_Butterworth));
Img_out = Img_out(1:size(Img_in,1),1:size(Img_in,2));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);plot3(u,v,H);title('巴特沃斯滤波器');
subplot(2,2,4);imshow(Img_out);title(['滤波后的图像,D0=',num2str(D0),',n=',num2str(n)]);
%计算功率谱比：
S = 0;S1 = 0;
[P,Q] = size(Img_fft_shift);
for a = 1:P
    for b=1:Q
        S1 = S1+(abs(Img_Butterworth(a,b)))^2;
        Img_out = (abs(Img_fft_shift(a,b)))^2;
        S=S+Img_out;
    end
end
p = S1/S;
end

结果与分析：

功率谱比：$p=0.9984$

功率谱比：$p=0.9942$

分析： 通过参数为 $D_0=200,n=5$ 的巴特沃斯低通滤波器分别对两幅图片进行滤波，可以看出，两幅图像均变得更加光滑。对第一幅图像而言，功率谱比要稍大一些。

2. 高斯低通滤波器：

Matlab代码：

function p=Gaussian(Img_in,D0)
[M,N]=size(Img_in);    M=2*M;N=2*N;
u = -M/2:(M/2-1);   v = -N/2:(N/2-1);
[u,v] = meshgrid(u,v);
D = sqrt(u.^2+v.^2);
%设计滤波器：
H = exp(-D.^2/(2.*(D0.^2))); 
Img_fft= fft2(Img_in,size(H,1),size(H,2));
Img_fft_shift = fftshift(Img_fft);
Img_Gaussian = Img_fft_shift.*H;
Img_out = ifft2(ifftshift(Img_Gaussian));
Img_out = Img_out(1:size(Img_in,1),1:size(Img_in,2));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);plot3(u,v,H);title('高斯滤波器');
subplot(2,2,4);imshow(Img_out);title(['滤波后的图像,D0=',num2str(D0)]);
%计算功率谱比：
S = 0;S1 = 0;
[P,Q] = size(Img_fft_shift);
for a = 1:P
    for b=1:Q
        S1 = S1+(abs(Img_Gaussian(a,b)))^2;
        Img_out = (abs(Img_fft_shift(a,b)))^2;
        S=S+Img_out;
    end
end
p = S1/S;
end

结果与分析：

功率谱比：$p=0.9839$

功率谱比 $p=0.9839$

分析： 通过参数为 $D_0=100$ 的高斯低通滤波器分别对两幅图像进行滤波，可以看出，两幅图像均变得光滑，其滤波效果和巴特我司滤波器基本相当，且高斯滤波略胜一筹。但是两者在截止频率$D_0$相同时，由于在过渡带处存在差异，所以功率谱比有所不同，当巴特沃斯低通滤波器的n越大时，过渡带越陡峭，功率谱比与高斯滤波器的差异也会越大。同时，当将两个滤波器的$D_0$减小时，图像会越来越模糊，功率谱也会越来越小，即滤波后包含的低频分量越来越少。

---

第二题

频域高通滤波器：设计高通滤波器包括butterworth and Gaussian，在频域增强边缘。选择半径和计算功率谱比，测试图像test3,4；分析各自优缺点；

1. 巴特沃斯高通滤波器：

MATLAB代码：

function p=Butterworth(Img_in,D0,n)
[M,N]=size(Img_in);    M=2*M;N=2*N;
u = -M/2:(M/2-1);   v = -N/2:(N/2-1);
[u,v] = meshgrid(u,v);
D = sqrt(u.^2+v.^2);
%设计滤波器：
H = 1./(1+(D0./D).^(2*n)); 
Img_fft= fft2(Img_in,size(H,1),size(H,2));
Img_fft_shift = fftshift(Img_fft);
Img_Butterworth = Img_fft_shift.*H;
Img_out = ifft2(ifftshift(Img_Butterworth));
Img_out = Img_out(1:size(Img_in,1),1:size(Img_in,2));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);plot3(u,v,H);title('巴特沃斯滤波器');
subplot(2,2,4);imshow(Img_out);title(['滤波后的图像,D0=',num2str(D0), ',n=',num2str(n)]);
%计算功率谱比：
S = 0;S1 = 0;
[P,Q] = size(Img_fft_shift);
for a = 1:P
    for b=1:Q
        S1 = S1+(abs(Img_Butterworth(a,b)))^2;
        Img_out = (abs(Img_fft_shift(a,b)))^2;
        S=S+Img_out;
    end
end
p = S1/S;
end

结果与分析：

功率谱比：$p=0.0296$

功率谱比：$p=0.0157$

分析： 通过参数为 $D_0=50,n=5$ 的巴特沃斯高通滤波器分别对两幅图像进行滤波，可以看到，通过高通滤波将图像的低频分量滤去，将高频分量保留提取出来，所以可以用来做图像的边缘提取。

2. 高斯高通滤波器：

MATLAB代码：

function p=Gaussian(Img_in,D0)
[M,N]=size(Img_in);    M=2*M;N=2*N;
u = -M/2:(M/2-1);   v = -N/2:(N/2-1);
[u,v] = meshgrid(u,v);
D = sqrt(u.^2+v.^2);
%设计滤波器：
H = 1-exp(-D.^2/(2.*(D0.^2))); 
Img_fft= fft2(Img_in,size(H,1),size(H,2));
Img_fft_shift = fftshift(Img_fft);
Img_Gaussian = Img_fft_shift.*H;
Img_out = ifft2(ifftshift(Img_Gaussian));
Img_out = Img_out(1:size(Img_in,1),1:size(Img_in,2));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);plot3(u,v,H);title('高斯滤波器');
subplot(2,2,4);imshow(Img_out);title(['滤波后的图像,D0=',num2str(D0)]);
%计算功率谱比：
S = 0;S1 = 0;
[P,Q] = size(Img_fft_shift);
for a = 1:P
    for b=1:Q
        S1 = S1+(abs(Img_Gaussian(a,b)))^2;
        Img_out = (abs(Img_fft_shift(a,b)))^2;
        S=S+Img_out;
    end
end
p = S1/S;
end

结果与分析：

功率谱比：$p=0.0208$

功率谱比：$p=0.0365$

分析： 通过参数为 $D_0=20$ 的高斯高通滤波器，分别对两幅图像进行高通滤波，由于滤去了低频分量，保留了高频分量，所以提取出了图像的边缘。对比两种滤波器可以看出，两者效果基本相当，高斯滤波略优一些。同时，如果调整程序中的参数可以发现，当截止频率 $D_0$ 增加时，两种滤波器所得到的图像边缘越来越清晰，但是当 $D_0$ 到一定程度时，由于滤去能量过多，所以图像会整体呈现黑色。且高通滤波器在滤波时，会将直流分量一起滤除，导致图像变暗。

---

第三题

其他高通滤波器：拉普拉斯和Unmask，对测试图像test3,4滤波；分析各自优缺点；比较并讨论空域低通高通滤波（Project3）与频域低通和高通的关系；

1. 拉普拉斯高通滤波器：

MATLAB代码：

function Laplacian(Img_in)
[M,N]=size(Img_in);    M=2*M;N=2*N;
u = -M/2:(M/2-1);   v = -N/2:(N/2-1);
[u,v] = meshgrid(u,v);
D = sqrt(u.^2+v.^2);
%设计滤波器：
H = -4.*pi.^2*(u.^2+v.^2); 
Img_fft= fft2(Img_in,size(H,1),size(H,2));
Img_fft_shift = fftshift(Img_fft);
Img_Laplacian = Img_fft_shift.*H;
Img_out = ifft2(ifftshift(Img_Laplacian));
Img_out = Img_out(1:size(Img_in,1),1:size(Img_in,2));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);plot3(u,v,H);title('拉普拉斯滤波器');
subplot(2,2,4);imshow(Img_out);title('滤波后的图像');
end

结果与分析：

分析： 通过拉普拉斯高通滤波分别对两幅图像进行处理，可以看出滤波器的边缘增强效果，但这个效果与之前的两个高通滤波器所得到的效果并不相同，拉普拉斯得到的滤波后图像有明显的线条感，对肉眼而言不是太“友好”。

2. Unmask高通滤波器：

MATLAB代码：

function Unmask(Img_in,D0,n)
Img_fft_shift=fftshift(fft2(Img_in));
[M,N]=size(Img_fft_shift);
for i=1:M
    for j=1:N
        d=sqrt((i-fix(M/2))^2+(j-fix(N/2))^2);
        if d==0
            H(i,j)=0;
        else
            H(i,j)=1/(1+0.414*(D0/d)^(2*n));
        end
        Img_out(i,j)=(1+H(i,j))*Img_fft_shift(i,j);
    end
end
Img_out=real(ifft2(ifftshift(Img_out)));
%输出图像：
figure;
subplot(2,2,1);imshow(Img_in);title("原始图像");
subplot(2,2,2);imshow(log(1+abs(Img_fft_shift)),[]);title('频域图像');
subplot(2,2,3);
u=-M/2:(M/2-1);v=-N/2:(N/2-1);
[u,v]=meshgrid(u,v);plot3(u,v,H);title('Unmask滤波器');
subplot(2,2,4);imshow(Img_out);title(['滤波后的图像,D0=',num2str(D0),',n=',num2str(n)]);
end

结果与分析：

分析： 通过Unmask高通滤波器分别对两幅图像进行高通滤波，经过仔细观察还是可以发现其边缘增强效果的。这个效果又不同与前述的几种高通滤波器的效果，让图像边缘的对比度有了一定的提升，总体来说图像更加清晰了。

空域滤波与频域滤波：

空域滤波定义为滤波函数与输入图像进行卷积，而频域滤波则定义为滤波函数与输入图像的傅里叶变换进行相乘。两者的纽带则是卷积定理和傅里叶变换。频率域增强技术和空间域增强技术有密切的联系，一方面，许多空域增强技术可借助频域概念来分析和帮助设计；另一方面，许多空域增强技术也可以通过频域来实现，而频域增强也可以通过空域实现。空域主要包括平滑滤波和锐化滤波。平滑滤波是要滤除不规则的噪声或干扰的影响，从频域的角度来看，不规则的噪声具有较高的频率，所以可用具有低通能力的频域滤波器来滤除。

 
 
 
 
 
 
 
 

---

参考资料

[1]冈萨雷斯.数字图像处理(第三版)北京:电子工业出版社,2011

---

 
 
 

数字图像处理作业链接：
数字图像处理第一次作业——Bmp格式与基本变换
数字图像处理第二次作业——图像仿射变换
数字图像处理第三次作业——基于直方图的图像空域操作
数字图像处理第四次作业——空域滤波与边缘检测
数字图像处理第五次作业——频域滤波器
数字图像处理第六次作业——图像噪声和恢复
数字图像处理第七次作业——图像直线检测