实验讲义

《基于MATLAB的信号与系统及数字信号处理仿真实验》

实验教学大纲

实验课名称：基于MATLAB的信号与系统及数字信号处理仿真实验 开课学期：第6学期 课程性质：实践教学课程 学 分：0.5 实验总学时：18 实验周学时：3

开课单位：电子工程系、网络与通信工程系

开设实验项目数：6（其中4个验证性、1个设计性,、1个综合性） 考核方法：实验过程及实验报告成绩综合

实验目录

实验一 MATLAB基本编程实验 .......................................................................................... 3 实验二 连续时间系统时域分析实验 ................................................................................... 14 实验三 离散时间系统时域分析及稳定性实验 .................................................................... 20 实验四 信号的频谱分析实验.............................................................................................. 22 实验五 数字滤波器设计与仿真实验 ................................................................................... 26 实验六 数字信号处理在双音多频拨号系统中的应用 ......................................................... 29 附录：实验报告模版 .......................................................................................................... 34

2

实验一 MATLAB基本编程实验

一、实验说明

1、 实验类型：验证性实验 2、 实验课时：3学时

二、实验目的

1、 熟悉MATLAB软件的基本编程方法。 2、 掌握常用基本运算方法。 3、 掌握简单的绘图命令。

4、 用MATLAB编程并会创建M文件、函数。 5、 熟悉常见连续信号和离散信号的生成方法。

三、实验原理

1、 MATLAB软件介绍

MATLAB是美国Mathworks公司于1984年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。MATLAB软件包括五大通用功能：数值计算功能（Numeric），符号运算功能（Symbolic）；数据可视化功能（Graphic），数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks公司已推出30多个应用工具箱。MATLAB在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。

MATLAB 6.5的工作桌面由标题栏、菜单栏、工具栏、命令窗口（Command Window）、工作空间窗口（Workspace）、当前目录窗口（Current Directory）、历史命令窗口（Command History）及状态栏组成，从而为用户使用MATLAB提供了集成的交互式图形界面，如图1.1所示。

图1.1 MATLAB的工作桌面

3

MATLAB的命令窗口是接收用户输入命令及输出数据显示的窗口，几乎所有的MATLAB行为都是在命令窗口进行的。当启动MATLAB软件时，命令窗口就做好了接收指令和输入的准备，并出现命令提示符（）。在命令提示符后输入指令通常会创建一个或多个变量。变量可以是多种类型的，包括函数和字符串，但通常的变量只是数据。这些变量被放置在MATLAB的工作空间中，工作空间窗口提供了变量的一些重要信息，包括变量的名称、维数大小、占用内存大小、以及数据类型等信息。查看工作空间的另一种方法是使用whos命令。在命令提示符后输入whos命令，工作空间中的内容概要将作为输出显示在命令窗口中。

有的命令可以用来清除不必要的数据，同时释放部分系统资源。clear命令可以用来清除工作空间的所有变量，如果要清除某一特定变量则需要在clear命令后加上该变量的名称。另外，clc命令用来清除命令窗口的内容。

如果希望将MATLAB所创建的变量及重要数据保留下来，则使用save命令，并在其后加上文件名，即可将整个工作空间保存为一个扩展名为.mat的文件。使用load命令，并在其后加上文件名，则可将MATLAB数据文件（.mat文件）中的数据加载到工作空间中。MATLAB历史命令窗口记录了每次输入的命令。在该窗口中可以对以前的历史命令进行查看、复制或者直接运行。

对于初学者而言，需要掌握的最重要且最有用的命令应为help命令。MATLAB命令和函数有数千个，而且许多命令的功能非常强大，调用形式多样。要想了解一个命令或函数，只需在命令提示符后输入help，并加上该命令或函数的名称，则MATLAB会给出其详细帮助信息。另外，MATLAB还精心设计了演示程序系统（Demo），内容包括MATLAB的内部主要函数和各个工具箱（Toolbox）的使用。初学者可以方便地通过这些演示程序及其给出的程序源代码进行直观的感受和学习。用户可以通过两种途径打开演示程序系统。一是在命令窗口中输入demo或demos命令并按Enter键；二是选择“help”→“Demos”菜单命令。 2、 MATLAB中常用基本运算

① 算术运算

MATLAB可以像一个简单的计算器一样使用，不论是实数运算还是复数运算都能轻松完成。标量的加法、减法、除法和幂运算均可通过常规符号“+”、“?”、“*”、“/”、以及“^”来完成。对于复数中的虚数单位，MATLAB采用预定义变量i或j表示，即i=j=?1。因此，一个复常量可用直角坐标形式来表示，例如，

A=-3-i*4 A=

-3.0000 - 4.0000i

将复常量-3-i4赋予了变量A。

一个复常量还可用极坐标的形式来表示，例如，

B=2*exp(i*pi/6) B=

1.7321 + 1.0000i

其中，pi是MATLAB预定义变量，pi=?。

复数的实部和虚部可以通过real和imag运算符来实现，而复数的模和辐角可以通过abs

4

和angle运算符来实现。但应注意辐角的单位为弧度。例如，复数A的模和辐角、复数B的实部和虚部的计算分别为

A_mag=abs(A) A_mag=

5

A_rad=angle(A) A_rad= -2.2143 B_real=real(B) B_real= 1.7321 B_imag=imag(B) B_imag= 1.0000

如果将弧度值用“度”来表示，则可进行转换，即

A_deg=angle(A)*180/pi A_deg= -126.8699

② 向量运算

向量是组成矩阵的基本元素之一，MATLAB具有关于向量运算的强大功能。一般地，向量被分为行向量和列向量。生成向量的方法有很多，我们主要介绍两种。

a) 直接输入向量：即把向量中的每个元素列举出来。向量元素要用“[ ]”括起来，元素

之间可用空格、逗号分隔生成行向量，用分号分隔生成列向量。例如，

A=[1,3,5,21] A=

1 3 5 21 B=[1;3;5;21] B= 1 3 5 21

b) 利用冒号表达式生成向量：这种方法用于生成等步长或均匀等分的行向量，其表达

式为x=x0:step:xn。其中，x0为初始值；step表示步长或增量；xn为结束值。如果step值缺省，则步长默认为1。例如，

C=0:2:10 C=

0 2 4 6 8 10 D=0:10 D=

0 1 2 3 4 5 6 7 8 9 10

在连续时间信号和离散时间信号的表示过程中，我们经常要用到冒号表达式。例如，对于0?t?1范围内的连续信号，可用冒号表达式“t=0:0.001:1;”来近似表达该区间，此时，向量t表示该区间以0.001为间隔的1001个点。

5

