Matlab课程设计任务书
学生姓名：专业班级：
指导教师：工作单位：武汉理工大学
题目: 基于Matlab的二维信号滤波设计
基础强化训练目的
1．理论目的
《数字信号处理》课程设计的目的之一是为了巩固课堂理论学习，并能用所学理论知识正确分析数字信号处理的基本问题和解释数字信号处理的基本现象。
2．实践目的
《数字信号处理》课程设计的目的之二是通过设计具体的各种滤波器掌握滤波器设计方法和步骤。
训练内容和要求
利用MATLAB－图像处理工具箱，（1）设计中值滤波器（medfilt2）对二维信号（数字图像）滤波器，并且显示原图像和滤波结果。（2）对二维信号（数字图像）进行FFT，然后进行低通滤波（即去掉高频分量），然后再反变换，得到滤波图像，然后反过来，进行高通滤波。
初始条件
1．MATLAB软件
2．数字信号处理与图像处理基础知识
时间安排：
第18周，安排任务
第18周，仿真设计
第18周，完成
指导教师签名：年月日
系主任（或责任教师）签名：年月日
目录
摘要................................................................ III Abstract............................................................... IV 1二维信号数字图像数据采集. (1)
1.1图像的选取 (1)
1.2 MATLAB读取方法 (1)
1.2.1 编辑M文件 (1)
1.2.2 图像的读取 (2)
1.2.3查看图像的格式 (3)
1.2.4 灰度值的获取 (4)
2 中值滤波处理 (5)
2.1 中值滤波 (5)
2.1.1中值滤波定义及原理 (5)
2.1.2 中值滤波的实现方法 (5)
2.2 中值滤波器（medfilt2）对二维信号（数字图像）滤波 (5)
3快速傅立叶变化（FFT） (7)
3.1快速傅立叶变化（FFT）的概念 (7)
3.1.1快速傅立叶变化（FFT）原理及计算公式 (7)
3.1.2函数应用 (8)
3.1.3离散傅立叶逆变换（IFFT） (9)
3.2二维信号（数字图像）的FFT及其逆变换 (10)
3.3对变换后的图像的说明 (11)
4 低通滤波处理 (12)
4.1 低通滤波器介绍 (12)
4.1.1 巴特沃斯滤波器 (12)
4.1.2 切比雪夫滤波器 (13)
4.2 二维信号（数字图像）的低通滤波处理 (13)
4.3 结果分析 (16)
5 高通滤波处理 (17)
5.1 高通滤波器介绍 (17)
5.2 二维信号（数字图像）的高通滤波处理 (17)
5.3 结果分析 (19)
6 总结与体会 (20)
7 参考文献 (21)
摘要
MATLAB 是一种直观、高效的计算机语言,同时又是一个科学计算平台。它为数据分析和数据可视化、算法和应用程序开发提供了最核心的数学和高级图形工具。根据它提供的500 多个数学和工程函数，工程技术人员和科学工作者可以在它的集成环境中交互或编程以完成各自的计算。
本次课程设计主要是利用MATLAB－图像处理工具箱，对二维信号（数字图像）进行快速傅立叶变换与反变换，并且设计中值滤波器、低通滤波器和高通滤波器对二维信号（数字图像）进行滤波。
关键词：MATLAB 二维数字图像处理数字滤波器中值滤波器低通滤波高通滤波
Abstract
MATLAB is an intuitive, efficient computer language, is also a scientific computing platform. It is a data analysis and data visualization, algorithm and application development provides the core mathematics and advanced graphics tools. According to it, providing more than 500 mathematics and engineering functions, engineering and technical personnel and scientific workers in its integrated environment interaction or programming to complete their calculation The curriculum design, mainly by using MATLAB-Image Processing Toolbox, for two-dimensional signal (digital images) for fast Fourier transform and inverse transform, and the design of median filter, low-pass filters and high-pass filter on the two-dimensional signal (digital image) for filtering.
Key words:MATLAB two-dimensional digital image processing digital filter Median Filter Low-pass filter high-pass filter
1二维信号数字图像数据采集
1.1图像的选取
选取一张图片，如图1所示，满足像素要求，但图像为三维彩图，不是灰度图，在下面的过程中会将彩图转换为灰度图。
图 1 选取的需要处理的原图
1.2 MATLAB读取方法
1.2.1 编辑M文件
MATLAB是解释型语言，就是说MATLAB命令行中敲入的命令在当前MATLAB进程中被解释运行。有两种方法可以使MATLAB的力量得到扩展——脚本和函数。这两种方法都用像emacs一样的文本编辑器中编写的m文件（因为扩展名是.m所以这样命名，m 文件还称点m文件）。m文件的好处在于它可以保存命令，还可以轻易地修改命令而无需重新敲入整个命令行。
新建m文件的方法为File->New->M-File
图 2 新建M文件的方法图
1.2.2 图像的读取
利用函数imread()可完成图形图像文件的读取，语法：
A=imread(filename,fmt)
[X,map]=imread(filename,fmt)
[...]=imread(filename)
[...]=imread(filename,idx) （只对TIF格式的文件）
[...]=imread(filename,ref) （只对HDF格式的文件）
通常，读取的大多数图像均为8bit，当这些图像加载到内存中时，Matlab就将其存放在类uint8中。此为Matlab还支持16bit的PNG和TIF图像，当读取这类文件时，Matlab 就将其存贮在uint16中。对于索引图像，即使图像阵列的本身为类uint8或类uint16，imread 函数仍将颜色映象表读取并存贮到一个双精度的浮点类型的阵列中。
MATLAB代码为：
P=imread('F:\猫咪.jpg'); %读取图像路径
imshow(P);%显示图像
title('原图');
图 3 MATLAB中读取图像的效果图
1.2.3查看图像的格式
使用MATLAB 命令isgray（）判断图像是否为灰度图，如果是灰度图将返回为1，不是将返回0。
我所选图的图片用isgray（）判断后返回值为0，证明不是灰度图，不符合题目的要求，所以要使用rgb2gray（）命令实现RGB图像到灰度图像的转化，转化后在使用命令isgray （）判断返回值为1，证明已经转换成灰度图了。
MATLAB命令：f=rgb2gray(P);
isgray(f)
再用MATLAB命令：figure
imshow(f) ;
title('灰白图像')
显示原图变成的灰度图像如图4：
图 4 原图变为灰度图像
1.2.4 灰度值的获取
直接在MATLAB命令行中输入存储图像矩阵的名字P即可得到灰度矩阵：由于矩阵较大，因此只选取前25行如图5所示：
图 5 灰度矩阵图
2 中值滤波处理
2.1 中值滤波
2.1.1中值滤波定义及原理
1）中值滤波法是一种非线性平滑技术，它将每一象素点的灰度值设置为该点某邻域窗口内的所有象素点灰度值的中值.
2）中值滤波原理
中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个拎域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。方法是去某种结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为g（x,y）=med{f(x-k,y-l),(k,l∈W)} ，其中,f(x,y)，g(x,y)分别为原始图像和处理后图像。W为二维模板，通常为2*2，3*3区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。
2.1.2 中值滤波的实现方法
1）通过从图像中的某个采样窗口取出奇数个数据进行排序；
2）用排序后的中值取代要处理的数据即可。
中值滤波法对消除椒盐噪音非常有效，在光学测量条纹图象的相位分析处理方法中有特殊作用，但在条纹中心分析方法中作用不大.
中值滤波在图像处理中,常用于用来保护边缘信息,是经典的平滑噪声的方法
2.2 中值滤波器（medfilt2）对二维信号（数字图像）滤波
利用函数medfilt2可以完成对二维信号（数字图像）的滤波
MATLAB代码为：
A=imread('F:\猫咪.jpg');
B=rgb2gray(A);
C=medfilt2(B);
figure
imshow(B);
figure
imshow(C);
title('中值滤波后');
显示图像为灰白图像和中值滤波后图像，分别如图6和图7所示。
图6 灰白图像显示
图7 中值滤波后图像显示
3快速傅立叶变化（FFT ）
在信号处理中，DFT 的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要
通过DFT 来实现。然而，当N 很大的时候，求一个N 点的DFT 要完成N*N 次复数乘法和N （N —1）次复数加法，其计算量相当大。1965年J.W.Cooley 和J.W.Tukey 巧妙地利用WN 因子的周期性和对称性，构造了一个DFT 快速算法，即快速傅立叶变换(FFT)。
3.1快速傅立叶变化（FFT ）的概念
3.1.1快速傅立叶变化（FFT ）原理及计算公式
有限序列长为N 的序列x(n)的DFT 变换为
nk
N n N
W n x k X ∑-==10)()( k=0，1，…，N-1
其逆变换为
nk N k N W
n x N n x --=∑=10)(1)( n=0，1，…，N-1
由于MATLAB 软件本身的特点，序列或向量元素下标从1开始记录，而不是从0开始。
因此，上述两式在MATLAB 中相对应的表达式为
nk N n N W
n x k X ∑-==10
)()( k=1，2，…，N-1 nk
N k N W
k X N k X -=∑=1)(1)( n=1，2，…，N
而下面所讨论使用的快速傅立叶变换（FFT ）并不是与DFT 不同的另外一种变换，而
是为减少DFT 计算次数的一种快速有效的算法。这种快速有效算法，主要是利用了nk N W 下面两个特性使长序列的DFT 分解为更小点数的DFT 所实现的。
(1)利用nk
N W 的对称性使DFT 运算中有些项合并
*)()(kn N kn N k n N k N W W W ==-- (2)利用nk
N W 的周期性和对称性使长序列的DFT 分解为更小点数的DFT
n N k N nk N nk N W W W )(+== 快速傅立叶变换（FFT ）算法正是基于这一基本思想而发展起来的。快速傅立叶变换算法
形式很多，但是基本上可以分为两大类，即按时间抽取（Decimation-IN-Time,简称DIT）法和按频率抽样（Decimation-IN-Frequency）法。
3.1.2函数应用
MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有Fft、Ifft，Fft2、Ifft2，Fftn、Ifftn和Fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方法来添补数据使之成为2的幂次。
1）fft和Ifft函数
调用方式
a．Y=fft（X）
参数说明
如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换；
如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；
如果X是（N*D）维数组，则是对每一个非单元素的维进行离散傅立叶变换。
b．Y=fft（X，N）
参数说明
N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。
c．Y=fft（X，[]，dim）或Y=fft（X，N,dim）
参数说明
在参数dim指定的维上进行离散傅立叶变换；
当X为矩阵时，dim用来指定变换的实施方向：dim=1，表明变换按列进行；dim=2，表明变换按行进行。
函数Ifft的参数应用于函数fft完全相同。
2）fft2和Ifft2函数
调用方式
a．Y=fft2（X）
参数说明
如果X是向量，则此傅立叶变换即变成一维傅立叶变换fft；
如果X是矩阵，则是计算该矩阵的二维快速傅立叶变换；
数据二维傅立叶变换fft2（X ）相当于fft （fft （X ）`）`,即先对X 的列做一维傅立叶变换，然后对变换结果的行做一维傅立叶变换。
b ．fft2（X,M,N ）
通过对X 进行补零或截取,使得X 成为（M*N ）的矩阵。
函数Ifft2的参数应用于函数Fft2完全相同。
fftn 、Ifftn 是对数据进行多维快速傅立叶变换，其应用与fft2、Ifft2类似。
3） fftshift 和Ifftshift 函数
调用方式
a ．Z=fftshift （Y ）
此函数可用于将傅立叶变换结果Y （频域数据）中的直流成分（即频率为0处的值）移到频谱的中间位置。
参数说明
如果Y 是向量，则交换Y 的左右两边；
如果Y 是矩阵，则交换Y 的一、三象限和二、四象限；
如果Y 是多维数组，则在数组的每一维交换其“半空间” 。
b ．Z=fftshift （Y ，dim ）
此函数将对矩阵的dim 维进行fftshift 操作。
函数Ifftshift 的参数应用与函数fftshift 完全相同。
3.1.3离散傅立叶逆变换（IFFT ）
逆变换IFFT 为
nk N k N W
n x N n x --=∑=10)(1)( n=0，1，…，N-1
通过下列修改，就可以用FFT 算法来实现逆FFT 运算。
增加一个归一化因子1/N ；
将nk N W 用其附属共轭kn N W -代替；
但是又因为，
***10*)]}([{1])([1)(k X FFT N
W k X N n x N k nk N ==∑-=
所以，求X(k)的逆FFT可以分为以下3个步骤：
a 取X(k)的共轭，得到)
X；
(*k
b 求)
X的FFT得N*x(k)；
(*k
c 取*x(n)的共轭，并除以N，即得x(n)。
采用这种方法可以直接用FFT程序来计算逆FFT,有关IFFT的具体应用，与FFT一一对应。
3.2二维信号（数字图像）的FFT及其逆变换
MATLAB程序：
A=imread('F:\猫咪.jpg');
I=rgb2gray(A);
C= fft2(double(I)); %对图像进行傅立叶变换
D=fftshift(fft2(double(I))); %将直流分量移到频谱图的中心
E=ifft2(D); %傅立叶反变换
figure, %画图像原始图
imshow(I);
title('傅里叶变换原图');
figure, %画直接变换频谱图
imshow(log(abs(D)+1),[]);
title('直接变换频谱图');
figure, %显示反傅立叶变换后的图像
imshow(abs(E),[]) ;
title('反傅里叶变换图');
图8 傅里叶变换原图 图9 直接变换频谱图
图10 反傅里叶变换图
3.3对变换后的图像的说明
在时域处理图像并不能得到很好的结果，我们需要把图像转化到频域空间去处理。在
频域空间里，Fourier 变换具有很好的性质，他使我们能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪声、显示点等的作用。本图经过傅里叶变换之后，能够很方便的去除噪音，使图像更清晰。当然傅里叶变换还有很多其他功能，要在后面更加深入的运用时才能体现出来。
4 低通滤波处理
4.1 低通滤波器介绍
低通滤波器是容许低于截至频率的信号通过，但高于截止频率的信号不能通过的
电子滤波装置。对于不同滤波器而言，每个频率的信号的减弱程度不同。当使用在音频应用时，它有时被称为高频剪切滤波器, 或高音消除滤波器。
低通滤波器概念有许多不同的形式，其中包括电子线路（如音频设备中使用的hiss 滤波器、平滑数据的数字算法、音障（acoustic barriers ）、图像模糊处理等等，这两个工具都通过剔除短期波动、保留长期发展趋势提供了信号的平滑形式。
低通滤波器有很多种，其中，最通用的就是巴特沃斯滤波器和切比雪夫滤波器。
4.1.1 巴特沃斯滤波器
巴特沃斯滤波器是电子滤波器的一种。巴特沃斯滤波器的特点是通频带的频率响
应曲线最平滑。这种滤波器最先由英国工程师斯替芬·巴特沃斯（Stephen Butterworth ）在1930年发表在英国《无线电工程》期刊的一篇论文中提出的。
巴特沃斯滤波器，又被称为“最平”的幅频响应滤波器。这是因为，该滤波器在
通带内具有最大平坦的幅度特性，而且随着频率升高呈现出出单调减小的特点。
N 阶低通巴特沃斯滤波器的特性函数为
N c
j j j G 22)(11|)(|ΩΩ+=Ω 其中，c Ω为通带宽度，即截止频率。当阶次N 增大时，滤波器的特性曲线变得
更加陡峭，其特性就是越接近于理想的矩形幅频特性2|)(|s G 。巴特沃斯滤波器属于全极点设计，它的极点由下式来决定。
N c N c
j s j s j G 2222)(11)(11|)(|Ω+=Ω-+=Ω 其中，N c P j S 21
)1(-Ω=。所以，在s 平面上有2N 个极点。这些极点是等间隔地分
布在半径为c Ω的圆周上，并且这些极点都是成复共轭对出现的，极点位置与虚轴是对称的，但不在虚轴上。
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波得图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。
一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。
4.1.2 切比雪夫滤波器
切比雪夫滤波器，是在通带或阻带上频率响应幅度等波纹波动的滤波器。切比雪夫滤波器也是滤波器的一种设计分类，其采用的是切比雪夫传递函数，也有高通、低通、带通、高通、带阻等多种滤波器类型。
巴特沃斯滤波器的频率特性曲线，在通带和阻带内都是单调的。因此当在通带的边界处满足性能指标时，通频带内肯定会有余量。所以更有效的设计方法是将精度均匀地分布在整个通带或阻带内，或者同时分布在两者之间。这可通过选择具有等波纹特性的逼近函数来实现，切比雪夫滤波器的振幅特性就具有这种等波纹特性。
根据频率响应曲线波动位置的不同，切比雪夫滤波器可以分为以下两种：一种是在通带（或称“通频带”）上频率响应幅度等波纹波动的滤波器称为“I型切比雪夫滤波器”，它在阻带内是单调的；另一种是在阻带（或称“阻频带”）上频率响应幅度等波纹波动的滤波器称为“II型切比雪夫滤波器”，它在通带内是单调的。
切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。
4.2 二维信号（数字图像）的低通滤波处理
先对二维信号（数字图像）进行FFT，然后进行低通滤波（即去掉高频分量），然后再反变换，得到滤波图像。
MATLAB代码为：
pic1=imread('F:\猫咪.jpg'); %读取原始图像
whos %原始图像信息
imshow(pic1); %显示原始图像
title('原始图像');
M=rgb2gray(pic1); %转换为灰白图像
figure,imshow(M)
title('灰白图像')
F=fft2(M); %对灰白图像进行傅里叶变换
F=fftshift(F);
figure;imshow(log(abs(F)),[]); %显示傅里叶变换图像
title('快速傅立叶变换（FFT）图像')
[M,N]=size(F);
nn=2; %二阶巴特沃斯低通
d0=100;
m=fix(M/2);n=fix(N/2);
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
h=1/(1+0.414*(d/d0)^(2*nn)); %过滤器传递函数
result(i,j)=h*F(i,j);
end
end
result=ifftshift(result);
J1=ifft2(result); %傅里叶反变换
J2=uint8(real(J1));
figure,imshow(J2); %显示滤波图像
title('低通滤波图像100')
显示图像分别为原始图像、灰白图像、快速傅立叶变换（FFT）图像、低通滤波图像。
其中快速傅立叶变换（FFT）图像如图11所示，低通滤波图像（d0=100）如图12所示，低通滤波图像（d0=50）如图13所示，低通滤波图像（d0=10）如图14所示。
图11 快速傅立叶变换（FFT）图像
图12 低通滤波图像（d0=100）
图13 低通滤波图像（d0=50）
图14 低通滤波图像（d0=10）
4.3 结果分析
根据低通滤波后图像显示可知：当改变d0的取值时，滤波结果不同。当d0=100时，滤波图像如图12所示，当d0=50时，滤波图像如图13所示，当d0=10时，滤波图像如图14所示，随着d0逐渐减小，图像越来越模糊，滤波带宽越窄，效果越来越差。可见，d0的取值与滤波效果有着密切的关系，d0越大，滤波带宽越宽，能通过低通滤波器的频率范
围越大，图像越清晰。