1、 课程设计任务书学生姓名: 专业班级: 指导教师: 工作单位: 题 目: 基于MATLAB的图像滤波设计 初始条件:1.MATLAB软件 2.滤波器处理相关函数要求完成的主要任务: (1)读入图像并分别加入高斯噪声、椒盐噪声和乘性噪声,并比较结果。(2)设计巴特沃斯低通滤波对图像进行低通滤波处理,显示结果。(3)设计高斯高通滤波器对图像进行处理,显示结果。(4)采用维纳滤波和中值滤波对图像进行处理,显示结果参考书:1. 信号与系统第一版 刘泉 江雪梅 主编 高等教育出版社2. 数字图像处理MATLAB版 冈萨雷斯 主编 电子工业出版社时间安排:第15周:任务安排、分组第16周:理论设计及仿真第
2、18周:撰写设计报告及答辩指导教师签名: 年 月 日系主任(或责任教师)签名: 年 月 日39摘 要31.MATLAB简介51.1 MATLAB的概况51.2 MATLAB产生的历史背景52.编程及运行结果72.1常见基本运算72.1.1极限的计算72.1.2微分的计算72.1.3积分的计算82.1.4级数的计算92.1.5求解代数方程102.1.6求解常微分方程102.2 矩阵基本计算112.2.1矩阵的最大值112.2.2矩阵的最小值112.2.3矩阵的均值122.2.4矩阵的方差132.2.5矩阵的转置132.2.6矩阵的逆142.2.7矩阵的行列式152.2.8矩阵的特征值计算152.
3、2.9矩阵的相乘162.2.10矩阵的右除和左除172.2.11矩阵的幂运算182.3 多项式基本计算182.3.1多项式加减运算182.3.2多项式乘除运算192.3.3多项式求导202.3.4求根和求值运算202.3.5多项式的部分分式展开212.3.6多项式的拟合222.3.7插值运算233.基于MATLAB的图像滤波设计253.1读入图像并分别加入高斯噪声、椒盐噪声和乘性噪声,并比较结果253.2设计巴特沃斯低通滤波对图像进行低通滤波处理,显示结果293.2.1叠加椒盐噪声的巴特沃斯低通滤波293.2.2叠加高斯噪声的巴特沃斯低通滤波313.2.3叠加乘性噪声的巴特沃斯低通滤波323.
4、3用MATLAB实现高斯高通滤波器对图像的处理333.4维纳滤波和中值滤波对图像进行处理354.总结38参考文献39摘 要 现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因,使得具有线性相位的FIR数字滤波器得到大力发展和广泛应用。在实际进行数字信号处理时,往往需要把信号的观察时间限制在一定的时间间隔内,只需要选择一段时间信号对其进行分析。取用有限个数据,即将信号数据截断的过程,就等于将信号进行加窗函数操作。这样操作以后,常常会发生频谱分量从其正常频谱扩展开来的现象,即所谓的“频谱泄漏”。当进行离散傅立叶变换时,时域中的截断是必需的,因此泄漏效应也是离散傅立叶变换所固有的,必须进行
5、抑制。而要对频谱泄漏进行抑制,可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣,或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。数字带通滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位重启响应函数的时域特性可分为两类:无限冲击响应滤波器(IIR),有限冲击响应滤波器(FIR)。与IIR滤波器相比,FIR的实现是递归的,总是稳定的;更重要的是,FIR滤波器在满足幅频响应要求的同时,可以获得严格的线性相位特性。因此,它在高保真的信号处理,如信号音频,图像处理,数据传输等领域得到广泛的应用。数字fir滤波
6、器的设计方法有很多种。如窗函数法设计,频率采样设计法和最优化设计法等。AbstractModern images, sounds, data communication of linear phase requirement is common. It is this reason, make with linear phase FIR digital filter to get strong development and extensive application.In the practical digital signal processing, often need to signa
7、l the observation time limit in certain interval of time, only need to choose a time signal to analyze it. So, take a finite number of data, forthcoming truncated signal data process, equals will signal is added window function operation. And such operation later, often happen spectrum component fro
8、m its normal phenomenon of spread spectrum, the so-called frequency spectrum leakage. When performing discrete Fourier transform, the time-domain truncated is necessary, therefore leakage effect is discrete Fourier transform inherent, must undertake inhibition. And in the behind of the FIR filters,
9、in the design for access limited long unit sampling response, need to use the window function truncation infinite long unit sampling response sequence. In addition, the power spectrum estimation also to meet a window function weighted problem. Thus, window function weighted technology in digital sig
10、nal processing the important position. Digital bandpass filter is used as a filtering time discrete signal digital system, based on sample data, mathematical treatment to achieve the purpose of frequency domain filtering. According to its unit restart response function of time domain properties can
11、be divided into two classes: infinite shock response filter (IIR), limited shock response filter (FIR). Therefore, it in high fidelity signal processing, such as signal audio, image processing, data transmission and other areas to be widely application.Digital fir filters design in many ways. Such a
12、s window function method design, frequency sampling design method and the optimum design method, etc.1.MATLAB简介1.1 MATLAB的概况 MATLAB是矩阵实验室(MatrixLaboratory)之意。除具备卓越的数值计算能力外,它还提供了专业水平的符号计算,文字处理,可视化建模仿真和实时控制等功能。 MATLAB的基本数据单位是矩阵,它的指令表达式与数学工程中常用的形式十分相似,故用MATLAB来解算问题要比用C、FORTRAN等语言完相同的事情简捷得多。当前流行的MATLAB
13、5.3/Simulink 3.0包括拥有数百个内部函数的主包和三十几种工具包(Toolbox).工具包又可以分为功能性工具包和学科工具包.功能工具包用来扩充MATLAB的符号计算,可视化建模仿真,文字处理及实时控制等功能.学科工具包是专业性比较强的工具包,控制工具包,信号处理工具包,通信工具包等都属于此类。开放性使MATLAB广受用户欢迎.除内部函数外,所有MATLAB主包文件和各种工具包都是可读可修改的文件,用户通过对源程序的修改或加入自己编写程序构造新的专用工具包。1.2 MATLAB产生的历史背景 在70年代中期,Cleve Moler博士和其同事在美国国家科学基金的资助下开发了调用EI
14、SPACK和LINPACK的FORTRAN子程序库.EISPACK是特征值求解的FOETRAN程序库,LINPACK是解线性方程的程序库.在当时,这两个程序库代表矩阵运算的最高水平。 到70年代后期,身为美国New Mexico大学计算机系系主任的Cleve Moler,在给学生讲授线性代数课程时,想教学生使用EISPACK和LINPACK程序库,但他发现学生用FORTRAN编写接口程序很费时间,于是他开始自己动手,利用业余时间为学生编写EISPACK和LINPACK的接口程序。Cleve Moler给这个接口程序取名为MATLAB,该名为矩阵(matrix)和实验室(labotatory)两
15、个英文单词的前三个字母的组合.在以后的数年里,MATLAB在多所大学里作为教学辅助软件使用,并作为面向大众的免费软件广为流传。 在当今30多个数学类科技应用软件中,就软件数学处理的原始内核而言,可分为两大类.一类是数值计算型软件,如MATLAB,Xmath,Gauss等,这类软件长于数值计算,对处理大批数据效率高;另一类是数学分析型软件,Mathematica,Maple等,这类软件以符号计算见长,能给出解析解和任意精确解,其缺点是处理大量数据时效率较低。MathWorks公司顺应多功能需求之潮流,在其卓越数值计算和图示能力的基础上,又率先在专业水平上开拓了其符号计算,文字处理,可视化建模和实
16、时控制能力,开发了适合多学科,多部门要求的新一代科技应用软件MATLAB.经过多年的国际竞争,MATLAB以经占据了数值软件市场的主导地位。在MATLAB进入市场前,国际上的许多软件包都是直接以FORTRANC语言等编程语言开发的。这种软件的缺点是使用面窄,接口简陋,程序结构不开放以及没有标准的基库,很难适应各学科的最新发展,因而很难推广。MATLAB的出现,为各国科学家开发学科软件提供了新的基础。在MATLAB问世不久的80年代中期,原先控制领域里的一些软件包纷纷被淘汰或在MATLAB上重建。 时至今日,经过MathWorks公司的不断完善,MATLAB已经发展成为适合多学科,多种工作平台的
17、功能强大大大型软件。在国外,MATLAB已经经受了多年考验。在欧美等高校,MATLAB已经成为线性代数,自动控制理论,数理统计,数字信号处理,时间序列分析,动态系统仿真等高级课程的基本教学工具;成为攻读学位的大学生,硕士生,博士生必须掌握的基本技能。在设计研究单位和工业部门,MATLAB被广泛用于科学研究和解决各种具体问题。在国内,特别是工程界,MATLAB一定会盛行起来。可以说,无论你从事工程方面的哪个学科,都能在MATLAB里找到合适的功能。2.编程及运行结果2.1常见基本运算2.1.1极限的计算MATLAB提供的命令limit()可以完成极限运算,其调用格式如下: limit(F,x,a
18、,left)该命令对表达式F求极限,独立变量x从左边趋近于a,函数中除F外的参数均可省略,left可换成right。例:求极限S= 1+1/x代码如下:clear;F=sym(1+1/x)limit(F,x,inf,left)运行结果:2.1.2微分的计算MATLAB提供的函数diff()可以完成对给定函数求导函数的运算,其调用格式如下:diff(fun,x,n),其意义是求函数fun关于变量x的n阶导数,n为1时可省略。这里的fun用上例的后一种方式来定义较为妥当。例:求函数y=log的一阶导。代码如下:clear;syms xy=log(1/x);dy=diff(y,x)运行结果: 2.1
19、.3积分的计算MATLAB中主要用int进行符号积分,用trapz,dblquad,quad,quad8等进行数值积分。R=int(s,v)%对符号表达式中指定的符号变量v计算不定积分,表达式R只是表达式函数s的一个原函数,后面没有带任何常数C。R=int(s)%对符号表达式s的确定的符号变量计算不定积分。R=int(s,a,b)%符号表达式s的定积分,a,b分别为积分的上,下限。trapz(x,y)梯形积分法,x时表示积分区间的离散化向量,y是与x同维数的向量,表示被积函数,z返回积分值。fblquad(fun,a,b,c,d)矩形区域二重积分,fun表示被积函数的M函数名,a,b分别为x的
20、上下限,c,d分别为y的上下线。例1:(不定积分)用符号积分命令int计算代码如下:clear;syms xint(x2*sin(x)运行结果: 例2:(定积分)计算数值积分用符号积分命令int,代码如下:clear;syms x;int(x4,x,-2,2)运行结果: 改用梯形积分法命令trapz计算积分,代码如下:clear;x=-2:0.01:1;y=x.2;trapz(x,y)运行结果: 2.1.4级数的计算MATLAB中主要用symsun,taylor求级数的和及Taylor展开。其中Symsum(s,v,a,b)表达式s关于变量v从a到b求和Taylor(f,a,n)将函数f在a点
21、展开为n-1阶Taylor多项式例:用Symsum计算代码如下:clear;clc;syms x yz=1/(x3);symsum(z,x,1,9)运行结果: 2.1.5求解代数方程MATLAB命令输人格式:solve(equ1,equ2,.equN),其中eqni表示第i个方程。例1:求解方程 x2+b*x+c=0代码如下:solve(a*x2+b*x+c)运行结果: 2.1.6求解常微分方程MATLAB中主要用dsolve求符号解析解,ode45,ode23,ode15s求数值解。s=dsolve(方程1, 方程2,初始条件1,初始条件2 ,自变量) 用字符串方程表示,自变量缺省值为t。导
22、数用D表示,2阶导数用D2表示,以此类推。S返回解析解。在方程组情形,s为一个符号结构。tout,yout=ode45(yprime,t0,tf,y0) 采用变步长四阶Runge-Kutta法和五阶Runge-Kutta-Felhberg法求数值解,yprime是用以表示f(t,y)的M文件名,t0表示自变量的初始值,tf表示自变量的终值,y0表示初始向量值。输出向量tout表示节点(t0,t1, ,tn)T,输出矩阵yout表示数值解,每一列对应y的一个分量。若无输出参数,则自动作出图形。例:求解微分方程的特解在初始条件下的特解:y(0)=1代码如下:dsolve(x2-1)*Dy+2*x*
23、y-cos(x)=0,y(0)=1,x)运行结果如下: 2.2 矩阵基本计算2.2.1矩阵的最大值求一个向量X的最大值的函数有两种调用格式,分别是:(1) y=max(X):返回向量X的最大值存入y,如果X中包含复数元素,则按模取最大值。(2) y,I=max(X):返回向量X的最大值存入y,最大值的序号存入I,如果X中包含复数元素,则按模取最大值。例:求出矩阵12,52,-25,63,45,1,0,12中的最大值及其所在位置代码如下:x=12,52,-25,63,45,1,0,12;y,l=max(x) %求向量x中的最大值及其该元素的位置运行结果: 2.2.2矩阵的最小值求一个向量X的最小
24、值的函数有两种调用格式,分别是:(1) y=min(X):返回向量X的最小值存入y,如果X中包含复数元素,则按模取最大值。(2) y,I=min (X):返回向量X的最小值存入y,最小值的序号存入I,如果X中包含复数元素,则按模取最小值。例:求出矩阵-5,85,74,42,-2,3,0,-568中的最小值及其所在位置代码如下:x=-5,85,74,42,-2,3,0,-568;y,l=min(x) %求向量x中的最小值及其该元素的位置运行结果: 2.2.3矩阵的均值mean(A,1)表示对列取平均,mean(A,2)表示对行取平均,mean(A)则默认为mean(A,1)例:求45,47,41
25、,42,46,48,42,45,44均值代码如下:A=45,47,41;42,46,48;42,45,44mean(A,1)mean(A,2)运行结果如下:2.2.4矩阵的方差var(a); % 默认来求var(a, 0); % 默认的公式(用N-1)var(a, 1); % 另外的公式(用N)var(a, 0, 1); % 对每列操作var(a, 0, 2); % 对每行操作var(a); % 检验var(a(:); % 通过直接访问矩阵的存储,来对矩阵进行操作例:求方差代码如下:a=rand(3,4) %(0,1)之间产生随机数b=var(a)c=var(a(:)运行结果:2.2.5矩阵的
26、转置 矩阵转置用符号“”来表示和实现, 如故Z是复数矩阵,则Z为它们的复数共轭转置矩阵,非共轭转置矩阵使用Z.或conj(Z)。 例:求33,34,33,36;32,33,31,33;38,37,30,37;34,36,37,34的转置。 代码如下: A=33,34,33,36;32,33,31,33;38,37,30,37;34,36,37,34 B=A 运行结果: 2.2.6矩阵的逆 MATLAB中求矩阵的逆用函数inv(A) 例:求矩阵A=13,14,15,16;12,13,11,13;18,17,10,17;14,16,17,14的逆 代码如下: A=13,14,15,16;12,13
27、,11,13;18,17,10,17;14,16,17,14 B=inv(A) 运行结果: 2.2.7矩阵的行列式 调用函数det格式:d = det(X) 例:求矩阵13,4,1,2;1,2,1,3;8,7,2,7;14,6,7,1的行列式 代码如下: A=13,4,1,2;1,2,1,3;8,7,2,7;14,6,7,1 B=det(A) 运行结果: 2.2.8矩阵的特征值计算 MATLAB中用函数eig()来进行求特征值和特征向量,其调用格式如下:d=eig(A),返回矩阵A的所有特征值。 V,D = eig(A),返回矩阵A 的特征值和特征向量,它们满足如下关系:A*V = V*D。
28、V,D = eig( A ,nobalance) ,在求解特征值和特征向量时不采用初期的平衡步骤。一般来说平衡步骤对输入矩阵进行调整,这使得计算出的特征值和特征向量更加准确。然而如果输入矩阵中确实含有值很小的元素(可能会导致截断误差),平衡步骤有可能加大这种误差,从而得到错误的特征值和特征向量。 d = eig(A,B),返回矩阵A 和B 的广义特征值。 V,D = eig(A,B),返回矩阵A 和B 的广义特征值和广义特征向量。 V,D = eig( A,B,flag),flag 有chol和qz 两种值。当flag =chol 时,计算广义特征值采用B 的Cholesky 分解来实现。当f
29、lag =qz时,无论矩阵的对称性如何,都采用QZ算法来求解广义特征值。 例:求矩阵33,24,15,26;12,23,31,33;38,27,20,37;14,16,22,14的特征值和特征向量 代码如下: A=33,24,15,26;12,23,31,33;38,27,20,37;14,16,22,14 V D=eig(A) eig(A) 运行结果: 2.2.9矩阵的相乘 例:求1,2,3,4;5,6,7,8;9,1,2,3和2;3;6;7的乘积。 代码如下: A=1,2,3,4;5,6,7,8;9,1,2,3 B=2;3;6;7 C=A*B 运行结果: 2.2.10矩阵的右除和左除 MA
30、TLAB中左除和右除分别用“”和“/”表示 例:求11,12,13;14,15,16;17,18,19和15,14,16;17,16,18;12,14,19的左除和右除 代码如下: A=11,12,13;14,15,16;17,18,19; B=15,14,16;17,16,18;12,14,19; C=AB D=A/B 运行结果: 2.2.11矩阵的幂运算 例:求矩阵1,2,3,4;5,6,7,8;7,5,4,3的3次幂 代码如下: A=1,2,3,4;5,6,7,8;7,5,4,3;4,5,6,7 B=3; C=AB 运行结果: 2.3 多项式基本计算2.3.1多项式加减运算 例:计算31
31、,12,13,16,87和5,9,6,6,3的加减乘除 编程如下: A=31,12,13,16,87; B=5,9,6,6,3; C=A+B D=A-B 运行结果: 2.3.2多项式乘除运算 conv(A,B)用于多项式的乘法; Decnov(A,B)用于多项式相除求余式。 例: P1=3,4,4,2,3,4,5 P2=6,7,8 P=conv(P1,P2) P3=deconv(P,P1) 运行结果: 2.3.3多项式求导 polyder(P) %对多项式P进行微分运算 例: 代码如下 P=1,2,3,4,5 P1=poly2sym(P) Q=polyder(P) 运行结果如下: 2.3.4求
32、根和求值运算 polyval(P) %按数组计算规则计算多项式P的值 polyvalm(P) %按矩阵计算规则计算多项式P的值 例:求多项式在x=2.5时的值。 代码如下: P=1,4,4 Pv=polyval(P,0.5) 运行结果: 2.3.5多项式的部分分式展开 调用函数residue格式:r,p,k = residue(b,a),b,a = residue(r,p,k)其中,b与a分别为两个多项式。其中a为分式之分母,b该分式之分子,两者均为向量型式;而r,p,k等则为行向量,分别代表展开后之分子、分母及余数之系数,若能除尽则k项应为零。r与p分别为余数与极数,其个数应相等,但应比a之
33、个数少一。 代码如下: b=3,4,5; a=8,-6,3; r,p,k=residue(b,a) 运行结果: 2.3.6多项式的拟合 MATLAB软件提供了基本的曲线拟合函数的命令多项式函数拟合:a=polyfit(xdata,ydata,n)其中n表示多项式的最高阶数,xdata,ydata为将要拟合的数据,它是用数组的方式输入输出参数a为拟合多项式 y=a1xn+.+anx+a n+1的系数 例:已知五个数据点:1,5.5,2,4.3,3,12.8,4,2.9,5,4.98,试画出这五个点拟合的三次曲线。 代码如下: x=1,2,3,4,5;y=5.5,4.3,1.28,0.29,0.4
34、98; p=polyfit(x,y,3) x2=1:0.1:5;y2=polyval(p,x2); plot(x,y,o,x2,y2) ,grid on 运行结果: 2.3.7插值运算 yi = interp1(x,Y,xi,method) %用指定的算法计算插值: nearest:最近邻点插值,直接完成计算; linear:线性插值(缺省方式),直接完成计算; spline:三次样条函数插值。 cubic: 分段三次Hermite插值。 对于超出x范围的xi的分量,使用方法nearest、linear、v5cubic的插值算法,相应地将返回NaN。对其他的方法,interp1将对超出的分量执
35、行外插值算法。 yi = interp1(x,Y,xi,method,extrap) yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围的xi中的分量的外插值extrapval,其值通常取NaN或0。 ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值: linear:双线性插值算法(缺省算法); nearest:最临近插值; spline:三次样条插值; cubic:双三次插值。 例:用二维三次插值函数模拟函数z(x,y)=sin(x)+cos(y)。 代码如下: clear; clc; x=li
36、nspace(0,2*pi,6); y=linspace(0,2*pi,6); X,Y=meshgrid(x,y); Z=cos(X)+sin(Y); surf(X,Y,Z) x1=linspace(0,2*pi,36); y1=x1; X1,Y1=meshgrid(x1,y1); Z1=interp2(X,Y,Z,X1,Y1,cubic); figure(2) surf(X1,Y1,Z1) 运行结果:3.基于MATLAB的图像滤波设计3.1读入图像并分别加入高斯噪声、椒盐噪声和乘性噪声,并比较结果 函数使用imnoise函数进行图片的噪声加入,其调用格式如下:J=imnoise(I,type
37、,parameters)其中,type是噪声的类型,有高斯噪声,椒盐噪声,乘性噪声。类型名分别为:guassian,salt & pepper,speckle,parameters可表示噪声密度D。 代码如下:clear;I=imread(C:Documents and SettingsAdministrator桌面花.jpg);subplot(2,2,1);imshow(I);title(原始图像);J2=imnoise(I,gaussian,0.1);subplot(2,2,3)imshow(J2)title(高斯噪声图像);J1=imnoise(I,salt & pepper,0.1);
38、subplot(2,2,2)imshow(J1)title(椒盐噪声图像);J3=imnoise(I,speckle,0.1);subplot(2,2,4)imshow(J3) title(乘性噪声图像); 运行结果,对比如下: 噪声密度D不同对图像干扰程度对比,程序如下: 高斯噪声:clear;I=imread(C:Documents and SettingsAdministrator桌面花.jpg);J1=imnoise(I,gaussian,0.03);subplot(1,3,1);imshow(J1);title(d=0.03时的图像);J2=imnoise(I,gaussian,0.
39、15);subplot(1,3,2)imshow(J2)title(d=0.15时的图像);J3=imnoise(I,gaussian,0.5);subplot(1,3,3)imshow(J3)title(d=0.5时的图像); 运行结果对比如下:椒盐噪声:clear;I=imread(C:Documents and SettingsAdministrator桌面花.jpg);J1=imnoise(I,salt & pepper,0.03);subplot(1,3,1);imshow(J1);title(d=0.03时的图像);J2=imnoise(I,salt & pepper,0.15);
40、subplot(1,3,2)imshow(J2)title(d=0.15时的图像);J3=imnoise(I,salt & pepper,0.5);subplot(1,3,3)imshow(J3)title(d=0.5时的图像); 运行结果对比如下:乘法噪声:I=imread(C:Documents and SettingsAdministrator桌面花.jpg);J1=imnoise(I,speckle,0.03);subplot(1,3,1);imshow(J1);title(d=0.03时的图像);J2=imnoise(I,speckle,0.15);subplot(1,3,2)ims
41、how(J2)title(d=0.15时的图像);J3=imnoise(I,speckle,0.5);subplot(1,3,3)imshow(J3)title(d=0.5时的图像); 运行结果对比如下:3.2设计巴特沃斯低通滤波对图像进行低通滤波处理,显示结果巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波得图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。3.2.1叠加椒盐噪声的巴特沃斯低通滤波 i=imread(C:Documents and SettingsAdministrator桌面花.
42、jpg); I=rgb2gray(i);I1=imnoise(I,salt & pepper,0.02);f=double(I1);g=fft2(f);g=fftshift(g);N1,N2=size(g);n=3; %阶次设为3d0=100; %此处d0为截止频率n1=fix(N1/2);n2=fix(N2/2);for i=1:N1 for j=1:N2 d=sqrt(i-n1)2+(j-n2)2); h=1/(1+0.414*(d/d0)(2*n); result(i,j)=h*g(i,j);end endresult=ifftshift(result);X2=ifft2(result);J1=uint8(real(X2);subplot(1,2,1),imshow(I1);title(受椒盐噪声污染的图像);subplot(1,2,2),imshow(J1);title(截止频率为100HZ的巴特沃斯低通滤波处