《MATLAB应用图像处理》课件第4章

4.1线性滤波器的空间域设计法4.2滤波器频域设计法

4.3图像滤波4.4典型滤波器介绍与设计习题

第4章MATLAB图像滤波

4.1线性滤波器的空间域设计法

4.1.1简介滤波是一种用来改变或增强图像的技术。例如，通过滤波可以强化或去除图像的某种特征，它能够处理图像的操作包括平滑、锐化、边缘增强等。滤波是一种邻域操作，在这种操作中，输出图像的像素点的值取决于输入图像某个邻域内的像素点的值。线性滤波中输出像素点的值是输入图像某个邻域内像素点的值的线性组合。4.1.2卷积线性滤波通过卷积来完成。卷积是一种邻域操作，操作中输出图像的像素点的值是对应输入像素点邻域内像素点值的加权和，其中权矩阵称为卷积核矩阵，也称为滤波器。卷积核是相关核旋转180°得到的。例如，图4.1所示图像A和卷积核h分别为(4.1)以输入图像像素点(2，4)为例说明卷积的计算过程。第一步：将卷积核绕其中心元素5旋转180°；第二步：将旋转后的卷积核平移，使其中心元素与像素点(2，4)重合；第三步：将卷积核元素与输入图像的对应像素相乘；第四步：将以上各项相乘的结果相加，得到卷积结果575。图4.1卷积计算方法4.1.3相关相关操作与卷积操作类似，也是求邻域内点的加权和，权矩阵为相关核。不同的是在计算过程中不需对相关核进行旋转操作。图4.2给出了图像A的像素点(2，4)处相关的计算示意图。图4.2相关的计算方法以图像像素点(2，4)为例说明相关的计算过程。第一步：将相关核平移，使其中心元素与像素点(2，4)重合；第二步：将相关核元素与图像的对应像素值相乘；第三步：将以上各项相乘的结果相加，得到卷积结果585。以上介绍表明，不管是在空间域内还是在频域内定义滤波器，关键在于给出滤波器的核矩阵，对于空间域中核矩阵的定义比较容易获得，我们将在第7章图像空间增强部分具体介绍。而滤波器频域内的描述我们接触得比较多，滤波器种类也比较丰富，设计并获取核矩阵比空间滤波器的复杂，MATLAB图像处理工具箱对频域滤波器给出了几个函数。本章将对此进行较为详细的介绍。另外，MATLAB图像处理工具箱还给出了几个对二维甚至高维图像数据进行滤波处理的函数，如imfilter等。4.2滤波器频域设计法

图像处理工具箱支持有限冲激响应(FIR)线性滤波器，即对单点或者单个脉冲具有有限范围的响应的滤波器。FIR滤波器具有如下几个特性，这些特性使得它便于在MATLAB图像处理环境中实现：

(1) FIR滤波器的系数便于用矩阵表示；

(2)二维FIR滤波器是一维FIR滤波器的自然拓展；

(3) FIR滤波器有多种可靠的众所周知的设计方法；

(4) FIR滤波器易于实现；

(5) FIR滤波器可以设计成具有线性相位以消除图像的扭曲。无限冲激响应(IIR)滤波器由于不稳定且相对于FIR滤波器来讲不便于设计，故而在MATLAB图像处理工具箱中没有提供对IIR滤波器的支持。MATLAB所有图像处理滤波器设计函数都返回一个FIR滤波器。4.2.1频率变换设计法频率变换法通过频率变换，将一个一维FIR滤波器变换成二维FIR滤波器。它设计出的滤波器保持了一维滤波器的大多数特性，特别是传输带宽、纹波特征等。由于设计一个具有特殊性质的一维滤波器比二维滤波器要简单得多，所以频域变换法可以产生很理想的效果。

MATLAB图像处理工具箱中的一个ftrans2函数实现了频域变换设计算法。ftrans的基本调用方法为

h=ftrans2(b,t)

h=ftrans2(b)第一种调用方法由一个一维数字滤波器的冲激响应b通过频域变换t生成二维滤波器h，t为变换矩阵；第二种调用方法调用t的默认值产生具有循环对称的变换矩阵，默认情况下，t采用McClellan变换矩阵：

(4.2)可以通过自己定义变换矩阵，得到不同的对称性。冲激响应b必须是一个一维的、长度为奇数的数字滤波器，可以通过信号处理工具箱中的fir1、fir2和firpm等函数产生。如果t是m×n维的，b的长度为Q，则所设计的滤波器h的尺度为((m-1)×(Q-1)/2+1)×((n-1)×(Q-1)/2+1)。函数ftrans2返回的二维滤波器的频率响应由以下变换定义：

(4.3)其中，B(ω)是一维滤波器b的傅立叶变换，即

