数字信号处理实验指导书

《数字信号处理》实验指导书

王莉

南京工业大学

自动化与电气工程学院

2008-04-17

目 录

实验一 信号、系统及系统响应…………………………………3 实验二 用双线性变换法设计IIR数字滤波器 ………………6 实验三 用窗函数法设计FIR数字滤波器 ………………… 10 附录 MATLAB信号处理工具箱函数……………………………14

实验一 信号、系统及系统响应

一．实验目的

1． 熟悉时域离散系统的时域特性。 2． 验证时域的卷积定理。

3． 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对离散信号及系

统响应进行频域分析。

二．实验内容

1． 观察信号xa(n)和系统h(n)的时域和频域特性，并绘出相应的曲线。①单位脉冲

序列：

xa(n)(n)；②系统单位脉冲响应序列：

h(n)(n)2.5(n1)2.5(n2)(n3)。

2． 利用线性卷积求信号xa(n)通过系统h(n)的响应ya(n)，比较所求响应ya(n)和

h(n)的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解

释所得结果。

3． 卷积定理的验证。将2中的信号换成xb(n)xb(nT)AeanTsin(0nT)u(n)，使

a=0.4，0=2.0734，A=1，T=1，重复实验2，求出yb(n)，绘出其频率特性Yb(ej)曲线；利用公式Y(ej)Xb(ej)H(ej)，并绘出Y(e)的幅频和相频特性曲线，与前面直接对yb(n)进行傅里叶变换所得频率特性曲线进行比较，验证时域卷积定理。

j三．实验设备及仪器

1． 计算机。 2． Matlab软件。

四. 实验线路及原理

1

离散信号和系统在时域均可用序列来表示。序列图形给人以形象直观的印象，它可加深我们对信号和系统的时域特征的理解。本实验将观察分析几种信号及系统的时域特性。

序列和信号的傅立叶变换是ω的连续函数，而计算机只能计算出有限个离散频率点的函数值。因此在取得频谱函数后，应该在0～2π之间取许多点，计算这些点的频谱函数的值，并取它们的包络，该包络才是需要的频率特性。当然，点数取得多一些，该包络才

jk能接近真正得频率特性。通常对X(ej)在［0，2π］上取模X(e)，绘出幅频特性曲

线进行观察分析。

一个时域离散线性非时变系统的输入/输出关系为

y(n)x(n)h(n)mx(m)h(nm) (1.1)

这里，y(n)为系统的输出序列，x(n)为输入序列。h(n)和x(n)可以是无限长，也可以是有限长。为了计算机绘图方便，主要讨论有限长情况。

上述卷积运算也可以在频域实现

Y(ej)X(ej)H(ej) (1.2)

(1.2)式右边的相乘是各频点的频谱值相乘。

五. 实验方法与步骤

1． 编制信号xa(n)和系统h(n)产生程序，并绘出时域特性曲线。

2． 求信号xa(n)和系统h(n)的傅立叶变换，并画出xa(n)和h(n)的幅频特性和相

频特性曲线。

3． 利用卷积函数conv（）求信号xa(n)通过系统h(n)的响应ya(n)，并绘出ya(n)的时域特性曲线和幅频特性曲线。 4． 求信号xb(n)xb(nT)AeanTa=0.4，0=2.0734，A=1，T=1sin(0nT)u(n)，

的时域特性曲线和频域特性曲线。

5． 利用卷积函数conv（）求信号xb(n)通过系统h(n)的响应yb(n)，并求yb(n)的

傅立叶变换，并绘出其幅频特性和相频特性曲线。 6． 利用公式Y(ej)Xb(ej)H(ej)计算Y(e)，并绘出Y(e 2

jjk)幅频特性和相频特性

曲线。

七. 实验报告内容与要求

1． 简述实验目的、实验原理及实验步骤。

2． 按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时

域和幅频特性曲线，并对所得结果进行分析和解释。 3． 总结实验中的主要结论。 4． 简要回答思考题。

八. 思考

1．信号的频率特性即信号的傅立叶变换利用MATLAB程序如何实现？

2．在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得Y(ej)Xb(ej)H(ej)，k=0,1,…,M-1，所得结果之

kkk间有无差异？为什么？

3

实验二 用双线性变换法设计IIR数字滤波器

一．实验目的

1． 熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 2． 掌握数字滤波器的计算机仿真方法。

3． 通过观察对实际心电图信号的滤波作用， 获得数字滤波的感性知识。

二．实验内容

1． 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通

带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π，π]频率区间上，最大衰减大于15dB。

2． 以0.02π为采样间隔，打印出数字滤波器在频率区间［0,π］上的幅频响应特性

曲线。

3． 用所设计的滤波器对实际心电图信号采样序列x(n)进行仿真滤波处理，并分别打

