Matlab实践报告模板

计算机科学与技术学院

课程设计报告

（ 2015— 2016学年 第 一 学期）

课程名称： MATLAB应用实践

2016 年 1月

实 习 报 告

实验一 基础练习

实验内容

练习1 基础操作和简单语句输入 1.用函数sym定义符号矩阵：

>> sym_matrix=('[a b c:Jack Help_Me NO_WAY]')

sym_matrix =

[a b c:Jack Help_Me NO_WAY] 2.用函数syms定义符号矩阵 >> syms a b c;

>> M1=sym('Classical'); >> M2=sym('Jazz'); >> M3=sym('Blues');

>> A=[a b c;M1,M2,M3;sym([2 3 5])] A =

[ a, b, c] [ Classical, Jazz, Blues] [ 2, 3, 5] 练习要求：

1.练习数据和符号的输入方式，自定义数据和变量将前面的命令在命令窗口中执行通过 2.输入 A=[7 1 5；2 5 6；3 1 5]，B=[1 1 1; 2 2 2; 3 3 3]，在命令窗口中执行下列表达式，掌握其含义：A(2, 3) A(:,2) A(3,:) A(:,1:2:3) A(:,3).*B(:,2)

- 1 -

实 习 报 告

A(:,3)*B(2,:) A*B A.*B A^2 A.^2 B/A B./A

>> A=[7 1 5;2 5 6;3 1 5]; >> B=[1 1 1;2 2 2;3 3 3]; >> A(2,3)

ans =

6

>> A(:,2)

ans =

1 5 1

>> A(3,:)

ans =

3 1 5

>> A(:,1:2:3)

ans =

7 5 2 6 3 5

>> A(:,3).*B(:,2)

ans =

5 12 15

- 2 -

实 习 报 告

>> A(:,3)*B(2,:)

ans =

10 10 10 12 12 12 10 10 10

>> A*B

ans =

24 24 24 30 30 30 20 20 20

>> A.*B

ans =

7 1 5 4 10 12 9 3 15

>> A^2

ans =

66 17 66 42 33 70 38 13 46

>> A.^2

ans =

49 1 25 4 25 36 9 1 25

>> B/A

- 3 -

实 习 报 告

ans =

0.1842 0.2105 -0.2368 0.3684 0.4211 -0.4737 0.5526 0.6316 -0.7105

>> B./A

ans =

0.1429 1.0000 0.2000 1.0000 0.4000 0.3333

1.0000 3.0000 0.6000

3.输入C=1：2：20，则C（i）表示什么？其中i=1，2，3，……，10； >> C=1:2:20, C =

1 3 5 7 9 11 13 15 17 19 则C(i)表示 1 3 5 7 9 4.查找已创建变量的信息，删除无用的变量；

练习2 编程

2.1 无条件循环

当需要无条件重复执行某些命令时,可以使用for循环:for 循环变量t=表达式1 : 达式2 :表达式 3 语句体 end

说明:表达式 1 为循环初值,表达式 2 为步长,表达式 3 为循环终值;当表达式 2 省略时则

默认步长为 1;for 语句允许嵌套.

例：矩阵输入程序生成 3×4 阶的 Hiltber 矩阵.

m=input('矩阵行数:m=');for i=1 : 3 n= input('矩阵列数:n='); for j=1 : 4 for i=1:m

H(i,j)=1/(i+j-1); for j=1:n

end disp(['输入第',num2str(i),'行,第', num2str(j),'列元素']) end A(i, j) = input (' ') end end

2.2 条件循环

1) if-else-then 语句

If-else-then 语句的常使用三种形式为:(1) if 逻辑表达式 (3) if 逻辑表达式 1语句体 语句体 1end elseif 逻辑表达式 2 语句体 2(2) if 逻辑表达式 1 elseif 逻辑表达

- 4 -

实 习 报 告

式 3 语句体 …else else 语句体 2 语句体 nend end

2) while 循环语句 while 循环的一般使用形式为:while 表达式 语句体 end 2.3 分支结构

若需要对不同的情形执行不同的操作,可用 switch 分支语句: switch 表达式(标量或字符串) case 值 1 语句体 1 case 值 2 语句体 2 Otherwise 语句体 n end

说明:当表达式不是\"case\"所列值时,执行 otherwise 语句体。 ① 建立 M 文件

将多个可执行的系统命令,用文本编辑器编辑后并存放在后缀为 .m 的文件中,若在 ATLAB 命令窗口中输入该 m-文件的文件名(不跟后缀.m!),即可依次执行该文件中的多 个命令.这个后缀为.m 的文件,也称为 Matlab 的脚本文件(Script File)。 注意:文件存放路径必须在 Matlab 能搜索的范围内. ② 建立函数文件

对于一些特殊用户函数,系统提供了一个用于创建用户函数的命令 function,以备用户随时调用 1. 格式:

function [输出变量列表]=fun_name(输入变量列表) 用户自定义的函数体

2. 函数文件名为:fun_name,注意:保存时文件名与函数名最好相同; 3. 存储路径:最好在系统的搜索路径上。 4. 调用方法:输出参量=fun_name (输入变量) 例:

计算 s = n!,在文本编辑器中输入: function s=pp(n); s=1;

for i=1:n s=s*i; end s;

在 MATLAB 命令窗口中输入:s=pp(5) 结果为 s = 120 练习要求：

1． 编写程序,计算 1+3+5+7+…+(2n+1)的值(用 input 语句输入 n 值). n=input('input: n='); s=1;

for i=1:n s=s+2*i+1; end s

- 5 -

实 习 报 告

>> Untitled11 input: n=7 s =

