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

基于计算机的数字信号处理

实验指导书

党向东 王海涛

沈阳理工大学信息科学与工程学院

二〇〇七年三月

目录

前言……………………………………………………………………...…………….1 第一章 第二章 第三章 第四章 第五章 第六章

MATLAB基础知识 ................................................................................. 1 MATLAB基本数值运算 ......................................................................... 4 MATLAB的图形处理功能 ..................................................................... 8 MATLAB的程序设计 ........................................................................... 11 常用数字信号处理函数 ......................................................................... 16 MATLAB在数字信号处理中的应用 ................................................... 23

实验一 常见离散信号的MATLAB产生和图形显示 ............................................. 33 实验二 离散系统的频率响应分析和零、极点分布 ........................................... 37 实验三 序列线性卷积、圆周卷积的计算及其关系的研究 ............................... 39 实验四 利用DFT分析信号的频谱 ..................................................................... 41 实验五 信号时间尺度变换的研究 ....................................................................... 43 实验六 快速傅里叶变换及其应用 ....................................................................... 47 实验七 IIR滤波器的实现与应用 ........................................................................ 56 实验八 FIR滤波器的实现与应用 ........................................................................ 61

前言

数字信号处理研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息。近年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已经日趋成熟。本课程以计算机为工具，通过一定量的实验项目，以验证所学的概念和算法。由于MATLAB软件的功能十分强大，使用起来也非常方便，在工程技术中尤其是信号处理领域得到了广泛的应用，因此以MATLAB作为本计算机实验课的计算机语言工具。希望大家通过本教材的学习及上机实践，能基本掌握MATLAB程序设计知识，能利用MATLAB进行简单的数字信号处理问题，利用其提供的工具箱能进行滤波器的设计，为理论知识的实用化而奠定基础。

第一章 MATLAB基础知识

§1-1 MATLAB软件简介

MATLAB，Matrix Laboratory的缩写，是由Mathworks公司开发的一套用于科学工程计算的可视化高性能语言，具有强大的矩阵运算能力。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个界面友好的用户环境，在这个环境中，问题与求解都能方便地以数学的语言（主要是矩阵形式）或图形方式表示出来。与大家常用的Fortran 和C等高级语言相比，MATLAB的语法规则更简单，更贴近人的思维方式，被称为“草稿纸式的语言”。 §1-2 MATLAB应用入门

1． MATLAB的安装与卸载

MATLAB软件在用户接口时具有较强的亲和力，其安装过程比较典型，直接运行光盘中的安装向导支撑程序SETUP.exe，按其提示一步步选择即可。MATLAB自身带有卸载程序，在其安装目录下有uninstall子目录，运行该目录下uninstall.exe的即可；也可以通过Windows系统的安装卸载程序进行卸载。 2． MATLAB的启动与退出

MATLAB安装完成后，会自动在Windows桌面上生成一个MATLAB图标，它是指向安装目录下\\bin\\win32\\matlab.exe的链接，双击这个图标即可来到MATLAB集成环境的基本窗口；也可以在开始菜单的程序选项中选择MATLAB快捷方式；还可以在MATLAB的安装路径的bin子目录中双击可执行文件matlab.exe。MATLAB的退出与普通WIN32的程序一样，值得一提的是它有一个自身专有的快捷键Ctrl+Q。初次启动MATLAB后，将进入MATLAB默认设置下的桌面平台如图1-1所示。

1

图1-1 MATLAB默认设置下的桌面平台

3． MATLAB的桌面平台

默认设置下的桌面平台包括6个窗口，分别是MATLAB主窗口、命令窗口(Command Window)、历史窗口Command History)、当前目录窗口Current Directory)、发行说明书窗口(Launch Pad)和工作间管理窗口(Workspace)。 3．1 MATLAB主窗口

MATLAB的其它几个窗口都包含在这个大的主窗口中，主窗口不能进行任何计算任务的操作，只用来进行一些整体的环境参数的设置。主要包括菜单栏（File、Edit、View、Web、Window和Help共6个下拉菜单）、工具栏（10个按钮控件）等。

工具栏各按钮控件及说明如下所示：

3．2 命令窗口(Command Window)

MATLAB的命令窗口如图1-2所示。其中“”为运算提示符，表示MATLAB正处在准备状态。当在提示符后输入一段运算式或命令并按Enter键后，MATLAB

2

将给出计算结果，然后再进入准备状态。

图1-2 MATLAB的命令窗口

3．3 MATLAB常用命令

MATLAB有一些嵌入函数，有时应用这些函数可以起到事半功倍的效果。MATLAB常用的控制命令见表1-1。

表1-1：MATLAB常用命令

命令 功能 显示或改变当前工作目录，与工具栏中 cd 同效 dir clc、home clf clear disp type 列出当前目录或指定目录下的文件和子目录清单，类似于DOS命令DIR 的所有显示内容，并把光标移到命令窗口的左上角 清除MATLAB当前图形窗口中的图形 清除内存中的变量和函数 显示变量的内容 列出指定文件的全部内容，类似于DOS命令TYPE exit、quit 退出MATLAB who whos what which lookfor 列出当前工作空间中的变量 列出当前工作空间中的变量的更多信息 列出当前目录或指定目录下的 .M文件、..MAT文件和 ..MEX文件 显示指定函数或文件的路径 按照指定的关键字查找所有相关的 .M文件

3

第二章 MATLAB基本数值运算

§2-1 MATLAB内部特殊变量和常数

MATLAB内部有很多变量和常数，用以表达特殊含义。常用的有：

⑴ 变量ans：指当前未定义变量名的答案。

⑵ 常数eps：表示浮点相对精度，其值是从1.0到下一个最大浮点数之间的差值。该变量值作为一些MATLAB函数计算的相对浮点相对精度，按IEEE标准，eps=2-52，近似为2.2204e-016。

⑶ 常数Inf：表示无穷大。当输入或计算中有除以0时产生Inf。 ⑷ 虚数单位i、j：表示复数虚部单位，相当于1。 ⑸ NaN：表示不定型值，是由0/0运算产生的。

⑹ 常数pi：表示圆周率π，其值为3.141 592 653 5 7„。 §2-2 变量类型

1．变量命名规则

MATLAB中对变量的命名应遵循以下规则：

1） 变量名可以由字母、数字和下划线混合组成，但必须以字母开头。 2） 字符长度不能大于31。 3） 变量命名区分大小写。

2． 局部变量和全局变量

局部变量是指那些每个函数体内自己定义的，不能从其它函数和MATLAB工作空间访问的变量。

全局变量是指用关键字“global”声明的变量。全局变量名应尽量大写，并能反映它本身的含义。如果需要在工作空间和几个函数中都能访问一个全局变量，必须在工作空间和这几个函数中都声明该变量是全局的。 §2-3 向量及其运算

向量运算是矢量运算的基础，向量也是组成矩阵的基本元素之一。 1． 向量的生成 1．1 直接输入向量

生成向量最直接的方法就是在命令窗口中直接输入。格式上的要求是，向量元素需要用“[ ]”括起来，元素之间可以用空格、逗号或分号分隔；需要注意的是，用空格和逗号分隔生成行向量，用分号分隔生成列向量。 例 A=[1,2,3] 或A=[1 2 3] %生成行向量 A=[1;2;3] %生成列向量 1．2 利用冒号表达式生成向量

4

冒号表达式的基本形式为x=x0：step：xn，其中x0、step、xn分别为给定数值，x0表示向量的首元素数值，xn表示向量尾元素数值限(只有当xn- x0恰为step值的整数倍时, xn才能成为尾数)，step表示从第二个元素开始，元素数值大小与前一个元素数值大小的差值。

例 在命令窗口，给向量a、b、c赋值。

>> a=1:2:12

a =

1 3 5 7 9 11

>> b=12:-2:1 b =

12 10 8 6 4 2 >> c=1:2:13 c =

1 3 5 7 9 11 13

1．3 特殊向量的生成

对于特殊的向量可直接调用MATLAB的函数生成。如y=linsoace(x1,x2,n)用于生成线性等分的n维行向量，使得y(1)=x1，y(n)=x2。

另外，向量还可以从矩阵中提取，还可以把向量看成是1×n阶（行向量）或n×1阶（列向量）的矩阵，以矩阵形式生成。

2． 向量的基本运算

2．1 加（减）与数乘计算

例 >> a=[1,2,3,4];b=[0,1,2,3];c=a-b

c =

1 1 1 1 >> d=a-3 d =

-2 -1 0 1 >> 4*a ans =

4 8 12 16

2．2 对位乘、点积计算

同维向量a与b的对位乘用c=a.*b实现，即c的每一个元素之值是a与b对应元素的乘积。同维向量a与b的点积，一个方法是利用函数dot来实现；另一种方法是先生成a与b的对位乘向量c， 再取c的各元素和即为a与b的点积。 例 >> a.*b

ans =

0 2 6 12

5

>> dot(a,b) %或sum(a.*b) ans =

20 §2-4 矩阵及其运算

MATLAB具有强大的矩阵运算和数据处理功能，对矩阵的处理必须遵从代数规则。

1． 矩阵的生成

（1）一般矩阵的生成

对于一般的矩阵，MATLAB的生成方法有很多种。最简单的方法是从键盘直接输入矩阵元素。直接输入矩阵元素时应注意：各行的元素之间用空格或逗号隔开，行与行之间用分号或回车隔开，用中括号把矩阵所有的元素括起来。 例 在工作空间产生一个3×3矩阵A可用MATLAB语言描述如下：

A=[1，2，3；4，5，6；7，8，9]

或 A=[ 1 2 3

4 5 6 7 8 9]

运行结果为

A=

1 2 3

4 5 6 7 8 9

（2）特殊矩阵的生成

对于特殊的矩阵可直接调用MATLAB的函数生成

用函数zeros生成全0矩阵：格式为B=zeros(m,n)生成m×n的全0阵。 用函数ones生成全1矩阵：格式为B=ones(m,n)生成m×n的全1阵。 用函数eye生成单位阵：格式为B=eye(m,n)生成m×n矩阵，其中对角线元素全为1，其它元素为0。

2．矩阵的运算

矩阵的运算有基本运算和函数运算两种类型。基本运算包括矩阵的加、减、乘、除、幂、转置、逆等，其主要特点是通过MATLAB提供的基本运算符+、－、*、/、^等即可完成。函数运算主要是通过调用MATLAB系统内置的运算函数来求取矩阵的，求秩，求特征值和特征相量，等等。需要时可以参阅联机帮助和相关参考书。

例 矩阵的基本运算

>> a=[1,2,3;4,5,6]; >> b=[6,5,4;3,2,1];

6

>> c=a+b %计算两个矩阵的和 c =

7 7 7 7 7 7

>> d=b' %计算矩阵b的转置 d =

6 3 5 2 4 1

>> e=a*d %做矩阵的乘法，必须满足矩阵乘法的基本要求 e =

28 10 73 28

>> f=det(e) %求矩阵e的行列式 f =

>> g=e^(-1) %求矩阵e的逆 g =

0.5185 -0.1852 -1.3519 0.5185

7

第三章 MATLAB的图形处理功能

从最原始版本的MATLAB开始，图形功能就已经成为基本的功能之一。随着MATLAB版本的逐步升级，MATLAB的图形工具箱从简单的点、线、面处理发展到了集二维图形、三维图形甚至四维表现图和对图形进行着色、消隐、光照处理、渲染及多视角处理等多项功能于一身的强大功能包。这里只简单讨论二维基本绘图命令及图形修饰命令。常用的绘图语句有figure、plot、subplot、stem等，图形修饰语句有title、axis、text等。

1．figure命令

figure有两种用法，只用一句figure命令，会创建一个新的图形窗口，并返回一个整数型的窗口编号。figure(n)表示将第n号图形窗口作为当前的图形窗口，并将其显示在所有窗口的最前面；如果该图形窗口不存在，则新建一个窗口，并赋以编号n。

2． plot命令

线性绘图函数。用法为plot(x, y,’s’)。参数x为横轴变量，y为纵轴变量，s用以控制图形的基本特征如颜色、粗细等的图形设置选项，通常可以省略。MATLAB语言中提供的图形控制符如表3-1所示。

表3-1：MATLAB语言中的图形设置选项