如果一个向量或一个标量与一个数进行运算，即“+”、“?”、“*”、“/”、以及“^”运算，则运算结果是将该向量的每一个元素与这个数逐一进行相应的运算所得到的新的向量。例如，

C=0:2:10; E=C/4 E=

0 0.5000 1.0000 1.5000 2.0000 2.5000

其中，第一行语句结束的分号是为了不显示C的结果；第二句没有分号则显示出E的结果。

一个向量中元素的个数可以通过命令“length”获得，例如， t=0:0.001:1;

L=length(t) L= 1001

③ 矩阵运算

MATLAB又称矩阵实验室，因此，MATLAB中矩阵的表示十分方便。例如，输入矩阵

?111213??212223?在MATLAB命令窗口中可输入下列命令得到，即 ????313233??a=[11 12 13;21 22 23;31 32 33] a=

11 12 13 21 22 23 31 32 33

其中，命令中整个矩阵用括号“[ ]”括起来；矩阵每一行的各个元素必须用逗号“，”或空格分开；矩阵的不同行之间必须用分号“；”或者按Enter键分开。

在矩阵的加减运算中，矩阵维数相同才能实行加减运算。矩阵的加法或减法运算是将矩阵的对应元素分别进行加法或减法运算。在矩阵的乘法运算中，要求两矩阵必须维数相容，即第一个矩阵的列数必须等于第二个矩阵的行数。例如，

a=[1 2 3;4 5 6] a=

1 2 3 4 5 6 b=[1 2; 3 4;5 6] b= 1 2 3 4 5 6 c=a*b c=

22 28 49 64

6

MATLAB中矩阵的点运算指维数相同的矩阵位置对应元素进行的算术运算，标量常数可以和矩阵进行任何点运算。常用的点运算包括“.*”、“./”、“.\\”、“.^”等。矩阵的加法和减法是在对应元素之间进行的，所以不存在点加法或点减法。

点乘运算，又称Hadamard乘积，是指两维数相同的矩阵或向量对应元素相乘，表示为

C=A.*B。点除运算是指两维数相同的矩阵或向量中各元素独立的除运算，包括点右除和点左除。其中，点右除表示为C＝A./B，意思是A对应元素除以B对应元素；点左除表示为C＝A.\\B，意思是B对应元素除以A对应元素。点幂运算指两维数相同的矩阵或向量各元素独立的幂运算，表达式为C＝A.^B。 ④ 符号运算

MATLAB符号运算工具箱提供的函数命令是专门研究符号运算功能的。符号运算是指符号之间的运算，其运算结果仍以标准的符号形式表达。符号运算是MATLAB的一个极其重要的组成部分，符号表示的解析式比数值解具有更好的通用性。在使用符号运算之前必须定义符号变量，并创建符号表达式。定义符号变量的语句格式为

syms 变量名

其中，各个变量名须用空格隔开。例如，定义x、y、z三个符号变量的语句格式为

syms x y z

我们可以用whos命令来查看所定义的符号变量，即

clear syms x y z whos

Name Size Bytes Class x 1x1 126 sym object y 1x1 126 sym object z 1x1 126 sym object Grand total is 6 elements using 378 bytes

可见，变量x、y、z必须通过符号对象定义，即sym object，才能参与符号运算。

另一种定义符号变量的语句格式为 sym('变量名')

例如，x、y、z三个符号变量定义的语句格式为

x=sym('x'); y=sym('y'); z=sym('z');

sym语句还可以用来定义符号表达式，语句格式为

sym('表达式')

例如，定义表达式x+1为符号表达式对象，语句为

sym('x+1');

另一种创建符号表达式的方法是先定义符号变量，然后直接写出符号表达式。例如，在

sin(t)e?2t?5MATLAB中创建符号表达式y?，其MATLAB源程序为 2cos(t)?t?1syms t

7

y=(sin(t).*exp(-2*t)+5)./(cos(t)+t.^2+1) y=

(sin(t)*exp(-2*t)+5)/(cos(t)+t^2+1)

3、 MATLAB中常用绘图方法

MATLAB的plot命令是绘制二维曲线的基本函数，它为数据的可视化提供了方便的途径。例如，函数 关于变量x的曲线绘制的语句格式为

plot(x,y)

其中，输出以向量x为横坐标，向量y为纵坐标，且按照向量x、y中元素的排列顺序有序绘制图形。但向量x与y必须拥有相同的长度。

绘制多幅图形的语句格式为 plot(x1,y1,'str1',x2,y2,'str2',...)

其中，用str1制定的方式，输出以x1为横坐标、y1为纵坐标的图形。用str2制定的方式，输出以x2为横坐标、y2为纵坐标的图形。若省略str，则MATLAB自动为每条曲线选择颜色与线型。

图形完成后，可以通过几个命令来调整显示结果。如grid on用来显示格线；axis([xmin,xmax,ymin,ymax])函数调整坐标轴的显示范围。其中，括号内的“，”可用空格代替；xlabel和ylabel命令可为横坐标和纵坐标加标注，标注的字符串必须用单引号引起来；title命令可在图形顶部加注标题。

用subplot命令可在一个图形窗口中按照规定的排列方式同时显示多个图形，方便图形的比较。其语句格式为

subplot(m,n,p)

或者

subplot(mnp)

其中，m和n表示在一个图形窗口中显示m行n列个图像，p表示第p个图像区域，即在第p个区域作图。

除了plot命令外，MATLAB提供了ezplot命令绘制符号表达式的曲线，其语句格式为

ezplot(y, [a,b])

其中，[a，b]参数表示符号表达式的自变量取值范围，默认值为[0,2?] 4、 MATLAB中的M文件

MATLAB是解释型语言，也就是说在MATLAB命令行中输入的命令在当前MATLAB进程中被解释运行，无需编译和链接等。MATLAB文件分为两类：M脚本文件（M-Script）和M函数（M-function），它们均为由ASCII码构成的文件，该文件可直接在文本编辑器中编写，称为M文件，保存的文件扩展名是.m。

8

M脚本文件包含一族由MATLAB语言所支持的语句，并保存为M文件。它类似于DOS下的批处理文件，不需要在其中输入参数，也不需要给出输出变量来接受处理结果。脚本仅是若干命令或函数的集合，用于执行特定的功能。其执行方式很简单，用户只需在MATLAB的提示符下键入该M文件的文件名，这样MATLAB就会自动执行该M文件中的各条语句，并将结果直接返回到MATLAB工作空间中。脚本M文件实际上是一系列MATLAB命令的集合，它的作用与在MATLAB命令窗口输入的一系列命令等效。

