摘 要
1.MATLAB简介
1.1 MATLAB的概况
1.2 MATLAB产生的历史背景
2.编程及运行结果
2.1常见基本运算
2.1.1极限的计算
2.1.2微分的计算
2.1.3积分的计算
2.1.4级数的计算
2.1.5求解代数方程
2.1.6求解常微分方程
2.2 矩阵基本计算
2.2.1矩阵的最大值
2.2.2矩阵的最小值
2.2.3矩阵的均值
2.2.4矩阵的方差
2.2.5矩阵的转置
2.2.6矩阵的逆
2.2.7矩阵的行列式
2.2.8矩阵的特征值计算
2.2.9矩阵的相乘
2.2.10矩阵的右除和左除
2.2.11矩阵的幂运算
2.3 多项式基本计算
2.3.1多项式加减运算
2.3.2多项式乘除运算
2.3.3多项式求导
2.3.4求根和求值运算
2.3.5多项式的部分分式展开
2.3.6多项式的拟合
2.3.7插值运算
3.基于MATLAB的图像滤波设计
3.1读入图像并分别加入高斯噪声、椒盐噪声和乘性噪声，并比较结果
3.2设计巴特沃斯低通滤波对图像进行低通滤波处理，显示结果
3.2.1叠加椒盐噪声的巴特沃斯低通滤波
3.2.2叠加高斯噪声的巴特沃斯低通滤波
3.2.3叠加乘性噪声的巴特沃斯低通滤波
3.3用MATLAB实现高斯高通滤波器对图像的处理
3.4维纳滤波和中值滤波对图像进行处理
4.总结
参考文献
摘 要
现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因，使得具有线性相位的FIR数字滤波器得到大力发展和广泛应用。
在实际进行数字信号处理时，往往需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进行分析。取用有限个数据，即将信号数据截断的过程，就等于将信号进行加窗函数操作。这样操作以后，常常会发生频谱分量从其正常频谱扩展开来的现象，即所谓的“频谱泄漏”。当进行离散傅立叶变换时，时域中的截断是必需的，因此泄漏效应也是离散傅立叶变换所固有的，必须进行抑制。而要对频谱泄漏进行抑制，可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。
数字带通滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位重启响应函数的时域特性可分为两类：无限冲击响应滤波器(IIR)，有限冲击响应滤波器(FIR)。与IIR滤波器相比，FIR的实现是递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如信号音频，图像处理，数据传输等领域得到广泛的应用。
数字fir滤波器的设计方法有很多种。如窗函数法设计，频率采样设计法和最优化设计法等。
Abstract
Modern 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 signal 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 from 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, 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 signal 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 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 as window function method design, frequency sampling design method and the optimum design method, etc.
1.MATLAB简介
1.1 MATLAB的概况
MATLAB是矩阵实验室（Matrix Laboratory）之意。除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学工程中常用的形式十分相似,故用MATLAB来解算问题要比用C、FORTRAN等语言完相同的事情简捷得多。
当前流行的MATLAB 5.3/Simulink 3.0包括拥有数百个内部函数的主包和三十几种工具包(Toolbox).工具包又可以分为功能性工具包和学科工具包.功能工具包用来扩充MATLAB的符号计算,可视化建模仿真,文字处理及实时控制等功能.学科工具包是专业性比较强的工具包,控制工具包,信号处理工具包,通信工具包等都属于此类。
开放性使MATLAB广受用户欢迎.除内部函数外,所有MATLAB主包文件和各种工具包都是可读可修改的文件,用户通过对源程序的修改或加入自己编写程序构造新的专用工具包。
1.2 MATLAB产生的历史背景
在70年代中期,Cleve Moler博士和其同事在美国国家科学基金的资助下开发了调用EISPACK和LINPACK的FORTRAN子程序库.EISPACK是特征值求解的FOETRAN程序库,LINPACK是解线性方程的程序库.在当时,这两个程序库代表矩阵运算的最高水平。
到70年代后期,身为美国New Mexico大学计算机系系主任的Cleve Moler,在给学生讲授线性代数课程时,想教学生使用EISPACK和LINPACK程序库,但他发现学生用FORTRAN编写接口程序很费时间,于是他开始自己动手,利用业余时间为学生编写EISPACK和LINPACK的接口程序。Cleve Moler给这个接口程序取名为MATLAB,该名为矩阵(matrix)和实验室(labotatory)两个英文单词的前三个字母的组合.在以后的数年里,MATLAB在多所大学里作为教学辅助软件使用,并作为面向大众的免费软件广为流传。
在当今30多个数学类科技应用软件中,就软件数学处理的原始内核而言,可分为两大类.一类是数值计算型软件,如MATLAB,Xmath,Gauss等,这类软件长于数值计算,对处理大批数据效率高；另一类是数学分析型软件,Mathematica,Maple等,这类软件以符号计算见长,能给出解析解和任意精确解,其缺点是处理大量数据时效率较低。MathWorks公司顺应多功能需求之潮流,在其卓越数值计算和图示能力的基础上,又率先在专业水平上开拓了其符号计算,文字处理,可视化建模和实时控制能力,开发了适合多学科,多部门要求的新一代科技应用软件MATLAB.经过多年的国际竞争,MATLAB以经占据了数值软件市场的主导地位。
在MATLAB进入市场前，国际上的许多软件包都是直接以FORTRANC语言等编程语言开发的。这种软件的缺点是使用面窄，接口简陋，程序结构不开放以及没有标准的基库，很难适应各学科的最新发展，因而很难推广。MATLAB的出现，为各国科学家开发学科软件提供了新的基础。在MATLAB问世不久的80年代中期，原先控制领域里的一些软件包纷纷被淘汰或在MATLAB上重建。
时至今日，经过MathWorks公司的不断完善，MATLAB已经发展成为适合多学科，多种工作平台的功能强大大大型软件。在国外，MATLAB已经经受了多年考验。在欧美等高校，MATLAB已经成为线性代数，自动控制理论，数理统计，数字信号处理，时间序列分析，动态系统仿真等高级课程的基本教学工具；成为攻读学位的大学生，硕士生，博士生必须掌握的基本技能。在设计研究单位和工业部门，MATLAB被广泛用于科学研究和解决各种具体问题。在国内，特别是工程界，MATLAB一定会盛行起来。可以说，无论你从事工程方面的哪个学科，都能在MATLAB里找到合适的功能。
2.编程及运行结果
2.1常见基本运算
2.1.1极限的计算
MATLAB提供的命令limit()可以完成极限运算，其调用格式如下: limit(F,x,a,'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 x
y=log(1/x);
dy=diff(y,x)
运行结果：
2.1.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的上下限，c,d分别为y的上下线。
例1：（不定积分）用符号积分命令int计算
代码如下：
clear;
syms x
int(x^2*sin(x))
运行结果：
例2：（定积分）计算数值积分
①用符号积分命令int，代码如下：
clear;
syms x;
int(x^4,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点展开为n-1阶Taylor多项式
例：用Symsum计算
代码如下：
clear;clc;
syms x y
z=1/(x^3);
symsum(z,x,1,9)
运行结果：
2.1.5求解代数方程
MATLAB命令输人格式：solve('equ1','equ2',...'equN')，其中eqni表示第i个方程。
例1：求解方程 x^2+b*x+c=0
代码如下：
solve('a*x^2+b*x+c')
运行结果：
2.1.6求解常微分方程
MATLAB中主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。s=dsolve(‘方程1', ‘方程2',…,'初始条件1'，'初始条件2' …,'自变量') 用字符串方程表示，自变量缺省值为t。导数用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('(x^2-1)*Dy+2*x*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的最小值的函数有两种调用格式，分别是：
(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,42,46,48,42,45,44]均值
代码如下：
A=[45,47,41;42,46,48;42,45,44]
mean(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矩阵的转置
矩阵转置用符号“`”来表示和实现, 如故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,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。
[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 分解来实现。当flag ='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矩阵的右除和左除
MATLAB中左除和右除分别用“\”和“/”表示
例：求[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=A\B
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=A^B
运行结果：
2.3 多项式基本计算
2.3.1多项式加减运算
例：计算[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求根和求值运算
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之个数少一。
代码如下：
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.498];
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将对超出的分量执行外插值算法。
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=linspace(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,parameters)其中，type是噪声的类型，有高斯噪声，椒盐噪声，乘性噪声。类型名分别为：guassian,salt & pepper,speckle，parameters可表示噪声密度D。
代码如下：
clear;
I=imread('C:\Documents and Settings\Administrator\桌面\花.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);
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 Settings\Administrator\桌面\花.jpg');
J1=imnoise(I,'gaussian',0.03);
subplot(1,3,1);
imshow(J1);
title('d=0.03时的图像');
J2=imnoise(I,'gaussian',0.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 Settings\Administrator\桌面\花.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);
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 Settings\Administrator\桌面\花.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)
imshow(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 Settings\Administrator\桌面\花.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; %阶次设为3
d0=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
end
result=ifftshift(result);
X2=ifft2(result);
J1=uint8(real(X2));
subplot(1,2,1),imshow(I1);
title('受椒盐噪声污染的图像');
subplot(1,2,2),imshow(J1);
title('截止频率为100HZ的巴特沃斯低通滤波处理后');
运行结果：
3.2.2叠加高斯噪声的巴特沃斯低通滤波
clear;
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg'); I=rgb2gray(i);
I1=imnoise(I,'gaussian',0.02);
f=double(I1);
g=fft2(f);
g=fftshift(g);
[N1,N2]=size(g);
n=3; %阶次设为3
d0=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
end
result=ifftshift(result);
X2=ifft2(result);
J1=uint8(real(X2));
subplot(1,2,1),imshow(I1);
title('受高斯噪声污染的图像');
subplot(1,2,2),imshow(J1);
title('截止频率为100HZ的巴特沃斯低通滤波处理后');
运行结果：
3.2.3叠加乘性噪声的巴特沃斯低通滤波
clear;
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg'); I=rgb2gray(i);
I1=imnoise(I,'speckle',0.02);
f=double(I1);
g=fft2(f);
g=fftshift(g);
[N1,N2]=size(g);
n=3; %阶次设为3
d0=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
end
result=ifftshift(result);
X2=ifft2(result);
J1=uint8(real(X2));
subplot(1,2,1),imshow(I1);
title('受乘性噪声污染的图像');
subplot(1,2,2),imshow(J1);
title('截止频率为100HZ的巴特沃斯低通滤波处理后');
运行结果：
3.3用MATLAB实现高斯高通滤波器对图像的处理
I1=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
subplot(1,2,1);
imshow(I1);
title('原图');
f=double(I1);
g=fft2(f);
g=fftshift(g);
[M,N]=size(g);
d0=100;
m=fix(M/2); n=fix(N/2);
result=g;
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
h=exp(-(d.^2)./(2*(d0^2)));
result(i,j)=(1-h)*g(i,j);
end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(1,2,2);
imshow(J2);
title('高斯高通滤波(d0=100)');
运行结果：
3.4维纳滤波和中值滤波对图像进行处理
随机信号或随机过程(random process)是普遍存在的。一方面，任何确定性信号经过测量后往往就会引入随机性误差而使该信号随机化；另一方面，任何信号本身都存在随机干扰，通常把对信号或系统功能起干扰作用的随机信号称之为噪声。噪声按功率谱密度划分可以分为白噪声(white noise)和色噪声(color noise)，我们把均值为0的白噪声叫纯随机信号(pure random signal)。因此，任何其它随机信号都可看成是纯随机信号与确定性信号并存的混合随机信号或简称为随机信号。要区别干扰(interference)和噪声( noise)两种事实和两个概念。非目标信号(nonobjective signal)都可叫干扰。干扰可以是确定信号，如国内的50Hz工频干扰。干扰也可以是噪声，纯随机信号(白噪声)加上一个直流成分(确定性信号)，就成了最简单的混合随机信号。医学数字信号处理的目的是要提取包含在随机信号中的确定成分，并探求它与生理、病理过程的关系，为医学决策提供一定的依据。例如从自发脑电中提取诱发脑电信号，就是把自发脑电看成是干扰信号，从中提取出需要的信息成分。因此我们需要寻找一种最佳线性滤波器，当信号和干扰以及随机噪声同时输入该滤波器时，在输出端能将信号尽可能精确地表现出来。
维纳滤波和卡尔曼滤波就是用来解决这样一类问题的方法：从噪声中提取出有用的信号。实际上，这种线性滤波方法也被看成是一种估计问题或者线性预测问题。对于在C中的图片而言，其进行均值滤波和维纳滤波的MATLAB的命令如下：
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
I=rgb2gray(i);
J1=imnoise(I,'gaussian',0,0.1); %高斯噪声 均值0 方差为0.1
J2=imnoise(I,'salt & pepper',0.1); %椒盐噪声 均值0 方差为0.1
J3=imnoise(I,'speckle',0.1); %乘性噪声 均值0 方差为0.1
K1=medfilt2(J1,[3,3]); %3*3的滤波窗
subplot(2,3,1),imshow(J1),title('高斯噪声');
subplot(2,3,4),imshow(K1,[]),title('对高斯噪声图像中值滤波');
K2=medfilt2(J2,[3,3]);
subplot(2,3,2),imshow(J2),title('椒盐噪声');
subplot(2,3,5),imshow(K2,[]),title('对椒盐噪声图像中值滤波');
K3=medfilt2(J3,[3,3]);
subplot(2,3,3),imshow(J3),title('乘性噪声');
subplot(2,3,6),imshow(K3,[]),title('对乘性噪声图像中值滤波');
运行结果：
以前面经常用的图片为例，再进行维纳滤波：
代码如下：
i=imread('C:\Documents and Settings\Administrator\桌面\花.jpg');
I=rgb2gray(i);
subplot(2,4,1),imshow(I); title('原始图像');
J1=imnoise(I,'gaussian',0,0.005);
subplot(2,4,2),imshow(J1);
title('高斯噪声图像');
K1=wiener2(J1);
subplot(2,4,6),imshow(K1);
title('维纳滤波');
J2=imnoise(I,'salt & pepper',0.005);
subplot(2,4,3),imshow(J2);
title('椒盐噪声图像');
K2=wiener2(J2);
subplot(2,4,7),imshow(K2);
title('维纳滤波');
J3=imnoise(I,'speckle',0.005);
subplot(2,4,4),imshow(J3);
title('乘性噪声图像');
K3=wiener2(J3);
subplot(2,4,8),imshow(K3);
title('维纳滤波');
运行结果：
参考文献
[1] 刘泉 《数字信号处理原理与实现》 电子工业出版社
[2] 刘泉 《信号与系统》 高等教育出版社
[3] 周博等编《MATLAB科学计算》 机械工业出版社
[4] 贺兴华等著.《MATLAB7.X图像处理》 人民邮电出版社
[5] 王洪元主编.《MATLAB语言以及在电子信息工程中的应用》 清华出版社.
[6] 张汗灵. 《MATLAB在图像处理中的应用》 清华大学出版社
[7] 董霖编 《MATLAB使用详解——基础、开发及工程应用》电子工业出版社
[8] Edward W.Kamen.《应用Web和MATLAB的信号与系统基础》.电子工业出版社
精彩内容
发表评论