参数 y m c g b w k r 含义 黄色 紫色 青色 绿色 蓝色 白色 黑色 红色 参数 · O × + * s d v 含义 点 圆 打叉 加号 星号 正方形 菱形 向下的三角形 参数 － ： －· -- ^ < > p 含义 实线 虚线 点划线 破折线 向上的三角形 向左的三角形 向右的三角形 五角星形 3．stem命令

绘制离散序列图，常用格式为stem(x)和stem(x，y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。

4．subplot命令

subplot（m，n，i）图形显示时分割窗口命令，把一个图形窗口分为m行，n列，m×n个小窗口，并指定第i个小窗口为当前窗口。 5．坐标轴标注

8

MATLAB中提供了许多关于坐标轴标注的函数，常用的函数有title、xlabel、ylabel等。其中，函数title是为图形添加标题，并将标题置于图形的顶部；而xlabel、ylabel是为x、y坐标轴添加标注，并分别将标注置于相应的坐标轴的边上。这三个函数的调用格式是大同小异的，这里仅以title为例介绍它们的调用格式。

title（‘标记’）

其中，标记可以文字，也可以是数学表达式。

6．文本标注

MATLAB语言对图形进行文本注释所提供的为text。它的调用格式为

text (x, y,‘文本’)

其中，(x, y)给定标注文本在图中添加的位置。

7．图例标注

在数值计算结果的绘图中，经常会出现在同一张图形中绘制多条曲线的情况，为了能更好地区分各条曲线，MATLAB语言提供了图例标注函数legend。它能为图形中所有的曲线进行自动标注，以其输入变量作为标注文本，具体的调用格式为

legend（‘标注1’，‘标注2’，„）

这里的标注1、标注2等分别对应绘图过程中按先后顺序所生成的曲线。 8．坐标轴的控制函数axis

函数axis用来控制坐标轴的刻度范围及显示形式。其最简单的调用格式为

axis（[xmin，xmax，ymin，ymax]）

即绘图时变量x∈[xmin，xmax]，变量y∈[ymin，ymax]。 例 绘制信号图形。

% fun0.m 定义文件名 x=0:0.1*pi:2*pi; %定义x向量

figure(1); %创建一个新的图形窗口,编号为1

subplot(2,2,1); %将窗口划分为2行,2列,在第一个窗口中作图 plot(x,sin(x)); %画图

title('正弦线'); %給图形加标题

subplot(2,2,2); %在第二个窗口中作图 plot(x,sin(x),'r') ; %画一正弦波红色 xlabel('X'); %给轴加说明 ylabel('SIN(X)'); %给轴加说明

subplot(2,2,3); %在第三个窗口中作图

plot(x,sin(x),'--',x,cos(x),'-.m*');%画一正弦波破折线

%画一余弦波紫色点划星线

9

legend('SIN(X)','COS(X)'); %给图形加图例标注 subplot(2,2,4); %在第四个窗口中作图 plot(x,sin(x),'r+--'); %画一正弦波红色破折线 text(4,0,'注记');

以上内容以文件名fun0.m存盘，在MATLAB命令窗口中执行

>> fun0↙

或将以上内容直接在MATLAB命令窗口中键入并↙，得到结果如下：

Fig.1

10

第四章 MATLAB的程序设计

MATLAB作为一种高级计算语言，它不仅可以如前面所介绍的那样，以一种人机交互式的命令行的方式工作，还可以像BASIC、FORTRAN、C等其它高级计算机语言一样进行控制流的程序设计，即编制一种以m为扩展名的文件，简称为M文件。而且， M文件的编写具有语法简单、可读性强、调试容易、调用方便等许多优点。 §4-1 M文件介绍

1． M文件的特点与形式

MATLAB实质上是一种解释性语言，就MATLAB本身来说，它并不能做任何事情，它就像DOS操作系统的一样，本身没有实现功能而只对用户发出的指令起解释执行的作用，即命令先送到MATLAB系统内解释，再运行得到结果。因此用户可以把所要实现的指令罗列编制成文件，再统一送入MATLAB系统中解释运行，这就是以m为扩展名的M文件。用户可以使用任何字处理软件对其进行编写或修改。正是M文件的这个特点造就了MATLAB强大的可开发性和可扩展性，Mathworks公司推出的一系列工具箱就是明证。

由于商用的MATLAB软件用C语言编写而成，因此，M文件的语法与C语言的十分相似。对广大的C语言爱好者来说，M文件的编写是相当容易的。

M文件有两种形式，即命令式和函数式。并且要注意 ·文件扩展名一定为m。

·标点符号的运用要恰到好处，建立好的书写风格，保持程序的可读性。 ·以符号%引导的行是注释行，不可执行，可供help命令查询； ·不需要用end语句作为M文件的结束标志；

·在运行此文件之前，需要把它所在目录加到MATLAB的搜索路径上去，或将文件所在目录设为当前目录。

2．命令式文件

命令式文件就是命令行的简单叠加，MATLAB会自动按顺序执行文件中的命令，其运行相当于在命令窗口中逐行输入并运行命令，因此，用户在编制此类文件时，只需把所要执行的命令按行编辑到指定的文件中，且变量不需预先定义，也不存在文件名对应问题，也可以访问存在于整个工作空间内的数据。但要注意

·命令式文件在运行中所产生的所有变量均为全局变量。也就是说，这些变量一旦生成，就一直保存在内存空间中，直到用户执行clear或quit时为止。 例 打开Medit窗口，编写如下程序 x=[1,3,4,7];

y=sum(x)/length(x);

以aver.m为文件名存盘。在命令窗口键入aver并回车即可运行该程序。

11

3．函数式文件

为了实现计算中的参数传递，需要用到函数式文件。函数式文件在MALTAB中应用十分广泛，MALTAB所提供的绝大多数功能函数都是由函数式文件实现的，这足以说明函数式文件的重要性。函数式文件的结构为：

function 输出参数=函数名（输入参数） 函数体 %注释

例 编制用于计算几个数的平均值的函数。

打开Medit窗口，编写如下程序

function y=lianxi(x)

n=length(x);%确定向量x的维数 if n==1 y=x;

end if n>=2

y=sum(x)/n;

end

编写完毕后，以lianxi.m存盘。

函数式文件的标志是第一行必为function语句。函数式文件可以有返回值如上例中的y，也可以只执行操作而无返回值，大多数函数式文件有返回值。函数体是函数的主体部分，它包括进行运算和赋值的所有MALTAB程序代码。函数体中可以包括流程控制、输入/输出、计算、赋值、注释以及函数调用和命令文件调用等。在函数体中完成对输出参数的计算。函数式文件执行后，只保留最后结果，不保留任何中间过程，所定义的变量也仅在函数内部起作用，并随调用的结果而被清除。函数调用的过程实际上就是参数的传递过程。

例 调用函数lianxi.m计算20，50，90，100，40的平均值。在命令窗口键入如下命令并回车即可。

x =[20,50,90,100,40]; y=lianxi(x);

该调用过程把变量x传递给了输入参数，然后把函数运算的返回值传给了输出参数y。

在编写函数式文件时，要特别注意函数中的

·只有文件名与函数名一一对应，才能保证调用成功。

·function后的语句定义函数名和输入、输出参数，在函数被调用过程中将按此输入、输出格式执行。

·要养成良好的注释习惯，以方便自己或其它用户的调用。 §4-2 程序控制语句

MALTAB的程序控制语句有循环语句、条件转移语句两种类型。

12

1．循环语句

MALTAB的循环语句包括for循环和while循环两种类型。 （1）for循环

语法格式：

for 循环变量=起始值：步长：终止值 循环体

end

起始值和终止值为一整型数，步长可以为整数或小数，省略步长时，默认步长为1。执行for循环时，判定循环变量的值是否大于（步长为负时则判定是否小于）终止值，不大于（步长为负时则不小于）则执行循环体，执行完毕后，再加上步长，这时循环变量的值若大于（步长为负时则小于）终止值后退出循环，否则继续。for循环最大的特点是，在一般情况下，此循环语句的循环次数是预先设定好的。

例 给矩阵a赋值。

function a=chuzhi k=5;

a=zeros(k,k); for m=1:k

for n=1:k

a(m,n)=1/(m+n-1); end end

以chuzhi.m存盘。然后在MALTAB命令窗口中执行以下命令

>> a=chuzhi 得到结果为：

a =

1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 说明：

·for语句一定要有end作为结束标志，否则下面的输入都被认为是for循环之内的内容。

·循环语句中的分号‘；’可防止中间结果的输出。 （2）while循环

语法格式：

while 表达式

13

循环体

end

其执行方式为：若表达式为真（运算值非0），则执行循环体；若表达式为假（运算结果为0），则退出循环体，执行end后的语句。while循环的特点是循环次数无法预先确定，因此在使用中要谨慎，以防止陷入死循环。 例 >> a=3;b=0;

while a

a=a-1; b=b+a; end >> b b = 3

>> a a = 0

同for循环比起来，while语句的判断控制可以是一个逻辑判断语句，因此它的适用范围会更广一些。 例 >> a=7;b=1; c=0

while (a>=1)&( b<=4) % a≥1与b≤4同时满足时执行循环体。

a=a-1; b=b+2;

c=3*a+2*b +c; end

>> [a,b,c] ans =

5 5 49

2．条件转移语句

MATLAB中的条件转移语句主要由if语句实现，其用法与其它高级语言相类似，基本语法格式如下：

if 表达式

执行语句1 else

执行语句2 end

在执行时，首先计算紧跟在关键字if后的表达式，若结果为真，则执行语句；若结果为假，则执行语句2。且else子句中可以嵌套if语句而形成elseif结

14

构，从而实现多路选择。

x0,x20例 欲实现一分段函数的计算。f(x)x 0x1 。

2x1x2 function f=pdbsline(x)

if x<0

f=0；

elseif x<1

f=x；

elseif x<2

f=2-x；

else

f=0；

end

以pdbsline.m存盘。然后在MALTAB命令窗口中执行以下命令

>> pdbsline (1.36) 得到结果为：

ans =

0. 00

15

第五章 常用数字信号处理函数

§5-1 波形发生器

1．基本信号序列的表示

离散时间信号是时间上不连续的“序列”，一般直接用x(n)表示。一个长度为N的序列x(n)在MATLAB中可以用一个N维行向量或列向量表示。常用的基本信号序列产生如下：

1n0

1） 单位抽样序列 (n)

0n0这一序列可用MATLAB中的zeros函数实现

x=[1，zeros(1,N-1)]；

1n0

2） 单位阶跃序列 u(n)

0n0

这一序列可用MATLAB中的ones函数实现

x=ones(1，N)； 3） 实指数序列 x(n)an n,aR

MATLAB实现：

n=0：N-1； x=a.^n；

4） 复指数序列 x(n)e(cj)n n

MATLAB实现：

n=0：N-1；

x=exp((c+j*w)*n)；

5） 随机序列

rand(1，N)产生在[0,1]上均匀分布的长度为N的随机序列。

randn(1，N)产生长度为N且均值为0方差为1的高斯随机序列，即白噪声序列。

2．基本波形的产生

1）三角波或锯齿波发生函数：sawtooth()

语法格式：sawtooth(t，width)。产生以2π为周期幅值范围在[-1，+1]之间的三角波或锯齿波。参数t为时间向量；width是[0，1]之间的数，它决定于函数在一个周期内上升部分和下降部分的比例。width=0.5产生三角波，width=1产生锯齿波，此时函数可简写为sawtooth(t)。 2）方波发生函数square()

语法格式：square(t)。产生以2π为周期幅值范围在[-1，+1]之间的方波，参数t为时间向量。

16

3）sinc发生函数：sinc()

语法格式：sinc(t)

sinc(t)sin(pi*t)/(pi*t);t0

t01例 信号产生举例

function fun1 clear all

t=0:0.0001:0.1; figure(1);

x1=sawtooth(2*pi*50*t); %在[0,0.1]时间内产生5个周期的锯齿波 subplot(2,2,1) plot(t,x1)

x2=sawtooth(2*pi*50*t,0.5);%在[0,0.1]时间内产生5个周期的三角波 subplot(2,2,2) plot(t,x2)

x3=square(2*pi*50*t); %在[0,0.1]时间内产生5个周期的方波 subplot(2,2,3) plot(t,x3)

axis([0,0.1,-1.2,1.2]) t=-4:0.1:4;

x4=sinc(t); %产生抽样函数 subplot(2,2,4) plot(t,x4)

运行结果如下：

§5-2 序列的操作

在数字信号处理中基本的序列运算及其MATLAB实现如下表5-1所示：

17

表5-1一些常用的序列运算及其matlab实现

信号加 信号乘 数乘 折叠 抽样和 抽样积 信号能量 序列运算 x(n)=x1(n)+x2(n) x(n)=x1(n)×x2(n) y(n)=k×x(n) y(n)=x(-n) n2其MATLAB实现 x=x1+x2 x=x1.*x2 y=k*x y=fliplr(x) y=sum(x(n1:n2)) y=prod(x(n1:n2)) 2ynn1n2x(n) 2yx(n) nn1Exnx(n)Ex=sum(abs(x).^2)) Px=sum(abs(x).^2))/n 信号功率 1N1Pxx(n) Nn0§5-3 常用数字信号处理函数