M函数文件不同于M脚本文件，是一种封装结构，通过外界提供输入量而得到函数文件的输出结果。函数是接受入口参数返回出口参数的M文件，程序在自己的工作空间中操作变量，与工作空间分开，无法访问。M函数文件和M脚本文件都是在编辑器中生成，通常以关键字function引导“函数声明行”，并罗列出函数与外界联系的全部“标称”输入输出宗量。它的一般形式为

function [output 1, output 2,…] = functionname(input1, input2,…) %[output 1, output 2,…] = functionname(input1, input2,…) Functionname %Some comments that explain what the function does go here. MATLAB command 1; MATLAB command 2; MATLAB command 3; ……

该函数的M文件名是functionname.m，在MATLAB命令窗口中可被其他M文件调用，例如，

[output1, output2] = functionname(input1, input2)

注意，MATLAB忽略了“%”后面的所有文字，因此，可以利用该符号写注释。以“；”结束一行可以停止输出打印，在一行的最后输入“…”可以续行，以便在下一行继续输入指令。M函数格式是MATLAB程序设计的主流，在一般情况下，不建议使用M脚本文件格式编程。

5、 MATLAB程序流程控制

MATLAB与其他高级编程语言一样，是一种结构化的编程语言。MATLAB程序流程控制结构一般可分为顺序结构、循环结构、以及条件分支结构。MATLAB中实现顺序结构的方法非常简单，只需将程序语句按顺序排列即可。在MATLAB中，循环结构可以由for语句循环结构和while语句循环结构两种方式来实现。条件分支结构可以由if语句分支结构和

9

switch语句分支结构两种方式来实现。下面我们主要介绍这几种程序流程控制。

① for循环结构

for循环结构用于在一定条件下多次循环执行处理某段指令，其语法格式为 for 循环变量＝初值：增量：终值

循环体 end

循环变量一般被定义为一个向量，这样循环变量从初值开始，循环体中的语句每被执行一次，变量值就增加一个增量，直到变量等于终值为止。增量可以根据需要设定，默认时为1。end代表循环体的结束部分。

② while循环结构

while循环结构也用于循环执行处理某段指令，但是与for循环结构不同的是在执行循环体之前先要判断循环执行的条件是否成立，即逻辑表达式为“真”还是“假”，如果条件成立，则执行；如果条件不成立，则终止循环。其语法格式为

while 逻辑表达式

循环体

end

③ if分支结构

if条件分支结构是通过判断逻辑表达式是否成立来决定是否执行制定的程序模块。其语法格式有两种，一种是单分支结构；另一种为多分支结构。其中，单分支结构语法格式为

if 逻辑表达式

程序模块 end

单分支结构语法格式的含义是，如果逻辑表达式为“真”，则执行程序模块，否则跳过该分支结构，按顺序结构执行下面的程序。

多分支结构的语法格式为 if 逻辑表达式1

程序模块1

else if 逻辑表达式2 （可选）

程序模块2

…… else

程序模块n end

多分支结构语法格式可理解为：首先判断if条件分支结构中的逻辑表达式1是否成立，如果成立则执行程序模块1；否则继续判断else if条件分支结构中的逻辑表达式2，如果成立则执行程序模块2；依次下去，如果结构中所有条件都不成立，则执行程序模块n。

④ switch分支结构

switch分支结构是根据表达式的取值结果不同来选择执行的程序模块，其语法格式为 switch 表达式

10

case 常量1

程序模块1 case 常量2

程序模块2

……

otherwise

程序模块n

end

其中，switch后面的表达式可以是任何类型，如数字、字符串等。当表达式的值与case后面的常量相等时，就执行对应的程序模块；如果所有常量都与表达式的值不等时，则执行otherwise后面的程序模块。

除了上述介绍的几种程序流程控制结构外，MATLAB为实现交互控制程序流程还提供了continue、break、pause、input、error、disp等命令。读者可通过doc或者help命令查看它们的具体使用。 6、 MATLAB中信号生成

信号按照自变量的取值是否连续可分为连续时间信号和离散时间信号。对信号进行时域分析，首先需要将信号随时间变化的规律用二维曲线表示出来。对于简单信号可以通过手工绘制其波形，但对于复杂的信号，手工绘制信号波形显得十分困难，且难以绘制精确的曲线。 在MATLAB中通常用三种方法来产生并表示信号，即（1）用MATLAB软件的funtool符合计算方法（图示化函数计算器）来产生并表示信号；（2）用MATLAB软件的信号处理工具箱（Signal Processing Toolbox）来产生并表示信号；（3）用MATLAB软件的仿真工具箱Simulink中的信号源模块。 ① funtool工具生成方法

在MATLAB环境下输入指令funtool，则回产生三个视窗。即

figure No.1：可轮流激活，显示figure No.3的计算结果。 figure No.2：可轮流激活，显示figure No.3的计算结果。

figure No.3：函数运算器，其功能有：f，g可输入函数表达式；x是自变量，在缺省时在[－2pi，2pi]的范围内；自由参数是a；在分别输入完毕后，按下面四排的任一运算操作键，则可在figure No.1或figure No.2产生相应的波形。 ② 信号处理工具箱函数生成方法

一种是用向量来表示信号，另一种则是用符合运算的方法来表示信号。用适当的MATLAB语句表示信号后，可以利用MATLAB的绘图命令绘制出直观的信号波形。常用预定义的信号生成函数：

a) 正余弦函数：sin、cos b) 随机函数： randn

11

c) 周期方波：square d) 周期锯齿波：sawtooth

e) Sinc函数： sinc

f)

Dirichlet函数： diric

以上函数的具体使用调用方法请参考MATLAB帮助文档。 ③ Simulink工具生成方法

在MATLAB的命令视窗下输入simulink指令，则会打开untitled和library simulink两个视窗。library simulink有7个子库，其中source是信号源子库，sinks

是显示器子库。子库中的任何模块都可以拖动到untitled视窗中，用鼠标把模块用连线按输入输出关系连接起来，就构成了仿真系统。在untitled视窗的菜单选simulation中的start，开始进行仿真，仿真执行完毕后，示波器上会显示出信号波形。

三、实验内容及步骤

在熟悉了MATLAB基本命令的基础上，完成以下实验。

（1）、 组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，

D=A-B，E=A .*B，F=A./B，G=A.^B并用画出A、B、C、D、E、F、G图形。 （2）、 用MATLAB实现信号x(t)?Sa(t)?sint(t?[?3?,3?])，并绘出函数曲线。 t（3）、 绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注

?a?x?t??sin?2?t?0?t?10s

?b?x?t??cos?100?t?sin(?t)（a）x?n??0.8（b）x?n??en0?t?4s（4）、 用MATLAB实现下列序列:

0?n?15

?0.2?3j?n0?n?15

0?n?15

