基于matlab的图像处理滤波器设计.docx

基于matlab的图像处理滤波器设计.docx

《基于matlab的图像处理滤波器设计.docx》由会员分享，可在线阅读，更多相关《基于matlab的图像处理滤波器设计.docx（56页珍藏版）》请在冰点文库上搜索。

基于matlab的图像处理滤波器设计

数字信号处理课程设计任务书

2011－2012学年第一学期　第15周－19周

题目

基于matlab的图像处理滤波器设计

内容及要求：

1、设计一个低通FIR滤波器和一个低通IIR滤波器；

2、分析比较上述两种滤波器的优劣；

3、分析上述两个模型的幅频特性、相频特性、相延迟、群延迟。

进度安排：

1、任务分配、查阅资料2天

2、方案论证3天

3、分析、设计、调试程序5天

4、书写、整理实验报告和小结3天

成员组成：

09044106苏青文08陈舒龙（组长）09方雪松

指导时间：

指导地点：

F618

任务下达

2011年12月16日

任务完成

2012年1月6日

考核方式

1.评阅□2.答辩□3.实际操作□4.其它□

指导教师

汪传忠

系（部）主任

王长坤

注：

1、此表一组一表二份，课程设计小组组长一份；任课教师授课时自带一份备查。

2、课程设计结束后与“课程设计小结”、“学生成绩单”一并交院教务存档。

摘要…………………………………………………………3

一、课程设计目的及要求………………………………………4

二、课程设计内容及任务安排…………………………………4

三、设计原理及设计方法………………………………………5

3.1FIR数字滤波器………………………………………5

3.2IIR数字滤波器………………………………………7

四、与设计相关的知识…………………………………………8

五、设计过程……………………………………………………14

5.1图像加噪处理及功率密度谱对比……………………14

5.2FIR滤波器的设计……………………………………22

（1）布莱克曼窗………………………………………..22

（2）海明窗………………………………………………26

（3）汉宁窗…………………………………………….28

（4）多尔夫-切比雪夫窗........................32

（5）巴特利特窗……………………………………….35

5.3IIR数字滤波器设计..........................38

（1）巴特沃兹………………………………………….38

（2）切比雪夫1................................41

（3）切比雪夫2……………………………………….43

（4）椭圆滤波.................................45

六、FIR和IIR的比较………………………………………….47

七、个人设计总结………………………………………………49

附录………………………………………………………………50

摘要

MATLAB软件具有使用简单、方便，易编程，语言简练，函数库可任意扩充，采用全新数据类型和面向对象编程技术等特点，具有数值分析容易、语法结构简单、矩阵数值运算高效、图形功能完备、图像处理方便等强大的功能，因此已被广泛应用在教学、科研和工程设计的各个领域。

传统的滤波器分析与设计均使用繁琐的公式计算，改变参数后需要重新运算，从而在分析与设计滤波器尤其是高阶滤波器时工作量特别大，利用MATLAB信号处理工具箱可以快速而有效地实现滤波器的分析、设计与仿真。

本课设介绍基于MATLAB环境的IIR数字滤波器和FIR数字滤波器的设计方法，原理和步骤。

IIR数字滤波器的常用设计方法有脉冲响应不变法和双线性变换法等，FIR数字滤波器的设计方法只要有窗函数设计法，频率抽样设计法和最优化设计法等，最后利用设计的滤波器对含噪声的图像进行滤波处理并验证所设计滤波器的性能，同时比较两种设计方法的最优性。

关键字：

数字信号处理MatlabFIR滤波器IIR滤波器

图像加噪图像滤波FIR与IIR滤波器比较

1、课程设计目的及要求

1.1目的

综合运用数字信号处理的理论知识进行频谱分析和滤波器设计FIR和IIR滤波器，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。

1、巩固和加深理解数字信号处理基本理论知识，提高综合运用这些知识进行程序设计的能力。

2、能根据需要查阅相关手册和资料，提高独立分析问题和解决问题的能力。

3、能正确使用仪器设备，掌握测试原理，掌握综合性程序设计的基本过程。

4、掌握时域离散时间信号和时域离散系统基本概念和基本性质。

5、掌握时域离散信号、系统的频域分析工具。

6、掌握基于Matlab的数字信号处理领域的计算机辅助分析与设计的基本方法。

7、培养团队协作精神、创新意识、严肃认真的治学态度和严谨求实的工作作风

1.2基本要求

1、学会MATLAB的使用，掌握MATLAB的程序设计方法。

2、掌握在Windows环境下数字图像处理的方法。

3、掌握数字信号处理的基本概念、基本理论和基本方法。

4、掌握MATLAB设计FIR和IIR数字滤波器的方法。

5、学会用MATLAB对数字图像进行分析和处理

二、课程设计内容及任务安排