对于线性离散系统,其表达方式有多种：

b0b1z1bMzM传递函数法 H(z)

a0a1z1aNzN(zy1)(zy2)(zyM)零极点增益法 H(z)K

(zx1)(zx2)(zxN)带余式的部分分式展开法 H(z)rNr1r2k1k2z1kMN1z(MN) 1111p1z1p2z1pNz差分方程法 y(n)ay(nk)bx(nm)

kmk1m0NM数字信号处理即研究信号通过系统后的变化，其中经常用到的函数有： 1． x=roots(a)

欲求多项式y(z)=a0+a1z –1+a2z –2+„ 的根，可用命令x=roots(a)来实现，其中a=[a0，a1，a2，„]是它的系数向量，x便是y(z)的根向量。

例 求多项式y(z)=1-1.6z –1+0.65z –2-0.05z –3的根。编程运行如下： >> a=[1,-1.6,0.65,-0.05]; x=roots(a)

x =

1.0000 0.5000

18

0.1000

从而原多项式可表达为y(z)=(1-z –1)(1-0.5z –1)(1-0.1z –1)。

2． b=poly(x)

欲求出y(z)=(z-x1)(z-x2)„的多项式表达，可用命令b=poly(x)来实现，其中x=[x1,x2,„]是它的根向量，b是欲求多项式的系数向量。 例 某滤波器的零极点增益法表达式为H(z)(z0.1)(z0.5)(z1)，求其

(z0.3)(z2)(z5)传递函数表达式。MATLAB程序及运行结果如下：

function [a,b]=fun2(x,y,k) b=k*poly(y); a=poly(x);

>> y=[0.1,0.5,1]; x=[0.3,2,5]; k=1; [ b, a]=fun2(y,x,k) b =

1.0000 -1.6000 0.6500 -0.0500 a =

1.0000 -7.3000 12.1000 -3.0000

123从而滤波器的传递函数表达式为H(z)11.6z0.65z0.05z 。

17.3z112.1z23.0z33．[r,p,k]=residuez(b,a)

由传递函数表达式求出带余式的部分分式展开式时可用命令[r,p,k]=residuez(b,a) 来实现，其中b、a分别是原系统传递函数表达式中的分子、分母系数向量，而r、p、k分别是该系统的带余式的部分分式展开式中的系数向量。

例 对 H(z)48z1作部分分式展开。程序和显示输出如下： 16z18z2>> b=[-4,8]; a=[1,6,8]; [r,p,k]=residuez(b,a) r = -12 8 p = -4 -2 k =

[ ]

这就是说，系统的部分分式展开形式为：

H(z)12181。

14z12z4．y(n)=conv(x,h)

19

求两个序列的线性卷积y(n)x(n)*h(n)kx(k)h(nk)可直接采用

命令y(n)=conv(x,h)来实现。其中特别要注意的是，函数conv默认序列从n=0开始。若{x(n)：nx1≤n≤nx2}，{h(n)：nh1≤n≤nh2}，则卷积结果{y(n)：ny1≤n≤ny2}，其中ny1= nx1+ nh1，ny2= nx2+ ny2。 例 已知x(n)={3，4，5}，h(n)={2，6，7，8}，求它们的线性卷积。 ．．>> x=[3,4,5]; h=[2,6,7,8]; n1=-1:1; n2=-1:2; y=conv(x,h)

y =

6 26 55 82 67 40 >>nb3=n1(1)+n2(1);

>>ne3=n1(length(x))+n2(length(h)); >>n3=nb3:ne3 n3 =

-2 -1 0 1 2 3 即序列x(n)与h(n)的线性卷积为y(n)={6，26，55，82，67，40}。 ．5．X=fft(x,N) 与x=ifft(X,N)

采用FFT算法计算一维序列x(n)的N点离散Fourier变换(DFT)

nk可以直接用命令X=fft(x,N)来实现，其逆运算用命令x X(k)x(n)WNn0N1=ifft(X,N)来实现。计算中若序列x（或X）的长度大于N时，截断x（或X）；

否则补零。

例 模拟信号x(t)=2sin(4πt)+5cos(8πt)，求其点的DFT的幅值谱及相位谱。

MATLAB实现：

function y=fun3(N) n=0:N-1; t=0.01*n; q=n*2*pi/N; x=fun30(t); y=fft(x,N); subplot(3,1,1);

plot(t,x); title('原信号'); grid subplot(3,1,2);

plot(q,abs(y)); title('幅值'); grid subplot(3,1,3);

plot(q,angle(y)); title('相位'); grid function x=fun30(t)

20

x=2*sin(4*pi*t)+5*cos(8*pi*t); >>y=fun3();

6．[H,w]=freqz(b,a,n,Fs)

函数freqz基于FFT算法计算以传递函数形式表达的数字滤波器的频率响应。命令freqz(b,a,N,Fs)能够自动绘制以b、a为系数的数字滤波器的N点幅频和相频曲线(Fs是抽样频率，此项可缺省)，命令[H,w]=freqz(b,a,n)计算并返回数字滤波器的N点频率响应H（abs(H)是幅频响应，angle(H)是相频响应），相应的N个分布在(0~π)范围内的频率点记录在w中。 例 绘制系统H(z)1的幅频和相频特性曲线。 10.9z1>> b=1; a=[1,-0.9]; freqz(b,a,512)

7．y=filter(b,a,x)或[y,zf]=filter(b,a,x,xic)

系统以传递函数的形式表达，b、a为其分子、分母多项式的系数向量。信号x通过该系统后的零状态响应（即该系统对信号x的滤波），可用命令y=filter(b,a,x)或[y,zf]= filter(b,a,x)来完成，其中zf是系统的最后状态。

由于存储空间常常需要进行分布滤波，估算初始和最后条件对分布滤波来说是非常有用的。假设现有两部分数据，每部分有5000个点，如下：

x1=randn(5000,1);

21

x2=randn(5000,1);

可能第一个序列x1对应前10分钟采集的数据，第二个序列x2对应后10分钟采集的数据。那么整个序列为x=[x1;x2]。如果现在没有足够的空间存储连接起来的序列x。对x1、x2分别进行滤波，为了保证滤波后序列的连续性，采用滤波x1的最后条件作为滤波x2的初始条件。

[y1,zf]=filter(b,a,x1); y2=filter(b,a,x2,zf);

若计算系统在初始状态（Y,X）下的全响应，可用命令zf=filtic(b,a,Y,X) 将（Y,X）等效成输入xic，再用命令[y,zf]=filter(b,a,x,xic)来完成。

例 求解差分方程y(n)=[x(n)+x(n-1)+x(n-2)]/3+0.95y(n-1)-0.9025y(n-2)，n≥0。其中x(n)=cos(πn/3)，且y(-1)=-2，y(-2)=-3。 用MATLAN语句实现如下： function y=chafen

b=[1,1,1]/3;

a=[1,-0.95,0.9025]; Y=[-2,-3];X=[0,0];

xic=filtic(b,a,Y,X); %初始状态等效成输入xic n=0:50;

x=cos(pi*n/3);

y=filter(b,a,x,xic); %全响应

plot(n,y) %或stem(n,y)

或者：

function chafen %用递推法求解 n=0:50; x0=cos(pi*n/3); x=[0,0,x0];

y=[-3,-2,zeros(1,51)]; for i=3:53

y(i)=(x(i)+x(i-1)+x(i-2))/3+0.95*y(i-1)-0.9025*y(i-2); end

z(1:51)=y(3:53); plot(n,z)

22

第六章 MATLAB在数字信号处理中的应用

§6-1圆周卷积的计算

对于无限长序列不能用MATLAB直接计算线性卷积，在MATLAB内部只提供了一个conv函数计算两个有限长序列的线性卷积。对于圆周卷积MATLAB内部没有提供现成的函数,我们可以按照定义式直接编程计算。 例6-1 已知两序列： x(n)0

0.8n0n11 h(n)1 0n5

0其它其它求它们的线性卷积yl(n)=h(n)*x(n)和N点的圆周卷积y(n)=h(n)x(n)，并

研究两者之间的关系。 实现程序： （1）计算圆周卷积的函数

function yc=circonv(x1,x2,N) %实现两序列x1和x2的圆周卷积 if length(x1)>N

err0r('N must not be less than length of x1'); end

if length(x2)>N

err0r('N must not be less than length of x2'); end

x1=[x1,zeros(1,N-length(x1))]; %填充序列x1(n)使其长度为N x2=[x2,zeros(1,N-length(x2))]; %填充序列使x2(n)其长度为N m=[0:1:N-1];

x2=x2(mod(-m,N)+1); %生成x2的圆周反转序列 H=zeros(N,N);

for n=1:1:N %生成计算圆周卷积的矩阵

H(n,:)=cirshiftd(x2,n-1,N); % x2圆周右移n位 end

yc=x1*H'; %计算圆周卷积

function y=cirshiftd(x,m,N) %序列的圆周移位 if length(x)>N

error('The length of x must be less than N'); end

x=[x,zeros(1,N-length(x))]; ％补零，长度变为N n=[0:1:N-1];

y=x(mod(n-m,N)+1); %得到输出

23

（2）研究两者之间的关系 function fun5 clear all;

[N,xn,hn]=fun50;

yln=conv(xn,hn); %直接用函数conv计算线性卷积

ycn=circonv(xn,hn,N); %用函数circonv计算N点的圆周卷积 ny1=[0:1:length(yln)-1]; ny2=[0:1:length(ycn)-1]; n=0:length(xn)-1; m=0:length(hn)-1; subplot(2,2,1);

stem(n,xn); xlabel('xn') subplot(2,2,2);

stem(m,hn); xlabel('hn') axis([0,16,0,4]); subplot(2,2,3);

stem(ny1,yln); xlabel('线性卷积') subplot(2,2,4);

stem(ny2,ycn); xlabel('圆周卷积') function [N,xn,hn]=fun50 n=[0:1:11];

xn=0.8.^n; %生成x(n) hn=ones(1,6); %生成h(n) N=12;

运行结果如图所示：

24

§6-2 利用FFT实现线性卷积

若序列x1(n)、x2(n)是长度分别为N1、N2的有限长序列，N点的圆周卷积yc(n)=x1(n)x2(n)，长度为N1+N2-1的线性卷积yl(n)=x1(n)*x2(n)。由DFT的性质可知：当N≥N1+N2-1时有yc(n)=yl(n)=IDFT{DFT [x1(n)]·DFT[x2(n)]}。序列较长时DFT运算通常用快速算法FFT实现。

例6-2 用FFT实现上节例中的两序列的线性卷积。 实现程序如下：

%fun6.m

n=[0:1:11];m=[0:1:5];

N1=length(n); N2=length(m);

xn=0.8.^n; %生成x(n) hn=ones(1,N2); %生成h(n) N=N1+N2-1; xk=fft(xn,N); hk=fft(hn,N); yk=xk.*hk; yn=ifft(yk,N);

if all(imag(xn)==0)&(all(imag(hn)==0)) yn=real(yn);

end

x=0:N-1;

stem(x,yn,'.'); 运行结果为：

§6-3 系统的单位冲激响应