（c）x?n??3cos?0.125?n?0.2???2sin?0.25?n?0.1??(d) 将(c)中的x(n)扩展为以16为周期的函数 x16?n??x?n?16? , 绘出四个周期。 (e) 将(c)中的x(n)扩展为以10为周期的函数x10?n??x?n?10?，绘出四个周期。 （5）、 利用MATLAB信号处理工具箱函数及Simulink工具分别绘出 (a) x?t??sin??t?

12

(b) 叠加随机噪声的x?t??sin??t??rand(t)，其中rand(t)为随机噪声生成函数。 (c) 信号x(t)?sinc(t)?sin??t?

四、思考题

比较实验内容第(3)题中（d）和（e）两小题的结果，试说明对于周期性信号，应当如何采样，（d）能保证周期扩展后与原信号保持一致?

五、实验报告要求

1、 简述实验目的及原理。 2、 按实验步骤附上实验程序。 3、 按实验步骤附上有关波形曲线。 4、 简要回答思考题。

13

实验二 连续时间系统时域分析实验

一、 实验说明

1. 实验类型：验证性实验 2. 实验课时：3学时

二、 实验目的

1. 2. 3. 4. 5. 6.

深刻理解卷积运算，利用离散卷积实现连续卷积运算； 深刻理解信号与系统的关系；

加深理解连续系统的冲激响应过程及原理。

掌握用MATLAB对连续系统的冲激响应进行仿真，观测响应波形； 学会用MATLAB求解连续系统冲激响应及阶跃响应； 初步学会利用SIMULINK的仿真。

三、实验原理

1. 利用离散卷积实现连续卷积 （1） 离散卷积和 调用函数：conv

S?conv(f1,f2)?i????f1(i)f(k?i)为离散卷积和，

?其中，f1(k), f2 (k) 为离散序列，K=…-2, -1, 0 , 1, 2, …。但是，conv函数只给出纵轴的序列值的大小，而不能给出卷积的X轴序号。为得到该值，进行以下分析：

对任意输入，设f1(k)非零区间n1~n2，长度L1=n2-n1+1；f2(k)非零区间m1~m2，长度L2=m2-m1+1。则：s(k)?f1(k)*f2(k)非零区间从n1+m1开始，长度为L=L1+L2-1，所以S（K）的非零区间为：n1+m1~ n1+m1+L-1。 （2） 连续卷积和离散卷积的关系

计算机本身不能直接处理连续信号，只能由离散信号进行近似： 设一系统（LTI）输入为P?(t)，输出为h?(t)，如图所示。

P?(t) P?(t) h?(t) 1 ?LTI t ?

14

P?(t)?h?(t)

?(t)?limP?(t)?limh?(t)?h(t)

??0??0若输入为f(t):

f(t)?f?(t)?得输出：

k????f(?k)P(t?k?)?

??y?(t)?k????f(?k)h??(t?k?)?

??当??0时：f(t)?limf?(t)?lim??0??0k????f(?k)P(t?k?)???f(?)?(t??)d?

????

y(t)?limy?(t)?lim?f(?k)h?(t?k?)????0??0k???????f(?)h(t??)d?

所以：

s(t)?f1(t)*f2(t)??f1(?)f2(t??)d??lim??0?f

1(k?)f2(t?k?)?如果只求离散点上的f 值f(?n)

f(n?)?k?????f?1(k?)f2(n??k?)?

???f1(k?)f2[(n?k)?]k???所以，可以用离散卷积和conv函数求连续卷积，只需?足够小以及在卷积和的基础上乘以

?。

（3） 连续卷积坐标的确定

设f1(t)非零值坐标范围：t1~t2，间隔P f2(t)非零值坐标范围：tt1~tt2，间隔P

s(t)?f1(t)*f2(t)非零值坐标：t1+tt1~t2+tt2+1

15

2. 连续系统的冲激响应

0?冲激函数是一种特殊的函数, 它的值在t≠0处均为零, 且有

0_??(t)dt?1, 因此以冲激

函数作为输入, 可以看成电路在t=0_到0+的时间区间内受到了激励, 从而使储能元件得到能量. 在t=0+以后, δ(t)为零, 此时电路中得电压、电流相当于零输入响应。

一个LTI系统，当其初始状态为零时，输入为单位电位冲激函数δ(t)所引起的响应称为单位冲激响应，简称冲激响应，用h(t)表示。如图2-1所示：

δ(tLTI系统 图2.1 冲激响应示意图

对LTI系统，设输入信号为f(t),冲激响应为h(t)，零状态响应为y(t)，则有：

y(t)=f(t)*h(t)

即h(t)属于连续系统的固有特性，与系统的输入无关。只要知道了系统的冲激响应，即可求得系统在不同输入时产生的输出。可见，求解系统的冲激响应对连续系统的分析具有非常重要的意义。

? MATLAB中利用信号处理工具箱函数实现连续系统响应方法 设微分方程：

h(t) ?ayii?0N(i)(t)??bjf(j)(t)

j?0Ma?[aNb?[bMaN?1bM?1aN?2bM?2??a1b1a0]b0] 均为降幂顺序，则：

① 冲激响应为：impulse(b,a) impulse(b,a,t) impulse(b,a,t1:p:t2) y=impulse( ) ② 阶跃响应为：step( ) ③ 零状态响应：lism(b,a,x,t)

? MATLAB中利用SIMULINK来实现连续系统冲激响应方法 ① 打开MATLAB应用程序，然后在命令窗口输入SIMULINK,回车；

16

② 打开FILE 菜单的NEW子菜单下的MODEL 如图并保存命名，在此记做exm 在SIMULINK中commonly used blocks模块下查找 scope,拖入exm; 在SIMULINK中source模块下查找 step,拖入exm;

在SIMULINK中continuous模块下查找deravitive 和integrator（两个） 分别拖入exm 在SIMULINK中math operator 模块下查找 add 和gain (两个) 分别拖入exm 然后用线条连接成下图的形式

17

注意需要双击模块设置如下：

Step模块：step time=0,final value=1,这样的信号才是u(t) Gain模块：分别设置为3 和2；

Add模块：参数Icon shape 设置为“rectangular” list of signs 设置为 “+--“ ；

Scope模块：双击Scope,出现的Scope窗口中选第二个按钮（Scope parameters），设置number

of axes 为3。

③ 双击exm 的运行按钮；在双击示波器得到如下图形

从上图可以很清晰的看到冲激响应的波形

18

冲激响应的方程如下：y??(t)?3y?(t)?2y(t)?f(t) 如图

f()t??32?y(t)

四、实验内容及步骤

1. 完成f1(t)与f2(t)两函数的卷积运算，f1(t)?e?2tu(t),f2(t)?u(t)?u(t?4)在一

个图形窗口中，画出f1(t)、f2(t)以及卷积结果。要求每个坐标系有标题、坐标轴名称。

2. 若系统模型为：

y??(t)?4y?(t)?4y(t)?f?(t)?3f(t) 其中 f(t)?e?tu(t) 求零状态响应，画出波形（函数本身画出一幅图，自己再画出一幅输入波形图）。 3. 若两个连续系统的微分方程分别为：

2y??(t)?7y?(t)?3y(t)?f(t) y??(t)?3y?(t)?2y(t)?f??(t)?f(t)

（1）调用MATLAB的impulse函数编程求解系统的冲激响应。通过不同的参数t绘制出不同时间范围内的波形。 （2）观察分析响应波形。

（3）调用step(b,a,t)对阶跃响应进行仿真观测。 4. 利用MATLAB中simulink仿真工具实现上题。

五．实验报告

1. 绘出给定系统冲激响应波形。 2. 分析说明冲激响应过程。

19

实验三 离散时间系统时域分析及稳定性实验

一、实验说明

1、 实验类型：验证性实验 2、 实验课时：2学时

二、实验目的

1、 掌握求系统响应的方法。

2、 掌握时域离散系统的时域特性。 3、 分析、观察及检验系统的稳定性。

三、实验原理与方法

在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。在计算机上适合用递推法求差分方程的解，最简单的方法是采用MATLAB语言的工具箱函数filter函数。也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。

系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。

系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。

实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的[19]。系统的稳态输出是指当n??时，系统的输出。如果系统稳定，信号加入系统后，系统输出的开始一段称为暂态效应，随n的加大，幅度趋于稳定，达到稳态输出。

注意在以下实验中均假设系统的初始状态为零。

四、实验内容及步骤

1、编制程序，包括产生输入信号、单位脉冲响应序列的子程序，用filter函数或conv函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。 2、给定一个低通滤波器的差分方程为

y(n)?0.05x(n)?0.05x(n?1)?0.9y(n?1) 输入信号 x1(n)?R8(n) x2(n)?u(n)

a) 分别求出系统对x1(n)?R8(n)和x2(n)?u(n)的响应序列，并画出其波形。 b) 求出系统的单位冲响应，画出其波形。 3、给定系统的单位脉冲响应为

