数字信号处理实验教学-电子教案

1、数字信号处理实验 电子教案讲义1 Matlab简介及其安装使用说明1讲义2 Matlab基本语句7讲义3 Matlab基本数值运算12讲义4 Matlab函数、及其调用方法15实验1 常见离散信号产生和实现19实验2 离散系统的时域分析21实验3 FFT算法的应用23实验4 离散系统的变换域分析26实验5 有限冲激响应数字滤波器设计31实验6 无限冲激响应数字滤波器设计35实验7 设计性和研究性实验40讲义1 Matlab简介及其安装使用说明一. MATLAB程序设计语言简介MATLAB，Matrix Laboratory的缩写，是由Mathworks公司开发的一套用于科学工程计算的可视化高性

2、能语言，具有强大的矩阵运算能力。 与大家常用的Fortran和C等高级语言相比，MATLAB的语法规则更简单，更贴近人的思维方式，被称之为“草稿纸式的语言”。截至目前，MATLAB已经发展到7.x版， 适用于所有32位的Windows操作系统， 按NTFS(NT文件系统)格式下完全安装约需 850 MB。MATLAB软件主要由主包、仿真系统和工具箱三大部分组成。二. MATLAB应用入门1. MATLAB的安装与卸载 MATLAB软件在用户接口设计上具有较强的亲和力，其安装过程比较典型， 直接运行光盘中的安装向导支撑程序SETUP.exe， 按其提示一步步选择即可。MATLAB自身带有卸载程序

3、，在其安装目录下有uninstall子目录，运行该目录下的uninstall.exe即可； 也可以通过Windows系统的安装卸载程序进行卸载。2. MATLAB的启动与退出 MATLAB安装完成后，会自动在Windows桌面上生成一个快捷方式， 它是指向安装目录下binwin32matlab.exe的链接， 双击它即可来到MATLAB集成环境的基本窗口，通常称之为命令窗口。 MATLAB的退出与普通WIN32的程序一样， 值得一提的是它有一个自身专有的快捷键Ctrl+Q。3. MATLAB界面简介 图 1 MATLAB基本界面命令窗口图2图3图4图5图6 指令历史图71)菜单栏菜单栏中包括F

4、ile、Edit、View、Web、Window和Help六个菜单项。这里着重介绍File、help项。File项：数据输入/输出的接口， 包括10个子项， 这里重点介绍其中的5个子项：New：新建文件项。 有四个选择：M File(*.M，文本格式的MATLAB程序文件，可以直接通过文件名的方式在MATLAB环境下解释运行；Figure(图形); Model(仿真模型文件)和GUI(可视化界面文件)。 Open：打开所有MATLAB支持的文件格式，系统将自动识别并采用相应的程序对文件进行处理。例如, 打开一个.m文件， 系统将自动打开M文件编辑器对它进行编辑。 Set Path.：设置工作路

5、径。可以打开路径设置(Set Path)对话框(图2)，将用户自己建立的目录加入MATLAB的目录系统中， 以便所编制的文件能够在MATLAB环境中直接调用。图8 路径设置对话框单击Add Folder.按钮可以将你的一个文件夹加入到系统路径中；Add with Subfolders.允许把一个文件夹包括其所有的子文件夹加入到系统路径中。这两种操作均可以直观地在右侧的路径栏内看到结果。选中一个加入的文件夹，你可以利用Move to Top(移至所有路径的最前面),Move Up(上移一个),Move Down(下移一个),Move to Bottom(移至所有路径的最后面)等四个按钮将改变文件

6、在系统路径中的排列位置以利于对文件的搜索使用， 也可以利用Remove按钮将其删除。对路径操作完毕后，按Save按钮予以保存; Help 项：Matlab Help：打开以html超文本形式存储的的帮助文件主页；Demos：打开matlab演示窗主页；About Matlab：Matlab注册图标、版本、制造商和用户信息；图9 Help选项图10 Help窗口2)命令行区进行通用操作，数值计算，编程和数据类型，输入输出，绘图，用户界面等命令，例如，命令：help函数名（*.m文件）；edit编辑函数、文件对输入命令的解释MATLAB按以下顺序进行： 检查它是否是工作空间中的变量， 是则显示变量