(4.4)T(ω1，ω2)是变换矩阵t的傅立叶变换，即

(4.5)返回的滤波器h是H(ω1，ω2)的逆傅立叶变换，即

(4.6)例如，以下代码设计了一个一维最优的波动FIR滤波器，并用频域变换方法设计一个具有同样特性的二维滤波器(结果如图4.3所示)：

b=firpm(10,[00.40.61],[1100]); %定义设计一维FIR滤波器

h=ftrans2(b); %频域变换法设计二维滤波器[H,w]=freqz(b,1,64,'whole'); %计算一维滤波器响应

colormap(jet(64))

plot(w/pi-1,fftshift(abs(H)))

figure,freqz2(h,[3232])

%计算二维滤波器响应并绘制图像图4.3波动低通滤波器频率变换设计法4.2.2频域采样法频域采样法基于一目标频率响应创建一个滤波器。给定一个频率响应矩阵，频域采样法创建一具有通过这些频响点的特性的滤波器，且对于给定点之间的行为没有限制。频域采样法所设计的滤波器的频响通过给定的点。

MATLAB图像处理工具箱中的fsamp2函数执行频域采样方法设计二维FIR滤波器。fsamp2返回的滤波器h通过给定的输入参数Hd。例如，创建一个11×11的滤波器，并且绘制所设计滤波器的频率响应，如图4.4所示。代码为Hd=zeros(11,11);Hd(4:8,4:8)=1;[f1,f2]=freqspace(11,'meshgrid');mesh(f1,f2,Hd),axis([-11-1101.2]),colormap(jet(64))h=fsamp2(Hd);figure,freqz2(h,[3232]),axis([-11-1101.2])图4.4低通滤波器频率采样法注意：对应于给定的频域响应采样点，所设计的滤波器的频响存在波动。这些波动是频域采样法设计滤波器不可避免的，只要在频率响应中存在尖的变化点或者边缘都会存在这个问题。可以通过一个大的滤波器来减弱这种波动，然而，大的滤波器并不减低波动的峰值，并且需要消耗更多的计算时间。为了获得一个对所要求频响更为平滑的滤波器，可以考虑使用频率变换法或者窗口法来设计滤波器。4.2.3窗口法二维窗函数设计法是二维滤波器设计的最早技术之一，它是一维窗函数的直接推广。二维窗函数设计法是用一个窗口函数与理想冲激响应相乘来产生响应的滤波器，这个方法逐渐减弱了理想脉冲响应。与频率采样法类似，窗口法产生一个频率响应接近于一个给定的频率响应的滤波器，但是相对于频率采样法，窗口法能够产生更好的滤波器。图像处理工具箱提供了两个窗口法滤波器设计函数，分别为fwind1和fwind2。其中fwind2函数用一个指定的二维窗口函数来设计一个二维的滤波器，是窗口法设计滤波器的基础，调用fwind1设计滤波器时也调用此函数。fwind1基于二维的窗口设计一个二维的滤波器，二维的窗口从用户指定的一个或者两个一维的窗口创建而来。常用的窗函数有矩形(Boxcar)窗、三角(Bartlett)窗、海明(Hamming)窗、汉宁(Hanning)窗、布莱克曼(Blackman)窗和凯撒(Kaiser)窗等。

fwind2的基本调用方法为

h=fwind2(Hd,win)

fwind2函数对期待频响Hd进行逆傅立叶变换，然后与窗函数win相乘得到所设计的滤波器h：

(4.7)(4.8)

fwind2函数返回值h是核矩阵，可以直接用于imfilter和filter2函数中，h的尺寸等于窗口win的尺寸，为了得到精确的结果，可以使用freqspace函数返回的频率点来生成Hd。

fwind2函数返回值h是核矩阵，可以直接用于imfilter和filter2函数中，h的尺寸等于窗口win的尺寸，为了得到精确的结果，可以使用freqspace函数返回的频率点来生成Hd。对上式进行傅立叶变换：

(4.9)此式表明，窗口函数相当于对指定的频响进行了平滑处理。此式表明，窗口函数相当于对指定的频响进行了平滑处理。

fwind1函数的基本调用方法如下：

h=fwin1(Hd,win)