h1(n)?R10(n)

h2(n)??(n)?2.5?(n?1)?2.5?(n?2)??(n?3)

20

用线性卷积法分别求系统h1(n)和h2(n)对x1(n)?R8(n)的输出响应，并画出波形。 4、给定一谐振器的差分方程为

y(n)?1.8237y(n?1)?0.9801y(n?2)?b0x(n)?b0x(n?2) 令 b0?1/100.49，谐振器的谐振频率为0.4rad。

a) 用实验方法检查系统是否稳定。输入信号为u(n)时，画出系统输出波形。 b) 给定输入信号为

x(n)?sin(0.014n)?sin(0.4n) 求出系统的输出响应，并画出其波形。

五、思考题

1、 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？

2、 如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面第一

个实验结果进行分析说明。

六、实验报告要求

1、 简述在时域求系统响应的方法。

2、 简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。 3、 对各实验所得结果进行简单分析和解释。 4、 简要回答思考题。

5、 打印程序清单和要求的各信号波形。

21

实验四 信号的频谱分析实验

一、实验说明

1、 实验类型：验证性实验 2、 实验课时：3学时

二、实验目的

1、 熟悉理解连续信号的时域采样和连续信号时域采样前后频谱变化关系； 2、 熟悉理解连续频谱的频域采样和频域采样点数的选择；

3、 学习掌握FFT对连续信号和离散信号进行频谱分析的方法；

4、 了解FFT对信号进行频谱分析过程中可能出现的分析误差及其原因；

三、实验原理

1、时域采样

?(j?)是① 对模拟信号xa(t)以间隔T进行时域等间隔理想采样，形成的采样信号的频谱X原模拟信号频谱Xa(j?)以采样角频率?s（?s?2?/T）为周期进行周期延拓。公式为：