印出滤波前后的心电图信号波形图， 观察总结滤波作用与效果。人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。 {x(n)} ={-4, -2, 0, -4, -6, -4, -2, -4, -6, -6,

-4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0,

0, -2, -2, -2, -2, 0}

三．实验设备及仪器

1． 计算机。 2． Matlab软件。

四. 实验线路及原理

4

IIR数字滤波器设计方法有两类，经常用到的一类设计方法是借助于模拟滤波器的设计方法进行的，其设计流程图如图2.1所示。

(1)确定数字滤波器指标参数(2)选择合适的数字滤波器设计方法将数字滤波器指标转换成相应的模拟滤波器指标(3)设计相应的数字滤波器系统函数Ha(s)(4)用所选设计方法将Ha(s)转换成数字滤波器系统函数H(z)图2.1 从模拟滤波器设计数字滤波器流程图

IIR数字滤波器经常用到的方法有：脉冲相应不变法和双线性变换法。本实验采用双线性变换法设计IIR数字滤波器。双线性变换法能够彻底消除频率混叠失真。其设计思想是算法逼近（脉冲相应不变法是波形逼近）， 用差分近似微分，用相邻两个采样的平均值近似微分方程中的瞬时值xa(t)和ya(t)，用模拟滤波器的微分方程变成差分方程，经Z变换得到H(z)。按这种思想导出的双线性映射公式为：

H(z)Ha(s)|T为采样间隔。ω与Ω的关系为：

21z1sT1z1 （2.1）

2tan （2.2） T2五. 实验方法与步骤

1． 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内

容，用双线性变换法设计数字滤波器系统函数H(z)。《数字信号处理》（第二版）（丁玉美 高西全编著）中例6.4.2已求出满足本实验要求的数字滤波器系统函数：

0.0007378(1z1)6 H(z)(11.268z10.705z2)(11.0106z10.3583z2)(10.904z10.215z2)Hk(z) （2.3）

K13 5

A(12z1z2)式中 Hk(z),k1,2,3 （2.4） 121BkzCkzA=0.09036

B1=1.2686, C1=-0.7051 B2=1.0106, C2=-0.3583 B3=0.9044, C3=-0.2155

根据设计指标，调用MATLAB信号处理工具箱函数buttord和butter，billinear，也可以得到H(z)。

由(2.3)式和(2.4)式可见，滤波器H(z)由三个二阶滤波器H1(z),H2(z)和H3(z)级联组成，如图2.2所示。

x(n)H1(z)y1(n)H2(z)H(z)y2(n)H3(z)y3(n)＝y(n)图2.2 滤波器H(z)的组成

2． 编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。设

yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk-1(n)为输入序列，如图2.2所示。由(2.4)式可得到差分方程：

yk(n)Ayk1(n)2Ayk1(n1)Ayk1(n2)Bkyk(n1)Ckyk(n2) （2.5）

当k=1时，yk1(n)x(n)。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即依次对k=1，2，3，求解差分方程(2.5)，最后得到y3(n)y(n)。可以直接调用MATLAB filter函数实现仿真。

3． 在通用计算机上运行仿真滤波程序，并调用通用绘图程序，完成实验内容(2)和(3)。

七. 实验报告内容与要求

1． 简述实验目的、实验原理及实验步骤。

2． 由所打印的|H(ej)|特性曲线及设计过程简述双线性变换法的特点。

ω

3． 对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。 4． 对比所设计的巴特沃斯模拟滤波器和数字滤波器的幅频特性。 5． 简要回答思考题。

6

八. 思考

21z11．用双线性变换法设计数字滤波器过程中，变换公式s中T的取值，对设

T1z1计结果有无影响?为什么?

7

实验三 用窗函数法设计FIR数字滤波器

一．实验目的

1． 掌握用窗函数法设计FIR数字滤波器的原理和方法。 2． 熟悉线性相位FIR数字滤波器特性。 3． 了解各种窗函数对滤波特性的影响。

二．实验内容

1． 用矩形窗设计一线性相位低通FIR数字滤波器，截至频率c4rad。窗口长

度N=15，33。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察过渡带带宽。总结窗口长度N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数

je,cjHd(e)0,c N12Hd(ej)，即

其中

11jjnhd(n)H(e)edd22sinc(na)(na)2． N=33，ccejaejnd

c4rad，用四种窗函数（矩型窗、汉宁窗、哈明窗和布莱克曼窗）

设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察过渡带带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。

三．实验设备及仪器

1． 计算机。 2． Matlab软件。

8

四. 实验线路及原理

如果所希望的滤波器的理想频率响应函数为Hd(ej)，则其对应的单位脉冲响应为

hd(n)12Hd(ej)ejnd (3.1)

窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼近hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数w(n)将hd(n)截断，并进行加权处理，得到：

