1、数字图像处理作业频域滤波器设计摘要在图像处理的过程中,消除图像的噪声干扰是一个非常重要的问题。本文利用 matlab 软件,采用频域滤波的方式,对图像进行低通和高通滤波处理。低通滤波是要保留图像中的低频分量而除去高频分量，由于图像中的边缘和噪声都对应图像傅里叶频谱中的高频部分，所以低通滤波可以除去或消弱噪声的影响并模糊边缘轮廓；高通滤波是要保留图像中的高频分量而除去低频分量，所以高通滤波可以保留较多的边缘轮廓信息。本文使用的低通滤波器有巴特沃斯滤波器和高斯滤波器，使用的高通滤波器有巴特沃斯滤波器、高斯滤波器、Laplacian 高通滤波器以及 Unmask 高通滤波器。实际应用中应该根据实际图
2、像中包含的噪声情况灵活地选取适当的滤波算法。1、 频域低通滤波器：设计低通滤波器包括 butterworth and Gaussian (选择合适的半径，计算功率谱比),平滑测试图像 test1 和 2。实验原理分析根据卷积定理，两个空间函数的卷积可以通过计算两个傅立叶变换函数的乘积的逆变换得到，如果 f(x, y)和 h(x, y)分别代表图像与空间滤波器，F(u, v)和 H(u, v)分别为响应的傅立叶变换（H(u, v)又称为传递函数） ，那么我们可以利用卷积定理来进行频域滤波。在频域空间，图像的信息表现为不同频率分量的组合。如果能让某个范围内的分量或某些频率的分量受到抑制，而让其他分
3、量不受影响，就可以改变输出图的频率分布，达到不同的增强目的。频域空间的增强方法的步骤：（1）将图像从图像空间转换到频域空间；（2）在频域空间对图像进行增强；（3）将增强后的图像再从频域空间转换到图像空间。低通滤波是要保留图像中的低频分量而除去高频分量。图像中的边缘和噪声都对应图像傅里叶频谱中的高频部分，所以低通滤波可以除去或消弱噪声的影响并模糊边缘轮廓。理想低通滤波器具有传递函数：其中 D0 为制定的非负数，D(u,v)为点(u,v)到滤波器中心的距离。功率谱比的定义： g(,)fPuvL其中， 为滤波前图像的功率谱， 为滤波后图像的功率谱。(,)fPuvg(,)uv频率计算公式为： ， 。2
4、(,)(,)fvFu2(,)(,)gPvG Butterworth 滤波器设计：理想低通滤波器在数学上定义得很清楚，在计算机模拟中也可实现，但在截断频率处直上直下的理想低通滤波器是不能用实际的电子器件实现的。n 阶 Butterworth 低通滤波器（BLPF）的传递函数（截止频率距原点的距离为 ）定义如下：（1）其中 ， 。 （2）22(uv)=+v-N/D， （ -M/） （ ）不同于 ILPF,BLPF 变换函数在通带与被滤除的频率之间没有明显的截断。对于有平滑传递函数的滤波器，定义一个截止频率的位置并使 H(u，v)幅度降到其最大值的一部分。在式（1）中，当 D（u，v）=D0 时，H
5、（u，v）=0.5（从最大值降到它的 50%） 。一阶的巴特沃斯滤波器没有振铃，在二阶中振铃通常很微小，这是因为与理想低通滤波器相比，它的通带与阻带之间没有明显的跳跃，在高低频率间的过渡比较光滑。巴特沃斯低通滤波器的处理结果比理想滤波器的要好，但阶数增高时振铃便成为一个重要因素。本次实验中设计实现了二阶巴特沃斯滤波器。根据以上原理设计 Butterworth 低通滤波器，其处理结果如下图示：低低低 test1.pgm Butterworth低低低低低 test1.pgm低低低 test2.tif Butterworth低低低低低 test2.tif理想低通滤波器有明显的振铃现象，而巴特沃斯滤波
6、器的效果较好。计算得 test1 的功率谱比 L=0.9939。test2 的功率谱比为 0.9902。Gaussian 滤波器设计：二维高斯低通滤波器，其传递函数的形式为：nDuH0),(1),(0（3）22(,)/(uv)=1-eDuvH，其中， 。 表示高斯曲线扩展的程度。使2(uv)=+N/D， （ -M/） （ ）=D0，可以将滤波器表示为：（4）220(,)/(uv)=eDuvH，其中，D0 是截止频率。当 D（u，v）=D0 时，滤波器下降到它最大值的 0.607倍处。由于高斯低通滤波器的傅里叶反变换也是高斯的，这就是说通过式（3）或式（4）的傅里叶反变换而得到的空间高斯滤波器将
7、没有振铃。根据以上分析，设计 Gaussian 低通滤波器，处理结果如下：低低低 test1.pgm Gaussian低低低低低低 test1.pgm(r=5)Gaussian低低低低低低test1.pgm(r=15) Gaussian低低低低低低test1.pgm(r=80)Gaussian低低低低低低test1.pgm(r=80) Gaussian低低低低低低test1.pgm(r=230)低低低 test2.tif Gaussian低低低低低低 test2.tif(r=5)Gaussian低低低低低低 test2.tif(r=15) Gaussian低低低低低低 test2.tif(r=
8、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：实验原理分析高通滤波是要保留图像中的高频分量而除去低频分量。理想高通滤波器传递函数表示为：
9、 Butterworth 滤波器设计：n 阶 Butterworth 高通滤波器（BLPF）的传递函数（截止频率距原点的距离为 ）定义如下：（5）其中 ， 。 （6）22(uv)=+v-N/D， （ -M/） （ ）不同于 ILPF,BLPF 变换函数在通带与被滤除的频率之间没有明显的截断。对于有平滑传递函数的滤波器，定义一个截止频率的位置并使 H(u，v)幅度降到其最大值的一部分。在式（1）中，当 D（u，v）=D0 时，H（u，v）=0.5（从最大值降到它的 50%） 。根据以上原理设计 Butterworth 高通滤波器，其处理结果如下图示：低低低 test3corrupt.pgm Bu
10、tterworth低低低低低 test3corrupt.pgm0 nuH0),(1,0),(10),(Dvuvu如如低低低 test4 copy.bmp Butterworth低低低低低 test4 copy.bmp计算得 test3 的功率谱比为 0.0851。test4 的功率谱比为 0.0547。Gaussian 滤波器设计：二维高斯高通滤波器，其传递函数的形式为：（7 ）22(,)/(uv)=1-eDuvH，其中， 。 表示高斯曲线扩展的程度。使2(uv)=+N/D， （ -M/） （ ）=D0，可以将滤波器表示为：（8）220(,)/(uv)=1-eDuvH，其中，D0 是截止频率。
11、当 D（u，v）=D0 时，滤波器下降到它最大值的 0.607倍处。由于高斯低通滤波器的傅里叶反变换也是高斯的，这就是说通过式（7）或式（8）的傅里叶反变换而得到的空间高斯滤波器将没有振铃。根据以上分析，设计 Gaussian 高通滤波器，处理结果如下：低低低 test3corrupt.pgm Gaussian低低低低低低 test3corrupt.pgm(r=5)Gaussian低低低低低低test3corrupt.pgm(r=15) Gaussian低低低低低低test3corrupt.pgm(r=30)低低低 test4 copy.bmp Gaussian低低低低低低 test4 cop
12、y.bmp(r=5)Gaussian低低低低低低test4 copy.bmp(r=15) Gaussian低低低低低低test4 copy.bmp(r=30)可见，当滤波器的半径不同时，对应的滤波效果也不同。半径越小，边缘效果越明显。一般图像中的大部分能量集中在低频分量里，高通滤波会将很多低频分量滤除，导致增强图中边缘得到加强但光滑区域灰度减弱变暗甚至接近黑色。为解决这个问题，可对频域里的高通滤波器的转移函数加一个常数以将一些低频分量加回去，获得既保持光滑区域又改善边缘区域对比度的效果。这样得到的滤波器称为高频增强滤波器。计算得 test3（r=5 ）的功率谱比 L= 0.0591，test4
13、（r=5）的功率谱比为 L= 0.4449。3、其他高通滤波器：拉普拉斯和 Unmask，对测试图像 test3,4 滤波；比较并讨论空域低通高通滤波（Project3）与频域低通和高通的关系；实验原理分析拉普拉斯高通滤波器公式表示如下：（9）()2)(nndfxjuF从这个简单的表达式可以得到：（10）222(,)(,)4()(,fxyfxyvu所以， （11）22(,),fuF即频域的拉普拉斯算子可以有如下滤波器实现：（12）2(,)4()Hvv前提是 F(u，v)的原点在进行图像变换之前已通过执行运算 中心化(,)1xyf了，使得变换中心（u,v）=(0,0)就是频率矩形的中点（ M/2
14、，N/2） 。因此。222(,)4(/)(/)uvuMvN根据以上分析，设计拉普拉斯算子高通滤波器，处理结果如下：低低低 test3corrupt.pgm Laplacian低低低低低 test3corrupt.pgm低低低 test4 copy.bmp Laplacian低低低低低 test4 copy.bmp由于拉普拉斯高通滤波器将原始图像完全加回到滤波后的结果中，因此解决了 Butterworth 滤波器和 Gaussian 滤波器除去了傅里叶变换的零频率成分的问题，从而使得滤波后的图像其背景的平均强度增加、变亮。但同时引入了噪声干扰，使得滤波后的图像有一定程度的失真。Unsharp m
15、asking 高通滤波器Unsharp masking 高通滤波器模板由以下公式确定：（13）g(,)(,)(,)maskLPxyffxy（14）1LPfHuvF（15）(,)(,)*g(,)maskxyxy当 K=1 时，为钝化模板；K1 时，为高频提升滤波器。由以上算法设计 Unsharp masking 高通滤波器，其中 使用 Butterworth 滤波算法实现，处理结(,)ask果如下图示：低低低 test3corrupt.pgm 低低低低低低低低 test3corrupt.pgm低低低 test4 copy.bmp 低低低低低低低低 test4 copy.bmp可见，反锐化掩膜后的
16、图像边缘信息更加清晰，但同时带来了过度锐化的问题，出现了多重轮廓。空域低通高通滤波与频域低通和高通的关系：空域滤波主要包括平滑滤波和锐化滤波。平滑滤波是要滤除不规则的噪声或干扰的影响，从频域角度看，不规则噪声具有较高的频率，所以可用具有低通能力的频域滤波器来滤除。由此可见，空域的平滑滤波对应频域的低通滤波。锐化滤波是要增强边缘和轮廓处的强度，从频域角度看，边缘和轮廓处都具有较高的频率，所以可用具有高通能力的频域滤波器来增强，由此可见空域的锐化滤波对应频域的高通滤波。附录一、 参考文献1 冈萨雷斯著.数字图像处理( 第三版).北京：电子工业出版社， 20102 杨杰 李庆著.数字图像处理及 MA
17、TLAB 实现学习与实验指导.北京：电子工业出版社，20103 苏金明 王永利著.MATLAB 图形图像. 北京：电子工业出版社， 20054 朱习军 隋思涟等著.MATLAB 在信号与图像处理中的应用. 北京：电子工业出版社，20095 百度文库.http:/ 百度文库.http:/ 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)
18、; % 傅立叶变换g=fftshift(g); % 转换数据矩阵M,N=size(g);nn=2; % 二阶巴特沃斯(Butterworth)低通滤波器d0=50;m=fix(M/2); n=fix(N/2);for i=1:Mfor j=1:Nd=sqrt(i-m)2+(j-n)2);h1=1/(1+0.414*(d/d0)(2*nn); % 计算低通滤波器传递函数result1(i,j)=h1*g(i,j);endendresult1=ifftshift(result1);J2=ifft2(result1);J3=uint8(real(J2);subplot(1,2,2);imshow(J
19、3);title(低通滤波图 test1.pgm); % 显示滤波处理后的图像S=0;S1=0;for i=1:Mfor j=1:NL=(abs(result1(I,j)2; %计算结果图像的功率谱S=S+L;endendfor i=1:Mfor j=1:NL1=(abs(g(I,j)2; %计算源图像的功率谱S1=S1+L1;endendL=S/S1 %计算功率谱比2、Gaussian 低通滤波器（以处理 test1.pgm 为例）I=imread(E:大三下图像处理英文课件 作业第五次作业test1.pgm);subplot(1,2,1);imshow(I);title(源图像 test1
20、.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:mfor v=1:nD=sqrt(u-M)2+(v-N)2);H=exp(-D2/(2*r2);F_result(u,v)=F_result(u,v)*H;endendG_result=ifftshift(F_result);g_result=ifft2(G_result);f=real(g_result);f=uint8(f);subplot(1,2,2);imsh
21、ow(f);title(Gaussian 低通滤波后的 test1.pgm(r=5);S=0;S1=0;for i=1:Mfor j=1:NL=(abs(F_result(i,j)2; %计算结果图像的功率谱S=S+L;endendfor i=1:Mfor j=1:NL1=(abs(g(i,j)2; %计算源图像的功率谱S1=S1+L1;endendL=S/S1 %计算功率谱比频域高通滤波器1、 Butterworth 滤波器（以处理 test3_corrupt.pgm 为例）I= imread(E:大三下图像处理英文课件作业 第五次作业test3_corrupt.pgm,pgm);figur
22、e;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:Mfor j=1:Nd=sqrt(i-m)2+(j-n)2); if (d=0)h2=0;elseh2=1/(1+0.414*(d0/d)(2*nn);% 计
23、算传递函数endresult2(i,j)=h2*g(i,j);endendresult3=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:Mfor j=1:NL=(abs(result2 (i,j)2; %计算结果图像的功率谱S=S+L;endendfor i=1:Mfor j=1:NL1=(abs(g(i,j)2; %计算源图像的功率谱S1=S1+L1;endend
24、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:mfor v=1:nD=sqrt(u-M)2+(v-N)2);H=1-exp(-
25、D2/(2*r2);F_result(u,v)=F_result(u,v)*H;endendG_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:Mfor j=1:NL=(abs(F_result (i,j)2; %计算结果图像的功率谱S=S+L;endendfor i=1:Mfor j=1:NL1=(abs(g
26、(i,j)2; %计算源图像的功率谱S1=S1+L1; endendL=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
27、:Mfor j=1:NH=1+4*pi2*(i-m)2+(j-n)2); % 计算高通滤波器传递函数result1(i,j)=H*g(i,j);endendresult1=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.bm
28、p,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:Mfor j=1:Nd=sqrt(i-m)2+(j-n)2); if (d=0)h2=0;elseh2=1/(1+0.414*(d0/d)(2*nn);% 计算传递函数endresult2(i,j)=(1+h2)*g(i,j);endendresult2=ifftshift(result2);J4=ifft2(result2);J5=uint8(real(J4);subplot(1,2,2);imshow(J5);title(反锐化掩膜滤波图 test4 copy.bmp); % 滤波后图像显示