1???a(t)]??Xa(j??jn?s) Xa(j?)?FT[xTn???② 采样频率?s必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的 频谱不产生频谱混叠。

利用计算机计算上式并不方便，下面我们导出另外一个公式，以便用计算机上进行实验。

?a(t)和模拟信号xa(t)之间的关系为： 理想采样信号x?a(t)?xa(t) xn?????(t?nT)

??对上式进行傅立叶变换，得到：

?(j?)??[x(t)?(t?nT)]e?j?tdt Xa?a???n??? ＝n??????????xa(t)?(t?nT)e?j?tdt

在上式的积分号内只有当t?nT时，才有非零值，因此：

?(j?)? Xan????j?nT x(nT)e?a上式中，在数值上xa(nT)＝x(n)，再将???T代入，得到：

22

?(j?)? Xan????x(n)e??j?n

上式的右边就是序列的傅立叶变换X(ej?)，即

?(j?)?X(ej?)Xa???T

上式说明理想采样信号的傅立叶变换可用相应的采样序列的傅立叶变换得到，只要将自变量

ω用?T代替即可。 2、频域采样

① 对信号x(n)的频谱函数X(ej?)在[0，2π]上等间隔采样N点，得到

XN(k)?X(ej?)??2?kN , k?0,1,2,?,N?1

则N点IDFT[XN(k)]得到的序列就是原序列x(n)以N为周期进行周期延拓后的主值区序列，公式为：

xN(n)?IDFT[XN(k)]N?[i????x(n?iN)]R?N(n)

② 由上式可知，频域采样点数N必须大于等于时域离散信号的长度M(即N≥M)，才能使

时域不产生混叠，则N点IDFT[XN(k)]得到的序列xN(n)就是原序列x(n),即

xN(n)=x(n)。如果N>M，xN(n)比原序列尾部多N-M个零点；如果N

短，因此。xN(n)与x(n)不相同。

3、FFT对信号进行频谱分析

用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2?/N，因此要求2?/N?D。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。 对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

四、实验内容及步骤

1、连续信号的FFT频谱分析

给定模拟信号，xa(t)?Ae??tsin(?0t)u(t) ，式中A=444.128，?=502π，

23

?0=502πrad/s，它的幅频特性曲线如图4.1

图4.1

xa(t)的幅频特性曲线

要求按照xa(t)的幅频特性曲线，选取三种采样频率，即Fs=1kHz，300Hz，200Hz。观测时间选Tp?50ms。

为使用FFT，首先用下面公式产生时域离散信号，对三种采样频率，采样序列按顺序用

x1(n)，x2(n)，x3(n)表示。

x(n)?xa(nT)?Ae??nTsin(?0nT)u(nT)

因为采样频率不同，得到的序列x1(n)，x2(n)，x3(n)的长度不同，长度（点数）用公式N?Tp?Fs计算。选FFT的变换点数为M=64，序列长度不够64的尾部加零。

X(k)=FFT[x(n)] ， k=0,1,2,3,-----,M-1 式中k代表的频率为 ?k?2?k。 M要求： 编写实验程序，计算x1(n)、x2(n)和x3(n)的幅度特性，并绘图显示。观察分析频谱混叠失真现象。

2、离散信号的FFT频谱分析

对以下非周期和周期序列进行谱分析。

?n?1,0?n?3?x1(n)??8?n,4?n?7

?0,othern??4?n,0?n?3?x2(n)??n?3,4?n?7

?0,othern?x3(n)?cos?

n4x4(n)?cos(?n/4)?cos(?n/8)

24

选择FFT的变换区间N为8和16 两种情况分别对以上序列进行频谱分析。分别打印其幅频特性曲线。并进行对比、分析和讨论。

五、思考题

1、对于周期序列，如果周期不知道，如何用FFT进行谱分析？ 2、如何选择FFT的变换区间？（包括非周期信号和周期信号）

3、当N=8时，x1(n)和x2(n)的幅频特性会相同吗？为什么？N=16 呢？

六、实验报告要求

1、 完成各个实验任务和要求。

2、 分析比较实验结果，简述由实验得到的主要结论 3、 简要回答思考题。

4、 附上程序清单和有关曲线图。

25

实验五 数字滤波器设计与仿真实验

一、实验说明

1、 实验类型：设计性实验 2、 实验课时：3学时

二、实验目的

1、 熟悉用数字滤波器设计的原理与方法；

2、 学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具

FDATOOL）设计各种类型数字滤波器，学会根据滤波需求确定滤波器指标参数。 3、 掌握数字滤波器的MATLAB实现方法。

4、 通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。

三、实验原理

1、 IIR数字滤波器设计原理

设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求调用如上函数直接设计IIR数字滤波器。 2、 FIR数字滤波器设计原理 设计FIR数字滤波器一般采用直接法（窗函数法、频率采样法和等波纹最佳逼近法），本实验采用的是窗函数法。基本设计过程是：①先根据给定的数字滤波器的指标选择合适的窗函数并计算窗口宽度N； ②调用MATLAB函数fir1设计一个FIR滤波器；③调用MATLAB快速卷积函数fftfilt实现对输入信号的滤波。MATLAB信号处理工具箱中的fir1函数可以实现窗函数法和频率采样法，remezord和remez函数可以实现FIR数字滤波器。本实验要求调用如上函数设计FIR数字滤波器。

四、实验内容及步骤

1、 调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数

还会自动绘图显示st的时域波形和幅频特性曲线，如图5.1所示。由图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。

图5.1 三路调幅信号st的时域波形和幅频特性曲线

26

2、 要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st

中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为60dB。

提示：抑制载波单频调幅信号的数学表示式为

1s(t)?cos(2?f0t)cos(2?fct)?[cos(2?(fc?f0)t)?cos(2?(fc?f0)t)]

2其中，cos(2?fct)称为载波，fc为载波频率，cos(2?f0t)称为单频调制信号，f0为调制正弦波信号频率，且满足fc?f0。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频fc?f0和差频fc?f0，这2个频率成分关于载波频率fc对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率fc对称的2根谱线，其中没有载频成分，故取名为抑制载波单频调幅信号。容易看出，图5.1中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。如果调制信号m(t)具有带限连续频谱，无直流成分，则s(t)?m(t)cos(2?fct)就是一般的抑制载波调幅信号。其频谱图是关于载波频率fc对称的2个边带（上下边带），在专业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号,简称双边带 (DSB) 信号。如果调制信号m(t)有直流成分，则s(t)?m(t)cos(2?fct)就是一般的双边带调幅信号。其频谱图是关于载波频率fc对称的2个边带（上下边带），并包含载频成分。

3、 编程序要求采用IIR的双线性变换法和FIR的窗函数法两类滤波器设计方法分别设计

低通、带通和高通三个滤波器，并绘图显示其幅频响应特性曲线。

4、 利用所设计的滤波器对信号产生函数mstg产生的信号st进行滤波处理，分离出st中

的三路不同载波频率的调幅信号y1(n)、y2(n)和y3(n)，并绘图显示y1(n)、y2(n)和y3(n)的时域波形，观察分离效果。

五、信号产生函数mstg清单

function st=mstg

%产生信号序列向量st,并显示st的时域波形和频谱

%st=mstg 返回三路调幅信号相加形成的混合信号，长度N=1600 N=1600 %N为信号st的长度。

Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp为采样时间 t=0:T:(N-1)*T;k=0:N-1;f=k/Tp; fc1=Fs/10; %第1路调幅信号的载波频率fc1=1000Hz, fm1=fc1/10; %第1路调幅信号的调制信号频率fm1=100Hz fc2=Fs/20; %第2路调幅信号的载波频率fc2=500Hz fm2=fc2/10; %第2路调幅信号的调制信号频率fm2=50Hz fc3=Fs/40; %第3路调幅信号的载波频率fc3=250Hz, fm3=fc3/10; %第3路调幅信号的调制信号频率fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号 xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号 xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号

27

st=xt1+xt2+xt3; %三路调幅信号相加 fxt=fft(st,N); %计算信号st的频谱

%====以下为绘图部分，绘制st的时域波形和幅频特性曲线====== subplot(3,1,1)

plot(t,st);grid;xlabel('t/s');ylabel('s(t)');

axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形') subplot(3,1,2)

stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b) s(t)的频谱') axis([0,Fs/5,0,1.2]);

xlabel('f/Hz');ylabel('幅度')

六、思考题

1、 请阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。

2、 信号产生函数mstg中采样点数N=800，对st进行N点FFT可以得到6根理想谱线。

如果取N=1000，可否得到6根理想谱线？为什么？N=2000呢？请改变函数mstg中采样点数N的值，观察频谱图验证您的判断是否正确。

3、 修改信号产生函数mstg，给每路调幅信号加入载波成分，产生调幅（AM）信号，重

复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别。

提示：AM信号表示式：s(t)?[1?cos(2?f0t)]cos(2?fct)。

4、 对比在同样的数字滤波器技术指标下，IIR类型数字滤波器和FIR数字滤波器的阶次

和信号分离效果。

七、实验报告要求

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

2、 画出实验主程序框图，打印程序清单。 3、 绘制三个分离滤波器的损耗函数曲线。

