数字图像处理三级项目
—高通、
低通、带通滤波器
摘要
在图像处理的过程中,消除图像的噪声干扰是一个非常重要的问题。利用matlab软件,采用频域滤波的方式,对图像进行低通和高通滤波处理。低通滤波是要保留图像中的低频分量而除去高频分量，由于图像中的边缘和噪声都对应图像傅里叶频谱中的高频部分，所以低通滤波可以除去或消弱噪声的影响并模糊边缘轮廓；高通滤波是要保留图像中的高频分量而除去低频分量，所以高通滤波可以保留较多的边缘轮廓信息。低通滤波器有巴特沃斯滤波器和高斯滤波器等等，本次设计使用的低通滤波器为****。高通滤波器有巴特沃斯滤波器、高斯滤波器、Laplacian高通滤波器以及Unmask高通滤波器等等，本次设计使用巴特沃斯高通滤波器。
1、 频域低通滤波器：设计低通滤波器包括 butterworth and Gaussian (选择合适的半径，计算功率谱比),平滑测试图像test1和2。
实验原理分析
根据卷积定理，两个空间函数的卷积可以通过计算两个傅立叶变换函数的乘积的逆变换得到，如果f(x, y)和h(x, y)分别代表图像与空间滤波器，F(u, v)和H(u, v)分别为响应的傅立叶变换（H(u, v)又称为传递函数），那么我们可以利用卷积定理来进行频域滤波。
在频域空间，图像的信息表现为不同频率分量的组合。如果能让某个范围内的分量或某些频率的分量受到抑制，而让其他分量不受影响，就可以改变输出图的频率分布，达到不同的增强目的。 频域空间的增强方法的步骤：
（1）将图像从图像空间转换到频域空间； （2）在频域空间对图像进行增强；
（3）将增强后的图像再从频域空间转换到图像空间。
低通滤波是要保留图像中的低频分量而除去高频分量。图像中的边缘和噪声都对应图像傅里叶频谱中的高频部分，所以低通滤波可以除去或消弱噪声的影响并模糊边缘轮廓。理想低通滤波器具有传递函数：
其中D0为制定的非负数，D(u,v)为点(u,v)到滤波器中心的距离。
功率谱比的定义：
L?Pg(u,v)Pf(u,v)
其中，Pf(u,v)为滤波前图像的功率谱，Pg(u,v)为滤波后图像的功率谱。
22频率计算公式为：Pf(u,v)?F(u,v)，Pg(u,v)?G(u,v)。
① Butterworth滤波器设计：
理想低通滤波器在数学上定义得很清楚，在计算机模拟中也可实现，但在截断频率处直上直下的理想低通滤波器是不能用实际的电子器件实现的。
n阶Butterworth低通滤波器（BLPF）的传递函数（截止频率距原点的距D0离为 ）定义如下：
1 H(u,v)?1??D(u,v)/D0?2n （1）
22其中，D(u，v)=（u-M/2）。 （2） +（v-N/2）不同于ILPF,BLPF变换函数在通带与被滤除的频率之间没有明显的截断。对
于有平滑传递函数的滤波器，定义一个截止频率的位置并使H(u，v)幅度降到其最大值的一部分。在式（1）中，当D（u，v）=D0时，H（u，v）=0.5（从最大值降到它的50%）。
一阶的巴特沃斯滤波器没有振铃，在二阶中振铃通常很微小，这是因为与理想低通滤波器相比，它的通带与阻带之间没有明显的跳跃，在高低频率间的过渡比较光滑。巴特沃斯低通滤波器的处理结果比理想滤波器的要好，但阶数增高时振铃便成为一个重要因素。本次实验中设计实现了二阶巴特沃斯滤波器。
根据以上原理设计Butterworth低通滤波器，其处理结果如下图示：
源图像test1.pgmButterworth低通滤波图test1.pgm
源图像test2.tifButterworth低通滤波图test2.tif
理想低通滤波器有明显的振铃现象，而巴特沃斯滤波器的效果较好。 计算得test1的功率谱比L=0.9939。test2的功率谱比为0.9902。 ②Gaussian滤波器设计：
二维高斯低通滤波器，其传递函数的形式为： H(u，v)=1-e?D2(u,v)/2?2 （3）
22其中，D(u，v)=（u-M/2）。?表示高斯曲线扩展的程度。使?=D0，+（v-N/2）可以将滤波器表示为：
H(u，v)=e?D2(u,v)/2D02 （4）
其中，D0是截止频率。当D（u，v）=D0时，滤波器下降到它最大值的0.607倍处。
由于高斯低通滤波器的傅里叶反变换也是高斯的，这就是说通过式（3）或式（4）的傅里叶反变换而得到的空间高斯滤波器将没有振铃。
根据以上分析，设计Gaussian低通滤波器，处理结果如下：
源图像test1.pgmGaussian低通滤波后的test1.pgm(r=5)
Gaussian低通滤波后的test1.pgm(r=15)Gaussian低通滤波后的test1.pgm(r=80)源图像test2.tifGaussian低通滤波后的test2.tif(r=15)Gaussian低通滤波后的test1.pgm(r=80)Gaussian低通滤波后的test1.pgm(r=230)Gaussian低通滤波后的test2.tif(r=5)Gaussian低通滤波后的test2.tif(r=30)
Gaussian低通滤波后的test2.tif(r=80)Gaussian低通滤波后的test2.tif(r=230)
可见，当滤波器的半径不同时，对应的滤波效果也不同。半径越小，平滑效果越明显，但半径过小，会使得图像变得模糊不清。
计算得test1（r=5）的功率谱比L= 0.4674。test2（r=5）的功率谱比为L=0.2930。
2、 频域高通滤波器：设计高通滤波器包括butterworth and Gaussian，在频域增强边缘。选择半径和计算功率谱比，测试图像test3,4：
实验原理分析
高通滤波是要保留图像中的高频分量而除去低频分量。理想高通滤波器传递函数表示为： 如D(u,v)?D0?0H(u,v)??
1如D(u,v)?D0?① Butterworth滤波器设计：
n阶Butterworth高通滤波器（BLPF）的传递函数（截止频率距原点的距D0离为 ）定义如下：
1 H(u,v)?2n1??D0D(u,v)? （5）
22其中， D(u，v)=（u-M/2）。 （6） +（v-N/2）不同于ILPF,BLPF变换函数在通带与被滤除的频率之间没有明显的截断。对
于有平滑传递函数的滤波器，定义一个截止频率的位置并使H(u，v)幅度降到其最大值的一部分。在式（1）中，当D（u，v）=D0时，H（u，v）=0.5（从最大值降到它的50%）。
根据以上原理设计Butterworth高通滤波器，其处理结果如下图示：
源图像test3corrupt.pgmButterworth高通滤波图test3corrupt.pgm
源图像test4 copy.bmpButterworth高通滤波图test4 copy.bmp
计算得test3的功率谱比为0.0851。test4的功率谱比为0.0547。
②Gaussian滤波器设计：
二维高斯高通滤波器，其传递函数的形式为：
H(u，v)=1-e?D2(u,v)/2?2 （7）
22其中，D(u，v)=（u-M/2）。?表示高斯曲线扩展的程度。使?=D0，+（v-N/2）可以将滤波器表示为：
H(u，v)=1-e?D2(u,v)/2D02 （8）
其中，D0是截止频率。当D（u，v）=D0时，滤波器下降到它最大值的0.607倍处。
由于高斯低通滤波器的傅里叶反变换也是高斯的，这就是说通过式（7）或式（8）的傅里叶反变换而得到的空间高斯滤波器将没有振铃。
根据以上分析，设计Gaussian高通滤波器，处理结果如下：
源图像test3corrupt.pgmGaussian高通滤波后的test3corrupt.pgm(r=5)
Gaussian高通滤波后的test3corrupt.pgm(r=15)Gaussian高通滤波后的test3corrupt.pgm(r=30)
源图像test4 copy.bmpGaussian高通滤波后的test4 copy.bmp(r=5)
Gaussian高通滤波后的test4 copy.bmp(r=15)Gaussian高通滤波后的test4 copy.bmp(r=30)
可见，当滤波器的半径不同时，对应的滤波效果也不同。半径越小，边缘效果越明显。
一般图像中的大部分能量集中在低频分量里，高通滤波会将很多低频分量滤除，导致增强图中边缘得到加强但光滑区域灰度减弱变暗甚至接近黑色。为解决这个问题，可对频域里的高通滤波器的转移函数加一个常数以将一些低频分量加回去，获得既保持光滑区域又改善边缘区域对比度的效果。这样得到的滤波器称为高频增强滤波器。 计算得test3（r=5）的功率谱比L= 0.0591，test4（r=5）的功率谱比为L= 0.4449。
3、其他高通滤波器：拉普拉斯和Unmask，对测试图像test3,4滤波；比较并讨论空域低通高通滤波（Project3）与频域低通和高通的关系；
实验原理分析
①拉普拉斯高通滤波器 公式表示如下：
dnf(x)?[]?(j2?u)nF(u) （9） ndx从这个简单的表达式可以得到：
?2f(x,y)?2f(x,y)?[?]??4?2(u2?v2)F(u,v) （10） 22?x?y所以， ?[?2f(x,y)]??4?2(u2?v2)F(u,v) （11） 即频域的拉普拉斯算子可以有如下滤波器实现：
H(u,v)??4?2(u2?v2) （12）
前提是F(u，v)的原点在进行图像变换之前已通过执行运算f(x,y)(?1)x?y中心化了，使得变换中心（u,v）=(0,0)就是频率矩形的中点（M/2，N/2）。因此
H(u,v)??4?2[(u?M/2)2?(v?N/2)2]。
根据以上分析，设计拉普拉斯算子高通滤波器，处理结果如下：
源图像test3corrupt.pgmLaplacian高通滤波图test3corrupt.pgm
源图像test4 copy.bmpLaplacian高通滤波图test4 copy.bmp
由于拉普拉斯高通滤波器将原始图像完全加回到滤波后的结果中，因此解决了Butterworth滤波器和Gaussian滤波器除去了傅里叶变换的零频率成分的问题，从而使得滤波后的图像其背景的平均强度增加、变亮。但同时引入了噪声干扰，使得滤波后的图像有一定程度的失真。
②Unsharp masking高通滤波器
Unsharp masking高通滤波器模板由以下公式确定：
x,y?) gmas(kf(x,?y)fx ,y （13） LP( fLP(x,y)???1[HLP(u,v)F(u,v)] （14）
x,y?)f(x,?y)kma*sg g(kx ( y （15）
当K=1时，为钝化模板；K>1时，为高频提升滤波器。由以上算法设计Unsharp masking高通滤波器，其中gmask(x,y)使用Butterworth滤波算法实现，处理结果如下图示：
源图像test3corrupt.pgm反锐化掩膜滤波图test3corrupt.pgm
源图像test4 copy.bmp反锐化掩膜滤波图test4 copy.bmp
可见，反锐化掩膜后的图像边缘信息更加清晰，但同时带来了过度锐化的
问题，出现了多重轮廓。
空域低通高通滤波与频域低通和高通的关系： 空域滤波主要包括平滑滤波和锐化滤波。平滑滤波是要滤除不规则的噪声或干扰的影响，从频域角度看，不规则噪声具有较高的频率，所以可用具有低通能力的频域滤波器来滤除。由此可见，空域的平滑滤波对应频域的低通滤波。锐化滤波是要增强边缘和轮廓处的强度，从频域角度看，边缘和轮廓处都具有较高的频率，所以可用具有高通能力的频域滤波器来增强，由此可见空域的锐化滤波对应频域的高通滤波。
附录
一、 参考文献
[1] 冈萨雷斯著.数字图像处理(第三版).北京：电子工业出版社，2010
[2] 杨杰 李庆著.数字图像处理及MATLAB实现——学习与实验指导.北京：电子工业出版社，2010
[3] 苏金明 王永利著.MATLAB图形图像. 北京：电子工业出版社，2005
[4] 朱习军 隋思涟等著.MATLAB在信号与图像处理中的应用. 北京：电子工业出版社，2009
[5]百度文库.http://wenku.http://www.32336.cn//view/1e4643380912a2161479292d.html
二、源代码：
频域低通滤波器
1. Butterworth低通滤波器（以处理test1.pgm为例）
I= imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test1.pgm','pgm'); figure;
subplot(1,2,1); imshow(I);
title('源图像test1.pgm'); f=double(I);
g=fft2(f); % 傅立叶变换 g=fftshift(g); % 转换数据矩阵 [M,N]=size(g);
nn=2; % 二阶巴特沃斯(Butterworth)低通滤波器 d0=50;
m=fix(M/2); n=fix(N/2); for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
h1=1/(1+0.414*(d/d0)^(2*nn)); % 计算低通滤波器传递函数 result1(i,j)=h1*g(i,j); end end
result1=ifftshift(result1); J2=ifft2(result1); J3=uint8(real(J2)); subplot(1,2,2); imshow(J3);
title('低通滤波图test1.pgm'); % 显示滤波处理后的图像 S=0; S1=0; for i=1:M
for j=1:N
L=(abs(result1(I,j)))^2; %计算结果图像的功率谱 S=S+L; end end
for i=1:M
for j=1:N
L1=(abs(g(I,j)))^2; %计算源图像的功率谱 S1=S1+L1; end end
L=S/S1 %计算功率谱比
2、Gaussian低通滤波器（以处理test1.pgm为例）
I=imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test1.pgm'); subplot(1,2,1); imshow(I);
title('源图像test1.pgm'); r=5;
Im=double(I); F=fft2(Im);
F_result=fftshift(F); g=F_result;
[m,n]=size(F_result); M=fix(m/2); N=fix(n/2); for u=1:m for v=1:n
D=sqrt((u-M)^2+(v-N)^2); H=exp(-D^2/(2*r^2));
F_result(u,v)=F_result(u,v)*H; end end
G_result=ifftshift(F_result); g_result=ifft2(G_result); f=real(g_result); f=uint8(f); subplot(1,2,2); imshow(f);
title('Gaussian低通滤波后的test1.pgm(r=5)'); S=0; S1=0; for i=1:M
for j=1:N
L=(abs(F_result(i,j)))^2; %计算结果图像的功率谱 S=S+L; end end
for i=1:M
for j=1:N
L1=(abs(g(i,j)))^2; %计算源图像的功率谱 S1=S1+L1; end end
L=S/S1 %计算功率谱比
频域高通滤波器
1、 Butterworth滤波器（以处理test3_corrupt.pgm为例）
I= imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test3_corrupt.pgm','pgm');
figure;
subplot(1,2,1); imshow(I);
title('源图像test3_corrupt.pgm');
f=double(I); % 数据类型转换，MATLAB不支持图像的无符号整型的计算 g=fft2(f); % 傅立叶变换 g=fftshift(g); % 转换数据矩阵 [M,N]=size(g);
nn=2; % 二阶巴特沃斯(Butterworth)高通滤波器 d0=5;
m=fix(M/2); n=fix(N/2); for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2); if (d==0) h2=0; else
h2=1/(1+0.414*(d0/d)^(2*nn));% 计算传递函数 end result2(i,j)=h2*g(i,j); end end
result3=ifftshift(result2); J4=ifft2(result3); J5=uint8(real(J4)); subplot(1,2,2); imshow(J5);
title('高通滤波图test3_corrupt.pgm'); % 滤波后图像显示 S=0; S1=0; for i=1:M
for j=1:N
L=(abs(result2 (i,j)))^2; %计算结果图像的功率谱 S=S+L; end end
for i=1:M
for j=1:N
L1=(abs(g(i,j)))^2; %计算源图像的功率谱 S1=S1+L1; end end
L=S/S1 %计算功率谱比
2、Gaussian滤波器（以处理test3_corrupt.pgm为例）
I=imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test3_corrupt.pgm'); subplot(1,2,1); imshow(I);
title('源图像test1.pgm'); r=5;
Im=double(I); F=fft2(Im);
F_result=fftshift(F);
g= F_result；
[m,n]=size(F_result); M=fix(m/2); N=fix(n/2); for u=1:m for v=1:n
D=sqrt((u-M)^2+(v-N)^2); H=1-exp(-D^2/(2*r^2));
F_result(u,v)=F_result(u,v)*H; end end
G_result=ifftshift(F_result); g_result=ifft2(G_result); f=real(g_result); f=uint8(f); subplot(1,2,2); imshow(f);
title('Gaussian低通滤波后的test3_corrupt.pgm (r=5)'); S=0; S1=0; for i=1:M
for j=1:N
L=(abs(F_result (i,j)))^2; %计算结果图像的功率谱 S=S+L; end end
for i=1:M
for j=1:N
L1=(abs(g(i,j)))^2; %计算源图像的功率谱 S1=S1+L1; end end
L=S/S1 %计算功率谱比
3、 Laplacian滤波器（以处理test3_corrupt.pgm为例）
I= imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test3_corrupt.pgm','pgm'); figure;
subplot(1,2,1); imshow(I);
title('源图像test3_corrupt.pgm'); f=double(I);
g=fft2(f); % 傅立叶变换 g=fftshift(g); % 转换数据矩阵 [M,N]=size(g);
m=fix(M/2); n=fix(N/2); for i=1:M
for j=1:N
H=1+4*pi^2*((i-m)^2+(j-n)^2); % 计算高通滤波器传递函数 result1(i,j)=H*g(i,j); end end
result1=ifftshift(result1); J2=ifft2(result1); J3=uint8(real(J2)); subplot(1,2,2); imshow(J3);
title('Laplacian高通滤波图test3_corrupt.pgm'); % 显示滤波处理后的图像
4、Unmask滤波器（以处理test4 copy.bmp为例）
I= imread('E:\\大三下\\图像处理英文课件\\作业\\第五次作业\\test4 copy.bmp','bmp'); figure;
subplot(1,2,1); imshow(I);
title('源图像test4 copy.bmp');
f=double(I); % 数据类型转换，MATLAB不支持图像的无符号整型的计算 g=fft2(f); % 傅立叶变换 g=fftshift(g); % 转换数据矩阵 [M,N]=size(g);
nn=2; d0=5;
m=fix(M/2); n=fix(N/2); for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
if (d==0) h2=0; else
h2=1/(1+0.414*(d0/d)^(2*nn));% 计算传递函数 end
result2(i,j)=(1+h2)*g(i,j); end end
result2=ifftshift(result2); J4=ifft2(result2); J5=uint8(real(J4)); subplot(1,2,2); imshow(J5);
title('反锐化掩膜滤波图test4 copy.bmp');
% 滤波后图像显示