7、内容。 检查它是否是嵌入函数，是则运行之。 检查它是否是子函数。 检查它是否是私有函数。 检查它是否是位于MATLAB搜索路径范围内的函数文件或脚本文件。请注意，如果有两个以上的方案与输入的命令相匹配， MATLAB将只执行第一个匹配。讲义2 Matlab基本语句一. 程序控制语句1. 循环语句 MATLAB的循环语句包括for循环和while循环两种类型。1)for循环语法格式：for 循环变量 = 起始值：步长：终止值循环体end起始值和终止值为一整形数，步长可以为整数或小数，省略步长时，默认步长为1。执行for循环时，判定循环变量的值是否大于(步长为负时则判定是否小于)终止值，不大于(步

8、长为负时则小于)则执行循环体，执行完毕后加上步长， 大于(步长为负时则小于)终止值后退出循环。例1 给矩阵A、B赋值。MATLAB 语句及运行结果如下：k=5;a=zeros(k, k) %矩阵赋零初值for m=1 : kfor n=1: ka(m,n)=1/(m+n-1); endendfor i=m : -1 : 1b(i)=i; end运行结果： a=1.0000 0.5000 0.3333 0.2500?0.20000.5000 0.3333 0.2500 0.2000 0.16670.3333 0.2500 0.2000 0.1667 0.14290.2500 0.2000 0.1

9、667 0.1429 0.12500.2000 0.1667 0.1429 0.1250 0.1111b= 1 2?3 4 52) while循环语法格式： while 表达式循环体end 其执行方式为：若表达式为真(运算值非0)，则执行循环体； 若表达式为假(运算结果为0)，则退出循环体，执行end后的语句。例2a=3;while aa=a-1end输出:a=2a=1a=02. 条件转移语句条件转移语句有if和switch两种。1) if语句MATLAB中if语句的用法与其他高级语言相类似， 其基本语法格式有以下几种：格式一：if 逻辑表达式执行语句end格式二：if 逻辑表达式执行语句1e

10、lse执行语句2end格式三：if 逻辑表达式1执行语句1else? if 逻辑表达式2执行语句2end2) switch语句switch语句的用法与其他高级语言相类似，其基本语法格式为：switch表达式(标量或字符串)case 值1语句1case 值2语句2otherwise语句nend二. 绘图语句常用的MATLAB绘图语句有figure、plot、subplot、stem等， 图形修饰语句有title、axis、text等。1figure figure有两种用法，只用一句figure命令，会创建一个新的图形窗口，并返回一个整数型的窗口编号。figure(n)表示将第n号图形窗口作为当前

11、的图形窗口， 并将其显示在所有窗口的最前面; 如果该图形窗口不存在， 则新建一个窗口，并赋以编号n。 2plot线型绘图函数。用法为plot(x,y,s)。参数x为横轴变量，y为纵轴变量，s用以控制图形的基本特征如颜色、粗细等，通常可以省略，常用方法如表1所示。表13Stem绘制离散序列图，常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。4subplotsubplot(m,n,i)图形显示时分割窗口命令，把一个图形窗口分为m行，n列，mn个小窗口，并指定第i个小窗口为当前窗口。5绘图修饰命令在绘制图形时，我们通常需要为图形