4、 绘制经过滤波分理出的三路调幅信号的时域波形。 5、 分析总结实验结果。 6、 简要回答思考题。

28

实验六 数字信号处理在双音多频拨号系统中的应用

一、实验说明

1、 实验类型：综合性实验 2、 实验课时：4学时

二、实验原理和内容

1、 双音多频拨号系统简介

双音多频（Dual Tone Multi Frequency, DTMF）信号是音频电话中的拨号信号，由美国AT&T贝尔公司实验室研制，并用于电话网络中。这种信号制式具有很高的拨号速度，且容易自动监测识别，很快就代替了原有的用脉冲计数方式的拨号制式。这种双音多频信号制式不仅用在电话网络中，还可以用于传输十进制数据的其它通信系统中，用于电子邮件和银行系统中。这些系统中用户可以用电话发送DTMF信号选择语音菜单进行操作。

DTMF信号系统是一个典型的小型信号处理系统，它要用数字方法产生模拟信号并进行传输，其中还用到了D/A变换器；在接收端用A/D变换器将其转换成数字信号，并进行数字信号处理与识别。为了系统的检测速度并降低成本，还开发一种特殊的DFT算法，称为戈泽尔(Goertzel)算法，这种算法既可以用硬件（专用芯片）实现，也可以用软件实现。下面首先介绍双音多频信号的产生方法和检测方法，包括戈泽尔算法，最后进行模拟实验。下面先介绍电话中的DTMF信号的组成。

在电话中，数字0~9的中每一个都用两个不同的单音频传输，所用的8个频率分成高频带和低频带两组，低频带有四个频率：679Hz,770Hz,852Hz和941Hz；高频带也有四个频率：1209Hz,1336Hz,1477Hz和1633Hz.。每一个数字均由高、低频带中各一个频率构成，例如1用697Hz和1209Hz两个频率，信号用sin(2?f1t)?sin(2?f2t)表示，其中f1?679Hz，

f2?1209Hz。这样8个频率形成16种不同的双频信号。具体号码以及符号对应的频率如

表6.1所示。表中最后一列在电话中暂时未用。

表6.1 双频拨号的频率分配

列 行 697Hz 770Hz 852Hz 942Hz 1209Hz 1 4 7 * 1336Hz 2 5 8 0 1477Hz 3 6 9 # 633Hz A B C D DTMF信号在电话中有两种作用，一个是用拨号信号去控制交换机接通被叫的用户电话机，另一个作用是控制电话机的各种动作，如播放留言、语音信箱等。 2、 电话中的双音多频（DTMF）信号的产生与检测 ① 双音多频信号的产生

假设时间连续的 DTMF信号用x(t)?sin(2?f1t)?sin(2?f2t)表示，式中f1和f2是按照表10.10.1选择的两个频率，f1代表低频带中的一个频率，f2代表高频带中的一个频率。

29

显然采用数字方法产生DTMF信号，方便而且体积小。下面介绍采用数字方法产生DTMF信号。规定用8KHz对DTMF信号进行采样，采样后得到时域离散信号为 x(n)?sin(2?f1n/8000)?sin(2?f2n/8000)

形成上面序列的方法有两种，即计算法和查表法。用计算法求正弦波的序列值容易，但实际中要占用一些计算时间，影响运行速度。查表法是预先将正弦波的各序列值计算出来，寄存在存储器中，运行时只要按顺序和一定的速度取出便可。这种方法要占用一定的存储空间，但是速度快。

因为采样频率是8000Hz，因此要求每125ms输出一个样本，得到的序列再送到D/A变换器和平滑滤波器，输出便是连续时间的DTMF信号。DTMF信号通过电话线路送到交换机。

② 双音多频信号的检测

在接收端，要对收到的双音多频信号进行检测，检测两个正弦波的频率是多少，以判断所对应的十进制数字或者符号。显然这里仍然要用数字方法进行检测，因此要将收到的时间连续 DTMF信号经过A/D变换，变成数字信号进行检测。检测的方法有两种，一种是用一组滤波器提取所关心的频率，根据有输出信号的2个滤波器判断相应的数字或符号。另一种是用DFT（FFT）对双音多频信号进行频谱分析，由信号的幅度谱，判断信号的两个频率，最后确定相应的数字或符号。当检测的音频数目较少时，用滤波器组实现更合适。FFT是DFT的快速算法，但当DFT的变换区间较小时，FFT快速算法的效果并不明显，而且还要占用很多内存，因此不如直接用DFT合适。下面介绍Goertzel算法，这种算法的实质是直接计算DFT的一种线性滤波方法。这里略去Goertzel算法的介绍，可以直接调用MATLAB信号处理工具箱中戈泽尔算法的函数Goertzel，计算N点DFT的几个感兴趣的频点的值。 3、 检测DTMF信号的DFT参数选择

用DFT检测模拟DTMF信号所含有的两个音频频率，是一个用DFT对模拟信号进行频谱分析的问题。根据第三章用DFT对模拟信号进行谱分析的理论，确定三个参数：（1）采样频率Fs，（2）DFT的变换点数N，（3）需要对信号的观察时间的长度Tp。这三个参数不能随意选取，要根据对信号频谱分析的要求进行确定。这里对信号频谱分析也有三个要求： a. 频率分辨率，b. 谱分析的频谱范围，c. 检测频率的准确性。

① 频谱分析的分辨率。

观察要检测的8个频率，相邻间隔最小的是第一和第二个频率，间隔是73Hz，要求DFT最少能够分辨相隔73Hz的两个频率，即要求Fmin?73Hz。DFT的分辨率和对信号的观察时间Tp有关，Tpmin?1/F?1/73?13.7ms 。考虑到可靠性，留有富裕量，要求按键的时间大于40ms。

② 频谱分析的频率范围

要检测的信号频率范围是697～1633Hz，但考虑到存在语音干扰，除了检测这8个频率外，还要检测它们的二次倍频的幅度大小，波形正常且干扰小的正弦波的二次倍频是很小的，如果发现二次谐波很大，则不能确定这是DTMF信号。这样频谱分析的频率范围为697～3266Hz。按照采样定理，最高频率不能超过折叠频率，即0.5Fs?3622Hz，由此要求最小的采样频率应为7.24KHz。因为数字电话总系统已经规定Fs＝8KHz，因此对频谱分析范围

30

的要求是一定满足的。按照Tpmin?13.7ms，Fs＝8KHz，算出对信号最少的采样点数为

Nmin?Tpmin?Fs?110。

③ 检测频率的准确性