64 2． 编写分段函数

的函数文件，存放于文件ff.m中，计算计算出

值.

解：function s=f(x) if 0elseif 1s;

命令行输出结果为：>> f(-3) >>f(2^0.5) >>f(inf)

ans = ans= ans=

0 1393/985 0

的

练习3 矩阵计算 3.1 矩阵的创建 1.加、减运算

例：在Matlab编辑器中建立m文件：LX0701.m A=[1,1,1;1,2,3;1,3,6] B=[8,1,6;3,5,7;4,9,2] C=A+B D=A-B

>> LX0701 A =

1 1 1 1 2 3 1 3 6 B =

- 6 -

实 习 报 告

8 1 6 3 5 7 4 9 2 C =

9 2 7 4 7 10 5 12 8 D =

-7 0 -5 -2 -3 -4

-3 -6 4 2.乘法

（1）两个矩阵相乘

例：在Mtalab编辑器中建立M文件：LX0702.m X=[2 3 4 5 1 2 2 1]; Y=[0 1 1 1 1 0 0 0 1 1 0 0]; Z=X*Y

>> LX0702 Z =

8 5 6 3 3 3

（2）矩阵的数乘：数乘矩阵 上例中：a=2*X X=[2 3 4 5 1 2 2 1]; Y=[0 1 1 1 1 0 0 0 1 1 0 0]; a=2*X

>> LX0702

- 7 -

实 习 报 告

a =

4 6 8 10 2 4 4 2

（3）向量的点乘（内积）：维数相同的两个向量的点乘。 命令：dot 向量点乘函数 例：X=[-1 0 2]; Y=[-2 -1 1]; Z=dot(X,Y) >> Untitled5 Z =

4

（4）向量叉乘

命令 cross 向量叉乘函数 例：a=[1 2 3]; b=[4 5 6]; c=cross(a,b) >> LX0703 c =

-3 6 -3 （5）混合积

例：计算向量a=(1 2 3);b=(4 5 6);和c=(-3 6 -3)的混合积 a=[1 2 3];b=[4 5 6];c=[-3 6 -3]; x=dot(a,cross(b,c)) >> LX0704 x =

54

注意：先叉乘后点乘，顺序不可颠倒。 3.矩阵的除法

Matlab 提供了两种除法运算：左除（\\）和右除（/）。一般情况下，x=a\\b 是方程 a*x =b 的解，而 x=b/a 是方程 x*a=b 的解 例：a=[1 2 3; 4 2 6; 7 4 9] b=[4; 1; 2]; x=a\\b

则显示：x= -1.5000 2.0000 0.5000

如果 a 为非奇异矩阵，则 a\\b 和 b/a 可通过 a 的逆矩阵与 b 阵得到：

- 8 -

实 习 报 告

a\\b = inv(a)*b b/a = b*inv(a)

4.矩阵乘方 运算符：^ 运算规则：

（1）当 A 为方阵，p 为大于 0 的整数时，A^P 表示 A 的 P 次方，即 A 自乘 P 次；p 为小于 0 的整数时，A^P 表示 A-1 的 P 次方。

（2）当 A 为方阵，p 为非整数时，则其中 V 为 A 的特征向量，为特征值矩阵 5.矩阵的转置 运算符：′

运算规则：与线性代数中矩阵的转置相同。 6.矩阵的逆矩阵

例1：A=[1 2 3;2 2 1;3 4 3]; inv(A)

>> LX07051

ans =

1.0000 3.0000 -2.0000 -1.5000 -3.0000 2.5000 1.0000 1.0000 -1.0000

例2：B=[1,2,3,1,0,0;2,2,1,0,1,0;3,4,3,0,0,1]; C=rref(B) %化行最简形 X=C(:,4:6) >> LX07052 C =

1.0000 0 0 1.0000 3.0000 -2.0000 0 1.0000 0 -1.5000 -3.0000 2.5000 0 0 1.0000 1.0000 1.0000 -1.0000 X =

1.0000 3.0000 -2.0000 -1.5000 -3.0000 2.5000

1.0000 1.0000 -1.0000 7.方阵的行列式

命令： det 计算行列式的值 例 ：A=[1 2 3;2 2 1;3 4 3]; D=det(A)

- 9 -

实 习 报 告

>> LX0706 D =

2.0000

3.2 符号矩阵的运算 1.符号矩阵的四则运算

把符号矩阵的四则运算简化为与数值矩阵完全相同的运算方式，其运算符为：加（＋），减（－）、乘（×）、除（/、\\）等或：符号矩阵的和（symadd），差（symsub），乘 (symmul)。 例：A=sym('[1/x,1/(x+1);1/(x+2),1/(x+3)]'); B=sym('[x,1;x+2,0]') C=B-A D=A\\B

>> LX7 B =

[ x, 1] [ x + 2, 0] C =

[ x - 1/x, 1 - 1/(x + 1)] [ x - 1/(x + 2) + 2, -1/(x + 3)] D =

[ -x*(2*x^2 + 7*x + 6), (x*(x^2 + 3*x + 2))/2] [ 2*(x + 1)^2*(x + 3), -(x*(x + 1)*(x + 3))/2] LX7 B =

[ x, 1] [ x + 2, 0] C =

[ x - 1/x, 1 - 1/(x + 1)] [ x - 1/(x + 2) + 2, -1/(x + 3)]

- 10 -

实 习 报 告

D =

[ -x*(2*x^2 + 7*x + 6), (x*(x^2 + 3*x + 2))/2] [ 2*(x + 1)^2*(x + 3), -(x*(x + 1)*(x + 3))/2] 2.其他基本运算