h=fwin1(Hd,win1,win2)第一种调用方法采用黄氏算法，从一维窗函数win得到近似回转对称的二维窗函数。设一维窗函数为win(t)，则应用黄氏算法得到的二维窗函数win(n1，n2)为(4.10)第二种调用方法采用两个一维窗函数直接相乘的方法得到二维窗函数：win(n1，n2) = win(n1)win(n2)(4.11)两种调用方法在得到二维窗函数后调用fwind2生成待设计的二维滤波器。例如，用fwind1函数从指定的频率响应Hd创建一个11×11的滤波器。调用hamming函数创建一个一维的窗函数，fwind1函数将其扩展成一个二维窗函数，代码如下：Hd=zeros(11,11);Hd(4:8,4:8)=1;[f1,f2]=freqspace(11,'meshgrid');mesh(f1,f2,Hd),axis([-11-1101.2]),colormap(jet(64))h=fwind1(Hd,hamming(11));figure,freqz2(h,[3232]),axis([-11-1101.2])运行结果如图4.5所示。图4.5窗函数法设计滤波器结果以上我们讨论了二维滤波器的三种不同的设计方法，这三种方法返回的都是可以直接用于计算的核矩阵。fsamp2、fwind1和fwind2函数都是基于待设计滤波器的频响幅值矩阵设计滤波器。频响就是滤波器对不同频率信号响应增益值函数。可以调用freqspace函数创建一个适当的待设计的频率响应矩阵。freqspace函数对任意响应都返回正确、均衡频率分布值。例如，创建一个理想的低通回转对称频响函数，频响的截断值是0.5，可以采用如下代码：

[f1,f2]=freqspace(25,'meshgrid');

Hd=zeros(25,25);d=sqrt(f1.^2+f2.^2)<0.5;

Hd(d)=1;

mesh(f1,f2,Hd)运行结果如图4.6所示。图4.6理想低通滤波器频响图注意：对应于这个频响，由fsamp2、fwind1和fwind2函数创建的滤波器是实数，这个结果适用于大多数的图像处理情况。为了得到类似的实数滤波器，频响函数必须关于频率原点(0，0)对称。可以调用函数freqz2来计算某一个二维滤波器的频响函数。freqz2对频率f1、f2进行归一化处理，使得其值位于[0，1]。其中1对应着采样频率的一半，或者π。如果没有输出参数，则freq2函数绘制该频响函数；调用方法[H,f1,f2]=freqz2(h)返回频响矩阵H、对应的频响点f1和f2。4.3图像滤波以上介绍了滤波器的三种设计方法，这些设计方法得到的都是滤波器的卷积核矩阵。这节介绍如何利用这些设计结果对具体的待滤波图像进行处理。

MATLAB提供了两个采用指定的滤波器对图像进行滤波处理的函数，分别为filter2和imfilter。其中，filter2函数的功能比较单一，而imfilter函数则相对丰富得多，下面我们对imfilter函数进行介绍。

imfilter函数的调用方法有两种：

B=imfilter(A,H)

B=imfilter(A,H,option1,option2,...)第一种调用方法利用滤波器H对多维数组A进行滤波处理。数组A可以是逻辑类型的或非稀疏的任意数据类型的数值数组。返回数组B的大小和类型与输入数组A的一致。在函数的内部都是采用双精度浮点类型进行计算的。如果A为整型或逻辑数组，则需要对计算结果进行阶段或者四舍五入等处理。相比于第一种调用方法，第二种调用方法多了调用选项，这些选项包括指定边界类型、输出大小控制以及计算类型等。具体的选项内容见表4.1～表4.3。表4.1边界类型控制选项表4.2输出大小控制选项表4.3计算类型选项

imfilter函数可以对图像施行卷积滤波或者相关滤波，其中卷积滤波为缺省模式。所需要的卷积核函数可以由频域滤波设计函数产生或者由空间域滤波器设计函数fspeical函数产生。例如，调用imfilter函数对图像施行滤波，使用一个5×5的矩阵。矩阵的元素都相等，这样的滤波器通常称为平均滤波。代码如下：I=imread('coins.png');h=ones(5,5)/25;%区域平均滤波I2=imfilter(I,h);imshow(I)figure,imshow(I2)运行结果如图4.7所示。图4.7利用imfilter函数对图像进行区域平均滤波处理函数imfilter能够处理的数据类型与图像的代数运算能够处理的图像数据类型一致。imfilter函数的输出类型与输入类型一致。imfilter函数滤波计算过程中使用双精度浮点类型，计算后再转化成输入数据的类型。在转化过程中如果计算数据超出输入数据类型所能够允许的范围，则对其进行截断处理。除了imfiter滤波器外，MATLAB还有别的二维和多维滤波器。如函数filter2施行二维相关滤波，conv2施行二维卷积滤波。convn施行多维卷积滤波。这些滤波函数总是将输入数据转化成double类型，输出类型也是double类型，并且总是采用零延拓方法对输入数据进行延拓。4.4典型滤波器介绍与设计