线性移不变系统的特性可用它在输入信号是单位冲激序列时的输出即系统的单位抽样响应h(n)来表征。有了它，就可得到此线性移不变系统对任意输入x(n)的输出y(n)。

25

例6-3 一个因果线性移不变系统：y(n)=0.81y(n-2)+x(n)-x(n-2)。 求：（1）单位冲激响应；（2）单位阶跃响应；（3）绘制系统的幅频和相频响应。

由差分方程直接得出系统函数H(z)=(1-z –2)/(1-0.81z –2)。从而可计算其它问题如下：

function [h,y]=chongji

b=[1,0,-1]; a=[1,0,-0.81];

xh=[1,zeros(1,50)]; % 产生单位冲激序列 xn=ones(1,51); % 产生单位阶跃序列 h=filter(b,a,xh); % 单位冲激响应 y=filter(b,a,xn); % 单位阶跃响应 figure(1); n=0:50;

subplot(2,1,1); stem(n,h,'.r')

legend('冲激响应'); subplot(2,1,2); stem(n,y,'k');

legend('阶跃响应'); figure(2);

freqz(b,a); % 绘制系统的幅频和相频响应

Fig.1 Fig.2

26

§6-4 IIR滤波器的设计与实现

无限长单位冲击响应（IIR）数字滤波器的设计方法有多种，对于一般条件下使用的数字滤波器，其常用的设计方法是基于模拟滤波器变换原理的经典设计：首先，将数字滤波器的技术指标转换成相应的模拟滤波器技术指标；设计模拟低通滤波器；然后，将设计好的模拟滤波器转换成满足给定技术指标的数字滤波器（常用算法有脉冲响应不变法和双线性变换）。在MATLAB的数字信号处理工具箱中，提供的有关设计函数如下：

1．模拟低通滤波器设计

模拟低通滤波器的逼近有巴特沃思型、切比雪夫型、和考尔型，分别用如下的函数实现。

[Z,P,K]=buttap(n)；返回一个n阶、巴特沃思型归一化的模拟低通滤波器的零极点增益模型

[Z,P,K]=cheb1ap(n,Rp)；n阶、通带内的最大衰减Rp、切比雪夫Ⅰ型 [Z,P,K]=cheb2ap(n,Rs)；n阶、阻带内的最小衰减Rs、切比雪夫Ⅱ型

[Z,P,K]=elliap(n,Rp,Rs)；n阶、通带内的最大衰减Rp、阻带内的最小衰减Rs、考尔型

2．模拟低通滤波器阶数n的选择函数

滤波器阶数的选择在整个滤波器的设计中占有十分重要的地位和作用。根据需要选择合适的滤波器阶数，MATLAB工具箱中提供了对应于各类模拟低通滤波器的阶数选择函数，如巴特沃思型的buttord、切比雪夫型的cheb1ord、cheb2ord和考尔型的ellipord。这些函数的调用格式大同小异，仅以buttord为例加以说明。

[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)

输入参数：Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减。

输出参数：n为符合要求的滤波器最小阶数，Wn为巴特沃思型模拟低通滤波器3dB截止频率。 ’s’：表示模拟域。

3．零极点增益模型到传递函数模型的转换 [num,den]=zp2tf(Z,P,K)

输入参数：Z，P，K分别表示系统的零极点增益模型的零点、极点和增益； 输出参数：num，den分别为同一系统传递函数模型的分子和分母多项式系数。 4．模拟域的频率变换

将归一化的模拟低通滤波器转换成所需要类型（低通、高通、带通和带阻）的模拟滤波器，可分别用如下命令实现：

[b,a]=lp2lp(Bap,Aap,Wn)；把传递函数形式的归一化模拟低通滤波器

27

原型转换成3dB截止频率为Wn的同型低通滤波器。

[b,a]=lphp(Bap,Aap,Wn)；转换成高通

[b,a]=lp2bp(Bap,Aap,W0,Bw)；转换成带通，W0：中心频率，Bw：带宽

[b,a]=lp2bs(Bap,Aap,W0,Bw)；转换成带阻 5．模拟滤波器数字化

[bz,az]=bilinear(b,a,Fs)：采用双线性变换法的映射关系。其中，Fs是采样频率。

[bz,az]=impinvar(b,a,Fs)：采用冲击响应不变法的映射关系。

例6-4 用双线性变换法设计一个Butterworth低通滤波器，要求其通带截止频率100Hz，阻带截止频率200Hz，通带衰减Rp小于2dB，阻带衰减大于15dB，采样频率Fs=500Hz。 MATLAB实现：

function [bz,az]=fun8 %求模拟滤波器参数

wp=100*2*pi; %通带截止频（率弧度/秒） ws=200*2*pi; %阻带截止频率（率弧度/秒） Rp=2; Rs=15;

Fs=500; %采样频率 Ts=1/Fs;