h(n)hd(n)w(n) (3.2)

h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数

H(ej)为

H(e)h(n)ejn (3.3)

jn0N1式中，N为所选窗函数w(n)的长度。

用窗函数设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表3.1。

表3.1 六种窗函数的基本参数

窗函数矩形窗三角形窗汉宁窗哈明窗布莱克曼窗凯塞窗（-7.865）过渡带宽4/N/N/N/N/N/N阻带最小衰减/dB-21-25-44-53-74-80

这样选定窗函数类型和长度N后，求出单位脉冲响应h(n)hd(n)w(n)，并按照式(3.3)求出H(ej)。H(ej)是否满足要求，要进行验算。如果H(ej)不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。 如果要求线性相位特性，则h(n)还必须满足：

9

h(n)h(N1n)

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择

h(n)h(N1n)一类，而不能选h(n)h(N1n)一类。

五. 实验方法与步骤

1． 复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。 2． 编写能产生理想滤波器时域特性的函数。

3． 用Matlab的函数实现矩型窗、汉宁窗、哈明窗和布莱克曼窗的窗函数。 4． 编写主程序。主程序框图如图3.1所示，仅供参考。其中幅度特性要求用dB表

示。

5． 也可以采用其他方式完成该实验，例如采用fir1函数。

开始读入窗口长度N计算hd(n)调用窗函数子程序求w(n)计算h(n)=hd(n)w(n)调用子程序计算H(ejw)=DFT[h(n)]调用绘图子程序绘制H(ejw)幅度相位曲线结束

图3.1 用窗函数法设计滤波器主程序框图

10

j可以用fft()函数来求h(n)的频率特性，画图时用20lgH(e)打印幅度特性。第k

点对应的频率k区间要选大些。

2k。为使曲线包络更接近H(ej)的幅度特性曲线，DFT变换N七. 实验报告内容与要求

1． 简述实验目的、实验原理及实验步骤。

2． 按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类

型对滤波特性的影响。

3． 总结用窗函数法设计FIR滤波器的主要特点。 4． 简要回答思考题。

八. 思考

1． 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计

线性相位低通滤波器? 写出设计步骤。

2． 什么是吉布斯效应？为什么会产生吉布斯效应？

11

附录 MATLAB信号处理工具箱函数

MATLAB是MathWorks公司于1982年推出的一种功能强大、效率高、交互性好的数值计算和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。随着其自身版本的不断提高，MATLAB的功能越来越强大，应用范围也越来越广，如广泛应用于数字信号处理、数字图像处理、仿真、自动控制、小波分析、神经网络等领域。与C语言或FORTRON语言作科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。MATLAB有着功能强大、丰富的函数工具箱，这是整个MATLAB语言得以如此快速发展的重要因素之一。这些函数工具箱大致可以分为两类：功能型工具箱提供了对内核的支持，它主要用来扩充MATLAB的符号计算功能、图形建模仿真功能、文字处理功能以及与硬件实时交互功能，属于MATLAB自身系统，可以用于多种学科。而领域型工具箱是专业性很强的学科研发性工具箱，应用于不同的学科。MATLAB信号处理工具箱包含了许多信号处理函数，这些函数的格式和说明很容易通过help命令查到。调用这些函数很容易验证数字信号处理习题答案，使数字滤波器分析与设计的程序非常简单。为了便于使用，下面简要介绍MATLAB的启动、编程和运行程序的过程，最后给出实验中可能用到的工具箱函数。

a) 启动MATLAB系统及编程

在Windows桌面上双击MATLAB，启动MATLAB系统，Windows工作平台上自动弹出一个窗口，其标题为MATLAB Command Window，称其为MATLAB命令窗口。该窗口是用户与MATLAB解释器进行通信的工作环境，所有的MATLAB函数、命令及程序都要在该命令窗口下运行。

在命令窗口中，用户可以发出MATLAB命令。每条命令输入并按回车键后，MATLAB系统便解释并执行之，显示命令执行结果。应当注意，在命令窗口输入命令只适合一些简单的运算程序，不便于编辑修改较复杂的程序。M文件为这种编辑运行程序过程提供了方便。下面简单介绍在MATLAB命令窗口状态下M文件的编辑和运行。

1）新建M文件

在MATLAB命令窗口状态，选File-New-Mfile，MATLAB系统弹出文件编辑窗口。 2）运行与修改M文件

12

MATLAB采用的是命令模式，用户每输入一条命令，回车后MATLAB就解释并执行这条命令，再根据要求显示运算结果。此外，MATLAB也可以执行某个文件中的MATLAB语句序列。这两种模式一起构成了MATLAB的解释环境。

当一个M文件编辑存盘后，在MATLAB命令窗口输入M文件名并回车，则MATLAB系统就会逐行解释并执行该M文件中的命令序列。