每人用自己的一张照片，并对照片进行加噪声处理，通过设计滤波器滤除噪声。

每组三人，其中一人设计FIR滤波器进行噪声处理，一人设计IIR滤波器滤除噪声，另外一人分析FIR和IIR滤波器的优劣比较。

其中在进行滤波器设计中，要求每种模型都要求设计，FIR滤波器设计包括布莱克曼窗、海明窗、汉宁窗、多尔夫—切比雪夫窗、巴特利特窗，IIR滤波器设计包括巴特沃兹、切比雪夫1和2，椭圆模型。

并且分析每个模型的幅频特性，相频特性，相延时，群延时。

最后得出结论。

具体实施方案：

1、苏青文主要负责设计IIR无限脉冲响应数字滤波器。

2、方雪松主要负责设计FIR有限脉冲响应数字滤波器。

3、组长陈舒龙主要负责图像处理、整理报告及FIR和IIR两种数字滤波器的比较。

三、设计原理及设计方法

数字滤波器设计的一个重要步骤就是确定一个可实现的传递函数G（z）来逼近指定的频率响应。

在G（z）确定下来之后，接下来要设计用一种合适的滤波器结构来实现它。

滤波器设计过程的第一步是估计传输函数阶数。

然后根据各种模型用Matlab函数得出传输函数。

3.1FIR数字滤波器

3.1.1基本介绍

FIR是有限冲激响应（FiniteImpulseResponse）的简称。

由线性系统理论可知，在某种适度条件下，输入到线性系统的一个冲击完全可以表征系统。

当我们处理有限的离散数据时，线形系统的响应（包括对冲击的响应）也是有限的。

若线性系统仅是一个空间滤波器，则通过简单地观察它对冲击的响应，我们就可以完全确定该滤波器。

通过这种方式确定的滤波器称为有限冲击响应（FIR）滤波器.

3.1.2设计方法

FIR滤波器设计的最简单方法是，对期望的理想频率响应进行离散时间傅里叶逆变换，得到双无限长冲击响应，对所得的无限系数简单地截尾为有限系数。

滤波器类型确定之后，接下来是估计滤波器阶数N，matlab里提供了阶数估计的函数，我们可以直接用，阶数估计之后就是确定传输函数num和den，每种模型都有这样的函数求出传输函数。

本设计用了5个模型加窗函数：

布莱克曼窗、海明窗、汉宁窗、多尔夫—切比雪夫窗、巴特利特窗，来设计FIR滤波器。

几种窗函数的性能比较

窗函数

主瓣宽度△ML

相对旁瓣级别（dB）

最小阻带衰减（dB）

过渡带带宽△ω

矩形窗

4π/（2M+1）

13.3

20.9

0.92π/M

巴特利特窗

4π/（M+1）

26.5

汉宁窗

8π/（2M+1）

31.5

43.9

3.11π/M

海明窗

8π/（2M+1）

42.7

54.5

3.32π/M

布拉克曼窗

12π/（2M+1）

58.1

75.3

5.56π/M

3.2IIR数字滤波器

3.2.1IIR数字滤波器的特点

　　a.封闭函数

IIR数字滤波器的系统函数可以写成封闭函数的形式。

b.IIR数字滤波器采用递归型结构

IIR数字滤波器采用递归型结构，即结构上带有反馈环路。

IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。

由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

c.借助成熟的模拟滤波器的成果

IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。

在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

d.需加相位校准网络

　IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络

3.2.2IIR数字滤波器的设计方法

IIR数字滤波器的最通用的方法是借助于模拟滤波器的设计方法，基于原型模拟传输函数的双线性变换。

模拟传输函数通常有四种模型：

巴特沃兹、切比雪夫1型和2型，椭圆模型。

模拟滤波器设计已经有了相当成熟的技术和方法，有完整的设计公式，还有比较完整的图表可以查询，因此设计数字滤波器可以充分利用这些丰富的资源来进行。

对于IIR数字滤波器的设计具体步骤如下：

（1）按照一定的规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。

（2）根据转换后的技术指标设计模拟低通滤波器G（s）（G（s）是低通滤波器的传递函数）。

（3）再按照一定的规则将G（s）转换成H（z）（H（z）是数字滤波器的传递函数）。

若设计的数字滤波器是低通的，上述的过程可以结束，若设计的是高通、带通或者是带阻滤波器，那么还需要下面的步骤：

将高通、带通或带阻数字滤波器的技术指标转换为低通模拟滤波器的技术指标，然后设计出低通G（s），再将G（s）转换为H（z）。

本设计应用matlab进行设计，对模型直接应用Matlab提供的函数来估计阶数，然后确定传输函数。

本设计设计了四种模型：

巴特沃兹、切比雪夫1型和2型，椭圆模型。

四、与设计相关的背景知识