4.4.1低通滤波器边缘和其他尖锐变化(如噪声)在图像的傅立叶变换中主要处于高频部分。因此可以通过衰减指定图像傅立叶变换中高频成分的范围来实现一种图像平滑的效果，用于衰减高频成分的滤波器即为低通滤波器。本节介绍四种低通滤波器：理想低通滤波器、梯形低通滤波器、Butterworth滤波器和指数低通滤波器。这四种滤波器涵盖了从非常尖锐(理想)到非常平坦(高斯)范围的滤波器函数。其中Butterworth滤波器有一个称为滤波器阶数的参数，当此参数较高时，Butterworth滤波器接近理想滤波器，从而可以实现从理想滤波器到平坦的高斯滤波器之间的过渡。

1.理想低通滤波器(ILPF)理想低通滤波器的传递函数为≤＞

(4.12)图4.8为理想低通滤波器的频响曲线和曲面。图4.8理想低通滤波器频率响应图式(4.12)中，，D0为截止频率。此传递函数表明，以径向频率D0为半径的圆区域内所有频率分量无损地通过，区域外的所有分量都被滤除。截断频率D0决定着通过滤波器的能量比，一般来讲，D0的选择要使图像中我们感兴趣的细节能量大部分通过，且截断大部分不感兴趣的部分。随着截断频率的增大，通过滤波器后图像的变化越小。

2.梯形低通滤波器梯形低通滤波器是改良的理想低通滤波器，截断频率由原来的陡降变为斜降，其传递函数为其中，D0、D1预先给定，D0称为截止频率。D(u，v)与理想低通滤波器中定义的相同。图4.9为梯形低通滤波器频率响应图。(4.13)图4.9梯形低通滤波器频率响应图

3. Butterworth滤波器(BLPF)

N阶Butterworth低通滤波器的传递函数为

(4.14)其中，D(u，v)与理想低通滤波器中定义的相同，D0为截止频率。不同于理想低通滤波器，Butterworth低通滤波器在通带与被滤除的频率之间没有明显的截断，是连续光滑的。随着阶数n的增大，Butterworth滤波器边缘越来越陡峭，越接近理想低频滤波器；随着阶数n的减小，Butterworth滤波器则接近于指数低通滤波器。图4.10为n=1和n=3时Buttterworth低通滤波器频率响应图。图4.10Butterworth低通滤波器频率响应图

4.指数低通滤波器指数低通滤波器也称为高斯低通滤波器，可用来寻找空间域和频率域之间的重要联系，其传递函数为

(4.15)其中D(u，v)与理想低通滤波器中定义的相同，D0为截止频率。以上介绍了四种低通滤波器，应用中应该依据具体问题选择所需的滤波器。表4.4给出了各类滤波器的性能特性。表4.4不同低通滤波器性能比较4.4.2高通滤波器上一节介绍了利用低通滤波器衰减图像傅立叶变换中的高频成分的方法使图像平滑，本节将对前面的低通滤波器进行反操作，以得到图像中的高频成分，这个过程称为图像的锐化。图像锐化的目的是突出图像中凸显的边缘信息，加强图像的轮廓特征，以便于人眼和机器识别。

1.理想高通滤波器理想高通滤波器的传递函数为

(4.16)式中，D(u，v) = ，D0为截止频率。图4.11为理想高通滤波器的频响曲线和曲面。图4.11理想高通滤波器频率响应图理想高通滤波器将以截断频率D0为半径的圆周内所有频率置零，圆周外的任何频率毫无衰减地通过。经过理想高通滤波器得到了原图像的边缘信息，截断频率越大时，其得到的信息越少。与理想低通滤波器相类似，理想高通滤波器会产生振铃现象。

2.梯形高通滤波器梯形高通滤波器是改良的理想高通滤波器，截断频率由原来的陡降变为斜降，其传递函数为

(4.17)其中，D0、D1预先给定，D0成为截止频率。D(u，v)与理想低通滤波器中定义的相同。图4.12为梯形高通滤波器频率响应图。图4.12梯形高通滤波器频率响应图

3. Butterworth高通滤波器

N阶Butterworth高通滤波器的传递函数为

(4.18)图4.13为n = 3时Butterworth高通滤波器频率响应图图4.13Butterworth高通滤波器频率响应图(n=3)

Butterworth高通滤波器相对于理想高通滤波器更为平滑，在滤除与保留的频率之间没有明显的分界，所以用此滤波器得到的输出图像的振铃效应不明显。

4.指数高通滤波器截断频率为D0的指数滤波器传递函数为

(4.19)指数高通滤波器得到的结果比前面所介绍的高通滤波器更为平滑，即使是对微小物体和细条用指数滤波器也是比较清晰的。4.4.3滤波器的MATLAB实现以上介绍了滤波器的频率域设计方法以及各类低通和高通滤波器，本节利用这些设计方法设计相应的滤波器对图像进行处理。以上介绍的滤波器的频域设计方法得到的是滤波器的空间卷积核。可以将图像进行傅立叶变换后得到的频率