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

1、 数字信号处理 课程设计任务书20112012学年 第一学期第 15 周 19 周 题目基于matlab的图像处理滤波器设计内容及要求：1、设计一个低通FIR滤波器和一个低通IIR滤波器；2、分析比较上述两种滤波器的优劣；3、分析上述两个模型的幅频特性、相频特性、相延迟、群延迟。进度安排： 1、任务分配、查阅资料 2天 2、方案论证 3天 3、分析、设计、调试程序 5天 4、书写、整理实验报告和小结 3天成员组成：09044106苏青文 08陈舒龙（组长） 09方雪松 指导时间： 指导地点： F 618 任务下达2011年12 月 16 日任务完成2012年1月 6日考核方式1.评阅 2.答辩

2、 3.实际操作 4.其它 指导教师汪传忠系（部）主任王长坤注：1、此表一组一表二份，课程设计小组组长一份；任课教师授课时自带一份备查。2、课程设计结束后与“课程设计小结”、“学生成绩单”一并交院教务存档。目录摘要3一、课程设计目的及要求4二、课程设计内容及任务安排4 三、设计原理及设计方法53.1 FIR数字滤波器5 3.2 IIR数字滤波器7四、与设计相关的知识8五、设计过程14 5.1图像加噪处理及功率密度谱对比14 5.2 FIR滤波器的设计22 (1)布莱克曼窗.22 (2)海明窗26 （3）汉宁窗.28 （4）多尔夫-切比雪夫窗.32 （5）巴特利特窗.35 5.3 IIR 数字滤波

3、器设计.38 （1）巴特沃兹.38 （2）切比雪夫1.41 （3）切比雪夫2 .43 （4）椭圆滤波.45六、FIR和IIR的比较.47七、个人设计总结49附录50摘要MATLAB软件具有使用简单、方便，易编程，语言简练，函数库可任意扩充，采用全新数据类型和面向对象编程技术等特点，具有数值分析容易、语法结构简单、矩阵数值运算高效、图形功能完备、图像处理方便等强大的功能，因此已被广泛应用在教学、科研和工程设计的各个领域。传统的滤波器分析与设计均使用繁琐的公式计算，改变参数后需要重新运算，从而在分析与设计滤波器尤其是高阶滤波器时工作量特别大，利用MATLAB信号处理工具箱可以快速而有效地实现滤波器

4、的分析、设计与仿真。本课设介绍基于MATLAB环境的IIR数字滤波器和FIR数字滤波器的设计方法，原理和步骤。IIR数字滤波器的常用设计方法有脉冲响应不变法和双线性变换法等，FIR数字滤波器的设计方法只要有窗函数设计法，频率抽样设计法和最优化设计法等，最后利用设计的滤波器对含噪声的图像进行滤波处理并验证所设计滤波器的性能，同时比较两种设计方法的最优性。关键字：数字信号处理 Matlab FIR滤波器 IIR滤波器 图像加噪 图像滤波 FIR与IIR滤波器比较 1、 课程设计目的及要求1.1目的 综合运用数字信号处理的理论知识进行频谱分析和滤波器设计FIR和IIR滤波器，通过理论推导得出相应结论

5、，再利用 MATLAB 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。 1、巩固和加深理解数字信号处理基本理论知识，提高综合运用这些知识进行程序设计的能力。2、能根据需要查阅相关手册和资料，提高独立分析问题和解决问题的能力。3、能正确使用仪器设备，掌握测试原理，掌握综合性程序设计的基本过程。4、掌握时域离散时间信号和时域离散系统基本概念和基本性质。5、掌握时域离散信号、系统的频域分析工具。6、掌握基于Matlab的数字信号处理领域的计算机辅助分析与设计的基本方法。7、培养团队协作精神、创新意识、严肃认真的治学态度和严谨求实的工作作风1.2基本要求 1、 学会 MATLAB 的