1、功率谱密度

功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。

这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆（ohm）时的实际功率。

由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。

维纳-辛钦定理（Wiener-Khinchintheorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。

信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。

如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。

功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变化情况。

保留了频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。

有两点需要注意：

（1）.功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。

（随机的频域序列）

（2）.功率概念和幅度概念的差别。

此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶矩是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。

功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。

一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。

功率谱密度的定义是单位频带内的“功率”（均方值）

功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。

数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。

2、uigetfile命令

功能：

创建标准的对话框并通过交互式操作取得文件名

uigetfile：

显示一个模态对话框，对话框列出了当前目录下的文件和目录，用于可以选择一个将要打开的文件名。

如果文件名是有效的且该文件存在，则当用户点击Open时函数uigetfile返回该文件名。

若不存在，uigetfile显示一个控制返回对话框值的错误提示信息，此时用户可以输入另外的文件名或点击Cancel按钮。

如果用户点击Cancel按钮或关闭对话框，函数uigetfile将返回0。

[FileName,PathName,FilterIndex]=uigetfile（FilterSpec）：

创建并显示对话框,列出当前目录下的文件。

FilterSpec决定文件的初始显示，它可以为一个文件全名或者包含通配符*，如‘*.doc’表示列出所有的文档文件。

如果FilterSpec是一个包含文件名的字符串，则该文件名所在的文件名域将被选中并显示出来，且该文件的扩展名将作为过滤作用；如果FilterSpec是一个包含文件路径的字符串，例如‘.’,’..’或’/’,例如‘../*.m’表示列出当前目录上的所有M文件；如果FilterSpec是一个单元数组，则其第一列为扩展名列表，第二列为描述列表；如果FilterSpec没有指定路径对象，则uigetfile将使用默认路径下的所有文件类型。

此函数的用法为

[FileName,PathName,FilterIndex]=uigetfile（FilterSpec,DialogTitle,DefaultName）

FileName：

返回的文件名

PathName：

返回的文件的路径名

FilterIndex：

选择的文件类型

FilterSpec：

文件类型设置

DialogTitle：

打开对话框的标题

DefaultName：

默认指向的文件名

3、Imread读取图像文件

函数功能：

该函数用于读取图片文件中的数据。

在matlab的命令窗口中输入docimread或者helpimread即可获得该函数的帮助信息。

调用格式：

　A=imread（filename,fmt）

　[X,map]=imread（filename,fmt）

　　[...]=imread（filename）

　　[...]=imread（URL,...）

　　[...]=imread（...,idx）（CUR,GIF,ICO,andTIFFonly）

　　[...]=imread（...,'PixelRegion',{ROWS,COLS}）（TIFFonly）

　　[...]=imread（...,'frames',idx）（GIFonly）

　　[...]=imread（...,ref）（HDFonly）

　　[...]=imread（...,'BackgroundColor',BG）（PNGonly）

　　[A,map,alpha]=imread（...）（ICO,CUR,andPNGonly

4、Imwrite图像数据保存

函数功能：

将图像数据写入到图像文件中，存储在磁盘上。

在matlab命令窗口中键入docimwrite或helpimwrite可以获得更多关于该函数的帮助信息。

　　调用格式：

　　imwrite（A,filename,fmt）

　　A是图像数据，filename是目标图像名字，fmt是要生成的图片的格式。

图像格式有：

bmp（1-bit、8-bit和24-bit）、gif（8-bit）、hdf、jpg（或jpeg）（8-bit、12-bit和16-bit）、jp2或jpx、pbm、pcx（8-bit）、pgm、png、pnm、ppm、ras、tif（或tiff）、xwd。

各种格式支持的图像位数不一样，比如bmp格式不支持16-bit，而png格式支持，又如gif只支持8-bit格式。

5、Imnoise加噪声处理函数

函数imnoise来使用噪声污染一幅图像，该函数的基本语法为：

　　g=imnoise（f,type,parameters）

　　f为是输入图像。

函数imnoise在给图像添加噪声之前，将它转换为范围[0,1]内的double类图像。

指定噪声参数时必须考虑到这一点。

　　g=imnoise（f,'gaussian',m,var）将均值M,方差为var的高斯噪声加到图像f上，默认值为均值是0，方差是0.01的噪声。

　　g=imnoise（f,'localvar',V）将均值为0，局部方差为V的高斯噪声添加到图像f上，其中V是与f大小相同的一个数组，它包含了每一个点的理想方差值。

　　g=imnoise（f,'localvar',image_intensity,var）将均值为0的高斯噪声添加到图像f中，其中噪声的局部方差var是图像f的亮度值的函数。

参量image_intensity和var是大小相同的向量，plot（image_intensity,var）绘制出噪声方差和图像亮度的函数关系。

向量image_intensity必须包含范围在[0,1]内的归一化亮度值。

　　g=imnoise（f,'salt&pepper',d）用椒盐噪声污染图像f，其中d是噪声密度（即包括噪声值的图像区域的百分比）。

因此，大约有d*numel（f）个像素受到影响。

默认的噪声密度为0.05。

　　g=imnoise（f,'speckle',var）用方程g=f+n*f）将乘性噪声添加到图像f上，其中n是均值为0，方差为var的均匀分布的随机噪声，var的默认值是0.04。

　　g=imnoise（f,'poisson'）从数据中生成泊松噪声，而不是将人工的噪声添加到数据中，为了遵守泊松统计，unit8和unit16类图像的亮度必须和光子的数量相符合。

当每个像素的光子数量大于65535时，就要使用双精度图像。

亮度值在0到1之间变化，并且对应于光子的数量除以10e12。

五、设计过程

5.1图像加噪处理及功率密度谱对比

（1）陈舒龙加噪前

加噪后

（2）苏青文

加噪前

加噪后

（3）方雪松加噪前

加噪后

5.2FIR滤波器

1，布莱克曼'blackman'

经多次调试下面参数滤波效果最好。

Fp=2200;Fs=3600;Rp=0.6;Rs=15;FT=20000;

[num,den]=FIR（Fp,Fs,Rp,Rs,FT,'blackman'）;

2，海明窗'hamming'

经多次调试下面参数滤波效果最好。

Fp=2400;Fs=3200;Rp=0.6;Rs=15;FT=20000;

[num,den]=FIR（Fp,Fs,Rp,Rs,FT,'hamming'）;

3，汉宁窗

经多次调试下面参数滤波效果最好。

Fp=2400;Fs=2800;Rp=0.6;Rs=15;FT=20000;

[num,den]=FIR（Fp,Fs,Rp,Rs,FT,'hanning'）;

4，多尔夫—切比雪夫窗

经多次调试下面参数滤波效果最好。

Fp=2300;Fs=3100;Rp=0.6;Rs=15;FT=20000;

[num,den]=FIR（Fp,Fs,Rp,Rs,FT,'cheby'）;

5，巴特利特窗'battlet'

经多次调试下面参数滤波效果最好。

Fp=1960;Fs=5200;Rp=0.6;Rs=15;FT=20000;

[num,den]=FIR（Fp,Fs,Rp,Rs,FT,'battlet'）;

5.3IIR数字滤波器设计

关键是确定Wp，Ws，Rp，Rs参数

用fornn=1:

1Fp=1940+nn*100;Fs=5300+nn*500;循环5次，Fs每次先500甚至更大增大，观察图像滤波效果，选择两张图像较好的进行大定位，然后以在两个值之间进行小范围循环5次，观察滤波效果，多次尝试，最终选择一个最好的参数。

1、巴特沃兹滤波

Wp=0.5;Ws=0.6;Rp=1.5;Rs=36;此参数效果最好。

[numden]=iirlp（Wp,Ws,'butter'）;

2.切比雪夫1

Wp=0.3;Ws=0.6;Rp=1.5;Rs=38;此参数效果最好。

[numden]=iirlp（Wp,Ws,'cheby1'）;

3.切比雪夫2型

Wp=0.5;Ws=0.6;Rp=1.5;Rs=38;此参数效果最好。

[numden]=iirlp（Wp,Ws,'cheby2'）;

4、椭圆滤波

Wp=0.5;Ws=0.6;Rp=1.5;Rs=8;此参数效果最好

[numden]=iirlp（Wp,Ws,'ellip'）;

比较结果分析

利用巴特沃兹、切比雪夫1、切比雪夫2、椭圆滤波器都可以进行低通、高通、带通、带阻滤波器的设计，但是各有特点。

巴特沃兹滤波器通带内的幅频响应曲线能得到最大限度的平滑，但牺牲了截止频率的坡度；切比雪夫1型滤波器通带内等波纹，阻带内单调；切比雪夫2型滤波器通带内单调，阻带内等波纹；椭圆滤波器阻带内和通带内都是等波纹的，但下降的坡度更大，而且可以以更低的阶数实现和其他两类滤波器一样的性能指标。

六、FIR和IIR比较

a.从性能上进行比较

从性能上来说，IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。

但是这个高效率是以相位的非线性为代价的。

选择性越好，则相位非线性越严重。

相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高5~10倍，结果，成本较高，信号延时也较大；如果按相同的选择性和相同的线性要求来说，则IIR滤波器就必须加全通网络进行相位较正，同样要大增加滤波器的节数和复杂性。

b.从结构上看

IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。

另外，在这种结构中，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引入寄生振荡。

相反，FIR滤波器主要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，运算误差也较小。

此外，FIR滤波器可以采用快速付里叶变换算法，在