[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %选择低通滤波器的最小阶数N

[Z,P,K]=buttap(N); %创建零极点模型的巴特沃思型模拟低通滤波器 [Bap,Aap]=zp2tf(Z,P,K); %把滤波器零极点模型转化为传递函数模型 [b,a]=lp2lp(Bap,Aap,Wn); %把模拟滤波器原型转换成截止频率为的低通滤

波器

[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的

转换

[H,W]=freqz(bz,az);

plot(W*Fs/(2*pi),20*log10(abs(H))); %绘制频率响应曲线 grid

xlabel('频率/Hz')； ylabel('幅频响应(dB)') 运行结果为：

28

§6-5 FIR滤波器的设计与实现

相对于无限冲激响应数字滤波器，有限冲激响应数字滤波器可以实现具有精确的线性相位、总是稳定的、过渡过程有限、硬件易实现等优点，但所需滤波器的阶数较高，其延迟也较大。

有限长单位冲击响应（FIR）数字滤波器系统的传递函数模型为：

H(z)h(n)zn，其设计方法有窗函数法和频率采样法等，在MATLAB的

n0N1数字信号处理工具箱中提供了fir1等设计函数。fir1是采用经典窗函数法设计线性相位FIR数字滤波器的函数，且具有标准低通、带通、高通和带阻等类型。

语法格式：

B=fir1(n,Wn,’ftype’,window)

其中n为FIR滤波器的阶数即窗口长度为n+1，对于高通、带阻滤波器 n只取偶数。Wn为对Nyquist频率（即抽样频率的二分之一）归一化的滤波器截止频率，取值范围0~1，即当Wn=1时，滤波器实际截止频率等于Nyquist频率。对于带通、带阻滤波器，Wn=[W1,W2]，且W1例6-5 用窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率Wp=0.2π，阻带截止频率Ws=0.3π，阻带衰减不小于40dB，通带衰减不大于3dB。设抽样频率Fs=1 Hz.

由阻带衰减不小于40dB的情况下，选择汉宁窗将比较经济。用汉宁窗生成线性相位FIR低通滤波器的实现程序如下： function [b,k]=fun7

wp=0.2*pi;

29

ws=0.3*pi;

wdelta=ws-wp; %计算过渡带宽

N=ceil(8*pi/wdelta); %ceil:沿正无穷方向取整，估算窗口长度 Wn=(0.2+0.3)*pi/2; %估算窗口截止频率

Wn=Wn/pi; %对Nyquist频率(即πrad/s)归一化 b=fir1(N,Wn,hanning(N+1)); %用汉宁窗生成线性相位FIR低通滤波器 freqz(b,1,512); %绘出此系统的频响特性 n=0:511;

W1=exp(-j*wp*n); W2=wxp(-j*ws*n);

HW1=20*log10(abs(sum(b.*W1))) HW2=20*log10(abs(sum(b.*W2))) If (abs(HW1)<3)&(abs(HW2)>40)

k=1; %满足技术指标 else

k=-1 %不满足技术指标 end

运行结果为：

常用窗函数的MATLAN命令如下： 窗函数 命令 矩形窗 boxcar(n) 三角形窗 triang(n) 汉宁窗 海明窗 布拉克曼窗 凯泽窗 hanning(n) hamming(n) blackman(n) kaiser(n,β) 注：窗函数的使用方法相同。n均为窗的长度，β为凯泽窗的形状参数，它们的取值参考教材352~353页。

30

例6-6 设抽样频率为Fs=1000Hz，已知原信号为x=sin(2π×80t)+2sin(2π×140t)，由于某种原因，信号被白噪声污染，实际获得的信号为xn=x+rand(size(t))，要求设计一个FIR滤波器恢复出原始信号。 滤波器设计要求；

频带(Hz) [0,65] [75,85] [95,125] [135,145] [155,500] 幅度 0 1 0 1 0

function lvbo %生成实际信号xn

Fs=1000; t=0:1/Fs:2;

x=sin(2*pi*80*t)+sin(2*pi*140*t); xn=x+randn(size(t));

%生成一个阶数为100的多带滤波器

n=100;

w1=[75,145]/(Fs/2);w2=[85,135]/(Fs/2);

%生成以w1为通带的带通滤波器

b1=fir1(n,w1,hamming(n+1));

% 生成以w2为阻带的带阻滤波器

b2=fir1(n,w2,’stop’,hamming(n+1));

%计算两个系统的频率响应

[H1,W]=freqz(b1,1,512,2); [H2,W]=freqz(b2,1,512,2);

%绘系统的频响特性

figure(1)

plot(W,abs(H1),’r’, W,abs(H2)); grid

xlabel('归一化频率(Nyquist频率=1)') ylabel('幅度')

%对实际信号xn进行滤波

y1=filter(b1,1,xn)； y=1.5*filter(b2,1,y1)；

%检查滤波效果

figure(2)； nn=500:750； tt=nn/Fs；

subplot(2,1,1)；

31

plot(tt,x(nn),tt,y(nn),'r')；

legend('原始信号','滤波信号')； grid

subplot(2,1,2)； plot(tt,xn(nn))；

legend('污染信号')；grid; xlabel('时间(秒)')；

Fig.1 FIR滤波器的幅频特性

Fig.2 滤波效果对比

32

实验一 常见离散信号的MATLAB产生和图形显示

一、实验目的

加深对常用离散信号的理解； 二、实验内容

编制程序产生以下几种信号（长度可输入确定），并绘出其图形。 三、实验原理 1． 单位抽样序列

1 n0 (n) 0 n0在MATLAB中可以利用zeros()函数实现。

x[1,zeros(1,N1)]; 或

xzeros(1,N);x(1)1;如果(n)在时间轴上延迟了k个单位，得到(nk)即： (nk)0 nk

2． 单位阶跃序列

1 n0 u(n)0 n0

在MATLAB中可以利用ones()函数实现。

xones(1,N);

3． 正弦序列

(n) x Asin(2fn/fs)在MATLAB中

n0:N1

xA*sin(2*pi*f*n/fsfai)

4． 复正弦序列

jn(n) x Ae在MATLAB中

n0:N1

xA*exp(j*w*n) 5． 实指数序列

x(n)Aan

在MATLAB中 n0:N1 xA*a.^n四、实验步骤

1． 产生点的单位抽样序列δ(n)、点并移位20位的单位抽样序列δ(n-20)

33

1 nk

及任意序列f(n)=8.0δ(n)+3.4δ(n-1)+1.8δ(n-2) +5.6δ(n-3)+ 2.9δ(n-4) +0.7δ(n-5)。

2． 产生N＝32点的单位阶跃序列及斜率为3，n0=4,点数为32点的斜坡序列g(n)=B(n-n0)

3． 产生幅度A=3，频率f＝100，初始相位＝1.2，点数为32点的正弦序列。 4． 产生幅度A=3，角频率ω＝314，点数为32点的复正弦序列。 5． 产生幅度A=3，a＝0.7，点数为32点的实指数序列。 五、实验要求

（1）熟悉MATLAB软件，预先阅读前面章节中有关预备知识。

（2）根据实验内容，编写实验程序，上机运行；并进行结果验证，如有不同，查找原因，解决问题。

（3）讨论复指数序列的性质。 六、程序示例 function sy1 clear all; N=;

unit_pulse(N) k=20;

shift_unit_pulse(N,k) N=8;

arbitrary_pulse(N) N=32;k=4;B=3; unit_step(N) slope(N,k,B)

A=3;f=100;fai=1.2; sine(N,A,f,fai) w=314;

complex_sine(N,A,w) a=0.7;

real_exponent(N,A,a)

function unit_pulse(N) % unit_pulse.m x=zeros(1,N); x(1)=1; n=0:N-1; figure(11);

stem(n,x);xlabel('单位抽样序列')

34

axis([-1 65 0 1.1])

function shift_unit_pulse(N,k) % shift_unit_pulse.m x=zeros(1,N); x(k+1)=1; n=0:N-1; figure(12);

stem(n,x);xlabel('移位20位的单位抽样序列') axis([-1 65 0 1.1])

function arbitrary_pulse(N) % arbitrary_pulse.m x=zeros(1,N);

x(1)=8.0;x(2)=3.4;x(3)=1.8; x(4)=5.6;x(5)=2.9;x(6)=0.7; n=0:N-1; figure(13);

stem(n,x);xlabel('任意序列f(n)') axis([-1 8 0 8.2])

function unit_step(N) % unit_step.m x=ones(1,N); n=0:N-1; figure(21);

stem(n,x);xlabel('单位阶跃序列') axis([-1 32 0 1.1])

function slope(N,k,B) % slope.m

x=[zeros(1,k) ones(1,N-k)]; for i=1:N

x(i)=B*x(i)*(i-k); end

n=0:N-1; figure(22);

stem(n,x);xlabel('斜坡序列')

35

axis([-1 32 0 90])

function sine(N,A,f,fai) % sine.m n=0:N-1;

x=A*sin(2*pi*f*(n/N)+fai); figure(3);

stem(n,x);xlabel('正弦序列') axis([-1 32 -3.2 3.2])

function complex_sine(N,A,w) % complex_sine.m n=0:N-1;

x=A*exp(j*w*n); figure(4);

stem(n,x);xlabel('复正弦序列') axis([-1 32 -3.2 3.2])

function real_exponent(N,A,a) % real_exponent.m n=0:N-1; x=A*a.^n; figure(5);

stem(n,x);xlabel('实指数序列') axis([-1 32 0 3.2])

36

实验二 离散系统的频率响应分析和零、极点分布

一、实验目的

加深对离散系统的频率响应分析和零、极点分布的概念理解。 二、实验内容

求系统

0.05280.797z10.1295z20.1295z30.797z40.0528z5

H(z)11.8107z12.4947z21.8801z30.9537z40.2336z5

（离散系统传递函数）

的零、极点和幅度频率响应。 三、实验原理

离散系统的时域方程为

NM

aky(nk)bmx(nm) k0m0其变换域分析方法如下：

频域

y[n]x[n]h[n]x[m]h[nm]Y(ej)X(ej)H(ej) m系统的频率响应为 Y(ej)b0b1ej...bMejMjH(e)X(ej)a0a1ej...aNejN

Z域

y[n]x[n]h[n]Y(z)b0b1z1...bMzM系统的传递函数为 H(z)X(z)a0a1z1...aNzN MM01H(z)mKmN分解因式 ，其中ym和Ni1mx[m]h[nm]Y(z)X(z)H(z)bk0kziz(1ymz1))k1

称为零、极点。

在MATLAB中，可以用函数[z,p,K]=tf2zp(b,a)求得有理分式形式的系统传递函数的零、极点，用函数zplane(z,p)绘出零、极点分布图；也可以用函数zplane(b,a)直接绘出有理分式形式的系统传递函数的零、极点分布图。

另外，在MATLAB中，可以用函数 [r,p,k]=residuez(b,a）完成部分分式展开计算；可以用函数sos=zp2sos(z,p,K)完成将高阶系统分解为2阶系统的串联。 四、实验步骤

1. 由系统传递函数求出系统的零极点，并进行部分分式展开； 2. 求出系统频率响应，画出幅频及相频特性曲线。 3. 应用不同方法实现上述1，2的要求。 五、实验要求

编程实现系统参数输入，绘出幅度频率响应曲线和零、极点分布图。

ak(1xkzxk

37

六、程序示例

%sy2-1_zero_pole_plot.m clear all;

b=[0.0528 0.797 0.1295 0.1295 0.797 0.0528]; a=[1 -1.8107 2.4947 -1.8801 0.9537 -0.2336]; g=b(1)/a(1);

%求得有理分式形式的系统传递函数的零、极点 [z,p,K]=tf2zp(b,a) %绘出零、极点分布图 zplane(z,p); %部分分式展开

[r,p,k]=residuez(b,a)

%sy2-2_transfer_function_response.m clear all;

b=[0.0528 0.797 0.1295 0.1295 0.797 0.0528]; a=[1 -1.8107 2.4947 -1.8801 0.9537 -0.2336]; %求系统频率响应 fs=1000;

[h,f]=freqz(b,a,256,fs); mag=abs(h); ph=angle(h); ph=ph*180/pi; subplot(211) plot(f,mag) grid

xlabel('frequency(Hz)'); ylabel('magnitude'); subplot(212) plot(f,ph) grid

xlabel('frequency(Hz)'); ylabel('phase');

38

实验三 序列线性卷积、圆周卷积的计算及其关系的研究

一、实验目的

1．实现序列线性卷积、圆周卷积的计算； 2．研究序列线性卷积、圆周卷积的关系。 二、实验内容

设有一离散因果系统，其输入输出关系由以下差分方程确定

y(n)-0.5y(n-1)=x(n)+0.5x(n-1)

1．求该系统的单位抽样响应；

2．利用卷积和求系统在x(n)=e j3n激励下的零状态响应y(n)。 三、实验原理

1．当一个离散时间序列信号x(n)通过一个线性离散时间系统（其单位冲击响应为h(n)）后的输出y(n)，用数学语言表达为x(n)与h(n)的线性卷积，即 y(n)= x(n)* h(n)=

h(k)x(nk)

k02．两个有限长的序列x1（n）、N1和x2（n）、N2的N点圆周卷积z（n）

Nx2（n）=z（n）= x1（n）○

x(k)x(nk)R12k0N1N(n)

3． 若x(n)的长度为N1，h(n)的长度为N2，在N＞N1+ N2－1的条件下，令

N1nN1N2nN100，x2(n) x1(n)x(n)0nN11h(n)0nN21Nx2（n） 则 x1（n）*x2（n）=x1（n）○

即可用序列的圆周卷积计算序列的线性卷积。 四、实验步骤

编写MATLAB应用程序，以实现如下功能 1．描述系统；

2．计算系统的单位抽样响应，并绘出频率特性；

3．直接计算x(n)与h(n)的线性卷积，并绘出频率特性；

3．在N＞N1+ N2－1、N=N1+ N2－1、N＜N1+ N2－1三种情况下，计算x(n)与h(n)的圆周卷积，并绘出频率特性。 五、实验要求

1．了解MATLAB软件，认真学习并练习每一个例题。

2．针对实验内容，首先进行理论分析，给出预测结果；再编写程序，上机运行，给出仿真结果；进行结果验证，如有不同，查找原因，解决问题。

39

六、程序示例

function [hn,yn]=shiyan1(M,N)

b=[1,0.5];a=[1,-0.5]; %线性离散系统的传递函数法的描述 n=0:M-1; %令M= N2（=N1） xh=[1,zeros(1,M-1)]; %创建单位抽样激励

hn=filter(b,a,xh); %求该系统的单位抽样响应hn figure(1);

stem(n,hn); %绘制系统的单位抽样响应hn的离散序列图 axis([0,M,0,1.1]);

xn=exp(j*3*n); %创建激励x(n)=e j3n yn=conv(xn,hn); %直接计算线性卷积y(n)

ynn=circonv(xn,hn,N); %调用函数circonv（23页例6-1）计算xn与hn %N点的圆周卷积 nh=0:1:length(yn)-1; ny=0:1:length(ynn)-1; figure(2)

subplot(3,1,1);

stem(n,abs(xn)) %绘制激励x(n)=e j3n幅度的离散序列图 axis([0,M,0,1.1]); subplot(3,1,2);

stem(nh,abs(yn)); %绘制零状态响应yn幅度的离散序列图 axis([0,N,0,1.1]); subplot(3,1,3);

stem(ny,abs(ynn),'r'); %绘制零状态响应ynn幅度的离散序列图 axis([0,N,0,1.1]); grid

figure(3) subplot(3,1,1)

stem(n,angle(xn)*180/pi); %绘制激励x(n)=e j3n幅角的离散序列图 axis([0,M,-180,180]); subplot(3,1,2);

stem(nh,angle(yn)*180/pi); %绘制零状态响应yn幅角的离散序列图 axis([0,N,-180,180]); subplot(3,1,3);

stem(ny,angle(ynn)*180/pi,'r'); %绘制零状态响应ynn幅角的离散序列图 axis([0,N,-180,180]); grid

40

实验四 利用DFT分析信号的频谱

一、实验目的

1．实现有限长序列的点的计算；

2．研究记录时间、采样频率、频域采样点数对频谱的影响； 二、实验内容

有一调幅信号

xa(t)=[1+cos(2π×50t)] cos(2π×500t)

确定合理的采样频率fS及记录长度T0，用DFT做频谱分析。要求能分辨xa(t)的所有频率分量。 三、实验原理

欲用DFT分析连续时间信号的频谱，要做到 1．将连续时间信号xa(t)离散为序列xa(nT)

时域采样定理：要想采样后能够不失真地还原出原信号，则采样频率fS必须大于两倍的信号谱最高频率fh。即fS>2fh。

时域采样间隔T=1/fS

xa(nT)= xa(t)|t=nT n为整数

2．设记录T0，则意味着将xa(nT)截短为一有限长的序列x(n)，其长度N1= T0/T，即在时域内乘一个矩形窗函数，用数学语言表达为 x(n)= xa(nT)RN1(n)

对周期序列而言，T0应为序列周期的整数倍，否则由于有限长的矩形窗函数影响，造成原信号频谱的“扩散”—频谱泄露，使最高频率有可能超过折叠频率（fS/2），产生频率混叠现象。

3．对有限长序列的频谱X(ejω)进行采样，在一个周期内采样点数为N2。 频域采样定理：对于N1点的有限长序列，频域采样不失真的条件是 N2≥N1

在N2＞N1时，对补零到长度为N2。（通常取N2=N1）

N214．X(k)=DFT[x(n)]=

x(n)en0j2nkN2=X(ej)2N2

四、实验步骤

1．编写一计算有限长序列的DFT的程序； 2．确定信号x(t)的周期Ti及最高频率fh； 3．设定fS=2kHz →T=0.5mS，将xa(t)离散化；

4．令记录长度T0=0.02S →N1=40，编程计算N2=40的X(ejω)、X(k)并绘制幅频图；

5．令记录长度T0=0.016S →N1=32，编程计算N2=32或的X(ejω)、X(k)并绘制幅频图；

41

6．令记录长度T0=0.02S →N1=80，编程计算N2=80或128的X(ejω)、X(k)并绘制幅频图；

7．令记录长度T0=0.05S →N1=100，编程计算N2=100或128的X(ejω)、X(k)并绘制幅频图 五、实验要求

仔细观查并解释步骤4~7中的各种现象。 六、程序示例

function shiyan2(N1,N2,M)

n0=[0:1:N1-1]; %截取原始信号N1点 n=[0:1:N2-1]; %补零N2-N1个 y=(1+cos(0.05*pi*n0)).*cos(0.5*pi*n0);

%原始信号，其周期N=40，频率分量为0.45*pi，0.5*pi，0.55*pi x=[y(1:1:N1),zeros(1,N2-N1)]; %信号序列 w=[0:1:2000]*2*pi/2000;

x0=x*exp(-j*n'*w); %信号序列的傅里叶变换

magx0=abs(x0); %信号序列的傅里叶变换的幅度 a=max(magx0);

magx0=magx0/a; %幅度归一化

x1=dft(x,N2); %信号序列的N2点离散傅里叶变换 magx1=abs(x1(1:1:N2/2+1))/a; k1=0:1:N2/2; w1=2*pi/N2*k1; figure(M)

subplot(2,1,1)

stem(n,x); %绘制信号序列的离散序列图 axis([0,N2,-2.5,2.5]) line([0,N2],[0,0]) subplot(2,1,2)

plot(w/pi,magx0) %绘制信号序列的傅里叶变换的幅频响应图 hold on

stem(w1/pi,magx1,'r') %绘制信号序列傅里叶变换的幅频响应离散序列图 axis([0,1,0,1.1])

function [Xk]=dft(xn,N)

%计算序列xn的N点离散傅里叶变换 n=[0:1:N-1];k=[0:1:N-1]; WN=exp(-j*2*pi/N); nk=n'*k;

WNnk=WN.^nk;

42

Xk=xn*WNnk;

实验五 信号时间尺度变换的研究

一、实验目的

1．研究抽取因子为整数时，抽取过程对频谱产生的影响； 2．研究插值因子为整数时，插值过程对频谱产生的影响；

3．在不产生频率响应的混叠失真的条件下，实现比值为有理数的抽样率转换。 二、实验内容 已知信号

33

x(t)=0.5+0.93*cos(0.025×2πt×10)+0.86*cos(0.050×2πt×10)+

33

0.79*cos(0.075×2πt×10)+0.72*cos(0.100×2πt×10)+0.65*cos(0.125×2πt×333

10)+0.58*cos(0.150×2πt×10)+0.51*cos(0.175×2πt×10)+

33

0.44*cos(0.200×2πt×10)+0.37*cos(0.225×2πt×10)+0.30*cos(0.250×2πt×333

10)+0.23*cos(0.275×2πt×10)+0.16*cos(0.300×2πt×10)+

33

0.09*cos(0.325×2πt×10)+0.02*cos(0.350×2πt×10)

其最高频率为fh=0.35kHz，抽样频率为fs=1kHz。现欲将此信号的抽样率转换为fs’=1.4kHz，设计该转换过程，并画出频谱的变化。 三、实验原理

实现抽样率的转换，是数字信号处理的一个重要内容。比较简单的办法是直接在数字域内对抽样信号作抽样率的变换，从而得到新的抽样信号。 1．序列的插值（插值因子I为整数）

插值就是在信号序列x(n)中，每相邻两个数据之间插入I-1个零，从而得到新的抽样信号z(n) 。用数学形式表达为

x(n/I)niI,i0,,N1 z(n)其它0插值意味着抽样率由fs转换为Ifs，这时

IN1Z(ej)z(n)en0jn z(iI)ei0N1jiI 

x(i)ei0N1jiIX(ejI)N1i0IikNIikx(i)WNX((k))N i0N1IN1或：Z(k)z(n)Wn0nkNIz(iI)W显然，经过插值，频谱（在数字域内）被压缩I倍，即在2π范围内将有I个完整周期。因而再加接一个数字低通滤波截取一个周期即可。 2．序列的抽取（抽取因子D为整数）

43

抽取就是在信号序列x(n)中，每D个数据中取出一个，从而得到新的抽样信号y(n) 。用数学形式表达为y(n)= x(Dn)

抽取意味着抽样率由fs压缩为fs/D，这时频谱（数字域）被展宽D倍，即

Y(ej)N11X(ej/D) DnkN Di0n0n0其中k/D为整数。即经过抽取，频谱被扩展D倍。因此，若原信号频谱的最高频率为ωh，抽取后最高频率为Dωh，当Dωh＞π时将发生频率响应的混叠失真。

3．比值为有理数的抽样率转换

若希望将信号x(n)的抽样率转换为I/D倍，比较合理且有效的办法是先对信号做I倍的插值，经过数字低通滤波后，再做D倍的抽取，抽样率由fs转换为(I/D)fs。 四、实验步骤

1．将连续信号x(t)离散为序列x(n)。

Y(k)Y(n)Wx(nD)WN1nkN 1ND1ik/Dx(i)WND1X(k/D)Dx(n)x(t)|tnfS

2．观察在对序列进行插值、抽取后频谱的变化。 3．实现对信号x(t)（即序列x(n)）的抽样率的转换。 五、实验要求

仔细观察I、D分别为各种整数时的频谱现象，并解释之。 六、程序示例

function shiyan31(I,D)

%观察对序列插值、抽取后频谱的变化 n0=[0:1:39]; N=40; n1=2/N*n0;

xn=0.5+0.7*cos(0.05*pi*n0)+0.4*cos(0.1*pi*n0)+0.1*cos(0.15*pi*n0); xk=dft(xn,N); magxk=abs(xk); subplot(3,1,1); stem(n1,magxk); %信号序列的插值 yn=chazhi(xn,I); M=length(yn); yk=dft(yn,M); magyk=abs(yk);

m=[0:1:M-1]; m1=2/M*m; subplot(3,1,2)

44

stem(m1,magyk); %信号序列的抽取 zn=chouqu(xn,D);

M=length(zn); m=[0:1:M-1]; m1=2/M*m; zk=dft(zn,M); magzk=abs(zk); subplot(3,1,3) stem(m1,magzk);

function shiyan32(I,D,wp)

%抽样转换率为I/D，先插值(I)，经低通滤波（截止频率wp*pi），再抽取(D)。 %有限频宽序列xn n0=[0:1:39]; N=40;

x1=0.5+.93*cos(0.05*pi*n0)+.86*cos(0.1*pi*n0); x2=.79*cos(0.15*pi*n0)+.72*cos(0.20*pi*n0); x3=.65*cos(0.25*pi*n0)+.58*cos(0.30*pi*n0); x4=.51*cos(0.35*pi*n0)+.44*cos(0.40*pi*n0); x5=.37*cos(0.45*pi*n0)+.30*cos(0.50*pi*n0); x6=.23*cos(0.55*pi*n0)+.16*cos(0.60*pi*n0); x7=.09*cos(0.65*pi*n0)+.02*cos(0.70*pi*n0); xn=x1+x2+x3+x4+x5+x6+x7; %序列的频谱xk xk=dft(xn,N); a=max(abs(xk));

magxk=abs(xk)/a; %归一化处理 n1=2/N*n0; subplot(4,1,1); stem(n1,magxk);

%信号序列的插值序列yn, 插值因子I, yn=chazhi(xn,I); M=length(yn);

yk=dft(yn,M); %序列插值后的频谱yk b=max(abs(yk));

magyk=abs(yk)/b; %归一化处理 m=[0:1:M-1]; m1=2/M*m; subplot(4,1,2) stem(m1,magyk);

%生成滤除高频分量的滤波器hk

45

np=wp*M/2;

hk=[1,zeros(1,M-1)]; for i=2:1:np hk(i)=1;

hk(M-i+2)=1; end

ykk=yk.*hk; %经过滤波后的频谱ykk magykk=abs(ykk)/b; subplot(4,1,3)

stem(m1,magykk);

%计算经过滤波后的序列ynn

ykkk=conj(ykk); %取复共轭函数:conj ynn=conj(dft(ykkk,M))/M;

%ynn序列的抽取序列zn，抽取因子D zn=chouqu(ynn,D); M=length(zn);

zk=dft(zn,M); %序列抽取后的频谱zk c=max(abs(zk));

magzk=abs(zk)/c; %归一化处理 m=[0:1:M-1]; m1=2/M*m; subplot(4,1,4) stem(m1,magzk);

function yn=chazhi(xn,I) N=length(xn); M=I*N;

yn=zeros(1,M); for n=1:1:N

m=(n-1)*I+1; yn(m)=xn(n); end

function yn=chouqu(xn,D) N=length(xn); m=0;

for n=1:D:N m=m+1;

yn(m)=xn(n);

46

end

实验六 快速傅里叶变换及其应用

一、实验目的

1在理论学习的基础上，通过本实验加深对DFT算法原理和基本性质的理解；熟悉FFT算法原理和FFT子程序的应用。

2学习应用FFT对连续信号和典型时域离散信号进行频谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。 3熟悉应用FFT实现两个序列的线性卷积的方法，理解掌握循环卷积与线性卷积之间的关系。

4熟悉应用单次N点DFT计算一个2N点实序列的DFT方法。 二、实验内容

(1) 观察衰减正弦序列x1(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现的位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现的位置，有无混叠和泄漏现象？说明产生现象的原因。

(2) 观察三角波和反三角波的时域和幅频特性，用N=8点FFT分析信号序列x2(n)和x3(n)的幅频特性，观察两者的序列形状和幅频曲线有什么异同？绘出两序列及其幅频特性曲线。

在x2(n)和x3(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？

(3) 一个连续信号含有3种频率成分，经采样得：

x(n)sin(2nf1/fs)sin(2nf2/fs)sin(2nf3/fs) n=0,1,...,N-1已知f120Hz, f220.5Hz, f340Hz, 采样频率fs100Hz, ① 求x(n)在0～128之间的DFT的X(k)；

② 求把①中的x(n)以补零方式使其加长到0～512之间后的DFT的X(k)； ③ 求x(n)在0～512之间的DFT的X(k)； 观察以上①②③问题结果有何不同？为什么？

(4) 利用FFT计算循环卷积。设：xa(n)1,2,3,4,5, xb(n)6,7,8,9，计算其5，6，7，8，9点循环卷积，并分析线性卷积与循环卷积的关系。 (5) 2N点实数序列

212 7n)cos(19n), n0,1,2,...,2N1cos( )a(nN2N 0, 其它nN=。用一个点的复数FFT程序，一次算出 Xa(m)DFT[a(n)]2N，并绘出 Xa(m)。 三、实验原理

N点序列的DFT和IDFT变换定义式如下： N1N11knkn X[k]x[n]WNx[n]X[k]WN， n0Nk0 47

2 jnkknNWNe 利用旋转因子 具有周期性，可以得到快速算法（FFT）。

在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。

关于频谱分析及圆周卷积和线性卷积计算原理可参见前面实验章节中内容。 利用一个N点复数FFT程序计算一个2N点实序列原理如下：

1． 首先把一个2N点的实序列a(n)按奇偶分成两部分a1(n)、a2(n)； 2． 再组合成一个新的N点的复数序列x(n)，N点奇序列作为虚部，N点偶

序列作为实部，即x(n)＝a1(n)＋i*a2(n)；

3． 对于N点的复数序列x(n)求其FFT结果X(m)＝Xr(m)+Xi(m),

m=0,1,„,N-1；

4． 为了从X(m)提取出想要的2N点的实序列a(n)的FFT结果

Xa(m)=Xar(m)+i*Xai(m),定义如下关系式：

X(m)Xr(Nm)X(m)Xr(Nm) Xr(m)rXr(m)r 22 X(m)Xi(Nm)Xi(m)Xi(Nm)Xi(m)iX(m)i 225． 将以上关系式的值代入下面的式子就可以求出最后的结果Xa(m)的实部

和虚部前N个点： mmX(m)X(m)cosX(m)sinriXr(m) arNNm0,1,,N1

mm Xai(m)Xr(m)sinXi(m)cosXr(m)NN

6．对于2N点的Xa(m)的后N个点Xa(N)～Xa(2N-1)应是Xa(0)～Xa(N-1)的复共轭，不需重新计算。 四、实验步骤

1． 复习DFT的定义、性质和用DFT进行频谱分析的有关内容； 2． 复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。 3． 编写主程序；

4． 编写信号产生子程序，产生以下典型信号供频谱分析用： (1) 衰减正弦序列 -anesin(2fn) 0n15 x1(n) 0 其它n

(2) 三角波序列 n 0n3 x2(n)8n 4n7 0 其它n

(3) 反三角波序列 4n 0n3 x(n)n4 4n73

0 其它n48

(4) 编写FFT频谱分析子程序，计算信号的DFT，绘制信号及幅频特性曲线。 (5) 编写FFT计算循环卷积子程序，计算两个序列线性卷积及5，6，7，8，9点的循环卷积，并绘制序列及卷积结果曲线。

(6) 编写一个N点复数FFT程序计算一个2N点实序列的DFT。 (7) 按实验内容要求，上机实验，并写出实验报告。 五、实验要求

1． 熟悉MATLAB软件及其应用。

2． 针对实验内容，首先进行理论分析，进行结果预测；再利用MATLAB编程完

成计算，绘出相应图形。并与理论计算相比较，说明实验结果的原因。 3． 总结实验所得主要结论，简要回答思考问题。 六、程序示例

function [y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11]=sy6 clear all; a=0.1;

[y1,y2,y3]=fun2_1(a); [y4,y5,y6,y7]=fun2_2(a); [y8,y9,y10,y11]=fun2_3; xa=[1 2 3 4 5]; xb=[6 7 8 9];

[ykl,yk5,yk6,yk7,yk8,yk9]=fun3(xa,xb); [an,Xa]=fun4;

function [y1,y2,y3]=fun2_1(a) %频谱分析子程序_1

%a=0.1,f=0.0625衰减正弦序列 N=16; f=0.0625; n=0:N-1;

[x1,x2,x3]=fun1(a,f); %产生衰减正弦序列x1(n) figure(1); subplot(2,3,1) stem(n,x1)

title('x1(n),a=0.1,f=0.0625'); subplot(2,3,4) y1=myditfft(x1); q=n*2*pi/N; plot(q,abs(y1))

title('FFT频谱N=16,f=0.0625'); %a=0.1,f=0.4375衰减正弦序列

49

f=0.4375;

[x1,x2,x3]=fun1(a,f); %产生衰减正弦序列x1(n) subplot(2,3,2) stem(n,x1)

title('x1(n),a=0.1,f=0.4375'); subplot(2,3,5) y2=myditfft(x1); plot(q,abs(y2))

title('FFT频谱N=16,f=0.4375'); %a=0.1,f=0.5625衰减正弦序列 f=0.5625;

[x1,x2,x3]=fun1(a,f); %产生衰减正弦序列x1(n) subplot(2,3,3) stem(n,x1)

title('x1(n),a=0.1,f=0.5625'); subplot(2,3,6) y3=myditfft(x1); plot(q,abs(y3))

title('FFT频谱N=16,f=0.5625');

function [y4,y5,y6,y7]=fun2_2(a) %频谱分析子程序_2 f=0.0625;

[x1,x2,x3]=fun1(a,f); %产生三角波序列x2(n)和反三角波序列x3(n) %N=8点FFT程序 N=8;n=0:N-1; figure(2); subplot(2,2,1) stem(n,x2)

title('三角波序列x2(n)'); subplot(2,2,3) y4=myditfft(x2); q=n*2*pi/N; plot(q,abs(y4))

title('FFT频谱N=8'); subplot(2,2,2) stem(n,x3)

title('反三角波序列x3(n)');

50

subplot(2,2,4) y5=myditfft(x3); plot(q,abs(y5))

title('FFT频谱N=8'); %末尾补零的16点FFT程序 M=16;m=0:M-1;

x2=[x2 zeros(1,M-N)]; %对三角波序列x2(n)末尾补零 x3=[x3 zeros(1,M-N)]; %对反三角波序列x3(n)末尾补零 figure(3); subplot(2,2,1) stem(m,x2)

title('末尾补零的三角波序列x2(n)'); subplot(2,2,3) y6=myditfft(x2); q=m*2*pi/M; plot(q,abs(y6))

title('FFT频谱N=16'); subplot(2,2,2) stem(m,x3)

title('末尾补零的反三角波序列x3(n)'); subplot(2,2,4) y7=myditfft(x3); plot(q,abs(y7))

title('FFT频谱N=16');

function [y8,y9,y10,y11]=fun2_3 %频谱分析子程序_3 %问题①的程序 N=128;n=0:N-1;

fs=100;f1=20;f2=20.5;f3=40;

xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs); y8=myditfft(xn); ay8=abs(y8(1:N/2)); figure(4); subplot(2,1,1) plot(n,xn)

ylabel('原始信号x(n)'); q=(0:N/2-1)*fs/N;

51

subplot(2,1,2) plot(q,ay8)

ylabel('FFT N=128'); %问题②的程序 M=512;m=0:M-1;

xn=[xn zeros(1,M-N)]; y9=myditfft(xn); ay9=abs(y9(1:M/2)); figure(5); subplot(2,1,1) plot(m,xn)

ylabel('补零原始信号x(n)'); q=(0:M/2-1)*fs/M; subplot(2,1,2) plot(q,ay9)

ylabel('FFT N=512'); %问题③的程序 n=0:M-1;

xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs); y10=myditfft(xn); ay10=abs(y10(1:M/2)); figure(6); subplot(2,1,1) plot(n,xn)

ylabel('原始信号x(n)'); q=(0:M/2-1)*fs/M; subplot(2,1,2) plot(q,ay10)

ylabel('FFT N=512');

function [x1,x2,x3]=fun1(a,f) %信号产生子程序 N1=16; n1=0:N1-1;

x1=exp(-a.*n1).*sin(2*pi*f*n1); %产生衰减正弦序列 x2=[[0:3] 8-[4:7]]; %产生三角波序列 x3=[4-[0:3] [4:7]-4]; %产生反三角波序列

52

function [ykl,yk5,yk6,yk7,yk8,yk9]=fun3(xa,xb) %计算线性卷积子程序

N1=length(xa); %求序列xa(n)长度 N2=length(xb); %求序列xb(n)长度 N=N1+N2-1;n=0:N-1; n1=0:N1-1;n2=0:N2-1;

n5=0:N-4;n6=0:N-3;n7=0:N-2;n9=0:N;

yk5=circonvfft(xa,xb,N-3); %用FFT计算5点圆周卷积 yk6=circonvfft(xa,xb,N-2); %yk7=circonvfft(xa,xb,N-1); %yk8=circonvfft(xa,xb,N); %yk9=circonvfft(xa,xb,N+1); %ykl=conv(xa,xb); %figure(7); subplot(2,4,1) stem(n1,xa) title('xa(n)'); axis([0 10 0 10]) subplot(2,4,2) stem(n2,xb) title('xb(n)'); axis([0 10 0 10]) subplot(2,4,3) stem(n,ykl)

title('线性卷积'); subplot(2,4,4) stem(n5,yk5)

title('5点圆周卷积'); axis([0 10 0 100]) subplot(2,4,5) stem(n6,yk6)

title('6点圆周卷积'); axis([0 10 0 100]) subplot(2,4,6) stem(n7,yk7)

title('7点圆周卷积'); subplot(2,4,7) stem(n,yk8)

用FFT计算6点圆周卷积 用FFT计算7点圆周卷积 用FFT计算8点圆周卷积 用FFT计算9点圆周卷积

直接计算xa(n)和xb(n)的线性卷积53

title('8点圆周卷积'); subplot(2,4,8) stem(n9,yk9)

title('9点圆周卷积');

function yk=circonvfft(x1,x2,N) %用FFT计算圆周卷积子程序 xk1=fft(x1,N); xk2=fft(x2,N); xk=xk1.*xk2; yk=ifft(xk,N); yk=real(yk);

function [an,Xa]=fun4 %2N点实数序列FFT子程序 M=;n=0:2*M-1; n1=0:2:2*M-2; n2=1:2:2*M-1;

an=cos(2*pi*7*n/M)+1/2*cos(2*pi*19*n/M); %产生2N点实数序列a(n) a1=cos(2*pi*7*n1/M)+1/2*cos(2*pi*19*n1/M); %产生N点偶数点序列 a2=cos(2*pi*7*n2/M)+1/2*cos(2*pi*19*n2/M); %产生N点奇数点序列 x=a1+j*a2; %将两个奇数点和偶数点序列组合成一个复数序列 X=myditfft(x);

Xr=real(X); %求取复数序列FFT的实部 Xi=imag(X); %求取复数序列FFT的虚部 m=2:M;

Xra=Xr(m)+Xr(M+2-m)/2; %求取Xr+(m) Xra=[Xr(1) Xra];

Xrs=Xr(m)-Xr(M+2-m)/2; %求取Xr-(m) Xrs=[0 Xrs];

Xia=Xi(m)+Xi(M+2-m)/2; %求取Xi+(m) Xia=[Xi(1) Xia];

Xis=Xi(m)-Xi(M+2-m)/2; %求取Xi-(m) Xis=[0 Xis]; m=0:M-1;

co=cos(pi*m/M); si=sin(pi*m/M);

Xar1=Xra+co.*Xia-si.*Xrs; %求取2N点序列a(n)前N点的FFT实部

Xar(m)

Xai1=Xis-si.*Xia-co.*Xrs; %求取2N点序列a(n)前N点的FFT虚部Xai(m)

Xa1=Xar1+i*Xai1; Xa2=Xa1';

Xa=[Xa1,Xa2.']; %得到2N点序列a(n)的FFT Xa(m) figure(8); subplot(2,1,1) plot(n,an);

title('|a(n)|'); subplot(2,1,2) plot(n,abs(Xa)); title('|Xa(m)|');

function y=myditfft(x)

%本程序对输入序列x实现DIT-FFT基2算法，点数取大于等于x长度的 %2的幂次。y=myditfft(x)

%----------------------------------------------------------------------

m=nextpow2(x);N=2^m; %求x的长度对应的2的最低幂次m if length(x)x=[x,zeros(1,N-length(x))]; %若x的长度不是2的幂，补零到2的整数幂 end

nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1:2^m数列的倒序 y=x(nxd); %将x倒序排列作为y的初始值 for mm=1:m %将DFT作m次基2分解，从左到右，对每次分解作DFT运算 Nmr=2^mm;u=1; %旋转因子u初始化为WN^0=1

WN=exp(-i*2*pi/Nmr); %本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr) for j=1:Nmr/2 %本次跨越间隔内的各次蝶形运算

for k=j:Nmr:N %本次蝶形运算的跨越间隔为Nmr=2^mm kp=k+Nmr/2; %确定蝶形运算的对应单元下标 t=y(kp)*u; %蝶形运算的乘积项 y(kp)=y(k)-t; %蝶形运算的加法项 y(k)=y(k)+t; %蝶形运算的加法项 end

u=u*WN; %修改旋转因子，多乘一个基本DFT因子WN end end

55

实验七 IIR滤波器的实现与应用

一、实验目的

1．掌握用模拟滤波器设计IIR数字滤波器的方法。 2．应用IIR数字滤波器实现对信号的处理。 二、实验内容

设信号x(t)=cos(2π×50t)+2×cos(2π×400t)，试将它的两个频率分量分离，并绘制它们的时域波形及频谱图。 四、实验原理 滤波原理

（一）用模拟滤波器设计IIR数字滤波器

由模拟滤波器设计数字滤波器是先设计一个合适的模拟滤波器，然后将其数字化，即将S平面映射到Z平面，得到所需要的数字滤波器。由于模拟滤波器的设计方法已很成熟，故此法设计数字滤波器准确、简便，是目前最普遍采用的方法。

实现上述设计思想有多种途径，其中之一的设计步骤是

1．确定所要设计的数字滤波器的性能指标及各数字特征频率{ωk}；

2．由冲激不变法或双线性变换法的变换关系将{ωk}变换为模拟特征频率{Ωk}；

冲激不变法：Ωk=ωk/2π

双线性变换法：Ωk=c×tan（ωk/2）（预畸）

3．把要求的Lp、Hp、Bp（带通）或Bs（带阻）特征频率参数转化为模拟低通滤波器的设计参数{ΩLk }；

Hp→Lp：Ωp(通带截止频率)=Ωph；Ωst（阻带截止频率）=Ωph2/Ωsth Bp→Lp：Ωp=Ωp2－Ωp1

2222

Ωst=min{|(Ωs2－Ω0)/Ωs2|，|(Ωs1－Ω0)/Ωs1|} Ω02=Ωp2×Ωp1

Bs→Lp：Ωp=Ωp2×Ωp1/（Ωp2－Ωp1）

222222

Ωst=min{|Ω0×Ωs2/(Ω0－Ωs2) |，|Ω0×Ωs1/(Ω0－Ωs1) |}

56

Ω02=Ωp2×Ωp1

4．按指标{ΩLk }设计一个模拟低通滤波器原型，设其系统函数为Ha(S)； 5．将模拟低通滤波器原型Ha(s)转化为所需类型的模拟滤波器Ha(P)；

Lp→Hp：S=Ωph2/P

2

Lp→Bp：S=(P－Ω02)/P

2

Lp→Bs：S= PΩ02/(Ω02－P)

6．将模拟滤波器Ha(p)转化为同类型的数字滤波器H(z)。

冲激不变法：Ha(P)Ak Ppk1kN

H(z)TAkpkT1 1ezk1N双线性变换法：H(z)Ha(P)|Pc1z11z1

通常取c=2/T；T为采样周期。

（二）应用IIR数字滤波器实现对信号的处理

对信号的处理包括滤波、变换、识别、压缩等，此处主要是应用IIR数字滤波器实现对信号的分离，即利用IIR数字滤波器选出有用的信号。 四、实验步骤

1．分析实验内容，设计实验方案

例如：

2．编程实现每一个子功能并调试之

例如： 设定采样频率fS，将连续信号x(t)离散为长度为N的序列x(n)。

设：fS=1kHz，

则：T=1/ fS=0.001秒 x(n)= x(t)| t=nT

=cos(0.1πn)+cos(0.8πn)

由于1/50与1/400的最小公倍数为1/50=0.02，因而x(t)的周期T0为0.02秒，序列的长度N=0.02/0.001=20。 其Matlab程序为[xn,wk,N]=shiyan40(fs,T0)

3．选择合理的滤波器参数，编程完成信号的分离。 五、实验要求

1．针对实验内容，首先进行理论分析，给出预测结果；

57

2．熟悉实验教材§6-4的内容；

3．编写每一个子功能的程序，上机运行；

4．研究滤波器参数值、采样频率、记录长度的大小对实验结果的影响。 六、程序示例

function [yl,yh]=shiyan49 fs=1600; %采样频率 Tt=0.02; %信号周期 T0=4*Tt; %记录长度 [xn,wk,N]=shiyan40(fs,T0); M=length(wk); if M==2

rp=1;rs=80; f1=wk(1)*fs/N; f2=wk(2)*fs/N; f0=f2-f1;

fpl=f1+f0/20; fstl=f2-f0/10; %模拟低通滤波器的特征参数 [bzl,azl]=shiyan42(fpl,fstl,rp,rs,fs);

fph=f2-f0/2;fsth=f1+f0/20; %模拟高通滤波器的特征参数 [bzh,azh]=shiyan43(fph,fsth,rp,rs,fs); end

ynl=filter(bzl,azl,xn); %序列xn通过数字低通滤波器，输出为ynl ynh=filter(bzh,azh,xn); %序列xn通过数字高通滤波器，输出为ynh knl=abs(fft(ynl)); % ynl的频谱knl

kl=knl/max(knl); % ynl的幅度归一化频谱kl knh=abs(fft(ynh)); kh=knh/max(knh); T=1/fs; figure(8)

t=0:T:(T0-T);

w=0:2*pi/N:(2*pi-2*pi/N); subplot(2,2,1); plot(t,ynl); subplot(2,2,2); stem(w,kl); subplot(2,2,3); plot(t,ynh); subplot(2,2,4); stem(w,kh);

%去掉滤波输出序列的头一个周期 yl=ynl(N/4+1:N); yh=ynh(N/4+1:N);

knl=abs(fft(yl)); kl=knl/max(knl); knh=abs(fft(yh)); kh=knh/max(knh);

58

N=length(kl); t=Tt:T:(T0-T);

w=0:2*pi/N:(2*pi-2*pi/N); figure(9)

subplot(2,2,1) plot(t,yl);

Xlabel('秒'); title('低频分量时域波形') subplot(2,2,2); stem(w,kl)

Xlabel('数字频率/弧度 ');title('低频分量归一化频谱图') subplot(2,2,3) plot(t,yh);

Xlabel('秒'); title('高频分量时域波形') subplot(2,2,4) stem(w,kh)

Xlabel('数字频率/弧度 '); title('高频分量归一化频谱图')

function [xn,wk,N]=shiyan40(fs,T0) T=1/fs;

t=0:T:(T0-T);

xn=cos(100*pi*t)+2*cos(800*pi*t); xk=fft(xn); N=length(xn); i=1; wk=0;

for m=1:1:(N+1)/2

if (abs(xk(m))>0.0001) wk(i)=m-1; i=i+1; end end

n=0:1:N-1; figure(1)

subplot(2,1,1) plot(t,xn); subplot(2,1,2)

stem(n,abs(xk),'r')

59

function [Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs) %生成归一化频率的模拟低通滤波器

[N,Wn]=buttord(mp,ms,rp,rs,'s');%mp/ms:通带/阻带截止频率(弧

度%/秒)

[z,p,k]=buttap(N); [Ba,Aa]=zp2tf(z,p,k);

function [bz,az]=shiyan42(fp,fst,rp,rs,fs) %数字低通滤波器的生成

W0=[0,fp,fst,fs/2]; rr=[rp,rp,rs,rs];%设计指标 mp=fp*2*pi;

ms=fst*2*pi; %mp/ms:通带/阻带截止频率（弧度/秒） [Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);

[b,a]=lp2lp(Ba,Aa,Wn);%将归一化频率的低通滤波器转换成截止频

%率为Wn的数字低通滤波器

[bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(2);

plot(W*fs/(2*pi),20*log10(abs(H))); hold on

plot(W0,-rr,'r') Xlabel('频率/Hz')

ylabel('幅频特性/dB') grid

function [bz,az]=shiyan43(fp,fst,rp,rs,fs) %数字高通滤波器的生成

W0=[0,fst,fp,fs/2]; rr=[rs,rs,rp,rp]; %设计指标

wp=fp*2*pi/fs; ws=fst*2*pi/fs; %模拟指标数字化 C=2*fs; %频率预畸 hp=C*tan(wp/2); hs=C*tan(ws/2);

mp=hp; ms=hp^2/hs; %模拟高通指标转化为模拟低通指标 [Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);

[b,a]=lp2hp(Ba,Aa,mp); %将归一化频率的低通滤波器转换成截止

%频率为Wn的高通滤波器

[bz,az]=bilinear(b,a,fs); [H,W]=freqz(bz,az); figure(3)

60

plot(W*fs/(2*pi),20*log10(abs(H))); hold on

plot(W0,-rr,'r') Xlabel('频率/Hz')

ylabel('幅频特性/dB') grid

实验八 FIR滤波器的实现与应用

一、实验目的

1．掌握用窗函数法设计FIR数字滤波器的方法。 2．应用FIR数字滤波器实现对信号的处理。 二、实验内容

36

设调制信号x(t)=3(1+0.5cos(2π×10t)) cos(2π×10t)，试截取它的载波信号，并绘出它的时域波形及频谱图。 三、实验原理

（一）用窗函数法设计数字滤波器

一个截止频率为ω0的理想数字低通滤波器，其表达式如下：

1Hd()0故其冲激响应序列为

||00||

1hd(n)2jnH()ed0sinc(0n) 这个滤波器是物理不可实现的，加窗函数ω(n)将其截短，保留冲激响应hd(n)

的中心部分，就可获得一个线性相位的FIR滤波器h(n)。即

h(n)=hd(n) ω(n)

最简单的窗函数是矩形窗，其ω(n)=RN(n)，问题是在通带和阻带的边缘存在波动（吉布斯效应），该现象不能通过增加滤波器的长度来消除，只能采用非线性窗如汉宁窗、海明窗等。

利用窗函数法设计FIR滤波器的步骤如下：

1．确定所要设计的数字滤波器的性能指标，确定理想的频率响应函数Hd(ejω)并计算hd(n)；

2．由最小阻带衰减选择窗的形状，由过渡带的宽度确定N的大小（参见教材的

61

窗函数基本参数表）；

3．求得所设计的滤波器的单位抽样响应h(n)；

4．计算H(ejω)=DTFT[h(n)]，检验是否满足设计指标。

在Matlab环境下，函数fir1就是基于这种窗函数方法的。给定滤波器的阶数和期望的理想滤波器的描述，它返回理想滤波器的加窗逆富立叶变换。因而用窗函数法设计FIR滤波器的步骤简化如下： 1）确定所要设计的数字滤波器的性能指标； 2）由最小阻带衰减选择窗的形状； 3）由过渡带的宽度确定N的大小； 4）调用函数fir1生成FIR滤波器。

（二）应用FIR数字滤波器实现对信号的处理

此处主要是应用FIR数字滤波器实现对调制信号的解调，即利用FIR数字滤波器选出基带信号。 四、实验步骤

1．分析实验内容，设计实验方案 例如：

2．编程实现每一个子功能并调试之 例如：FIR数字带通滤波器的设计

FIR数字带通滤波器设计指标：通带截止频率ωp1、ωp2，阻带截止频率ωs1、ωs2，阻带最小衰减rs。

首先根据rs的大小选择窗形，计算过渡带ω=min(ωp1－ωs1，ωs2－ωp2)，由此确定N的大小，程序如

function [windowxing,jieshu]=shiyan51(rs,wguo) 在此基础上确定FIR数字带通滤波器的系数b，编写程序如 function b=shiyan52(wp,ws,rs)

3．选择合理的滤波器参数，编程完成信号的截取。 五、实验要求

1．针对实验内容，首先进行理论分析，给出预测结果； 2．熟悉实验教材§6-5的内容；

3．编写每一个子功能的程序，上机运行；

4．研究滤波器参数值、采样频率、记录长度的大小对实验结果的影响。 六、程序示例

function yn=shiyan59

k=10; %记录长度为信号周期的k倍 [xn,freq,N,fs,Tt]=shiyan50(k); T=1/fs;

fp=[9800,10200];

62

fst=[9200,10800]; rs=60;

wp=(2*pi/fs).*fp; ws=(2*pi/fs).*fst;

b=shiyan52(wp,ws,rs); %生成带通数字滤波器 yn=filter(b,1,xn); %对信号xn进行滤波 kk=5; %kkynn=yn(kk*N/k+1:N); %去掉yn的前kk个周期 knn=abs(fft(ynn)); %计算滤波后信号ynn的频谱 M=length(ynn);

t=kk*Tt:T:(k*Tt-T); m=[0:M-1]; f=m.*(fs/M); figure(9)

subplot(2,1,1) plot(t,ynn);

Xlabel('时间/S');

legend('输出信号的波形'); subplot(2,1,2)

stem(f(1:M/2),knn(1:M/2)); Xlabel('频率/Hz');

legend('输出信号的频谱');

function [xn,freq,N,fs,Tt]=shiyan50(k) %生成信号序列 fs=40000; %设定抽样频率fs/Hz Tt=0.0001; %信号周期Tt/S T=1/fs;

T0=k*Tt; %设定记录长度T0/S t=0:T:(T0-T);

xn=2.*(1+0.5.*cos(2000*pi*t)).*cos(20000*pi*t); %信号序列xn xk=fft(xn); %计算信号序列xn的频谱xk N=length(xn); %计算信号序列的长度 i=1;

m=[0:N-1];

f=m.*(fs/N); %将数字频率转换成模拟线性频率/Hz freq=0; % freq记录信号xn的频率分量 for m=1:1:(N+1)/2

if (abs(xk(m))>0.0001)

63

freq(i)=m-1; i=i+1; end end

figure(1)

subplot(2,1,1)

plot(t,xn);Xlabel('时间');legend('信号波形'); subplot(2,1,2)

stem(f(1:N/2),xk(1:N/2),'r');

Xlabel('频率/Hz');legend('信号频谱'); grid

function [windowxing,jieshu]=shiyan51(rs,wguo)

%根据rs、过渡带wguo的大小选择窗形windowxing，确定阶数jieshu if rs>0

N(1)=ceil(4*pi/wguo);as(1)=21;str(1).window=boxcar(N(1)+1); N(2)=ceil(8*pi/wguo);as(2)=25;str(2).window=triang(N(2)+1); N(3)=ceil(8*pi/wguo);as(3)=44;str(3).window=hanning(N(3)+1); N(4)=ceil(8*pi/wguo);as(4)=53;str(4).window=hamming(N(4)+1); N(5)=ceil(12*pi/wguo);as(5)=74;str(5).window=blackman(N(5)+1); N(6)=ceil(10*pi/wguo);as(6)=80;str(6).window=kaiser(N(6)+1,7.865); i=0;

while i<6 i=i+1; if as(i)N(i)=1000000; end end

[Nmin,m]=min(N)

windowxing=str(m).window; jieshu=Nmin; else

error('rs > 0') end

function b=shiyan52(wp,ws,rs) %生成滤波器

wguo=min(abs(wp-ws));

wn=(wp+ws)./2;

[windowxing,N]=shiyan51(rs,wguo); b=fir1(N,wn/pi,windowxing); figure(2)

[H,W]=freqz(b,1);

plot(W/pi,20*log10(abs(H))); grid

65