应当注意，如果已经编辑好的M文件不在MATLAB系统搜索路径中，用户应当用path命令在MATLAB的搜索路径中添加新的搜索路径。例如，M文件存在C:\\MYFILES路径，下面的命令：

path(‘C:MYFILES’，path)；

将搜索路径改为先搜索C:\\MYFILES，再在当前的搜索路径目录中搜索。

b) 实验用MATLAB函数简介 1）abs

功能：求绝对值（模值） 格式：y=abs(x)

说明：y=abs(x)用于计算实数x的绝对值。当x为复数时得到x的模（幅度值）。当x为矢量时，计算其每个元素的模，返回模向量y。

2） angle

功能：求相角。 格式：Ph=angle(x)

说明：Ph=angle(x)用于计算复矢量x的相角(rad)。Ph值介于-π和+π之间。 3）conv

功能：计算离散卷积、多项式相乘。 格式：yn=conv(hn, xn)

说明：yn=conv(hn, xn)用于计算序列hn和xn的卷积，卷积的结果序列为yn。 例如，程序：

hn=[1, 1, 1, 1]；xn=[1, 1, 1, 1]；yn=conv(hn, xn)； 运行结果为： yn=[1 2 3 4 5 3 2 1]

如果N和M分别表示hn和xn的长度，则yn的长度为N+M-1。 4）filter

功能：利用IIR滤波器或者FIR滤波器对数据进行滤波。

13

格式：y=filter(b, a, x) [y, zf]=filter(b, a, x) y=filter(b, a, x, zi)

说明：y=filter(b, a, x)利用给定系数矢量a和b对x中的数据进行滤波，结果放入y矢量中，y的长度取max(N, M)。

y=filter(b, a, x, zi)可在zi中指定x的初始状态。

[y, zf]=filter(b, a, x)除了得到矢量y外，还得到x的最终状态矢量zf。 5）freqz

功能：计算数字滤波器H(z)的频率响应。 格式：[h, w]=freqz(b, a, n) [h, f]=freqz(b, a, n, Fs) h=freqz(b, a, w) h=freqz(b, a, f, Fs) freqz(b, a)

说明：freqz用于计算数字滤波器H(z)的频率响应波器系统函数H(z)的分子和分母多项式系数。

B(z)b0b1z1bM1z(M1)bMzMH(z)A(z)1a1z1aN1z(N1)aNzNH(ej)。矢量a和b分别为数字滤

[h, w]=freqz(b, a, n)计算出数字滤波器n个频率点上的频率响应，存放在h矢量中，n个频率存放在矢量w中。Freqz函数自动将这n个频点均匀设置在频率范围[0,2π]上。缺省whole时，n个频点均匀设置在频率范围[0,π]上。缺省w和M时，freqz自动选取512个频率点计算。

[h, f ]=freqz(b, a, n, Fs)用于对

H(ej)在[o,Fs/2]上等间隔采样n点，采样点频率及相应

频率响应值分别记录在f和h中，由用户指定Fs（单位为Hz）值。

jh=freqz(b, a, w)用于对H(e)在[0,2π]上进行采样，采样频率点由矢量w指定。 jh=freqz(b, a, f, Fs)用于对对H(e)在[0,Fs]上进行采样，采样频率点由矢量f指定。

freqz(b, a)用于在当前图形窗口中绘制出幅频和相频特性曲线。 6）impz

功能：计算数字滤波器的单位脉冲响应。

14

格式：[h, t]=impz(b, a) [h, t]=impz(b, a, N) [h, t]=impz(b, a, n, Fs) impz(b, a)

说明：[h, t]=impz(b, a)计算出单位脉冲响应序列向量h，并返回时间列向量t。自动选择h的样点数。调用参数b和a分别为数字滤波器系统函数H(z)的分子和分母多项式系数矢量。

[h, t]=impz(b, a, N)计算出单位脉冲响应N个样值。 如果N是一个整数向量，则仅计算出这些整数点的单位脉冲响应样值。

impz(b, a)则完成相应的计算功能，并自动调用stem(n, hn)绘图。其他调用格式用help命令查看。

例如，系统的差分方程y1(n)0.9y1(n1)x(n)对应的系统函数为 H(z)110.9z1

计算出单位脉冲响应50个样值，并给出绘图的程序如下： B=1；A=[1，-0.9]；impz(B, A, 50)； 7）fft

功能：一维快速傅立叶变换。 格式：y=fft(x) y=fft(x, n)

说明：y=fft(x) 利用FFT算法计算矢量x的离散傅立叶变换，当x为矩阵时，y为矩阵x的每一列FFT。y=fft(x, n)采用n点FFT。当x长度小于n时，fft函数自动在x尾部补零，以构成n点数据；当x的长度大于n时，fft截取x前面的n点数据进行FFT。

15