这是一个用DFT检测正弦波频率是否准确的问题。序列的N点DFT是对序列频谱函数在0～2?区间的N点等间隔采样，如果是一个周期序列，截取周期序列的整数倍周期，进行DFT，其采样点刚好在周期信号的频率上，DFT的幅度最大处就是信号的准确频率。分析这些DTMF信号，不可能经过采样得到周期序列，因此存在检测频率的准确性问题。 DFT的频率采样点频率为?k?2?k/N（k=0,1,2,---,N-1），相应的模拟域采样点频率为，希望选择一个合适的N，使用该公式算出的fk能接近要检fk?Fsk/N（k=0,1,2,---,N-1）

测的频率，或者用8个频率中的任一个频率

fk?代入公式fk??Fsk/N中时，得到的k值最

接近整数值，这样虽然用幅度最大点检测的频率有误差，但可以准确判断所对应的DTMF频率，即可以准确判断所对应的数字或符号。经过分析研究认为N＝205是最好的。按照Fs＝8KHz，N＝205，算出8个频率及其二次谐波对应k值，和k取整数时的频率误差见表6.2。

表6.2

8个基频 最近的DFT的 二次谐波 对应的 最近的 绝对误差 绝对误差 Hz Hz 整数k值 k值 k值 整数k值 697 770 852 941 1209 1336 1477 1633 17.861 19.531 21.833 24.113 30.981 34.235 37.848 41.846 18 20 22 24 31 34 38 42 0.139 0.269 0.167 0.113 0.019 0.235 0.152 0.154 1394 1540 1704 1882 2418 2672 2954 3266 35.024 38.692 42.813 47.285 60.752 67.134 74.219 82.058 35 39 43 47 61 67 74 82 0.024 0.308 0.187 0.285 0.248 0.134 0.219 0.058 通过以上分析，确定Fs＝8KHz，N＝205，Tp?40ms。

4、 DTMF信号的产生与识别仿真实验

下面先介绍MATLAB工具箱函数goertzel，然后介绍DTMF信号的产生与识别仿真实验程序。Goerztel函数的调用格式额为

Xgk=goertzel(xn,K)

xn是被变换的时域序列，用于DTMF信号检测时，xn就是DTMF信号的205个采样值。 K是要求计算的DFT[xn]的频点序号向量，用N表示xn的长度，则要求1≤K≤N。由表10.2.2可知，如果只计算DTMF信号8个基频时，

K=[18，20，22，24，31，34，38，42]，

如果同时计算8个基频及其二次谐波时，

K=[18，20，22，24，31，34，35，38，39，42，43，47，61，67，74，82]。

Xgk是变换结果向量，其中存放的是由K指定的频率点的DFT[x(n)]的值。设X(k)=

31

DFT[x(n)]，则Xgk(i)?X(K(i)), i?1,2,?,length(K)。

DTMF信号的产生与识别仿真实验在MATLAB环境下进行，编写仿真程序，运行程序，送入6位电话号码，程序自动产生每一位号码数字相应的DTMF信号，并送出双频声音，再用DFT进行谱分析，显示每一位号码数字的DTMF信号的DFT幅度谱，安照幅度谱的最大值确定对应的频率，再安照频率确定每一位对应的号码数字，最后输出6位电话号码。 本实验程序较复杂，所以将仿真程序提供给读者，只要求读者读懂程序，直接运行程序仿真。程序名为exp6。程序分四段：第一段（2—7行）设置参数，并读入6位电话号码；第二段（9—20行）根据键入的6位电话号码产生时域离散DTMF信号，并连续发出6位号码对应的双音频声音；第三段（22—25行）对时域离散DTMF信号进行频率检测，画出幅度谱；第四段（26—33行）根据幅度谱的两个峰值，分别查找并确定输入6位电话号码。根据程序中的注释很容易分析编程思想和处理算法。程序清单如下： %实验6程序：exp6.m

% DTMF双频拨号信号的生成和检测程序 %clear all;clc;

tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; % DTMF信号代表的16个数 N=205;K=[18,20,22,24,31,34,38,42];

f1=[697,770,852,941]; % 行频率向量 f2=[1209,1336,1477,1633]; % 列频率向量 TN=input('键入6位电话号码= '); % 输入6位数字

TNr=0; %接收端电话号码初值为零 for l=1:6;

d=fix(TN/10^(6-l)); TN=TN-d*10^(6-l); for p=1:4; for q=1:4;

if tm(p,q)==abs(d); break,end % 检测码相符的列号q end

if tm(p,q)==abs(d); break,end % 检测码相符的行号p end

n=0:1023; % 为了发声，加长序列 x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号 sound(x,8000); % 发出声音 pause(0.1)

% 接收检测端的程序

X=goertzel(x(1:205),K+1); % 用Goertzel算法计算八点DFT样本 val = abs(X); % 列出八点DFT向量 subplot(3,2,l);

stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出DFT(k)幅度 axis([10 50 0 120])

limit = 80; % for s=5:8;

if val(s) > limit, break, end % 查找列号

32

end for r=1:4;

if val(r) > limit, break, end % 查找行号 end

TNr=TNr+tm(r,s-4)*10^(6-l); end

disp('接收端检测到的号码为：') % 显示接收到的字符

disp(TNr)

运行程序，根据提示键入6位电话号码123456，回车后可以听见6位电话号码对应的DTMF信号的声音，并输出相应的6幅频谱图如图7.1所示，左上角的第一个图在k=18和k=31两点出现峰值，所以对应第一位号码数字1。最后显示检测到的电话号码123456。

图6.1 6位电话号码123456的DTMF信号在8个近似基频点的DFT幅度

5、 实验主要内容 （1）、 详细阅读以上实验原理，熟悉基本检测步骤。 （2）、 运行仿真程序exp6.m，任意送入6位电话号码，打印出相应的幅度谱。观察程序

运行结果，对照表6.1和表6.2，判断程序谱分析的正确性。 （3）、 分析该仿真程序，将产生、检测和识别6位电话号码的程序改为能产生、检测和识

别8位电话号码的程序，并运行一次，打印出相应的幅度谱和8位电话号码。

五、实验报告

1、 分析程序exp6.m，画出仿真程序流程图。

2、 打印6位和8位电话号码DTMF信号的幅度谱。

3、 简述DTMF信号的参数：采样频率、DFT的变换点数以及观测时间的确定原则。

33

附录：实验报告模版

XXX系实验报告

课程名称：《基于MATLAB的信号与系统及数字信号处理仿真实验》 实验项目名称：XXX

指导教师： 成绩： 实验时间：XXXX.XX.XX 实验地点：信息学院4层机房 班级：电子信息工程/通信工程 姓名：XXX 学号：2005117001

一、实验目的 二、实验原理: 三、实验环境

四、实验结果及分析 五、参考资料 六、实验心得

34

本文来源：https://www.bwwdw.com/article/7ks7.html