12、添加各种注记以增加可读性。 在plot语句后使用title(标题)可以在图形上方添加标题， 使用xlabel(标记)或ylabel(标记)为X轴或Y轴添加说明，使用text(X值、Y值、想加的标示)可以在图形中任意位置添加标示。例3 画图基本语句如图1所示。MATLAB 语句及运行结果如下：x=0:0.1*pi:2*pi; %定义x向量figure(1); %创建一个新的图形窗口，编号为1subplot(2,2,1); %将窗口划分为2行，2列，在第1个窗口中作图plot(x,sin(x); %画图title(正弦线); %给图形加标题subplot(2,2,2); %在第2个窗口中作图plo

13、t(x,sin(x),r); %画一正弦波，红色xlabel(X); %给x轴加说明ylabel(SIN(X); %给y轴加说明subplot(2,2,3); %在第3个窗口中作图plot(x,sin(x),-); %画一正弦波，破折线subplot(2,2,4); %在第4个窗口中作图plot(x,sin(x),r+); %画一正弦波，红色+线text(4,0,注记);图1讲义3 Matlab基本数值运算一. MATLAB内部特殊变量和常数 MATLAB内部有很多变量和常数，用以表达特殊含义。常用的有： 1. 变量ans： 指示当前未定义变量名的答案； 2. 常数eps：表示浮点相对精度，其

14、值是从1.0到下一个最大浮点数之间的差值。该变量值作为一些MATLAB函数计算的相对浮点精度，按IEEE标准，如：，近似为2.2204e-016； 3. 常数Inf：表示无穷大。当输入或计算中有除以0时产生Inf； 4. 虚数单位i,j：表示复数虚部单位，相当于； 5. NaN：表示不定型值，是由0/0运算产生的。 6. 常数pi：表示圆周率，其值为3.1415926535897； MATLAB中可表示的数字的近似范围从, 1. 有效数字表示的典型例子如下：1234.56789,123456.789E-2,1.23456789e3（format指令可以控制显示格式） 2. 复数形式：3.5+4

15、*j, -2.1-7.4*j（i 也可以）取绝对值：abs()语法格式：abs(x)。当x为实数时计算x的绝对值；x为复数时得到的是复数的模值；x为字符串时得到各字符的ASCII码。取相角：angle()语法格式：angle(z)。求复矢量或复矩阵的相角，结果为一个以弧度为单位介于-和+之间的值。 二. 变量1. 变量命名规则 MATLAB中对变量的命名应遵循以下规则：（1）变量名可以由字母、数字和下划线混合组成，但必须以字母开头；（2）字符长度不能大于31；（3）变量命名区分大小写。2. 局部变量和全局变量 局部变量是指那些每个函数体内自己定义的，不能从其他函数和MATLAB工作空间访问的变

16、量。全局变量是指用关键字“global”声明的变量。全局变量名应尽量大写，并能反映它本身的含义。如果需要在工作空间和几个函数中都能访问一个全局变量，必须在工作空间和这几个函数中都声明该变量是全局的。 三. 矩阵及其运算MATLAB具有强大的矩阵运算和数据处理功能， 对矩阵的处理必须遵从代数规则。1矩阵生成 1）一般矩阵的生成对于一般的矩阵MATLAB的生成方法有多种。最简单的方法是从键盘直接输入矩阵元素。直接输入矩阵元素时应注意：各元素之间用空格或逗号隔开，用分号或回车结束矩阵行，用中括号把矩阵所有元素括起来。例1 在工作空间产生一个33矩阵A可用MATLAB语言描述如下：A=1 2 3; 4

17、 5 6; 7 8 9或A=1 2 34 5 67 8 9运行结果：A=1 2 34 5 67 8 9Size(A)得到矩阵的大小，ans =3 32）特殊矩阵的生成对于特殊的矩阵可直接调用MATLAB的函数生成。用函数zeros生成全0矩阵：格式 B=zeros(m,n)生成mn的全0阵。用函数ones生成全1矩阵：格式 B=ones(m,n)生成mn的全1阵。用函数eye生成单位阵：格式 B=eye(m,n)生成mn矩阵， 其中对角线元素全为1，其他元素为0。2矩阵的运算 矩阵的运算有基本运算和函数运算两种类型。基本运算包括矩阵的加、减、乘、除、乘方、求转置、求逆等，其主要特点是通过MAT

18、LAB提供的基本运算符、*、/（）、等即可完成。函数运算主要是通过调用MATLAB系统内置的运算函数来求取矩阵的行列式（det(A)）,求秩（rank(A)）,求逆(inv(A)，求特征值和特征向量（V， Deig(A)）, 求Jordan标准形（jordan(A)）和矩阵分解等。需要用时可以参阅联机帮助和相关参考书。（.*,.,表示逐个元素的乘积和相除；矩阵X/Y 相当于X*inv(Y), XY相当于inv(Y)*X）例2 矩阵的基本运算。 A=1, 2, 3; 4, 5, 6;B =6, 5, 4; 3, 2, 1;C =A+B %计算两个矩阵的和D =B %计算矩阵B的转置E=A*D %

19、做矩阵乘法，必须要满足矩阵乘法的基本要求E应该是2阶方阵F=det(E) %求E的行列式值G=E(-1) %求E的逆 输出结果：C=7 7 77 7 7D=6 35 24 1E=28 1073 28F=54G=0.5185 -0.1852-1.3519 0.5185讲义4 Matlab函数、及其调用方法在MATLAB语言中，M文件有两种形式：脚本和函数。脚本没有输入/输出参数，只是一些函数和命令的组合。它可以在MATLAB环境下直接执行，也可以访问存在于整个工作空间内的数据。由脚本建立的变量在脚本执行完后仍将保留在工作空间中可以继续对其进行操作，直到使用clear命令对其清除为止。函数是MAT

20、LAB语言的重要组成部分。MATLAB提供的各种工具箱中的M文件几乎都是以函数的形式给出的。函数接收输入参数，返回输出参数，且只能访问该函数本身工作空间中的变量，从命令窗或其他函数中不能对其工作空间的变量进行访问。1函数结构MATLAB语言中提供的函数通常由以下五个部分组成：（1）函数定义行: 以function开头，函数名(必须与文件名相同)及函数输入输出参数在此定义；（2）H1行：第一注释行，供lookfor和help在线帮助使用；（3）函数帮助文件；通常包括函数输入输出参数的含义，调用格式说明；（4）函数体：它包括进行运算和赋值的所有MATLAB程序代码。函数体中可以包括流程控制、输入/

21、输出、计算、赋值、注释以及函数调用和脚本文件调用等。 在函数体中完成对输出参数的计算；（5） 注释。这五个部分中最重要的是函数定义行和函数体。函数定义行是一个MATLAB函数所必需的，其他各部分的内容可以没有， 这种函数称为空函数。例1function sa=circle(r,s)% Circle plot a circle of radii r in the line specified by s% r raddi% s line color% sa area of the circle% circle(r) use blue line to draw a circle of radii r

22、% circle(r,s) use s to draw circle% sa=circle(r) compute the area of the circle and draw it in blue% sa=circle(r,s) compute the area of the circle and draw it in color sif nargin2 error(to many input );endclf; t=0:pi/100:2*pi;x=r*exp(i*t);if nargout=0plot(x,s);elsesa=pi*r*r;fill(real(x),imag(x),s)en

23、daxis(square); %makes the current axis box square in size.例2function y = imp_fun(n,n0)% IMP_FUN Unit impulse function.% IMP_FUN(N,N0), where N is a vector of sequential integers,% returns a vector the same length as N with zeros everywhere except % N = N0.m1,n1 = size(n);if (m1 = 1 | n1 = 1)error(Th

24、e sample vector must be one-dimensional.);endy = zeros(m1,n1);i = find(n = n0);if isempty(i)returnendy(i(1) = 1; % where n = n0, set output to 12函数调用函数调用的过程实际上就是参数传递的过程。例如，在一个脚本文件里调用函数“max”可采用如下方式：n=1:20;a=sin(2*pi*n/20);Y,I=max(a);该调用过程把变量“a”传给了函数中的输入参数“x”,然后把函数运算的返回值传给输出参数“Y”和“I”。其中，Y是a序列的最大值，I是最大

25、值Y对应的坐标值。例3构造：yn = n-3:调用函数： n = 0:6;y = imp_fun(n,3);stem(n,y)图1例4构造：yn = 5n - n-2:调用函数： y=5*imp_fun(n,0) - imp_fun(n,2);stem(n,y)；图2实验1 常见离散信号产生和实现一、实验目的1、加深对常用离散信号的理解；2、熟悉使用MATLAB在时域中产生一些基本的离散时间信号。 二、实验原理1、单位抽样序列 在MATLAB中可以利用函数实现。 2、单位阶越序列 在MATLAB中可以利用函数实现: 3、正弦序列 在MATLAB中实现过程如下： 4、复指数序列 在MATLAB中

26、实现过程如下： 5、指数序列 在MATLAB中实现过程如下： 三、预习要求1、预先阅读实验讲义（MATLAB基础介绍）；2、讨论正弦序列、复指数序列的性质。A绘出信号，当、时、时的信号实部和虚部图；当时呢？此时信号周期为多少？B绘出信号的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期。3、使用帮助功能学习square(方波),sawtooth(锯齿波)和sinc函数，并绘图。四、实验内容 编制程序产生上述5种信号，长度可输入确定，函数需要的参数可输入确定，并绘出其图形。实验2 离散系统的时域分析一、实验目的1、熟悉并掌握离散系统的差分方程表示法；2、加深对

27、冲激响应和卷积分析方法的理解。 二、实验原理在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号,记系统单位冲激响应，则系统响应为如下的卷积计算式：当时，hn是有限长度的（），称系统为FIR系统；反之，称系统为IIR系统。三、预习要求1、在MATLAB中，熟悉利用函数实现差分方程的仿真；2、在MATLAB中，熟悉用函数 计算卷积，用求系统冲激响应的过程。 四、实验内容1、以下程序中分别使用conv和filter函数计算h和x的卷积y和y1，运行程序，并分析y和y1是否有差别，为什么要使用

28、xn补零后的x1来产生y1；具体分析当hn有i个值，xn有j个值，使用filter完成卷积功能，需要如何补零？% Program P2_7clf;h = 3 2 1 -2 1 0 -4 0 3; %impulse responsex = 1 -2 3 -4 3 2 1; %input sequencey = conv(h,x);n = 0:14;subplot(2,1,1);stem(n,y);xlabel(Time index n); ylabel(Amplitude);title(Output Obtained by Convolution); grid;x1 = x zeros(1,8)

29、;y1 = filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel(Time index n); ylabel(Amplitude);title(Output Generated by Filtering); grid; 2、编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。要求分别用 filter、conv、impz三种函数完成。,，给出理论计算结果和程序计算结果并讨论。实验3 FFT算法的应用一、实验目的1、加深对离散信号的DFT的理解；2、在MATLAB中实现FFT算法。二、实验原理N点序列的DFT和IDFT变换定义式如下：,利用旋转

30、因子具有周期性，可以得到快速算法（FFT）。在MATLAB中，可以用函数和计算N点序列的DFT正、反变换。三、预习要求1、在MATLAB中，熟悉函数fft、ifft的使用；2、阅读扩展练习中的实例，学习在MATLAB中的实现FFT算法的实现；3、利用MATLAB编程完成计算，绘出相应图形。并与理论计算相比较，说明实验结果的原因。 四、实验内容1、2N点实数序列N=64。用一个64点的复数FFT程序，一次算出，并绘出 的图形。2、已知某序列在单位圆上的N=64等分样点的Z变换为:。用N点IFFT程序计算出和。五、扩展练习例1：对连续的单一频率周期信号 按采样频率采样，截取长度N分别选N =20和

31、N =16，观察其DFT结果的幅度谱。解：此时离散序列，即k=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下：k=8;n1=0:1:19;xa1=sin(2*pi*n1/k);subplot(2,2,1)plot(n1,xa1)xlabel(t/T);ylabel(x(n);xk1=fft(xa1);xk1=abs(xk1);subplot(2,2,2)stem(n1,xk1)xlabel(k);ylabel(X(k);n2=0:1:15;xa2=sin(2*pi*n2/k);subplot(2,2,3)plot(n2,xa2)xlabel(t/T);ylabe

32、l(x(n);xk2=fft(xa2);xk2=abs(xk2);subplot(2,2,4)stem(n2,xk2)xlabel(k);ylabel(X(k);图 不同的截取长度的正弦信号及其DFT结果计算结果示于图，(a)和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。实验4 离散系统的变换域分析一、实验目的1、熟悉对离散系统的频率响应分析方法；2、加深对零、极点分布的概念理解。 二、实验原

33、理离散系统的时域方程为其变换域分析方法如下：频域:系统的频率响应为:Z域:系统的转移函数为:分解因式:，其中和称为零、极点。三、预习要求1. 在MATLAB中，熟悉函数tf2zp、zplane、freqz、residuez、zp2sos的使用，其中：z，p，K=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点；zplane（z，p）绘制零、极点分布图；h=freqz(num,den,w)求系统的单位频率响应；r，p，k=residuez（num，den）完成部分分式展开计算；sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。 2. 阅读扩展练习中的实例，

34、学习频率分析法在MATLAB中的实现； 3. 编程实现系统参数输入，绘出幅度频率响应和相位响应曲线和零、极点分布图。 四、实验内容求系统的零、极点和幅度频率响应和相位响应。五、扩展练习例1: 求下列直接型系统函数的零、极点，并将它转换成二阶节形式解:用MATLAB计算程序如下：num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den);m=abs(p);disp(零点);disp(z);disp(极点);disp(p);disp(增益系数);disp(k);sos=zp2sos(z,p,k);disp(二阶节);di

35、sp(real(sos);zplane(num,den)输入到“num”和“den”的分别为分子和分母多项式的系数。计算求得零、极点增益系数和二阶节的系数：零点:0.9615-0.5730-0.1443 + 0.5850i-0.1443 - 0.5850i极点:0.5276+0.6997i0.5276-0.6997i-0.5776+0.5635i-0.5776-0.5635i增益系数:1二阶节:1.0000 -0.3885 -0.5509 1.0000 1.15520 0.65111.0000 0.28850 0.36300 1.0000 -1.0552 0.7679系统函数的二阶节形式为：极

36、点图见图:图 系统函数的零、极点图例2: 差分方程所对应的系统的频率响应。解:差分方程所对应的系统函数为：用MATLAB计算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle(实部)xlabel(omega/pi);ylabel(幅度)subplot(2,2,2);plot(w/pi,imag(h);gridtitle(虚部)xlabel(omega/pi);ylabel(Amplitu

37、de)subplot(2,2,3);plot(w/pi,abs(h);gridtitle(幅度谱)xlabel(omega/pi);ylabel(幅值)subplot(2,2,4);plot(w/pi,angle(h);gridtitle(相位谱)xlabel(omega/pi);ylabel(弧度) 图实验5 有限冲激响应数字滤波器设计一、实验目的： 1、加深对数字滤波器的常用指标理解。 2、学习数字滤波器的设计方法。 二、实验原理：低通滤波器的常用指标：（1）通带边缘频率；（2）阻带边缘频率；（3）通带起伏；（4）通带峰值起伏，（5）阻带起伏，最小阻带衰减。三、预习要求1. 在MATLAB

38、中，熟悉函数fir1、kaiserord 、remezord、remez的使用； B = fir1(n,Wn,high,noscale)设计滤波器； n,Wn,beta,ftype = kaiserord(f,a,dev)估计滤波器阶数； n,fo,ao,w = remezord (f,a,dev,fs)计算等波纹滤波器阶数n和加权函数w()； B=remez(n,f,a)进行等波纹滤波器的设计。 2. 阅读扩展练习中的实例，学习FIR滤波器的设计方法及其在MATLAB中的实现； 3. 给出FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 数字滤波器有I

39、IR和FIR两种类型，它们的特点和设计方法不同。四、实验内容：利用MATLAB编程,分别用窗函数法和等波纹滤波器法设计两种FIR数字滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：。阻带边缘频率：，最小阻带衰减：。 五、扩展练习例1: 用凯塞窗设计一FIR低通滤波器，通带边界频率，阻带边界频率，阻带衰减不小于50dB。解: 首先由过渡带宽和阻带衰减来决定凯塞窗的N和上图给出了以上设计的频率特性，(a) 为N=30直接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的MATLAB程序为：wn=kaiser(30,4.55);nn=0:1:29;alfa=(30-1)/2;hd=si

40、n(0.4*pi*(nn-alfa)./(pi*(nn-alfa);h=hd.*wn;h1,w1=freqz(h,1);或者：b = fir1(29,0.4,kaiser(30,4.55);h1,w1=freqz(b,1);plot(w1/pi,20*log10(abs(h1);axis(0,1,-80,10); grid;xlabel(归一化频率/p);ylabel(幅度/dB);还可以使用n,Wn,beta,ftype=kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器：fcuts = 0.3 0.5; %归一化频率omega/pimags = 1 0;devs

41、= 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs); %计算出凯塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),noscale);freqz(hh);实际中，一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W()，调用函数remez可进行等波纹滤波器的设计，直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界例2：利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=8

42、00Hz，阻带边界fr=1000Hz，通带波动阻带最小衰减At=40dB，采样频率fs=4000Hz。解:在MATLAB中可以用remezord 和remez两个函数设计，其结果如图2，MATLAB程序如下：fedge=800 1000; mval=1 0; dev=0.0559 0.01; fs=4000; N,fpts,mag,wt=remezord(fedge,mval,dev,fs);b=remez(N,fpts,mag,wt);h,w=freqz(b,1,256);plot(w*2000/pi,20*log10(abs(h);grid;xlabel(频率/Hz);ylabel(幅度/

43、dB);所得图像如下所示：实验6 无限冲激响应数字滤波器设计一、实验目的1、掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法;2、熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。二、实验原理在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1)利用buttord和cheb1ord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；2)num,den=butter（N,Wn）(巴特沃斯)和num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；3）lp2hp，

44、lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；4)使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；5)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。三、预习要求1. 在MATLAB中，熟悉函数butter、cheby1、cheby2的使用，其中：num,den=butter（N,Wn）巴特沃斯滤波器设计；num,den=cheby1（N,Wn）切比雪夫1型滤波器设计；num,den=cheby2（N,Wn）切比雪夫2型滤波器设计。2. 阅读扩展练习中的实例，学习在MATLAB中进行数字滤波器的设计； 3.

45、给出IIR数字滤波器参数和滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 四、实验内容利用MATLAB编程，用脉冲响应不变法和双线性变换法设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：；阻带边缘频率：，最小阻带衰减： 。五、扩展练习例1：设采样周期T=250s（采样频率fs =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc =1kHz。B,A=butter(3,2*pi*1000,s);num1,den1=impinvar(B,A,4000);h1,w=freqz(num1,den1);B,A=

46、butter(3,2/0.00025,s);num2,den2=bilinear(B,A,4000);h2,w=freqz(num2,den2);f=w/pi*2000; plot(f,abs(h1),-.,f,abs(h2),-);grid;xlabel(频率/Hz )ylabel(幅值/dB)程序中第一个butter的边界频率21000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应

47、，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。同时，也看到双线性变换法，在z=-1即=或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在=处的三阶传输零点通过映射形成的。下图给出了MATLAB计算的结果。 例2: 设计一数字高通滤波器，它的通带为400500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。wc=2*1000*tan(2*pi*400/(2*1000);wt=2*1000*tan(2*pi*317/(2*1000);N,wn=cheb1ord(wc,wt,0.5,19,s);B,A=cheb

48、y1(N,0.5,wn,high,s);num,den=bilinear(B,A,1000);h,w=freqz(num,den);f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel()ylabel(幅度/dB) 下图给出了MATLAB计算的结果。例3： 设计一巴特沃兹带通滤波器，其dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3 = 120kHz处的最小衰减大于dB，采样频率fs=400kHz。w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0 wr,3,10,s)