6、使用，掌握 MATLAB 的程序设计方法。 2、掌握在 Windows 环境下数字图像处理的方法。 3、 掌握数字信号处理的基本概念、基本理论和基本方法。 4、 掌握 MATLAB 设计 FIR 和 IIR 数字滤波器的方法。 5、 学会用 MATLAB 对数字图像进行分析和处理二、课程设计内容及任务安排 每人用自己的一张照片，并对照片进行加噪声处理，通过设计滤波器滤除噪声。每组三人，其中一人设计FIR滤波器进行噪声处理，一人设计IIR滤波器滤除噪声，另外一人分析FIR和IIR滤波器的优劣比较。其中在进行滤波器设计中，要求每种模型都要求设计，FIR滤波器设计包括布莱克曼窗、海明窗、汉宁窗、多尔

7、夫切比雪夫窗、巴特利特窗，IIR滤波器设计包括巴特沃兹、切比雪夫1和2，椭圆模型。并且分析每个模型的幅频特性，相频特性，相延时，群延时。最后得出结论。具体实施方案：1、苏青文主要负责设计IIR无限脉冲响应数字滤波器。2、方雪松主要负责设计FIR有限脉冲响应数字滤波器。3、组长陈舒龙主要负责图像处理、整理报告及FIR和IIR两种数字滤波器的比较。三、设计原理及设计方法数字滤波器设计的一个重要步骤就是确定一个可实现的传递函数G(z)来逼近指定的频率响应。在G(z)确定下来之后，接下来要设计用一种合适的滤波器结构来实现它。滤波器设计过程的第一步是估计传输函数阶数。然后根据各种模型用Matlab函数得

8、出传输函数。3.1 FIR数字滤波器3.1.1 基本介绍FIR是有限冲激响应（Finite Impulse Response）的简称。 由线性系统理论可知，在某种适度条件下，输入到线性系统的一个冲击完全可以表征系统。当我们处理有限的离散数据时，线形系统的响应（包括对冲击的响应）也是有限的。若线性系统仅是一个空间滤波器，则通过简单地观察它对冲击的响应，我们就可以完全确定该滤波器。通过这种方式确定的滤波器称为有限冲击响应（FIR)滤波器. 3.1.2 设计方法 FIR滤波器设计的最简单方法是，对期望的理想频率响应进行离散时间傅里叶逆变换，得到双无限长冲击响应，对所得的无限系数简单地截尾为有限系数。

9、滤波器类型确定之后，接下来是估计滤波器阶数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.543.93.11/M海明窗8/(2M+1)42.754.5 3.32/M布拉克曼窗 12/(2M+1)5

10、8.175.3 5.56/M3.2 IIR数字滤波器3.2.1 IIR数字滤波器的特点 a.封闭函数IIR数字滤波器的系统函数可以写成封闭函数的形式。 b.IIR数字滤波器采用递归型结构IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。 c.借助成熟的模拟滤波器的成果IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量

11、比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 d.需加相位校准网络IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络3.2.2 IIR数字滤波器的设计方法IIR数字滤波器的最通用的方法是借助于模拟滤波器的设计方法，基于原型模拟传输函数的双线性变换。模拟传输函数通常有四种模型：巴特沃兹、切比雪夫1型和2型，椭圆模型。模拟滤波器设计已经有了相当成熟的技术和方法，有完整的设计公式，还有比较完整的图表可以查询，因此设计数字滤波器可以充分利用这些丰富的资源来进行。 对于I

12、IR数字滤波器的设计具体步骤如下： (1)按照一定的规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。 (2)根据转换后的技术指标设计模拟低通滤波器G(s)(G(s)是低通滤波器的传递函数)。 (3)再按照一定的规则将G(s)转换成H(z)(H(z)是数字滤波器的传递函数)。若设计的数字滤波器是低通的，上述的过程可以结束，若设计的是高通、带通或者是带阻滤波器，那么还需要下面的步骤：将高通、带通或带阻数字滤波器的技术指标转换为低通模拟滤波器的技术指标，然后设计出低通G(s)，再将G(s)转换为H(z)。 本设计应用matlab进行设计，对模型直接应用Matlab提供的函数来估计阶数

13、，然后确定传输函数。本设计设计了四种模型：巴特沃兹、切比雪夫1型和2型，椭圆模型。四、与设计相关的背景知识1、功率谱密度 功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆(ohm)时的实际功率。由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。维纳-辛钦定理（Wiener-Khinchin theorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。 信号的功率谱密度当且仅当信号是广义