符号矩阵的其他一些基本运算包括转置（'）、行列式（det）、逆（inv）、秩（rank）、 幂（^）和指数（exp和expm）等都与数值矩阵相同 3.符号矩阵的简化

符号工具箱中提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作函数。（1）

(1)因式分解

命令：factor 符号表达式因式分解函数 格式：factor(s)

说明：s为符号矩阵或符号表达式。常用于多项式的因式分解。 例：将x-1分解因式 >> syms x

>> factor(x^9-1) ans =

(x - 1)*(x^2 + x + 1)*(x^6 + x^3 + 1) (2)符号矩阵的展开

命令：expand 符号表达式展开函数 格式：expand(s)

说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数

函数和对数函数的展开中

例3-10 将(x+1)^3、sin(x+y)展开 syms x y

p=expand((x+1)^3) q=expand(sin(x+y))

>> LX0710 p =

x^3 + 3*x^2 + 3*x + 1 q =

cos(x)*sin(y) + cos(y)*sin(x) (3)同类式合并

命令：Collect 合并系数函数

格式：Collect(s,v) 将s中的变量v的同幂项系数合并。

- 11 -

实 习 报 告

Collect(s) s — 矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合

并。

(4)符号简化

命令：simple或simplify 寻找符号矩阵或符号表达式的最简型 格式：Simple(s) s — 矩阵或表达式

说明：Simple(s)将表达式s的长度化到最短。若还想让表达式更加精美，可使用函数Pretty。

格式：Pretty(s) 使表达式s更加精美.

练习 4 、秩与线性相关性

4.1 矩阵和向量组的秩以及向量组的线性相关性。 矩阵A的秩是矩阵A中最高阶非零子式的阶数；向量组的秩通常由该向量组构成的矩阵来计

算。

命令：rank

格式：rank(A) A为矩阵式向量组构成的矩阵

例4-1 求向量组α1＝（1 -2 2 3）α2＝（-2 4 -1 3）α3＝（-1 2 0 3）α4＝（0 6 2 3）α5＝（2 -6 3 4）的秩，并判断其线性相关性。 A=[1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4]; B=rank(A)

>> Untitled7 B =

3

由于秩为3<向量个数，因此向量组线性相关。 4.2 向量组的最大无关组 命令：rref

格式：rref(A) A为矩阵

例4-2 求向量组a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)的一个最大无关组。 a1=[1 -2 2 3]'; a2=[-2 4 -1 3]'; a3=[-1 2 0 3]'; a4=[0 6 2 3]'; a5=[2 -6 3 4]'; A=[a1 a2 a3 a4 a5]

format rat %以有理格式输出 B=rref(A) %求A的行最简形

>> LX0715 A =

- 12 -

实 习 报 告

1 -2 -1 0 2 -2 4 2 6 -6 2 -1 0 2 3 3 3 3 3 4 B =

1 0 1/3 0 16/9 0 1 2/3 0 -1/9 0 0 0 1 -1/3 0 0 0 0 0 从B中可以得到：向量a1 a2 a4为其中一个最大无关组。

练习5 线性方程的组的求解

5.1求线性方程组的唯一解或特解（第一类问题） 这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 —— 直接法；另一类是解大型稀疏矩阵 —— 迭代法。

利用矩阵除法求线性方程组的特解（或一个解） 方程：AX=b 解法：X=A\\b

解：A=[5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5]; B=[1 0 0 0 1]';

R_A=rank(A) %求秩 X=A\\B %求解

>> LX0716

R_A =

5

- 13 -

实 习 报 告

X =

1507/665 -229/133 37/35 -79/133

212/665 这就是方程组的解

解：A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; B=[1 4 0]'; X=A\\B

>> Untitled10 X =

0 0 -8/15 3/5

- 14 -

实 习 报 告

实验二 离散时间信号分析

一、实验目的

1．掌握各种常用的序列，理解其数学表达式和波形表示。 2．掌握在计算机中生成及绘制数字信号波形的方法。

3．掌握序列的相加、相乘、移位、反褶等基本运算及计算机实现与作用。 4．掌握线性卷积软件实现的方法。

5．掌握计算机的使用方法和常用系统软件及应用软件的使用。

6．通过编程，上机调试程序，进一步增强使用计算机解决问题的能力。 二、实验原理 1．序列的基本概念

离散时间信号在数学上可用时间序列 {x(n)} 来表示，其中x(n) 代表序列的第 n个数字，n 代表时间的序列，n 的取值范围为 -常用序列有：单位脉冲序列（单位抽样）δ(n)、单位阶跃序列u(n)、矩形序列Rn(n)、 实指数序列、复指数序列、正弦型序列等。 3．序列的基本运算

序列的运算包括移位、反褶、和、积、标乘、累加、差分运算等。 4．序列的卷积运算 Y(n)=

m-x(m)h(nm)x(n)*h(n)   

上式的运算关系称为卷积运算，式中  代表两个序列卷积运算。两个序列的卷积是一个序列与另一个序列反褶后逐次移位乘积之和，故称为离散卷积，也称两序列的线性卷积。其计算的过程包括以下 4 个步骤。

（1）反褶：先将x(n)和h(n)的变量 n 换成 m ，变成x(m)和h(m)，再将h(-m) 以纵 轴为对称轴反褶成h(-m)。

（2）移位：将h(-m)移位 n ，得h(n-m)。当 n 为正数时，右移 n 位；当 n 为负数时，

- 15 -

实 习 报 告

左移 n 位。

（3）相乘：将h(n-m)和x(m)的对应点值相乘。

