频率域滤波的MATLAB设计与实现.docx
《频率域滤波的MATLAB设计与实现.docx》由会员分享，可在线阅读，更多相关《频率域滤波的MATLAB设计与实现.docx（27页珍藏版）》请在冰豆网上搜索。
频率域滤波的MATLAB设计与实现
综合课程设计
设计题目频率域滤波的MATLAB设计与实现
专业名称
班级学号
学生姓名
指导教师
设计时间
摘要
图像处理主要应用于对图像视觉效果的改善，如去噪处理、图像增强、几何校正等。
而本次主要是进行频率域增强技术，它是增强技术的重要组成部分。
本次设计重点用MATLAB对一幅图像作Fourier变换，然后对变换后得到的频谱图像进行分别进行理想低通滤波、二阶Butterworth低通滤波和高斯低通滤波，将原图和变换图放在一起进行比较，观察其图像，对其结果进行分
关键词：
Matlab，Fourier变换，Butterworth，高斯（Gaussian），低通滤波
1.数字图像处理
1.1发展概况：
数字图像处理最早出现于20世纪50年代，当时的电子计算机已经发展到一定水平，人们开始利用计算机来处理图形和图像信息。
数字图像处理作为一门学科大约形成于20世纪60年代初期。
早期的图像处理的目的是改善图像的质量，它以人为对象，以改善人的视觉效果为目的。
图像处理中，输入的是质量低的图像，输出的是改善质量后的图像，常用的图像处理方法有图像增强、复原、编码、压缩等。
首次获得实际成功应用的是美国喷气推进实验室（JPL）。
他们对航天探测器徘徊者7号在1964年发回的几千张月球照片使用了图像处理技术，如几何校正、灰度变换、去除噪声等方法进行处理，并考虑了太阳位置和月球环境的影响，由计算机成功地绘制出月球表面地图，获得了巨大的成功。
随后又对探测飞船发回的近十万张照片进行更为复杂的图像处理，以致获得了月球的地形图、彩色图及全景镶嵌图，获得了非凡的成果，为人类登月创举奠定了坚实的基础，也推动了数字图像处理这门学科的诞生。
1.2关键技术：
1）图像变换：
由于图像阵列很大，直接在空间域中进行处理，涉及计算量很大。
因此，往往采用各种图像变换的方法，如傅立叶变换、沃尔什变换、离散余弦变换等间接处理技术，将空间域的处理转换为变换域处理，不仅可减少计算量，而且可获得更有效的处理（如傅立叶变换可在频域中进行数字滤波处理）。
目前新兴研究的小波变换在时域和频域中都具有良好的局部化特性，它在图像处理中也有着广泛而有效的应用。
2）图像编码压缩：
图像编码压缩技术可减少描述图像的数据量（即比特数），以便节省图像传输、处理时间和减少所占用的存储器容量。
压缩可以在不失真的前提下获得，也可以在允许的失真条件下进行。
编码是压缩技术中最重要的方法，它在图像处理技术中是发展最早且比较成熟的技术。
3）图像增强和复原：
按照特定的需要突出一幅图像中的某些信息或强化某些感兴趣的特征，将原来不清晰的图片变得清晰，使之改善图像质量和丰富信息量，提高图像的视觉效果和图像成分的清晰度，加强图像判读和识别效果的图像处理的方法。
图像增强不考虑图像降质的原因，突出图像中所感兴趣的部分：
如强化图像高频分量，可使图像中物体轮廓清晰，细节明显；如强化低频分量可减少图像中噪声影响[1]。
4）图像分割：
是将图像分为若干个特定的、具有独特性质的区域，其中每一个区域都是像素的一个连续集合。
它是图像处理到图像分析的关键步骤。
常用的分割方法主要分一下几类：
基于阀值的分割方法、基于区域的分割方法、基于边缘的分割方法。
虽然近年来提出了很多新的分割方法，但并没有一种适用于所有图像的分割方法。
在实际应用中，通常将多种分割算法有效的结合在一起使用以获得更好的图像分割效果[2]。
5）图像描述：
将图像分割为区域后，接下来通常要将分割区域加以表示与描述，以方便计算机处理。
图像描述也是图像识别的必要前提。
作为最简单的二值图像可采用其几何特性描述物体的特性，一般图像的描述方法采用二维形状描述，它有边界描述和区域描述两类方法[1]。
对于特殊的纹理图像可采用二维纹理特征描述。
随着图像处理研究的深入发展，已经开始进行三维物体描述的研究，提出了体积描述、表面描述，广义圆柱体描述等方法[4]。
6）灰度变换：
其运算简单，可以达到增强对比度去除噪声的效果。
但是该方法对于原图中所有灰度级的变换是相同的。
然而在实际情况中，人们更关心图像中的目标，对背景不太关注，希望对目标的灰度级进行增强，对背景的灰度级进行压缩。
这样，不但能提高目标的对比度，还可以更清楚的显示目标内部的细节变化，并且忽略了人们不关心的背景的部分细节。
即使原灰度级的范围较大，该方法也可以得到满意的效果[3]。
7）空间滤波：
一种采用滤波处理的影响增强方法。
其理论基础是空间卷积。
目的是改善影像质量，包括出去高频噪声与干扰，及影像边缘增强、线性增强以及去模糊等。
2．频率域滤波的产生背景及意义
2.1傅立叶级数和变换简史：
法国数学家傅立叶指出任何周期函数都可以表示为不同频率的正弦和或余弦之和的形式，每个正弦项和/或余弦乘以不同的系数（现在称该和为傅立叶级数）。
无论函数多么复杂，只要它是周期的，并且满足某些适度的数学条件，都可以用这样的和来表示。
我们现在认为这是理所当然的，但在当时，这个概念第一次出现之后，一个复杂函数可以表示为简单的正弦和余弦之和的概念一点也不直观，所以傅立叶思想遭到怀疑是不足为奇的。
甚至非周期函数也可用正弦和/或余弦乘以加权函数的积分来表示。
在这种情况下的公式就是傅立叶变换，其作用在多数理论和应用学科中甚至远大于傅立叶级数。
用傅立叶级数或变换表示的函数特征完全可以通过傅立叶反变换来重建，而不会丢失任何信息。
这是这种表示方法的最重要特征之一，因为它可以使我们工作于“傅立叶域”，而且在返回到函数的原始域时不会丢失任何信息。
总之，傅立叶级数和变换是解决实际问题的工具，它作为基础工具被广泛的学习和使用。
傅立叶概念的最初应用是在热扩散领域，在该领域，人们考虑用微分方程来表示热流动，并且使用这种方法第一次获得了结论。
在过去一个世纪，特别是后50年，傅立叶的思想使整个工业和学术界都空前繁荣。
早在20世纪60年代，数字计算的出现和快速傅立叶变换算法的“发现”在信号处理领域产生了巨大变革。
这两种核心技术第一次允许人们对医学监视器和扫描仪到现代电子通信的异常重要的信号进行实际处理。
2.2频率域滤波的意义：
滤波法对于模糊图像的复原，边缘的强化和噪声的去除方面都有明显的效果。
熟练的应用不同的滤波图像处理方法对试听资料的检验有很大的帮助，达到改善图像质量的目的。
图像的平滑主要应用在去除图像的噪声上；图像的锐化则可以增强图像的细节和边缘；而图像的复原则可对运动模糊图像，离焦模糊图像等进行复原。
近年来，计算机技术的迅速发展，为数字图像处理提供了强大的软件和硬件支持，也促进了数字图像处理理论和方法的不断深入，使其迅速成为一门新兴技术，并被广泛应用于各个领域。
滤波法新理论的不断提出，新算法的不断应用都使这门技术不断前进成为图像处理技术的重要一员。
3.频率域滤波的常用方法
频率域处理法[4]是在图像的某种变换域内，对变换域的系数进行运算，然后在反变换到原来的空域得到增强的图像，这是一种间接处理方法。
比如，先对图像进行傅里叶变换，再对图像的频域进行滤波处理，最后将频域处理后的图像变换值反变换到空间域，从而得到增强后的图像[5]。
频率域数字图像滤波是通过不同的滤波器在频率域对图像进行处理的方法[6]。
它可分为频率域平滑，频率域锐化和同态滤波增强。
对于一幅图像来说，高频部分大致对应图像中的边缘细节，低频部分大致对应着图像中过渡比较平缓的部分。
每一类滤波法根据具体的算法不同都可以分为很多种，如中值滤波，现在就推广出许多种中值滤波器，如加权中值滤波器（weightedmedianfilter）,多级中值滤波器（multistagemedianfilerner）,组合滤波器（hybridfiler或称L1滤波器），堆滤波器（stackfilter）和置换滤波器（permutationfilter）等。
3.1低通滤波
◆理想低通滤波
在以原点为圆心、以D0为半径的圆内，无衰减地通过所有频率，而在该圆外“切断”所有频率的二维低通滤泼器，称为理想低通滤波器（ILPF）；它由一下函数决定：
其中D0是一个正常数，D（u，v）是频率域中点（u，v）与频率矩形中心的距离：
D（u，v）=[（u-P/2）2+（v-Q/2）2]1/2
3.1.1理想低通滤波器的截面图
如下图3.1所示；
图3.1理想低通滤波器截面图
理想低通滤波器具有一个特性——振铃效果，效果图如下3.2所示
图3.2理想低通滤波器的振铃效果
◆巴特沃斯低通滤波
物理上可实现（理想低通滤波器在数学上定义得很清楚，在计算机模拟中也可实现，但在截断频率处直上直下的理想低通滤波器是不能用实际的电子器件实现的）减少振铃效应，高低频率间的过渡比较光滑，
n阶Butterworth低通滤波器的传递函数为：
其中，D0为截止频率。
=1时，
=0.5，它的特性是传递函数比较平滑，连续衰减，而不像理想滤波器那样陡峭变化，即明显的不连续。
因此采用该滤波器滤波在抑制噪声的同时，图像边缘的模糊程度大大减小，没有振铃效应产生，滤波效果如图3.3所示。
图3.3Butterworth低通滤波效果
巴特沃斯滤波器中阶数对振铃现象的影响：
阶数越高，越明显，如下图3.4所示：
图3.4巴特沃斯滤波器阶数对振铃现象的影响
◆高斯低通滤波
高斯（Gaussian）低通滤波器的传递函数为：
其中，为标准偏差。
令＝D0,我们可以根据截止参数D0得到表达式：
当D（u,v）=D0时，滤波器H（u,v）由最大值1下降为0.607。
GLPF没有振铃现象，但与阶数为2的BLPF相比，其通带要宽些，这样对应的空间滤波器的灰度级轮廓更窄些，因而平滑效果要差些。
对于巴特沃斯低通滤波器和高斯低通滤波器，振铃现象从严重到无，但平滑效果从好到差，BLPF可以看成ILPF和GLPF的过渡，阶为1时与GLPF差不多，阶数越高越接近BPLG.
如下图3.5表示出了高斯低通滤波器对于不同D0值的滤波效果；
图3.5高斯（Gaussian）低通滤波器对于不同的D0值的滤波效果
3.2高通滤波
◆理想高通滤波
一个理想高通滤波器（IHPF）定义为：
其中，D0是截止频率，D（u，v）由下面公式给出：
D（u，v）=[（u-P/2）2+（v-Q/2）2]1/2
如同ILPF一样，IHPF在物理上也是无法实现的，但是IHPF可以用于解释空间域的振铃等现象。
下图3.6即为理想高通滤波器的滤波效果：
图3.6理想高通滤波器的滤波效果
◆巴特沃斯高通滤波
巴特沃斯n阶截止频率为DO的巴特沃斯高通滤波器的传递函数为：
其中：
D（u，v）=[（u-P/2）2+（v-Q/2）2]1/2
由频域滤波模型Q（U，v）=F（U，v）H（U，v）知，F（U，v）中的低频（小于D0）成分，因乘上一个远小于1的H（U，v）值而被衰减。
而高频成分却被乘以一个接近于1的H（U，v）值而保留，这即是所谓的高通滤波的原理。
当截止频率D0越大，滤掉的低频成分越多，同样损失的高频成份也越多。
如下图3.7即为巴特沃斯高通滤波器的滤波效果：
图3.7n阶巴特沃斯高通滤波器滤波效果
◆高斯高通滤波
其中，为标准偏差。
通过令＝D0,我们可以根据截止参数D0得到的值。
下图3.8即为高斯高通滤波器的滤波效果：
图3.8高斯高通滤波器滤波效果
3.3带阻滤波
带阻滤波器阻止一定频率范围内的信号通过而允许其它频率范围内的信号通过。
◆理想带阻滤波器
理想带阻滤波器的传递函数：
这里，W是频带的宽度，D0是频带的中心半径。
◆巴特沃斯带阻滤波器：
n阶的巴特沃思带阻滤波器的表达式为：
◆高斯带阻滤波器
高斯带阻滤波器的表达式为：
下图3.9是理想带阻滤波器、阶数为1的巴特沃斯带阻滤波器和高斯带阻滤波器的透视图：
图3.9理想滤波器、巴特沃思滤波器（阶数为1）和高斯带阻滤波器的透视图
3.4带通滤波
带通滤波器执行与带阻滤波器相反的操作，带通滤波器的传递函数据相应的带阻滤波器的传递函数并应用下式得到的：
4．原理及实现
4.1频率域增强基本理论
不对Fourier变换（FT）和图像的频率域处理技术有所了解，就不可能完全理解图像增强这个最基本的图像处理任务。
频域增强指在图像的频率域内，对图像的变换系数（频率成分）直接进行运算，然后通过Fourier逆变换以获得图像的增强效果。
一般来说，图像的边缘和噪声对应Fourier变换中的高频部分，所以低通滤波能够平滑图像、去除噪声。
图像灰度发生聚变的部分与频谱的高频分量对应，所以采用高频滤波器衰减或抑制低频分量，能够对图像进行锐化处理。
卷积理论是频域技术的基础，设函数f（x,y）与算子h（x,y）的卷积结果是g（x, y），即g（x,y）=h（x,y）*f（x,y），那么根据卷积定理在频域有：
其中G（u,v），H（u,v），F（u,v）分别是g（x,y），h（x,y），f（x,y）的傅立叶（或其它）变换，H（u,v）是转移函数。
在具体增强应用中，f（x,y）是给定的（所以F（u,v）可利用变换得到），需要确定的是H（u,v），这样具有所需特性的g（x,y）就可算出G（u,v）而得到：
4.2傅立叶变换
傅里叶变换是将时域信号分解为不同频率的正弦信号或余弦函数叠加之和。
傅立叶变换是数字图像处理技术的基础，其通过在时空域和频率域来回切换图像，对图像的信息特征进行提取和分析，简化了计算工作量，被喻为描述图像信息的第二种语言，广泛应用于图像变换，图像编码与压缩，图像分割，图像重建等。
因此，对涉及数字图像处理的工作者，深入研究和掌握傅立叶变换及其扩展形式的特性，是很有价值得。
把傅立叶变换的理论通其物理解释相结合，将有助于解决大多数图像处理问题。
傅里叶变换可分为连续傅里叶变换、离散傅里叶变换、快速傅里叶变换。
快速傅里叶变换（FFT）是计算离散傅里叶变换（DFT）的快速算法。
离散傅里叶变换运算量巨大，计算时间长，即运算时间很长。
而快速傅里叶变换的提出将傅里叶变换的复杂度由降到了，很大程度上减少了计算量。
令
，
，u=0，1，2，…，M-1
则
，
4.3频率域理想低通（ILPF）滤波器
一个二维的理想低通滤波器（ILPF）的转换函数满足（是一个分段函数）
其中：
D0为截止频率D（u,v）为距离函数D（u,v）=（u2+v2）1/2
4.3.1理想低通滤波器的截面图
理想低通滤波器截面图如图4.1所示：
图4.1理想低通滤波器截面图
振铃效果——理想低通滤波器的特性
理想低通滤波后图
原图
图4.2理想低通滤波器振铃效果
4.3频率域巴特沃兹（Butterworth）低通滤波器
物理上可实现（理想低通滤波器在数学上定义得很清楚，在计算机模拟中也可实现，但在截断频率处直上直下的理想低通滤波器是不能用实际的电子器件实现的）减少振铃效应，高低频率间的过渡比较光滑n阶Butterworth低通滤波器的传递函数为：
式中，D0为截止频率。
=1时，
=0.5，它的特性是传递函数比较平滑，连续衰减，而不像理想滤波器那样陡峭变化，即明显的不连续。
因此采用该滤波器滤波在抑制噪声的同时，图像边缘的模糊程度大大减小，没有振铃效应产生，如图2.1所示。
图4.3Butterworth低通滤波效果
阶数对振铃现象的影响：
阶数越高，越明显，如图2.2所示
图4.4阶数对振铃现象的影响
4.4频率域高斯（Gaussian）低通滤波器
高斯（Gaussian）低通滤波器的传递函数为：
其中，
为标准偏差。
通过令
＝D0,我们可以根据截止参数D0得到表达式：
当D（u,v）=D0时，滤波器H（u,v）由最大值1下降为0.607。
GLPF没有振铃现象，但与阶数为2的BLPF相比，其通带要宽些，这样对应的空间滤波器的灰度级轮廓更窄些，因而平滑效果要差些。
以上两种滤波器，振铃现象从严重到无，但平滑效果从好到差，BLPF可以看成ILPF和GLPF的过渡，阶为1时与GLPF差不多，阶越高越接近BPLG.。
图4.5高斯（Gaussian）低通滤波器对于不同的D0值的滤波效果
5．程序设计
5.1算法设计（程序设计流程图）
算法流程图如图3.1所示：
图3.1程序设计流程图
5.2对灰度图像进行Fourier变换的程序
I=imread（'Testpicture.jpg'）;
f=double（I）;%数据类型转换，MATLAB不支持图像的无符号整型的计算
g=fft2（f）;
figure,imshow（log（1+abs（fftshift（g）））,[]）%将频谱显示为一幅图像
title（'原图Fourier变换后频谱图'）;
5.3频率域理想低通滤波器
I=imread（'Test1.jpg'）;
I=rgb2gray（I）;
Do=input（'请输入Do值:
'）;
subplot（221）,imshow（I）;
title（'原图像'）;
s=fftshift（fft2（I））;
subplot（223）,
imshow（abs（s）,[]）;
title（'图像傅里叶变换所得频谱'）;
subplot（224）,
imshow（log（abs（s））,[]）;
title（'图像傅里叶变换取对数所得频谱'）;
[a,b]=size（s）;
a0=round（a/2）;
b0=round（b/2）;
fori=1:
a
forj=1:
b
distance=sqrt（（i-a0）^2+（j-b0）^2）;
ifdistance<=Doh=1;
elseh=0;
end;
s（i,j）=h*s（i,j）;
end;
end;
s=uint8（real（ifft2（ifftshift（s））））;
subplot（222）,
imshow（s）;
title（sprintf（'%s%d','理想低通滤波结果Do=',Do））
5.4二阶巴特沃斯（Butterworth）低通滤波程序
此函数的功能为：
使用巴特沃斯低通滤波原理实现图像柔和平滑，为便比较，通过设置不同的d0值运行程序,如果对各个d0一起运行，运行速度奇慢，所以对各个d0单独运行。
I=imread（'a.jpg'）;
imshow（I）;title（'原图'）;
s=fftshift（fft2（I））;
[M,N]=size（s）;%分别返回s的行数到M中，列数到N中
n=2;%对n赋初值
%BLPF滤波，d0=15，30，80
d0=input（'请输入d0值:
'）;%初始化d0,d0依次取15，30，80
n1=floor（M/2）;%对M/2进行取整
n2=floor（N/2）;%对N/2进行取整
fori=1:
M
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）;%点（i,j）到傅立叶变换中心的距离
h=1/（1+（d/d0）^（2*n））;%BLPF滤波函数
s（i,j）=h*s（i,j）;%BLPF滤波后的频域表示
end
end
s=ifftshift（s）;%对s进行反FFT移动
%对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数
s=uint8（real（ifft2（s）））;
figure;%创建图形图像对象
imshow（s）;%显示BLPF滤波处理后的图像ttitle（sprintf（'%s%d%s','二阶Butterworth低通滤波结果（D0=',d0,'）'））
5.5高斯（Gaussian）低通滤波程序
此函数的功能为：
使用高斯低通滤波原理实现图像越柔和平滑，为便比较，通过设置不同的d0值运行程序,如果对各个d0一起运行,运行速度奇慢，所以对各个d0单独运行。
%高斯低通滤波程序
ima=imread（'a.jpg'）;
ima=rgb2gray（ima）;
imshow（ima）;
title（'原图'）;
ima0=im2double（ima）;
Do=input（'请输入d0值:
'）;%初始化d0,d0依次取15，30，80
[Row,Col]=size（ima）;
ima=im2double（ima）;
forix=1:
Row%滤波中心化
foriy=1:
Col
ima（ix,iy）=（-1）^（ix+iy）*ima（ix,iy）;
end
end
FourIma=fft2（ima）;%中心变换后的傅立叶变换
[FRow,FCol]=size（FourIma）;
foru=1:
FRow%Do的高斯低通滤波器
forv=1:
FCol
D（u,v）=[（u-Row/2）^2+（v-Col/2）^2]^0.5;
factc=-D（u,v）^2/（2*Do^2）;
H（u,v）=exp（factc）;
G（u,v）=H（u,v）*FourIma（u,v）;
end
end
IFourIma=ifft2（G）;%高斯低通滤波后的反傅立叶变换
forix=1:
FRow%反傅立叶变换后*（-1）^（x+y）并取实部
foriy=1:
FCol
LowPass（ix,iy）=（-1）^（ix+iy）*IFourIma（ix,iy）;
end
end
LowPass=real（LowPass）;
figure;
imshow（LowPass）;%Do=15高斯低通滤波结果
title（sprintf（'%s%d%s','Do=',Do,'高斯低通滤波结果'））
6．结果与分析
我们还可以通过函数colorbar给一个坐标轴对象添加加一条色带，该色带可以指示出该图像中不同颜色的数据值，以显示出图像的灰度级。
输入以下代码：
f=imread（'Test1.jpg'）;
imshow（f）;
colorbar;
可得到添加色带后的图像如图6.1所示，从图6.3可知该图像是数据类型为uint8的灰度图像，其灰度级范围从0~255，级灰度级为256。
图6.1添加色带后的灰度图
6.1对灰度图像进行Fourier变换后的频谱图
在MATLAB命令输入窗中，输入Fourier变换程序，按回车即生成Fourier变换后频谱图如图6.2所示：
图6.2原图Fourier变换后频谱图
6.2二阶巴特沃斯（Butterworth）低通滤波结果与分析
在MATLAB命令输入窗中，输入巴特沃斯（Butterworth）低通滤波程序，按回车即生成原图如图6.3，Do=15的二阶Butterworth低通滤波结果如图6.4所示，Do=30的二阶Butterworth低通滤波结果如图6.5所示，Do=80的二阶Butterworth低通滤波结果如图6.6所示。
图6.3原图
图6.4Do=15的二阶Butterworth低通滤波结果
图6.5Do=30的二阶Butterworth低通滤波结果
图6.6Do=80的二阶Butterwort
升级会员 