14、的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变化情况。保留了频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。 有两点需要注意：（1）. 功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列）（2）. 功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存

15、在性取决于二阶矩是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。 功率谱密度的定义是单位频带内的“功率”（均方值） 功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。2、ui

16、getfile命令功能：创建标准的对话框并通过交互式操作取得文件名uigetfile：显示一个模态对话框，对话框列出了当前目录下的文件和目录，用于可以选择一个将要打开的文件名。如果文件名是有效的且该文件存在，则当用户点击Open时函数uigetfile返回该文件名。若不存在，uigetfile显示一个控制返回对话框值的错误提示信息，此时用户可以输入另外的文件名或点击Cancel按钮。如果用户点击Cancel按钮或关闭对话框，函数uigetfile将返回0。FileName,PathName,FilterIndex = uigetfile(FilterSpec)：创建并显示对话框,列出当前目录下

17、的文件。FilterSpec决定文件的初始显示，它可以为一个文件全名或者包含通配符*，如*.doc表示列出所有的文档文件。如果FilterSpec是一个包含文件名的字符串，则该文件名所在的文件名域将被选中并显示出来，且该文件的扩展名将作为过滤作用；如果FilterSpec是一个包含文件路径的字符串，例如.,.或/,例如./*.m表示列出当前目录上的所有M文件；如果FilterSpec是一个单元数组，则其第一列为扩展名列表，第二列为描述列表；如果FilterSpec没有指定路径对象，则uigetfile将使用默认路径下的所有文件类型。此函数的用法为FileName,PathName,Filter

18、Index= uigetfile(FilterSpec,DialogTitle,DefaultName)FileName：返回的文件名PathName：返回的文件的路径名FilterIndex：选择的文件类型FilterSpec：文件类型设置DialogTitle：打开对话框的标题DefaultName：默认指向的文件名3、Imread 读取图像文件 函数功能：该函数用于读取图片文件中的数据。在matlab的命令窗口中输入doc imread或者help imread即可获得该函数的帮助信息。 调用格式：A = imread(filename,fmt) X,map = imread(filen

19、ame,fmt) . = imread(filename) . = imread(URL,.) . = imread(.,idx) (CUR, GIF, ICO, and TIFF only) . = imread(.,'PixelRegion',ROWS, COLS) (TIFF only) . = imread(.,'frames',idx) (GIF only) . = imread(.,ref) (HDF only) . = imread(.,'BackgroundColor',BG) (PNG only) A,map,alpha = i

20、mread(.) (ICO, CUR, and PNG only4、Imwrite 图像数据保存 函数功能：将图像数据写入到图像文件中， 存储在磁盘上。在matlab命令窗口中键入doc imwrite或help imwrite可以获得更多关于该函数的帮助信息。 调用格式： 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）

21、、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上

22、，默认值为均值是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必须包含范

23、围在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和

24、unit16类图像的亮度必须和光子的数量相符合。当每个像素的光子数量大于65535时，就要使用双精度图像。亮度值在0到1之间变化，并且对应于光子的数量除以10e12。五、设计过程5.1 图像加噪处理及功率密度谱对比(1)陈舒龙 加噪前加噪后 (2)苏青文 加噪前加噪后(3)方雪松 加噪前加噪后5.2 FIR滤波器1，布莱克曼'blackman'经多次调试下面参数滤波效果最好。Fp=2200; Fs=3600; Rp=0.6; Rs=15; FT=20000; num,den=FIR(Fp,Fs,Rp,Rs,FT,'blackman'); 2，海明窗'ha

25、mming'经多次调试下面参数滤波效果最好。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,d

26、en=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.3 IIR数字滤波器设计关键是确定Wp，Ws，Rp，Rs参数用 for nn=1:1 Fp=1940+nn*100;Fs=5300+nn*500;循环5次，Fs每次先500甚至更大增大，观察图像滤波效果，选择两张图像较好的进行大定位，然后以在两个值之间

27、进行小范围循环5次，观察滤波效果，多次尝试，最终选择一个最好的参数。1、巴特沃兹滤波Wp=0.5; Ws=0.6;Rp=1.5; Rs=36;此参数效果最好。 num den=iirlp(Wp,Ws,'butter'); 2. 切比雪夫1 Wp=0.3; Ws=0.6;Rp=1.5; Rs=38;此参数效果最好。 num den=iirlp(Wp,Ws,' cheby1');3.切比雪夫2型Wp=0.5; Ws=0.6;Rp=1.5; Rs=38;此参数效果最好。num den=iirlp(Wp,Ws,' cheby2');4、椭圆滤波Wp=0.

28、5; Ws=0.6;Rp=1.5; Rs=8;此参数效果最好num den=iirlp(Wp,Ws,'ellip');比较结果分析利用巴特沃兹、切比雪夫1、切比雪夫2、椭圆滤波器都可以进行低通、高通、带通、带阻滤波器的设计，但是各有特点。巴特沃兹滤波器通带内的幅频响应曲线能得到最大限度的平滑，但牺牲了截止频率的坡度；切比雪夫1型滤波器通带内等波纹，阻带内单调；切比雪夫2型滤波器通带内单调，阻带内等波纹；椭圆滤波器阻带内和通带内都是等波纹的，但下降的坡度更大，而且可以以更低的阶数实现和其他两类滤波器一样的性能指标。六、FIR和IIR比较a.从性能上进行比较 从性能上来说，IIR滤

29、波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的。选择性越好，则相位非线性越严重。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高510倍，结果，成本较高，信号延时也较大；如果按相同的选择性和相同的线性要求来说，则IIR滤波器就必须加全通网络进行相位较正，同样要大增加滤波器的节数和复杂性。b.从结构上看IIR滤波器必须采用递归结构，极点位置必须在

30、单位圆内，否则系统将不稳定。另外，在这种结构中，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引入寄生振荡。相反，FIR滤波器主要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，运算误差也较小。此外，FIR滤波器可以采用快速付里叶变换算法，在相同阶数的条件下，运算速度可以快得多。两者的优缺点如下：FIR有限冲激响应滤波器原理：无反馈（无递归），极点在原点，卷积（或FFT卷积，部分也可用递归计算）优点：线性相位，频谱特性好缺点：速度较慢（由于阶数多延时长）应用场合：需线性相位的滤波（如图像信号等），线性相位是指不同频率的信号经过滤波器其延时一样。FIR的应用灵活

31、，易于适应某些特殊的应用。IIR无限冲激响应滤波器原理：有反馈（有递归），极点可位于单位圆内任何地方，递归计算优点：阶数少（延时少），借助成熟的模拟滤波器的成果（如巴特沃斯、契比雪夫和椭圆滤波器等）缺点：非线性相位，易出现震荡（不稳定，因为递归计算）应用场合：不需线性相位的滤波（如音频信号等）总之，可以将IIR滤波器和FIR滤波器的性能指标以及优缺点归纳如下：IIR 数字滤波器FIR 数字滤波器（1）相位一般是非线性的（1）相位可以做到严格线性（2）不一定稳定（2）一定是稳定的（3）对频率分量的选择性好（零极点可同时起作用）（3）选择性差（4）相同性能下阶次较低（4）相同性能下阶次高（5）有噪

32、声反馈，噪声大（5）噪声小（6）设计有封闭形式的公式，一次完成（6）没有封闭形式的设计公式，须靠经验与反复调试（7）对计算手段的要求较低（7）一般需用计算机计算七、个人设计总结1.陈舒龙本次课设我作为组长给咱组员分配任务，不过由于苏青文同学家里有事提前回家了，他之前负责做的IIR滤波器设计不太准确，所以我就做了较大修改。通过对巴特沃兹、切比雪夫1型和2型、椭圆型滤波器的设计，从图形直观分析验证理论知识。比较这四种滤波器发现：巴特沃兹滤波器有最宽的过渡带，有单调递减的增益响应。两种类型的切比雪夫滤波器有相等宽度的过渡带，该宽度小于巴特沃兹但大于椭圆。切比雪夫1型在过渡带中有一个比切比雪夫2型稍微

33、快一些的下降。在通带中，切比雪夫2型的幅度响应与巴特沃兹的几乎一样。椭圆有最窄的过渡带并具有等波纹通带和等波纹阻带响应。 另外我的本职是分析比较IIR和FIR，其中，FIR滤波器的最大好处是稳定、线性相位和广泛的适用范围，而它的最大缺点是阶数高，从而带来时延大、存储单元多等问题。但IIR与FIR滤波器各有所长，所以在实际应用时应该从多方面考虑来加以选择。从使用要求上来看，在对相位要求不敏感的场合，如语言通信等，选用IIR较为合适，这样可以充分发挥其经济高效的特点；对于图像信号处理，数据传输等以波形携带信息的系统，则对线性相位要求较高。如果有条件，采用FIR滤波器较好。当然，在实际应用中可能还要

34、考虑更多方面的因素。不论IIR和FIR，阶数越高，信号延迟越大；同时在IIR滤波器中，阶数越高，系数的精度要求越高，否则很容易造成有限字长的误差使极点移到单位圆外。因此在阶数选择上是综合考虑的。2.方雪松通过本次的课程设计，让我平时在课本上学习的知识得以运用。在这次课设中，我主要任务是负责FIR滤波器的设计，虽然设计过程中遇到许多问题，但通过与同学的交流以及在网上查阅相关资料都逐渐得以解答。通过对FIR滤波器的设计让我认识到它可进行准确的线性相位设计且结构具有稳定的量化滤波器系数，FIR滤波器可以得到严格的线性相位，然而由于FIR滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到高的选择

35、性，成本较高，信号延时也较大。经过本次课程设计后，让我对以前在课堂上学习的知识有了更深层次的理解以及FIR滤波器的设计和基本应用，最后和其他两位组员合作共同完成了本次课程设计。附录：%图像滤波function DE_NOISE_IMAGEfilename, pathname = uigetfile( .'*.jpg','MATLAB Files (*.jpg)' '*.m', 'M-files (*.m)' . '*.fig','Figures (*.fig)' . '*.mat',

36、'MAT-files (*.mat)' . '*.mdl','Models (*.mdl)' . '*.*', 'All Files (*.*)', . 'Pick a file');q,qq=size(pathname);S=''q1,qq1=size(filename);RNAME=pathname;RNAME(qq+1)=S(1);RNAME(qq+2:qq+1+qq1)=filename(1:qq1);sh=imread(RNAME,'jpg');%选被噪声污

37、染的图象；M,N,O=size(sh);%sh=double(sh)/255;R=rand*0.025;if R<0.015 R=0.015;enddh1=sh;for o=1:3 dh1(:,:,o)=imnoise(dh1(:,:,o),'gaussian',0.015);enddh=dh1;dhL=dh1;dhI=dh1;%参数设计% Wp=0.5; Ws=0.6; Rp=1; Rs=40; %'iir_type butter','cheby1','cheby2','ellip'.%fir-type &

38、#39;blackman'',''hamming'',''hanning'',''cheby'',''kaiser''%num den=firlp(Wp,Ws,'blackman');num den=iirlp(Wp,Ws,'ellip');%从此处选择滤波器类型for o=1:3 A(1:M,1:N)=dh(1:M,1:N,o); AL=A; AI=A; %用设计的滤波器对图象滤波 %按列处理； for k=1:N x

39、=AL(:,k); AL(:,k)=filter(num,den,x);%设计出的IIR滤波器的参数num,den end dhL(1:M,1:N,o)=AL(1:M,1:N); %行列终和 %DH(1:M,1:N,o)=(dhI(1:M,1:N,o)+dhL(1:M,1:N,o)/2; DH(1:M,1:N,o)=dhL(1:M,1:N,o); %有约束滤波 endfigure;%图片subplot(1,3,1);imshow(sh);title('原始图像');subplot(1,3,2);imshow(dh1);title('加噪图像');subplot(

40、1,3,3);imshow(DH);title('去噪图像');figure;%功率谱subplot(3,1,1);PGJ_INPUT(sh);grid on;axis(0,1,-50,40);title('原始图像功率谱');subplot(3,1,2);PGJ_INPUT(dh1);grid on;axis(0,1,-50,40);title('加噪功率谱');subplot(3,1,3);PGJ_INPUT(DH);grid on;%axis(0,1,-100,60);title('去噪功率谱');end%FIR滤波器设计f

41、unction num den=firlp(Wp,Ws,windows)%N num den=firlp(windows,ftype)%函数名：FIR低通滤波器的设计.%参数：窗函数windows,滤波器响应类型 ftype.%返回：滤波器的阶数N和系数num den.%滤波器指标 Rp=0.1;Rs=25; dp=1-10(-Rp/20);ds=10(-Rs/20); %估计阶数N Wn,beta,ftype=kaiserord(Wp Ws,1 0,dp ds);N=N+18;switch lower(windows) case 'blackman' taper=blackm

42、an(N+1); case 'hamming' N=N-10; taper=hamming(N+1); case 'hanning' N=N-13; taper=hanning(N+1); case 'cheby' N=N+6; taper=chebwin(N+1); case 'kaiser' N=N-15; taper=kaiser(N+1,beta); case 'battlet' N=N+120; taper=bartlett(N+1); otherwise error('Unknown windo

43、ws of ''%s'' ,please select one from ''blackman'',''hamming'',''hanning'',''cheby'',''kaiser''.',windows)endnum=fir1(N,Wn,ftype,taper);den=1; freqRespone1(num,den,windows,Wp,Ws);%频率响应(subf)%function f

44、reqRespone(num,den,FT,methods,Wp,Ws)function freqRespone(num,den,methods,Wp,Ws)h wa=freqz(num,den,512);%滤波器频率响应figure,subplot(2,2,1);plot(wa,(abs(h);gridylabel('振幅(dB)');xlabel('频率(Hz)');title('频率响应','(IIR-',methods,'逼近','Wp=',num2str(Wp),' Hz,'

45、;,'Ws=',num2str(Ws),' Hz.)');subplot(2,2,2);plot(wa,unwrap(angle(h)*180/pi);gridxlabel('频率(Hz)');ylabel('相位(degree)');subplot(2,2,3);plot(wa,grpdelay(h);xlabel('频率(Hz)');ylabel('群延时(grpdelay)');subplot(2,2,4);plot(wa,phasedelay(h);xlabel('频率(Hz)

46、9;);ylabel('相延时(phasedelay)')figure,zplane(num,den);title('零极点图','(IIR-',methods,'逼近)');%频率响应%function freqRespone1(num,den,windows,Wp,Ws)h,wa=freqz(num,den,512);%滤波器频率响应figure;subplot(2,2,1);xlabel('阶数(N)'),ylabel('h(N)');title('冲激响应','(N=

47、',. num2str(length(num),',FIR-',windows,'窗','Wp=',num2str(Wp),' Hz,','Ws=',num2str(Ws),' Hz.)');plot(wa,abs(h);gridxlabel('频率(Hz)');ylabel('振幅(dB)');subplot(2,2,2);plot(wa,unwrap(angle(h)*180/pi);gridxlabel('频率(Hz)');ylabel(

48、'相位(degree)');subplot(2,2,3);plot(wa,grpdelay(h);xlabel('频率(Hz)');ylabel('群延时(grpdelay)');subplot(2,2,4);plot(wa,phasedelay(h);xlabel('频率(Hz)');ylabel('相延时(phasedelay)')figure,zplane(num,den);%title('零极点图'(FIR-',windows,'窗')%双向性变换法设计IIR滤波器%

49、function num den=iirlp(Wp,Ws,methods)%N num den=iirlp(methods,ftype)%函数名：IIR低通滤波器的设计.%参数：滤波器的逼近类型 methods：'butter','cheby1','cheby2','ellip'.%返回：滤波器的阶数N和传输函数num den.ftype='low'%滤波器设计%滤波器指标Rp=1;Rs=20;%设计switch lower(methods) case 'butter'%巴特沃斯逼近 N,Wn=buttord(Wp,Ws,Rp,Rs); num den=butter(N,Wn,ftype); freqRespone(num,den,methods,Wp,Ws); case 'cheby1'%切比