（4）求和：将以上所有对应点的乘积累加起来，即得y(n)。 三、主要实验仪器及材料 微型计算机 、 Matlab6.5 教学版 四、实验要求

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

2．列出计算卷积的公式，画出程序框图，并列出实验程序清单（可略）（包括必要的程 序说明）。

3．记录调试运行情况及所遇问题的解决方法。 4．给出实验结果，并对结果作出分析。 5．简要回答思考题。 五、实验内容 1．知识准备

认真复习以上基础理论，理解本实验所用到的实验原理。 2．离散时间信号（序列）的产生

利用 MATLAB 语言编程产生和绘制下列有限长序列： （1）单位脉冲序列δ(n） 在matlab编辑器中建立M文件： n1 = -5:5; x1=[(n1-0)==0]; stem(n1,x1);

title('单位脉冲序列'); axis([-5,5,0,1]);

- 16 -

实 习 报 告

（2）单位阶跃序列u(n) 在matlab编辑器中建立M文件： n1=-5:5; x1=[(n1-0)>=0]; stem(n1,x1);

title('阶跃序列');

axis([-5,5,0,1.1]);

（3）矩形序列Rs(n)

在matlab编辑器中建立M文件： n1 = -5:5; x1=[(n1+1)>=0]; x2=[(n1-4>=0)]; x3=x1-x2;

- 17 -

实 习 报 告

stem(n1,x3);title('矩形序列') axis([-5,5,0,1.1])

（4）正弦型序列 x(n)=Asin（n/5/3） 在matlab编辑器中建立M文件： n = 0:40;

x1=sin(0.2*pi*n+(1/3)*pi); subplot(2,1,2);

stem(n,x1);title('正弦序列'）

（5）任意序列

X(n)=δ(n)+2δ(n-1)+3δ(n-2)+4δ(n-3)+5δ(n-4) 在matlab编辑器中建立M文件：

- 18 -

实 习 报 告

n=-5:10

x=[zeros(1,5),1,2,3,4,5,zeros(1,6)]; stem(n,x,'filled'); title('任意序列'); x1able('时间(n)'); yiable('序列值x(n)');

Y(n)=δ(n)+2δ(n-1)+δ(n-2)+2δ(n-3) 在matlab编辑器中建立M文件： n=-5:10

x=[zeros(1,5),1,2,1,2,zeros(1,7)]; stem(n,x,'filled'); title('任意序列'); x1able('时间(n)'); yiable('序列值h(n)');

- 19 -

实 习 报 告

   3．序列的运算

利用 MATLAB 语言编程完成上述两序列的移位、反褶、和、积、标乘、累加等运算， 并绘制运算后序列的波形。 （1）移位

在matlab编辑器中建立M文件： N=10; n=-N:N; x1=[n==0]; x2=[(n-1)==0]; x3=[(n-2)==0]; x4=[(n-3)==0]; x5=[(n-4)==0]; n1=n+2; n2=n-1; n3=-n;

xn=x1+2*x2+3*x3+4*x4+5*x5; hn=x1+2*x2+x3+2*x4; subplot(2,2,1); stem(n,xn,'filled') axis([-N,N,0,5]) xlabel('n')

- 20 -

实 习 报 告

ylabel('xn') subplot(2,2,3); stem(n1,xn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('x(n-2)') title('移位') subplot(2,2,2); stem(n,hn,'filled') axis([-N,N,0,2]) xlabel('n') ylabel('hn') subplot(2,2,4); stem(n2,hn,'filled') axis([-N,N,0,2]) xlabel('n') ylabel('h(n+1)') title('移位')

（2）反褶 N=10; n=-N:N; x1=[n==0];

x2=[(n-1)==0];

- 21 -

实 习 报 告

x3=[(n-2)==0]; x4=[(n-3)==0]; x5=[(n-4)==0]; n1=n+2; n2=n-1; n3=-n;

xn=x1+2*x2+3*x3+4*x4+5*x5; hn=x1+2*x2+x3+2*x4; subplot(2,2,1); stem(n,xn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('xn') subplot(2,2,2); stem(n3,xn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('x(-n)') title('反褶') subplot(2,2,3); stem(n,hn,'filled') axis([-N,N,0,2]) xlabel('n') ylabel('hn') subplot(2,2,4); stem(n3,hn,'filled') axis([-N,N,0,2]) xlabel('n') ylabel('h(-n)') title('反褶')

- 22 -

实 习 报 告

（3）和 N=10; n=-N:N;

x1=[n==0]; x2=[(n-1)==0]; x3=[(n-2)==0]; x4=[(n-3)==0]; x5=[(n-4)==0];

xn=x1+2*x2+3*x3+4*x4+5*x5; hn=x1+2*x2+x3+2*x4; sn=xn+hn; subplot(3,1,1); stem(n,xn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('xn') subplot(3,1,2); stem(n,hn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('hn') subplot(3,1,3); stem(n,sn,'filled') axis([-N,N,0,8]) xlabel('n')

- 23 -

实 习 报 告

ylabel('sn（和）')

（4）积 N=10; n=-N:N;

x1=[n==0]; x2=[(n-1)==0]; x3=[(n-2)==0]; x4=[(n-3)==0]; x5=[(n-4)==0];

xn=x1+2*x2+3*x3+4*x4+5*x5; hn=x1+2*x2+x3+2*x4; mn=xn.*hn; subplot(3,1,1); stem(n,xn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('xn') subplot(3,1,2); stem(n,hn,'filled') axis([-N,N,0,5]) xlabel('n') ylabel('hn') subplot(3,1,3); stem(n,mn,'filled') axis([-N,N,0,8])

- 24 -

实 习 报 告

xlabel('n') ylabel('mn')

（5）标乘 （6）累加 3.卷积 N=10; n=0:N;

x1=[n==0]; x2=[(n-1)==0]; x3=[(n-2)==0]; x4=[(n-3)==0]; x5=[(n-4)==0];

xn=x1+2*x2+3*x3+4*x4+5*x5; hn=x1+2*x2+x3+2*x4; k1=length(xn); k2=length(hn); L=k1+k2-1; y0=zeros(1,L); for i=1:L for m=1:k2 k=i-m+1; if(k>=1&&k<=k1)

y0(i)=y0(i)+xn(m)*hn(k); end end

- 25 -

实 习 报 告

end

yn=conv(xn,hn); n=0:k1-1; subplot(4,1,1); stem(n,xn,'filled') axis([0,N,0,5]) xlabel('n') ylabel('xn') n=0:k2-1; subplot(4,1,2); stem(n,hn,'filled') axis([0,N,0,2]) xlabel('n') ylabel('hn') n=0:L-1; subplot(4,1,3); stem(n,y0,'filled') axis([0,N,0,20]) xlabel('n') ylabel('y0') n=0:L-1; subplot(4,1,4); stem(n,y0,'filled') axis([0,N,0,20]) xlabel('n') ylabel('yn=xn*hn')

- 26 -

实 习 报 告

五、思考题

1．如何产生方波信号序列和锯齿波信号序列？ 方波信号序列 N=10;

A=zeros(1,2); B=ones(1,2); y=ones(1,N); for i=1:N

if rem(i,2)==0 y(i)=0; end end j=0; k=1; for i=1:N m=1;

for x=j:0.01:i

z=ones(1,(i-j)/0.01); if m<=100

f(k)=y(i)*z(1,m); k=k+1; m=m+1; end end j=j+1; end plot(f)

axis([0,1000,-1,2]) title('方波序列')

- 27 -

实 习 报 告

锯齿波信号 N=5;

x=-N:0.01:N; y=1-mod(x,1); plot(x,y) axis([-N,N,0,1]) title('锯齿波')

- 28 -

实 习 报 告

2．实验中所产生的正弦序列的频率是多少？是否是周期序列？ 答：F=0.1Hz,是周期序列

- 29 -

实 习 报 告

实验三 数字信号处理综合设计

一、实验目的

1．学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法； 2．掌握在 Windows 环境下语音信号采集的方法； 3．掌握数字信号处理的基本概念、基本理论和基本方法； 4．掌握 MATLAB 设计 FIR 和 IIR 数字滤波器的方法； 5．学会用 MATLAB 对信号进行分析和处理。 二、实验原理

参考《数字信号处理》教材。 三、主要实验仪器及材料 微型计算机 、 Matlab6.5 教学版 四、实验内容 1．语音信号的采集

要求利用 windows 下的录音机或其他软件，录制一段自己的话音，时间控制在 1 秒左 右。然后在 MATLAB 软件平台下，利用函数 wavread 对语音信号进行采样，记住采样频率

和采样点数。通过 wavread 函数的使用，要求理解采样频率、采样位数等概念。 wavread 函数调用格式：

y=wavread(file)，读取 file 所规定的 wav 文件，返回采样值放在向量 y 中。 [y,fs,nbits]=wavread(file)，采样值放在向量 y 中，fs 表示采样频率（Hz），nbits 表示采 样位数。

y=wavread(file,N)，读取前 N 点的采样值放在向量 y 中。

y=wavread(file,[N1,N2])，读取从 N1 点到 N2 点的采样值放在向量 y 中。 2．语音信号的频谱分析

要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在 MATLAB 中， 可以利用函数 fft 对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱

- 30 -

实 习 报 告

特性 的理解。

3．设计数字滤波器和画出频率响应

根据语音信号的特点给出有关滤波器的性能指标：

1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB； 2）高通滤波器性能指标，fc=4800 Hz，fp=5000 Hz As=100dB，Ap=1dB；

3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz， As=100dB，Ap=1dB。

要求学生首先用窗函数法设计上面要求的三种滤波器，在 MATLAB 中，可以利用函数 fir1 设计 FIR 滤波器；然后在用双线性变换法设计上面要求的三种滤波器，在 MATLAB 中，

可以利用函数 butte、cheby1 和 ellip 设计 IIR 滤波器；最后，利用 MATLAB 中的函数 freqz

画出各滤波器的频率响应。 4．用滤波器对信号进行滤波

比较两种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波，在 MATLAB 中，FIR 滤波器利用函数 fftfilt 对信号进行滤波，IIR 滤波器利用函数 filter 对信号 进行滤波。

5．比较滤波前后语音信号的波形及频谱

要求在一个窗口同时画出滤波前后的波形及频谱。 6．回放语音信号

在 MATLAB 中，函数 sound 可以对声音进行回放。其调用格式： sound(x,fs,bits)；

可以感觉滤波前后的声音有变化。 五、实验思考

1．双线性变换法中 Ω 和 ω 之间的关系是非线性的，在实验中你注意到这种非线性关系

了吗？从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 2．能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？

- 31 -

实 习 报 告

六、实验要求

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

2.按照实验步骤及要求，比较各种情况下的滤波性能。 3.总结实验所得主要结论。 4.简要回答思考题。 七、实验结果及实验分析 1、语音信号的采集及频谱分析 Matlab的程序：

function rr1011

[y, fs, bits]=wavread('pch.wav'); sound(y,fs,bits); n=length(y); Y=fft(y,n); subplot(2,1,1); plot(y);

title('语音信号波形'); xlabel('时间'); ylabel('幅值') subplot(2,1,2); plot(abs(Y)); title('语音信号频谱'); xlabel('时间'); ylabel('幅值');

2、信号的波形和频谱

- 32 -

实 习 报 告

3、设计数字滤波器和画出频率响应

1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB； （1）用kaiser窗函数设计FIR低通滤波器： Matlab的程序： %FIR低通滤波器 fp=1000; fc=1200; As=100; Ap=1; fs=8000;

wc=2*fc/fs;wp=2*fp/fs;

N=ceil((As-7.95)/(14.36*(wc-wp)/2))+1;

Window=kaiser(N+1,beta) b=fir1(N,wc,Window); freqz(b,1,512,fs)

- 33 -

实 习 报 告

%plot(w*8000*0.5/pi,abs(h)); title('FIR低通滤波器（凯撒窗口）'); xlabel('频率/Hz'); ylabel('幅值');

实验分析：由图像可以看出：窗函数低通滤波器有很大的阻带条纹，在通带基本上没有条纹。

（2）用双线性法设计IIR低通滤波器 Matlab程序：

fp=1000;fc=1200;As=100;Ap=1;ffs=22050 wp=2*fp/ffs;wc=2*fc/ffs; [n,wn]=ellipord(wp,wc,Ap,As); [num,den]=ellip(n,Ap,As,wn); freqz(num,den,512,ffs); title('IIR低通滤波器'); xlabel('频率/Hz'); ylabel('幅值');

- 34 -

实 习 报 告

2）高通滤波器性能指标，fc=4800 Hz，fp=5000 Hz As=100dB，Ap=1dB； （1）用kaiser窗函数设计FIR高通滤波器： Matlab的程序：

%高通滤波器 fp=1000; fc=1200; As=100; Ap=1; fs=8000;

wc=2*fc/fs;wp=2*fp/fs;

N=ceil((As-7.95)/(14.36*(wc-wp)/2))+1; Window=kaiser(N+1,beta) b=fir1(N,wc,Window); freqz(b,1,512,fs);

title('FIR高通滤波器（凯撒窗口）'); xlabel('频率/Hz'); ylabel('幅值');

- 35 -

实 习 报 告

实验分析：由图像可知：kaiser窗函数高通滤波器在通带基本没波纹，过渡带较窄，基本符合要求。

（2）用双线性法设计IIR高通滤波器： Matlab程序;

%高通滤波器

fp=5000;fc=4800;As=100;Ap=1;ffs=22050 wp=2*fp/ffs;wc=2*fc/ffs; [n,wn]=ellipord(wp,wc,Ap,As); [num,den]=ellip(n,Ap,As,wn); freqz(num,den,512,ffs); title('IIR高通滤波器'); xlabel('频率/Hz'); ylabel('幅值');

- 36 -

实 习 报 告

3)带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz， As=100dB，Ap=1dB。

(1)用kaiser窗函数设计FIR带通滤波器： Matlab程序：

%FIR带通滤波器 fp1=1200;fp2=3000; fc1=1000;fc2=3000 As=100;Ap=1; fs=22050; wc1=2*pi*fc1/fs;

wp1=2*pi*fp1/fs;wc2=2*pi*fc2/fs;wp2=2*pi*fp2/fs; deltaw1=(wc1+wp1)/(2*pi);

deltaw2=(wc2+wp2)/(2*pi);ws=wp1-wc1; deltaw=[deltaw1,deltaw2]; n0=ceil((As-7.95)/2.285/ws); beta=0.1102*(As-8.7);

- 37 -

实 习 报 告

wdkai=(kaiser(n0+1,beta)); b0=fir1(n0,deltaw,wdkai); freqz(b0,1);

plot(w*fs*0.5/pi,20*log10(abs(h))); title('FIR带通滤波器（凯撒窗口）'); xlabel('频率/Hz'); ylabel('幅值');

实验分析：由图像可知：kaiser窗函数带通滤波器具有较窄的过渡带，通带上基本没波纹，阻带上波纹较大。

（2）用双线性法设计巴特沃斯IIR带通滤波器 Matlab程序：

%带通滤波器

fp=[1200,3000];fc=[1000,3200];As=100;Ap=1;ffs=22050 wp=2*fp/ffs;wc=2*fc/ffs; [n,wn]=ellipord(wp,wc,Ap,As); [num,den]=ellip(n,Ap,As,wn);

- 38 -

实 习 报 告

freqz(num,den,256,ffs); title('IIR带通滤波器'); xlabel('频率/Hz'); ylabel('幅值');

（3）契比雪夫设计IIR带通滤波器 Matlab程序：

%用Chebyshev I设计IIR带通滤波器： fp1=1200; fp2=3000; fc1=1000; fc2=3200; Fs=22050; wc1=fc1*2/Fs; wc2=fc2*2/Fs; wp1=fp1*2/Fs; wp2=fp2*2/Fs; Rp=1;

- 39 -

实 习 报 告

Rs=10; wp=[wp1,wp2]; wc=[wc1,wc2];

[N,wn]=cheb1ord(wp,wc,Rp,Rs,'s'); [bz,az]=cheby1(N,Rp,wc,'bandpass'); [H,W]=freqz(bz,az); plot(W*Fs/(2*pi),abs(H)); grid;

xlabel('频率/HZ'); ylabel('幅度');

title('chebyshev1带通IIR数字滤波器');

4、用滤波器对信号进行滤波，比较滤波前后语音信号的波形及频谱，回放语音信号。 对以上的图像分析我们发现，FIR滤波器的更加接近理想矩形，过渡带窄，通带基本无波纹，衰减基本上能符合要求，但是在阻带的波纹值较大。巴特沃兹滤波器的过渡带很宽，

- 40 -

实 习 报 告

很难满足要求。综合考虑，选择基于Kaiser窗函数的FIR低通滤波器、基于Kaiser窗函数的FIR高通滤波器、基于Kaiser窗函数的FIR带通滤波器为最佳。 （1）a.使用kaiser窗函数FIR低通滤波器对语音信号进行滤波 Matlab程序：

%FIR低通滤波比较

[y,fs,bits]=wavread('pch.wav'); sound(y,fs,bits); n=length(y); Y=fft(y,n);

fp=1000;fc=1200;As=100;Ap=1;%给定指标 wc=2*pi*fc/fs;%归一化频率 wp=2*pi*fp/fs; ws=wc-wp;

deltaw=(wc+wp)/(2*pi);%归一化过渡带宽计算 n0=ceil((As-8)/2.285/ws);%计算所需滤波器长度 beta=0.1102*(As-8.7);%计算kaiser窗的值 wdkai=(kaiser(n0+1,beta));%kaiser窗函数 b0=fir1(n0,deltaw,wdkai); figure(1);

[h,w]=freqz(b0,1);

plot(w*fs*0.5/pi,20*log10(abs(h))); axis([0,4000,-200,0]);grid on x=fftfilt(b0,y); X=fft(x,fs); figure(2);

subplot(2,2,1);plot(abs(Y)); axis([0,5000,0,200]); title('低通滤波前信号频谱'); subplot(2,2,2);plot(abs(X)); axis([0,5000,0,200]);

- 41 -

实 习 报 告

title('低通滤波后信号频谱'); subplot(2,2,3); plot(y);

title('低通滤波前信号波形'); subplot(2,2,4);plot(x); title('低通滤波前信号波形'); sound(x,fs,bits);

b.使用butterIIR低通滤波器对语音信号进行滤波

%--------------------采用butter低通滤波器-----------------% clear;clc;

[y,Fs]=audioread('pch.wav');

fp=1000;fs=1200;As=100;Ap=1;%相应的参数指标 ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;Bt=ws-wp;%归一化频率 [N,wc]=buttord(wp,ws,Ap,As,'s'); [B,A]=butter(N,wc,'s'); [Bz,Az]=bilinear(B,A,Fs);

- 42 -

实 习 报 告

figure(1);

[h,w]=freqz(Bz,Az);

plot(w*44100*0.5/pi,20*log10(abs(h))); grid on;

%语音信号的长度,FFT变换 n=length(y); yn=fft(y,n); x=fftfilt(Bz,y); xn=fft(x,fs); figure(2);

subplot(2,2,1);plot(abs(yn));title('原始信号频谱'); subplot(2,2,2);plot(abs(xn));title('滤波后信号频谱'); subplot(2,2,3);plot(y);title('原始信号波形'); subplot(2,2,4);plot(x);title('滤波后信号波形'); sound(x,Fs);

（2）使用kaiser窗函数FIR高通滤波器对语音信号进行滤波

- 43 -

实 习 报 告

Matlab程序：

%高通滤波：

[y,fs,bits]=wavread('pch.wav'); sound(y,fs,bits); n=length(y); Y=fft(y,n);

fc=4800;fp=5000;As=100;Ap=1;%给定指标 wc=2*pi*fc/fs;%归一化频率 wp=2*pi*fp/fs;

deltaw=(wc+wp)/(2*pi);%截止频率计算 ws=wp-wc;%归一化带宽计算

n0=ceil((As-8)/2.285/ws);%计算所需滤波器长度 beta=0.1102*(As-8.7);%计算kaiser窗的值 wdkai=(kaiser(n0,beta));%kaiser窗函数 b0=fir1(n0-1,deltaw,'high',wdkai); figure(1);

[h,w]=freqz(b0,1);

plot(w*fs*0.5/pi,20*log10(abs(h))); axis([0,8000,-200,0]);grid on x=fftfilt(b0,y); X=fft(x,fs); figure(2); subplot(2,2,1); plot(abs(Y)); axis([0,5000,0,200]); title('高通滤波前信号频谱'); subplot(2,2,2); plot(abs(X)); axis([0,8000,0,50]); title('高通滤波后信号频谱');

- 44 -

实 习 报 告

subplot(2,2,3);plot(y); title('高通滤波前信号波形'); subplot(2,2,4);plot(x); title('高通滤波前信号波形'); sound(x,fs,bits);

b.使用butterIIR高通滤波器对语音信号进行滤波

%--------------------采用butter低通滤波器-----------------% clear;clc;

[y,Fs]=audioread('pch.wav');

fp=4800;fs=5000;As=100;Ap=1;%相应的参数指标 ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;Bt=ws-wp;%归一化频率 [N,wc]=buttord(wp,ws,Ap,As); [B,A]=butter(N,wc,'high'); [Bz,Az]=bilinear(B,A,Fs);

- 45 -

实 习 报 告

figure(1);

[h,w]=freqz(Bz,Az);

plot(w*44100*0.5/pi,20*log10(abs(h))); grid on;

%语音信号的长度,FFT变换 n=length(y); yn=fft(y,n); x=fftfilt(Bz,y); xn=fft(x,fs); figure(2);

subplot(2,2,1);plot(abs(yn));title('原始信号频谱'); subplot(2,2,2);plot(abs(xn));title('滤波后信号频谱'); subplot(2,2,3);plot(y);title('原始信号波形'); subplot(2,2,4);plot(x);title('滤波后信号波形'); sound(x,Fs);

- 46 -

实 习 报 告

（3）a.使用kaiser窗函数FIR高通滤波器对语音信号进行滤波 Matlab程序：

%FIR带通滤波器滤波比较 fp1=1200;fp2=3000; fc1=1000;fc2=3000 As=100;Ap=1; fs=22050;

x1=wavread('pch.wav'); wc1=2*pi*fc1/fs;

wp1=2*pi*fp1/fs;wc2=2*pi*fc2/fs;wp2=2*pi*fp2/fs; deltaw1=(wc1+wp1)/2;

deltaw2=(wc2+wp2)/2;ws=wp1-wc1; N=ceil(8*pi/ws);

deltaw=[deltaw1,deltaw2]; %wdkai=(kaiser(n0+1,beta)); %b0=fir1(n0,deltaw,wdkai); [b,a]=fir1(N,deltaw,'bandpass'); figure(1); freqz(b,a,512); title('FIR带通滤波器'); xlabel('频率/Hz'); ylabel('幅值'); figure(2);

f1=fftfilt(a,b,x1); y1=fft(f1,1024); y2=fft(x1,1024); subplot(2,1,1); plot(abs(y2));

title('FIR带通滤波器滤波前的频谱');

- 47 -

实 习 报 告

xlabel('频率/Hz'); ylabel('幅值'); subplot(2,1,2); plot(abs(x1));

title('FIR带通滤波器滤波后的频谱'); xlabel('频率/Hz'); ylabel('幅值');

b.使用butterIIR带通滤波器对语音信号进行滤波

%--------------------采用butter带通滤波器-----------------% clear;clc;

[y,Fs]=audioread('pch.wav');

fp1=1200;fp2=3000;fs1=1000;fs2=3200;As=100;Ap=1;%给定指标 wp=[2*fp1*pi/Fs,2*fp2*pi/Fs];

- 48 -

实 习 报 告

ws=[2*fs1*pi/Fs,2*fs2*pi/Fs]; [N,wpo]=ellipord(wp,ws,Ap,As);%

[B,A]=ellip(N,Ap,As,wpo);%调用ellip函数计算带通滤波器系统函数系数向量B和A [Bz,Az]=bilinear(B,A,Fs); figure(1);

[h,w]=freqz(Bz,Az);

plot(w*44100*0.5/pi,20*log10(abs(h))); grid on;

%语音信号的长度,FFT变换 n=length(y); yn=fft(y,n); x=fftfilt(Bz,y); xn=fft(x,Fs); figure(2);

subplot(2,2,1);plot(abs(yn));title('原始信号频谱'); subplot(2,2,2);plot(abs(xn));title('滤波后信号频谱'); subplot(2,2,3);plot(y);title('原始信号波形'); subplot(2,2,4);plot(x);title('滤波后信号波形'); sound(x,Fs);

- 49 -

实 习 报 告

由滤波前后的对比我们可以发现，基于Kaiser窗函数的滤波器基本上能让通带内的信号通过，滤掉阻带内的信号。 回放语音信号：

回放语音信号之后可以发现，声音发生了变化，证明了滤波的效果。经过低通滤波器之后的声音变得低沉了很多，经过高通滤波器之后声音消失不见。经过带通滤波器之后的声音也发生变化，跟原声也比较接近。这是由声音本身的频谱特性和滤波器的特性决定的。 八、实验思考

1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗?从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？

答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的butter数字滤波器，从图中可以看到这种非线性关系。 2．能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？

答：IIR数字滤波器的设计实际上是求解滤波器的系数和，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼

- 50 -

实 习 报 告

近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。但是它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。

- 51 -

实 习 报 告

实验四 simulink简单应运

一、实验目的

1．学会 SIMULINK 的使用，掌握 SIMULINK 设计方法； 2．SIMULINK 设计 FIR 和 IIR 数字滤波器的方法； 3．学会用 SIMULINK 对信号进行分析和处理。 二、实验内容

1．初步 SIMULINK 的使用方法，识别常用的模块，如信号源模块 Sources 接收模块 Sinks，以及数学运算模块 Math Operations，接口模块 Ports & Operations。 任务 1 ：构建如下图所示的 Bus Signal 系统。

任务 2 ：参照下面系统，实现 1-100 求和运算

- 52 -

实 习 报 告

任务 3 ：1.利用相关模块实现 4 个 D 触发器实现 16 进制计数器 2.理解 PCM 编解码原理，构建如下系统实现 PCM 便解码仿真

3.理解数字滤波设计方法，相关参数自拟，实现数字滤波器设计

- 53 -

实 习 报 告

三、实验要求

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

2.完成实验任务及要求，简述参数设定、实现方法和过程 3.总结实验所得主要结论。 四、实验分析和应运

1 ：构建如下图所示的 Bus Signal 系统。

- 54 -

实 习 报 告

2 ：参照下面系统，实现 1-100 求和运算

结果显示：

3.（1）利用相关模块实现 4 个 D 触发器实现 16 进制计数器

- 55 -

实 习 报 告

（2）理解 PCM 编解码原理，构建如下系统实现 PCM 便解码仿真

- 56 -

实 习 报 告

仿真显示：

（3）理解数字滤波设计方法，相关参数自拟，实现数字滤波器设计

显示：

- 57 -

实 习 报 告

- 58 -

实 习 报 告

- 59 -

实 习 报 告

五、实验小结

通过Simulink简单应用，掌握了simulink的使用以及simulink的设计方法。了解了simulink的元器件库里各种模块的作用。通过用simulink设计FIR和IIR数字滤波器更好的理解和运用simulink。

- 60 -

实 习 报 告

- 61 -