实验一 信号的时域描述与运算
一、实验目的
①掌握信号的MATLAB表示及其可视化方法。 ②掌握信号基本时域运算的MATLAB实现方法。
③利用MATLAB分析常用信号，加深对信号时域特性的理解。
二、实验原理与方法
1. 连续时间信号的MATLAB表示
连续时间信号指的是在连续时间范围内有定义的信号，即除了若干个不连续点外，在任何时刻信号都有定义。在MATLAB中连续时间信号可以用两种方法来表示，即向量表示法和符号对象表示法。
从严格意义上来说，MATLAB并不能处理连续时间信号，在MATLAB中连续时间信号是用等时间间隔采样后的采样值来近似表示的，当采样间隔足够小时，这些采样值就可以很好地近似表示出连续时间信号，这种表示方法称为向量表示法。表示一个连续时间信号需要使用两个向量，其中一个向量用于表示信号的时间范围，另一个向量表示连续时间信号在该时间范围内的采样值。例如一个正弦信号可以表示如下：
>> t=0:0.01:10; >> x=sin(t);
利用plot(t,x)命令可以绘制上述信号的时域波形，如图1所示。
如果连续时间信号可以用表达式来描述，则还可以采用符号表达式來表示信号。例如对于上述正弦信号，可以用符号对象表示如下：
>> x=sin(t); >> ezplot(X);
利用ezplot(x)命令可以绘制上述信号的时域波形
10.80.60.40.20-0.2-0.4-0.6-0.8-10124567Time(seconds)图1 利用向量表示连续时间信号38910
常用的信号产生函数 函数名 heaviside sin cos sinc exp 功能 单位阶跃函数 正弦函数 余弦函数 sinc函数 指数函数 sin(t) 函数名 rectpuls tripuls square sawtooth 功能 门函数 三角脉冲函数 周期方波 周期锯齿波或三角波 10.50-0.5-1-6-4024t图 2 利用符号对象表示连续时间信号-26
2.连续时间信号的时域运算
对连续时间信号的运算包括两信号相加、相乘、微分、积分，以及位移、反转、尺度变换（尺度伸缩）等。
1）相加和相乘
信号相加和相乘指两信号对应时刻的值相加和相乘，对于两个采用向量表示的可以直接使用算术运算的运算符“+”和“*”来计算，此时要求表示两信号的向量时间范围和采样间隔相同。采用符号对象表示的两个信号，可以直接根据符号对象的运算规则运算。
2）微分和积分
对于向量表示法表示的连续时间信号，可以通过数值计算的方法计算信号的微分和积分。这里微分使用差分来近似求取的，由时间向量[表示的连续时间信号，其微分可以通过下式求得
t1,t2,,tN]和采样值向量[x1,x2,,xN]
x'(t)|ttkxk1xk,k1,2,,N1t
xk1xk。
其中t表示采样间隔。MATLAB中用diff函数来计算差分
连续时间信号的定积分可以由MATLAB的qud函数实现，调用格式为
quad ('function_name',a,b)
其中，function_name为被积函数名，a、b为积分区间。
对于符号对象表示的连续时间信号，MATLAB提供了diff函数和quad函数分别用于求微分和积分。
3.离散时间信号的MATLAB表示
离散时间信号仅在一些离散时刻有定义。在MATLAB中离散时间信号需要使用两个向量来表示，其中一个向量用于表示离散的时间点，另一个向量表示在这些时间点上的值。例如对于如下时间信号
x(n){3,2,1,2,1,1,2,3} 采用MATLAB可以表示如下： >> n=-3:4;
>> x=[-3 2 -1 2 1 -1 2 3]; >> stem(n,x,'filled'); >> xlabel('n');
>> title('x(n)');
Stem函数用于绘制离散时间信号波形，为了与我们表示离散时间信号的习惯相同，在绘图时一般需要添加‘filled’选项，以绘制实心的杆状图形。上述命令绘制的信号时域波形如图3所示。
x(n)3210-1-2-3-3-2-11n图3 离散时间信号示例0234
4.离散时间信号的时域运算
离散时间信号的相加相乘是将两个信号对应的时间点上的值相加或相乘，可以直接使用算术运算的运算符“+”和“*”来计算。
离散时间信号的位移，则可看作是将表示时间的向量平移，而表示对应时间点上的值的向量不变。
离散时间信号的反转，则可以看作是将表示时间的向量和表示对应时间点上的值的向量以零点为基准点，一纵轴为对称轴反折，向量的反折可以利用MATLAB的fliplr函数实现。
三、实验内容
（1）、利用MATLAB绘制下列连续时间信号的波形：
0.5tx(t)(1e)u(t) 1、
MATLAB代码如下： syms t;
x=(1-exp(-0.5*t))*heaviside(t); ezplot(t,x);
xlabel('t')； ylabel('x'); title('x(t)');
程序运行结果如下：
x(t)2.521.510.50-0.5-1-1.50123t456错误！未找到引用源。
xt)[u(t)u(t2)] 2、x(t)cos(MATLAB代码如下： syms t;
x=cos(pi*t)*[heaviside(t)-heaviside(t-2)]; ezplot(x,[-1,3]); xlabel('t')； ylabel('x'); title('x(t)');
程序运行结果如下：
x(t)10.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.500.51t1.522.53错误！未找到引用源。
xx(t)3、
|t|cos(2t)[u(t2)u(t2)]2
MATLAB代码如下： syms t;
x=abs(t)/2*cos(pi*t)*(heaviside(t+2)-heaviside(t-2)); ezplot(x,[-3,3]); xlabel('t')； ylabel('x'); title('x(t)');
程序运行结果如下：
x(t)0.80.60.40.20-0.2-0.4-0.6-3x-2-10t123
4、错误！未找到引用源。 MATLAB代码如下： syms t;
x=exp(-t)*sin(2*pi*t)*(heaviside(t)-heaviside(t-2)); ezplot(x,[-1,4]); axis([-1 4 -2 2]); xlabel('t')； ylabel('x'); title('x(t)');
程序运行结果如下：
x(t)21.510.50-0.5-1-1.5-2-1x-0.500.511.5t22.533.54
（2）、利用MATLAB绘制下列离散时间信号的波形：
1、x(n)u(n3) MATLAB代码如下： n=-4:7; x=ones(size(n)); stem(n,x,'filled'); xlabel('n'); ylabel('x'); title('x[n]');
程序运行结果如下：
nx(n)(1/2)u(n) 2、
MATLAB代码如下： n=-1:10;
x=(-1/2).^n.*heaviside(n); stem(n,x,'filled'); xlabel('n'); ylabel('x'); title('x[n]');
程序运行结果如下：
3、x(n)n[u(n)u(n5)] MATLAB代码如下： n=-1:10;
x=n.*[heaviside(n)-heaviside(n-5)]; stem(n,x,'filled'); xlabel('n'); ylabel('x'); title('x[n]');
程序运行结果如下：
4、x(n)sin(n/2)u(n) MATLAB代码如下： n=-1:10;
x=sin(n*pi/2).*heaviside(n); stem(n,x,'filled'); xlabel('n'); ylabel('x'); title('x[n]');
程序运行结果如下：
（3）、利用MATLAB生成并绘制连续周期矩形波信号，要求周期为2，峰值为3，显示三个周期波形。
MATLAB代码如下： t=0:0.01:6; x=3.*square(t*pi); plot(t,x); axis([0 6 -4 4]); xlabel('n'); ylabel('x'); title('x[n]');
程序运行结果如下：
x[n]43210-1-2-3-4x0123n456
(4)、已知信号x1(t)4t,(0t4)，及信号错误！未找到引用源。,用MATLAB绘出下列信号的波形
1、
x3(t)x1(t)x2(t)
MATLAB代码如下： syms t;
x1=(4-t)*(heaviside(t)-heaviside(t-4)); x2=sin(2*pi*t); x3=x1+x2; ezplot(x3,[-1,5]); xlabel('t'); ylabel('x3'); title('x(t)');
程序运行结果如下：
x(t)543x3210-1-1012t345
2、x4(t)x1(t)x2(t) MATLAB代码如下： syms t;
x1=(4-t)*(heaviside(t)-heaviside(t-4)); x2=sin(2*pi*t); x4=x1*x2; ezplot(x4,[-1,5]); xlabel('t'); ylabel('x4'); title('x(t)');
程序运行结果如下：
x(t)4321x40-1-2-3-1012t345错误！未找到引用源。
3、
x5(t)x1(t)x1(t)
MATLAB代码如下： syms t;
x1=(4-t)*(heaviside(t)-heaviside(t-4)); x1=(4+t)*(heaviside(-t)-heaviside(t+4)); x5=x1+x2; ezplot(t,x5);
程序运行结果如下：
错误！未找到引用源。
4、
x6(t)x2(t)x3(t1)
MATLAB代码如下： syms t;
x1=(4-t)*(heaviside(t)-heaviside(t-4)); x2=sin(2*pi*t); x3=x1+x2; x4=subs(x3,t,t-1); x6=x2*x4; ezplot(t,x6);
程序运行结果如下：
x = t, y = sin(2 t) (sin(2 (t - 1)) - (heaviside(t - 1) - heaviside(t - 5)) (t - 5))432y10-1-2-10123x4567
（5）、已知离散时间信号x(n)波形（
x(n)[0,1,2,3,3,3,3]），用MATLAB
绘出x(n)、x(-n)、x(n+2)和x(n-2)的波形。
MATLAB代码如下： n=-3:4;
x=[0,1,2,3,3,3,3,0]; n1=-fliplr(n); x1=fliplr(x); n2=n+2; n3=n-2; subplot(221); stem(n,x,'filled'); subplot(222); stem(n1,x1,'filled'); subplot(223); stem(n2,x,'filled');
subplot(224); stem(n3,x,'filled'); 程序运行结果如下：
3322110-43-20240-43-202422110-202460-6-4-202
（6）、有MATLAB编程绘制下列信号的时域波形，观察信号是否为周期信号？若是周期信号，周期是多少？若不是周期信号，请说明原因。
x(t)1cos(t)2cos(t)cos(2t)43241、
MATLAB代码如下： syms t;
x=1+cos(pi/4*t-pi/3)+2*cos(pi/2*t-pi/4)+cos(2*pi*t); ezplot(x,[-20,20]); 程序运行结果：
cos(2 t) + 2 cos(( t)/2 - /4) +...+ 1 data143210-1-2-3 -20-15-10-50t5101520
信号分析：该信号为周期信号，周期为8
2、x(t)sin(t)2sin(t) MATLAB代码如下： syms t;
x=sin(t)+2*sin(pi*t); ezplot(x,[-20,20]); 程序运行结果如下：
2 sin( t) + sin(t) 3data2210-1-2-3 -20-15-10-50t5101520
信号分析：该信号为非周期信号错误！未找到引用源。
x(n)23sin(3、
MATLAB代码如下： n=1:20;
2n)38
x=2+3*sin(2*n/3*pi-pi/8); stem(n,x,'filled'); 程序运行结果如下：
5data34.543.532.521.510.50 024681012141618 20
信号分析：该信号为周期信号，周期为3错误！未找到引用源。
nnnx(n)cos()sin()cos()632 4、
MATLAB代码如下： n=1:20;
x=cos(n*pi/6)+sin(n*pi/3)+cos(n*pi/2); stem(n,x,'filled'); 程序运行结果如下：
信号分析：该信号为周期信号，周期为12
四、体会与建议
本次实验主要是针对MATLAB在信号与系统中应用的一些基本操作是熟悉，对于一种新的编程语言的初步接触，第一次实验感觉有些生涩，才深知课前预习的重要性，结合以往学过的知识，感受到MATLAB在信号分析中不可或缺的作用。
2data41.510.50-0.5-1-1.5-2-2.5 0 2468101214161820
实验二 LTI系统的时域分析
一、实验目的
1、掌握利用MATLAB对系统进行时域分析的方法。
2、掌握连续时间系统零状态响应、冲激响应和阶跃响应的求解方法。 3、掌握求解离散时间系统响应、单位抽样响应的方法。
4、加深对卷积积分和卷积和的理解。掌握利用计算机进行卷积积分和卷积和计算的方法。
二、实验原理与方法
1、连续时间系统时域分析的MATLAB实现
（1)、连续时间系统的MATLAB表示
LTI连续系统通常可以由系统微分方程描述，设LTI因果系统的微分方程一般式为：
则在MATLAB里，可以建立系统模型如下： b=[bM,bM-......b0]; a=[aN,aN-1......a0]; sys=tf(b,a);
其中，tf是用于创建系统模型的函数，向量a与b的元素是以微分方程求导的降幂次序来排列的，如果有缺项，应用0补齐。
(2)、连续时间系统的零状态响应
零状态响应是指系统的初始状态为零，仅由输入信号引起的响应。MATLAB提供了一个用于求解零状态响应的函数lism，其调用格式如下：
lsim(sys,x,t)绘出输入信号及响应的波形，x和t表示输入信号数值向量和时间向量。 y= lsim(sys,x,t)这种调用格式不绘出波形，而是返回响应的数值向量。
（3）、连续时间系统的冲激响应与阶跃响应。
MATLAB提供了函数impluse来求指定时间范围内，由模型sys描述的连续时间系统的单位冲激响应。其调用格式如下：
impulse(sys)在默认的时间范围内绘出系统冲激响应的时域波形。 impulse(sys，T)绘出系统在0-T范围内冲激响应的时域波形。
impulse(sys，ts：tp：te)绘出系统在ts-tp范围内，以tp为时间间隔取样的冲激响应的时域波形。
[y,t]= impulse(…)该调用格式不绘出单位冲激响应波形，而是返回单位冲激响应的数值向量及其对应的时间向量。
函数step用于求解单位阶跃响应，函数step同样也有如下几种调用格式： step(sys) step（sys，T）
step（sys，ts：tp：te） [y，t]=step（…）
2、离散时间系统时域分析的MATLAB实现
（1）、离散时间系统的MATLAB表示。
LTI离散系统通常可以由系统差分方程描述，设描述系统的差分方程为：
ak0Nky(nk)bx(nr)r0rM
则在MATLAB里，可以用如下两个向量来表示这个系统： b=[b0,b1……,bM] a=[a0,a1……,aN]
(2）、离散时间系统对任意输入的响应。
MATLAB提供求LTI离散系统响应的专用函数fliter，该函数用于求取差分方程描述的离散时间系统在指定时间范围内对输入序列所产生的响应，其调用格式如下：
y=filter（b,a,x）
其中，x为输入序列，y为输出序列，x，y所对应的时间区间必须相同。
（3）、离散时间系统的单位抽样响应。
MATLAB提供了函数impz来求指定时间范围内，由向量b和a描述的离散时间系统的单 位抽样响应。其调用格式如下：
impz(b,a)在默认的时间范围内绘出系统抽样响应的时域波形。。
impz(b,a，T)绘出系统在0-N范围内冲激响应的时域波形。 impz(b,a，ns:ne)绘出系统在ns-ne范围的冲激响应的时域波形。
[y,n]= impz(…)该调用格式不绘出单位冲激响应波形，而是返回单位冲激响应的数值向量及其对应的时间向量。
3、卷积和与卷积积分
（1）、离散时间序列的卷积和
卷积和是离散系统时域分析的基本方法之一，离散时间序列x1(n)和x2(n)的卷积和
x(n)定义如下：
x(n)x1(n)x2(n)x1(k)x2(nk)
MATLAB提供了函数conv来求两个离散序列的卷积和。其调用格式如下：
x=conv(x1,x2)
（2）、连续时间信号的卷积积分
卷积积分是连续系统时域分析的有效方法和工具之一，连续时间信号和的卷积积分定义如下:
用户可根据书上内容自定义一个用于计算卷积积分的通用函数sconv。
三、实验内容
1、已知描述模拟低通、高通、带通和带阻滤波器的微分方程如下，试采用MATLAB绘出各系统的单位冲激响应和单位阶跃响应波形。
MATLAB代码如下： b=[0 0 1]; a=[1 2.^0.5 1]; sys=tf(b,a);
subplot(121); impulse(sys); subplot(122); step(sys);
程序运行结果如下：
Impulse Response0.60.50.40.3udetiplmA0.20.10-0.10510Time (seconds)
MATLAB代码如下： b=[1 0 0]; a=[1 2.^0.5 1]; sys=tf(b,a); subplot(121); impulse(sys); subplot(122); step(sys);
程序运行结果如下：
Step Response1.41.210.8udetiplmA0.60.40.2002468Time (seconds)
Impulse Response0.20-0.2-0.4eud-0.6tiplmA-0.8-1-1.2-1.4-1.602468Time (seconds)
MATLAB代码如下： b=[0 1 0]; a=[1 1 1]; sys=tf(b,a); subplot(121); impulse(sys); subplot(122); step(sys);
程序运行结果如下：
Step Response10.80.6e0.4udtiplmA0.20-0.2-0.40246Time (seconds)
Impulse Response10.80.6de0.4utiplmA0.20-0.2-0.40510Time (seconds)
MATLAB代码如下： b=[1 0 1]; a=[1 1 1]; subplot(121); impulse(sys); subplot(122); step(sys);
程序运行结果如下：
Step Response0.60.50.40.3udetiplmA0.20.10-0.10510Time (seconds)
Impulse Response0.40.20-0.2-0.4-0.6-0.8-11.110.90.80.70.60.50.4Step ResponseAmplitude05Time (seconds)10Amplitude0510Time (seconds)
2、已知某系统可以由微分方程描述
①请利用MATLAB绘出该系统冲激响应和阶跃响应的时域波形； MATLAB代码如下： b=[0 0 1]; a=[1 1 6]; sys=tf(b,a); subplot(121); impulse(sys); subplot(122); step(sys); 程序运行结果：
Impulse Response0.50.40.30.20.10-0.1-0.20.350.30.250.20.150.10.050Step ResponseAmplitude0510Amplitude0510Time (seconds)Time (seconds)
②根据冲激响应的时域波形分析系统的稳定性； 信号分析：单位冲击响应绝对可积，系统稳定。 ③如果系统的输入为，求系统的零状态响应。 MATLAB代码如下： b=[0 0 1]; a=[1 1 6]; sys=tf(b,a);
t=[0:0.01:10];x=exp(-t); lsim(sys,x,t); 程序运行结果如下：
Linear Simulation Results1.210.80.60.40.20-0.2Amplitude012345678910Time (seconds)
3、已知描述离散系统的差分方程如下，试采用MATLAB绘出各系统的单位抽样响应，并根据单位抽样响应的时域波形分析系统的稳定性。
①y(n)+3y(n-1)+2y(n-2)=x(n) MATLAB代码如下： b=[1 0 0]; a=[1 3 2]; impz(b,a,10); 程序运行结果：
Impulse Response6004002000 Amplitude-200-400-600-800-1000-1200012345n (samples)6789
信号分析：单位抽样响应不绝度可积，系统不稳定。
②y(n)-0.5y(n-1)+0.8y(n-2)=x(n)-3x(n-1) MATLAB代码如下： b=[1 -3 0]; a=[1 -0.5 0.8]; impz(b,a,10); 程序运行结果如下：
Impulse Response2.521.51 Amplitude0.50-0.5-1-1.5-2-2.5012345n (samples)6789
信号分析：单位抽样响应绝对可积分，系统稳定。
4、已知系统可以由如下差分方程描述
y(n)+y(n-1)+0.25y(n-2)=x(n)
试采用MATLAB绘出该系统的单位抽样响应波形和阶跃响应波形。
MATLAB代码如下： b=[1 0 0]; a=[1 1 0.25]; subplot(121); impz(b,a,10); subplot(122); stepz(b,a,10); 程序运行结果如下：
Impulse ResponseStep Response110.80.90.60.80.40.7e0.20.6ududetipl0tipl0.5mmAA -0.2 0.4-0.40.3-0.60.2-0.80.1-102468002468n (samples)n (samples)5、采用MATLAB计算如下两个序列的卷积，并绘出图形x1(n)[1,2,1,1] x2(n)1,(2n2)
MATLAB代码如下： x1=[0 0 1 2 1 1 0]; x2=[0 1 1 1 1 1 0]; x x =
0 0 0 1 3 4 5 5 4 2 1 0
n=-3:4;
x=[1 3 4 5 5 4 2 1]; stem(n,x,'filled');
0
程序运行结果如下：
54.543.532.521.510.50-3-2-101234
6、已知某LTI离散系统，其单位抽样响应h(n)=sin(0.5n)，n>=0，系统的输入为x(n)=sin(0.2n)，n>=0，计算当n=0,1,2，…，40时系统的零状态响应y(n)，绘出x(n)，h(n)和y(n)时域波形。
MATLAB代码如下： n=0:40; h=sin(0.5*n); x=sin(0.2*n); y=conv(h,x); subplot(311); stem(n,h,'filled'); subplot(312); stem(n,x,'filled'); subplot(313); stem(y,'filled');
程序运行结果如下：
10-110-150-5051015202530354005101520253035400102030405060708090
7、已知两个连续信号，是采用MATLAB求这两个信号的卷积
MATLAB代码如下： Sconv.m
function [x,t]=sconv(x1,x2,t1,t2,dt) x=conv(x1,x2); x=x*dt;
t0=t1(1)+t2(1);
l=length(x1)+length(x2)-2; t=t0:dt:(t0+l*dt); end dt=0.001; t1=(-2):dt:2; t2=(-1):dt:1;
x1=2.*(heaviside(t1+2)-heaviside(t1-2)); x2=heaviside(t2+1)-heaviside(t2-1); subplot(221); plot(t1,x1); xlabel('t(s)'); title('x_1(t)'); subplot(222); plot(t2,x2); xlabel('t(s)'); title('x_2(t)');
[x,t]=sconv(x1,x2,t1,t2,dt); subplot(212); plot(t,x); xlabel('t(s)');
title('x(t)=x_1(t)*x_2(t) 程序运行结果如下：
x1(t)210.81.50.6x2(t)1-2-10t(s)120.4-1-0.5x(t)=x1(t)*x2(t) t=0.0010t(s)0.5143210-3-2-10t(s)123
四、体会与建议
虽然在反复操作中逐渐找到了些许诀窍，但感觉计算机绘图还是有些地方与自己的想法相左，不能用正确的细致的编程语言完成想要的波形，总感觉要绕一些弯子才能够正确分析，此外，对于一些绘图的细节依然关注程度不够，导致无法全面的精确的反应出信号的某些特性。在以后的实验中还需要加强锻炼。
实验三 信号频域分析
一、实验目的
①深入理解信号频谱的概念，掌握信号的频域分析方法。 ②观察典型周期信号和非周期信号的频谱，掌握其频谱特性。
二、实验原理与方法
1、连续周期信号的频谱分析
如果周期信号满足狄里赫利条件，就可以展开为傅里叶级数形式，即
x(t)ckejkw0tk （1）
1ckT0x(t)ejkw0tdtT0 （2）
式中，分。
式（1）和式（2）定义为周期信号复指数形式的傅里叶级数，系数
T0表示基波周期，w02/T0为基波频率，T0()表示任一个基波周期内的积
ck称为x(t)的傅里
叶系数。周期信号的傅里叶级数还可以由三角函数的线性组合来表示，即
x(t)a0akcoskw0tbksinkw0tk1k1 （3）
其中：
a0122T0x(t)dt,akTx(t)coskw0tdt,bkTx(t)sinkw0tdtT0T00T00 （4）
式（3）中同频率的正弦项和余弦项可以合并，从而得到三角函数形式的傅里叶级数，即
x(t)A0Akcos(kw0tk)k1 （5）
其中：
2A0a0,Akakbk2,karctanbkak （6）
可见，任何满足狄里赫利条件的周期信号都可以表示成一组谐波关系的复指数函数或三角函数的叠加。一般来说周期信号表示为傅里叶级数时需要无限多项才能完全逼近原信号，但在实际应用中经常采用有限项级数来替代，所选项数越多就越逼近原信号。
2、连续非周期信号的频谱分析
对于非周期连续时间信号，吸纳后的傅里叶变换和傅里叶逆变换定义为
X(w)x(t)12x(t)ejwtdt （7）
X(w)ejwtdw （8）
式（7）和式（8）把信号的时域特性和频域特性联系起来，确立了非周期信号x(t)和频谱X(w)之间的关系。
采用MATLAB可以方便地求取非周期连续时间信号的傅里叶变换，这里我们介绍常用的集中方法。
（1）、符号运算法
MATLAB的符号数学工具箱提供了直接求解傅里叶变换和反变换的函数，fourier函数和ifourier函数，基本调用格式为
X=fourier(x) X=ifourier(X)
默认的时域变量为t，频域变量为w。
（2）、数值积分法
除了采用符号运算的方法外，我们还可以利用MATLAB的quad函数，采用数值积分的方法来进行连续信号的频谱分析，quad函数是一个用来计算数值积分的函数。利用quad函数可以计算非周期连续时间信号的频谱。Quad函数的一般调用格式为：
y=quad(fun,a,b)
y=quad(fun,a,b,TOL,TRACE,p1,p2,…)
其中fun指定被积函数，可以采用inline命令来创建，也可以通过传递函数句柄的形式来指定，a、b表示定积分的下限和上限，TOL表示允许的相对或绝对积分误差，TRACE表
示以被积函数的点绘图形式来跟踪该函数的返回值，如果TOL和TRACE为空矩阵，则使用缺省值，“p1,p2，…”表示被积函数出时间t之外所需的其他额外输入参数。
（3）、数值近似法
我们还可以利用MATLAB的数值计算的方法近似计算连续时间傅里叶变换。傅里叶变换X(w)可以由式（9）近似计算
X(w)x(t)ejwtdtlimx(k)ejwk0 （9）
当x(t)为时限信号，且足够小，式（9）可以演变为
X(w)x(k)ejkwkab （10）
而式（10）中求和部分又可以表示成一个行向量和一个列向量的乘积
x(k)ekabjkwejawj(a1)we[x(a),x((a1)),,x(b)].....jbwe （11）
式（11）可以很方便地利用MATLAB实现。
3、离散周期时间信号的频域分析
基波周期为N的周期序列x(n)可以用N个成谐波关系的复指数序列的加权和表示，即
x(n)kNckejk(2/N)n （12）
这里k= ck则称为离散 ck可以由式（13）确定。 ck1NkNx(n)ejk(2/N)n （13） 傅里叶系数 ck也称为x(n)的频谱系数，而且可以证明 ck是以N为周期的离散频率序 列。这说明了周期的离散时间函数对应于频域为周期的离散频率。 这里，我们用周期N与傅里叶系数 ck的乘积来表示周期离散时间信号的频谱，即 （14） X(k)NckkNx(n)ejk(2/N)nX(k)可以利用MATLAB提供的函数fft用来计算，调用格式为 Xfft(x) 该函数返回X(k)一个周期内的值，其中x表示x(n)一个周期内的样本值。 4、离散非周期时间信号的频域分析 非周期序列x(n)可以表示成一组复指数序列的连续和 x(n)其中 122X(ej)ejnd （15） X(e)jnx(n)ejn （16） 式（16）称为x(n)的离散时间傅里叶变换，式（15）和式（16）确立了非周期离散时 jjX(e)X(e)是连续频率的函数，x(n)间信号及其离散时间傅里叶变换之间的关系。jX(e)是周期的连续频率函数，其周期为2。可见，非周期离散时间称为频谱函数，且 函数对应于频域中是一个连续的周期的频率函数。 对于有限长的时间序列，式（16）可以表示为 X(e)x(n)ejnjnn1nNejn1jn2e[x(n1),x(n2),,x(nN)].....jnNe (17) 式（17）可以方便地利用MATLAB实现。 三、实验内容 1、已知x(t)是周期矩形脉冲信号。 ①计算该信号的傅里叶级数； ②利用MATLAB绘出由前N次谐波合成的信号波形，观察随着N的变化合成信号波形的变化规律； MATLAB代码如下： t=-1.5:0.01:1.5; N=input('N='); T=input('T='); a=input('a='); A=input('A='); x=zeros(size(t)); x=A*a/T; for n=1:1:N x=x+2*A/(n*pi)*sin(n*pi*a/T)*cos(n*2*pi/T*t); end plot(t,x); 程序运行结果如下： N=7;T=1;a=0.5;A=1; 1.210.80.60.40.20-0.2-1.5-1-0.500.511.5 N=20;T=1;a=0.5;A=1; 1.210.80.60.40.20-0.2-1.5-1-0.500.511.5 ③利用MATLAB绘出周期矩形脉冲信号的频谱，观察参数T和变化时对频谱波形的影响。 MATLAB代码如下： N=10; T=input('T='); a=input('a='); A=input('A='); n1=-N:-1; c1=A./(n1.*pi).*sin(pi.*a.*n1./T); c0=A.*a./T; n2=1:N; c2=A./(n2.*pi).*sin(pi.*a.*n2./T); cn=[c1 c0 c2]; n=-N:N; subplot(211); stem(n,abs(cn),'filled'); subplot(212); stem(n,angle(cn),'filled'); 程序运行结果如下： a=0.5;A=1;T=1; 0.50.40.30.20.10-1043210-10-8-6-4-20246810-8-6-4-20246810 Q1-1、什么是吉伯斯现象？产生吉伯斯现象的原因是什么？ Answer：将具有不连续点的周期函数（如矩形脉冲）进行傅立叶级数展开后，选取有限项进行合成。当选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时，该峰起值趋于一个常数，大约等于总跳变值的9%。这种现象称为吉伯斯现象。原因是在不连续点附近所有的正弦信号均具有相同的变化趋势，该趋势在有限项内无法被消除。 Q1-2、以周期矩形脉冲信号为例，说明周期信号的频谱有什么特点。 Answer：周期信号的频谱是具有周期性的一系列的脉冲信号。 Q1-3、周期矩形脉冲信号的有效频带宽度与信号的时域宽度之间有什么关系？ Answer：时域宽度越大，有效频带宽度越小。 Q1-4、随着矩形脉冲信号参数的变化，其频谱结构如何变化？ Answer：频谱包络形状不变，过零点不变，普贤间隔随着T变大而缩小。 2、已知x(t)是如图所示的矩形脉冲信号。 ①求该信号的傅里叶变换； ②利用MATLAB绘出矩形脉冲信号的频谱，观察矩形脉冲信号宽度变化时对频谱波形的影响； ③让矩形脉冲信号的面积始终等于1，改变矩形脉冲宽度，观察矩形脉冲信号时域波形和频谱随矩形脉冲宽度的变化趋势。 MATLAB代码如下： syms t w; A=input('A='); a=input('a='); X=int(A*exp(-j*w*t),t,-a/2,a/2); subplot(211); axis([-6*pi 6*pi 0 6]); ezplot(abs(X)); xlabel('\\omega'); ylabel('Magtitude'); title('|X(\\omega)|'); subplot(212); ezplot(angle(X),[-6*pi,6*pi]); xlabel('\\omega'); ylabel('Angle'); title(' |X()|1Magtitude0.50-15-10-5051015 |X()|2Magtitude1.510.50-6-4-20246 |X()|1Magtitude0.50-6-4-20246 Answer：矩形脉冲信号频谱为连续函数，周期矩形脉冲信号的频谱为一系列的脉冲。 Q2-2、根据矩形脉冲宽度变化时频谱的变化规律，说明信号的有效频带宽度与其时域宽度之间有什么关系？ Answer：信号有效频带宽度越大，时域宽度越小。 3、已知x[n]是如图所示的周期方波序列。 利用MATLAB绘制周期方波序列的频谱波形，改变参数N和N1的大小，观察频谱波形的变化趋势。 MATLAB代码如下： N1=input('N1='); N=input('N='); n1=0:N1; x1=ones(size(n1)); n2=N1+1:N-N1-1; x2=zeros(size(n2)); n3=N-N1:N; x3=ones(size(n3)); n=0:N; x=[x1 x2 x3]; X=fft(x); subplot(211); stem(n,x,'filled'); xlabel('n'); title('x(n)'); subplot(212); stem(n,X,'filled'); xlabel('k'); title('X(k)'); 程序运行结果如下： N1=2;N=9; x(n)10.5001234nX(k)567896420-201234k56789Q3-1、以周期方波序列为例，说明周期序列与连续周期信号的频谱有何异同。 Answer：周期序列的频谱向外越来越大，连续周期信号频谱则是中间向两边越来越小。 Q3-2、随着周期方波序列占空比的变化，其频谱如何随之变化？ Answer：随着占空比越来越大，频谱密度也越来越大。 4、已知一矩形脉冲序列。利用MATLAB绘制周期方波序列的频谱波形，改变矩形脉冲序列的宽度，观察频谱波形的变化趋势。 MATLAB代码如下： N1=input('N1='); w=-pi:0.01:pi; n=-N1:N1; x=ones(size(n)); X=x*exp(-j*n'*w); subplot(211); stem(n,x,'filled') xlabel('n'); title('x(n)'); subplot(212); plot(w/pi,abs(X)); xlabel('\\Omega/\\pi'); title('|X(e^j^Omega)|'); 程序运行结果如下： N1=4 x(n)10.50-4-3-2-10n|X(ejOmega)|12341050-1-0.8-0.6-0.4-0.20/0.20.40.60.81 N1=5 x(n)10.50-5-4-3-2-10n|X(ejOmega)|12345151050-1-0.8-0.6-0.4-0.20/0.20.40.60.81 Q4-1、随着矩形脉冲序列宽度的变化，其频谱如何随之变化？其宽度与频谱的有效宽度有何关系？ Answer：宽度越大，频谱变化就越密集。矩形脉冲信号序列宽度越大，频谱有效宽度越大。 四、体会与建议 对与信号的频域分析，在信号与系统课程中我的学习比较细致，各种类型的题目也比较容易找到思路，在编程过程中也能够正确分析。 实验四 LTI系统的频域分析 一、实验目的 ①加深对LTI系统频域响应基本概念的掌握和理解 ②学习和掌握LTI系统频域特性的分析方法 二、实验原理和方法 1、连续时间系统的频域响应 系统的频域响应定义为系统单位冲激响应的傅里叶变换，即 HY/X H()jh()ed 若LTI连续时间系统的冲激响应为h(t)，输入为x(t)，根据时域分析可知 Ytxtht对本式求傅里叶变换得 YXH 所以，频率响应还可以由零状态相应和输入的傅里叶变换之比得到 HY/XH 反映的是系统的固有属性，与外部激励无关，又可以表示为 H()H()ej()其中 H()称为系统的幅度响应， ()被称为系统的相位响应。 对于由下述微分方程描述的LTI连续时间系统 aynn0N(n)(t)bmx(m)(t)m0M 其频率响应H（jw）可以表示为（8-34）所示的jw的有理多项式。 Y()bM(j)MbM1(j)M1bM2(j)M2...b0H()X()aN(j)NaN1(j)N1aN2(j)N2...a0 MATLAB的信号处理工具箱提供了专门的函数freqs，用来分析连续时间系统的频率响应，该函数有下列几种调用格式： [h，w]=freqs(b,a)计算默认频率范围内200个频率点上的频率响应的取样值，这200个频率点记录在w中。 h=freqs(b,a,w) b、a分别为表示H(jw)的有理多项式中分子和分母多项式的系数向量，w为频率取样点，返回值h就是频率响应在频率取样点上的数值向量。 [h,w]=freqs(b,a,n)计算默认频率范围内n个频率点上的频率响应的取样值，这n个频率点记录在w中。 Freqs(b,a…) 这种调用格式不返回频率响应的取样值，而是以对数坐标的方式绘出系统的幅频响应飞相频响应。 2、 离散时间系统的频率响应 LTI离散时间系统的频率响应定义为单位抽样响应h(n)的离散时间傅里叶变换。 H(e)jnh(n)ejn 对于任意的输入信号x(n)，输入与输出信号的离散时间傅里叶变换有如下关系 Y(ej)H(ej)X(ej) 因此，系统的频率响应还可以表示为 H(ej)Y(ej)/X(ej) 当系统输入信号为x(n)时，系统的输出 y(n)ejnh(n)kej(nk)h(k)ejnH(ej) 由式（8-38）可知，虚指数信号通过LTI离散时间系统后信号的频率不变，信号的幅度 jH(e))表示了系统对不同频率信号的衰减量。 由系统频率响应的幅度值确定，所以 jH(e)是复值函数，可用幅度和相位表示。 一般情况下离散系统的频率响应 H(ej)H(ej)ej()若LTI离散系统可以由如下差分方程描述。 ay(ni)bx(nj)iji0j0NM jH(e)可以表示为的有理多项式。 则由式（8-37）描述的离散时间系统的频率响应 Y(ej)b0b1ej.....bMejMH(e)jX(e)a0a1ej.....aNejN jMatlab中有专门的函数freqz，与freqs函数用法类似，不再赘述。 三、实验内容 1、已知一个RLC电路构造的二阶高通滤波器。 （1）、计算该电路系统的频率响应及高通截止频率； （2）、利用MATLAB绘制幅度响应和相位响应曲线，比较系统的频率特性与理论计算的结果是否一致。 MATLAB代码如下： b=[1 0 0]; a=[1 10 50]; [H,w]=freqs(b,a); subplot(211); plot(w,abs(H)); set(gca,'xtick',[0,10]); set(gca,'ytick',[0 0.4 0.707 1]); xlabel('\\omega(rad/s)'); ylabel('Magnitude'); title('|H(j\\omega)|'); grid on; subplot(212); plot(w,angle(H)); set(gca,'xtick',[0,10]); xlabel('\\omega(rad/s)'); ylabel('Phase'); title('\\phi(j\\omega)'); grid on; 程序运行结果如下： |H(j)|1Magnitude0.7070.40010(rad/s)(j)43Phase210010(rad/s) 2、已知一个RC电路 （1）、对不同的RC值，用MATLAB绘出系统的幅度响应曲线，观察实验结果，分析RC电路具有什么样的频率特性？系统的频率特性随着RC值的改变，有何变化规律？ MATLAB代码如下： C=input('C='); R=input('R=') A=C*R; b=[1]; a=[A 1]; [H,w]=freqs(b,a); subplot(211); plot(w,abs(H)); set(gca,'ytick',[0 0.4 0.707 1]); xlabel('\\omega(rad/s)'); ylabel('Magnitude'); title('|H(j\\omega)|'); grid on; subplot(212); plot(w,angle(H)); xlabel('\\omega(rad/s)'); ylabel('Phase'); title('\\phi(j\\omega)'); grid on; 程序运行结果如下： C=0.05F R=10欧姆 |H(j)|1Magnitude0.7070.40012345(rad/s)(j)6789100-0.5-1-1.5Phase012345(rad/s)678910 C=1F R=5欧姆 |H(j)|1Magnitude0.7070.4000.10.20.30.40.50.6(rad/s)(j)0.70.80.910-0.5-1-1.5Phase00.10.20.30.40.50.6(rad/s)0.70.80.91 （2）、系统的输入信号x（t）=cos（100t）+cos（3000t），t=0-0.2s，该信号包含了一个低频分量和一个高频分量。试确定适当的RC值，滤掉信号中的高频分量。并绘出滤波前后的时域信号波形及系统的频率响应曲线。 MATLAB代码如下： b=[1]; a=[0.001 1]; sys=tf(b,a); t=0:0.001:0.2;x=cos(100*t)+cos(3000*t); lsim(sys,x,t); 程序运行结果如下： Linear Simulation Results21.510.5Amplitude0-0.5-1-1.5-200.020.040.060.080.10.120.140.160.180.2Time (seconds) 3、已知离散系统的系统框图 （1）、写出M=8时系统的差分方程及系统函数； （2）、利用MATLAB计算系统的单位抽样响应； MATLAB代码如下： a=[1]; b=[1 1 1 1 1 1 1 1 1]; impz(b,a,0:20); 程序运行结果如下： Impulse Response10.90.80.7 Amplitude0.60.50.40.30.20.10024681012n (samples)14161820 （3）、试利用MATLAB绘出其系统的零极点分布图、幅频和相频特性曲线，并分析该系统具有怎样的频率特性。 MATLAB代码如下： b=[1 1 1 1 1 1 1 1 1]; a=[1]; zplane(b,a) 程序运行结果如下： 10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.518 MATLAB代码如下： b=[1 1 1 1 1 1 1 1 1]; a=[1]; [H,w]=freqz(b,a); subplot(211); plot(w/pi,abs(H)); xlabel('\\omega(\\pi)'); ylabel('Magnitude'); title('|H(e^j^\\omega)|'); grid on; subplot(212); plot(w/pi,angle(H)/pi); set(gca,'xtick',[0,10]); xlabel('\\omega(\\pi)'); ylabel('Phase(\\pi)'); title('\heta(\\omega)'); grid on; 程序运行结果如下： |H(ej)|10Magnitude5000.10.20.30.40.5()()0.60.70.80.910.50-0.5-1Phase()0() 4、已知一离散时间LTI系统的频率响应，输入信号为x（n）=cos（0.3∏n）+0.5cos(0.8∏n)，计算系统对于x（n）的响应y（n）。 MATLAB代码如下： n=-20:20; x=cos(0.3*pi*n)+0.5*cos(0.8*pi*n); subplot(211); stem(n,x,'filled'); y=2*cos(0.3*pi*n); subplot(212); stem(n,y,'filled'); 程序运行结果如下： 210-1-2-20210-1-2-20-15-10-505101520-15-10-505101520 信号分析：为低通网络。 四、体会与建议 在编程过程中出了一些问题，题目表述不甚清晰，造成了一些误解和歧义。对信号的周期性与否说明太过于模糊，而且编程语言有一定的局限性，无法对理想信号进行分析，自身也需要加强对实际信号的理解还有误差分析的能力。 实验五 连续时间系统的复频域分析 一、实验目的 ①掌握拉普拉斯变换及其反变换的定义，并掌握MATLAB实现方法。 ②学习和掌握连续时间系统系统函数的定义及复频域分析方法。 ③掌握系统零极点的定义，加深理解系统零极点分布与系统特性的关系。 二、实验原理与方法 1、拉普拉斯变换 连续时间信号x(t)的拉普拉斯变换定义为 X(s)x(t)estdt （1） 拉普拉斯反变换定义为 x(t)2j1jjX(s)estds （2） 在MATLAB中，可以采用符号数学工具箱的laplace函数和ilaplace函数进行拉氏变换和反拉氏变换。 L=laplace(F)符号表达式F的拉氏变换，F中时间变量为t，返回变量为s的结果表达式。 L=laplace(F,t)用t替换结果中的变量s。 F=ilaplace(L)以s为变量的符号表达式L的拉氏反变换，返回时间变量为t的结果表达式。 F=ilaplace(L,x)用x替换结果中的变量t。 除了上述ilaplace 函数，还可以采用部分分式法，求解拉普拉斯逆变换，具体原理如下： 当 X (s)为有理分式时，它可以表示为两个多项式之比： N(s)bMsMbM1sM1b0X(s)D(s)aNsNaN1sN1a0 （3） 式（3）可以用部分分式法展成一下形式 X(s)rr1r2...Nsp1sp2spN （4） 通过查常用拉普拉斯变换对，可以由式（1-2）求得拉普拉斯逆变换。 利用 MATLAB 的residue 函数可以将 X (s)展成式（1-2）所示的部分分式展开式，该函数的调用格式为：[r,p,k] = residue(b,a) 其中b、a 为分子和分母多项式系数向量，r、p、k 分别为上述展开式中的部分分式系数、极点和直项多项式系数。 2、连续时间系统的系统函数 连续时间系统的系统函数是系统单位冲激响应的拉氏变换 H(s)h(t)estdt （5） 此外，连续时间系统的系统函数还可以由系统输入和系统输出信号的拉氏变换之比得到 H(s)Y(s)/X(s) （6） 单位冲激响应h(t)反映了系统的固有性质，而H(s)从复频域反映了系统的固有性质。由式（6）描述的连续时间系统，其系统函数为s的有理函数 bMsMbM1sM1b0H(s)aNsNaN1sN1a0 （7） 3、连续时间系统的零极点分析 系统的零点是指式（7）的分子多项式为零的点，极点指使分母多项式为零的点，零点使系统的值为零，极点使系统函数的值无穷大。通常将系统函数的零极点绘在s平面上，零点用表示，极点用表示，这样得到的图形称为零极点的分布图。 由零极点的定义可知，零点和极点分别指式（7）的分子多项式和分母多项式的根。利用MATLAB求多项式的根可以通过函数roots来实现，该函数的调用格式为： r=roots(c) c为多项式的系数向量，返回值r为多项式的根向量。 分别对式（7）的分子多项式和分母多项式求根即可得到零极点。 此外，在MATLAB中还提供了更简便的方法来求取零极点和绘制系统函数的零极点分布图，即利用pzmap函数，该函数的调用格式为： pzmap(sys)绘出由系统模型sys描述的系统的零极点分布图。 [p,z]=pzmap(sys) 这种调用方法返回极点和零点，而不绘出零极点分布图。其中sys为系统传函模型，由t命令sys=tf(b,a)实现，b、a为传递函数的分子多项式和分母多项式的系数向量。 MATLAB还为用户提供了两个专用函数tf2zp和zp2tf来实现系统传递函数模型和零极点增益模型的转换，其调用格式为： [z,p,k]=tf2zp(b,a) [b,a]=`zp2tf(z,p,k) 其中b、a为传递函数的分子多项式和分母多项式的系数向量，返回值z为零点列向量，p为极点列向量，k为系统函数零极点形式的增益。 研究系统函数的零极点分步不仅可以了解系统冲激响应的形式，还可以了解系统的频率特性以及判断系统的稳定性。 （1）零极点分布于冲激响应的关系 系统的极点位置决定着系统冲激响应h(t)的波形，冲激响应的幅值是由系统函数的零点和极点共同决定的，系统的零点位置只影响冲激响应的幅度和相位，不影响波形。 （2）零极点分步与系统频率响应的关系 系统函数的零极点分布不仅决定了系统函数H(s)，也决定了系统的频率响应H()，根据系统的零极点分布情况，可以由几何矢量法分析系统的频率响应。 （3）零极点分布与系统稳定性的关系 稳定性是系统的固有性质，与激励信号无关，由于系统函数H(s)包含了系统的所有固有的性质，因而可以根据系统函数的零极点分布判断系统的稳定性。因果系统稳定的充要条件是H(s)的全部极点位于s平面的左半平面。 三、实验内容 （1）、已知系统的冲激响应h(t)u(t)u(t2),输入信号x(t)u(t)，是采用复频域的方法求解系统的响应，编写MATLAB程序实现。 使用卷积定理求解，先分别求h(t)和x(t)的拉氏变换H(s)和X(s)然后根据式（6）求出输出Y(s)H(s)X(s)，最后对Y(s)进行拉普拉斯反变换即可得到系统的响应。 MATLAB程序如下： f1=sym('heaviside(t)-heaviside(t-2)'); f2=sym('heaviside(t)'); F1=laplace(f1); F2=laplace(f2); F=F1*F2; f=ilaplace(F); F f 结果： F = -(exp(-2*s)/s - 1/s)/s f = t - heaviside(t - 2)*(t - 2) 故系统响应为 y(t)t(t2)u(t2) （2）、已知因果连续时间系统的系统函数分别如下，试采用MATLAB画出其零极点分布图，求解系统的冲激响应 h(t)和频率响应H(),并判断系统是否稳定。 H(s)○1 1s32s22s1 MATLAB程序及运行结果如下: b=[1]; a=[1 2 2 1]; sys=tf(b,a); pzmap(sys); [p,z]=pzmap(sys); p z 结果： p = -1.0000 + 0.0000i -0.5000 + 0.8660i -0.5000 - 0.8660i z = Empty matrix: 0-by-1 系统没有零点，极点为p=-1，-0.5±0.866i； 另：pzmap(sys) Pole-Zero Map10.80.6Imaginary Axis (seconds-1)0.40.20-0.2-0.4-0.6-0.8-1-1.4-1.2-1-0.8-0.6-0.4-0.20Real Axis (seconds-1) b=[1]; a=[1 2 2 1]; sys=tf(b,a); pzmap(sys); [r,p,k]=residue(b,a); r p k 结果： r = 1.0000 + 0.0000i -0.5000 - 0.2887i -0.5000 + 0.2887i p = -1.0000 + 0.0000i -0.5000 + 0.8660i -0.5000 - 0.8660i k = [] 故系统冲激响应为： h(t)=(exp(-1)+(-0.5-0.2887*i).*exp(-0.5+0.866*i)+(-0.5+0.2887*i).*exp(-0.5-0.866*i)).*heaviside(t) Impulse Response0.450.40.350.3Amplitude0.250.20.150.10.050-0.0502468101214Time (seconds)频率响应： b=[1]; a=[1 2 2 1]; [H w]=freqs(b,a); subplot(211); plot(w,abs(H)); xlabel('\\omega(rad/s)'); ylabel('Magnitude'); title('|H(j\\omega)|'); grid on; subplot(212); plot(w,angle(H)); xlabel('\\omega(rad/s)'); ylable('Phase'); title('\\phi(\\omgea)'); grid on; |H(j)|1Magnitude0.50012345(rad/s)678910420-2-4012345(rad/s)678910由于该因果系统的所有极点都位于S 平面的左半平面，所以系统是稳定的。 s21H(s)5s2s43s33s23s2 ○2 MATLAB程序及运行结果如下: b=[1 0 1]; a=[1 2 -3 3 3 2]; sys=tf(b,a); pzmap(sys); [p,z]=pzmap(sys); p z 结果： p = -3.1704 + 0.0000i 0.9669 + 0.9540i 0.9669 - 0.9540i -0.3817 + 0.4430i -0.3817 - 0.4430i z = 0.0000 + 1.0000i 0.0000 - 1.0000i 另：pzmap(sys) Pole-Zero Map10.80.6Imaginary Axis (seconds-1)0.40.20-0.2-0.4-0.6-0.8-1-3.5-3-2.5-2-1.5-1-0.500.51Real Axis (seconds-1) 则有系统的零点为z=±i，极点为p=-3.1704，0.9669±0.9540i，-0.3817±0.4430i， b=[1 0 1]; a=[1 2 -3 3 3 2]; [r p k]=residue(b,a); r p k 结果： r = 0.0769 + 0.0000i -0.0300 - 0.0881i -0.0300 + 0.0881i -0.0085 - 0.1436i -0.0085 + 0.1436i p = -3.1704 + 0.0000i 0.9669 + 0.9540i 0.9669 - 0.9540i -0.3817 + 0.4430i -0.3817 - 0.4430i k = [] 故系统冲激响应为： h(t)=(0.0769.*exp(-3.1704.*t)+(-0.03-0.0881.*i).*exp((0.9669+0.954.*i).t)+(-0.03+ 0.0881.*i).*exp((0.9669 - 0.9540.*i).*t)+(-0.0085 - 0.1436.*i).*exp((-0.3817 + 0.4430.*i).*t)+(-0.0085+0.1436i).*exp((-0.3817-0.4430.*i).t)).*heaviside(t); 10.50-0.5-1-1.5-2-2.5x 1028Impulse ResponseAmplitude010203040506070Time (seconds) 频率响应： b=[1 0 1]; a=[1 2 -3 3 3 2]; [H,w]=freqs(b,a); subplot(211); plot(w,abs(H)); xlabel('\\omega(rad/s)'); ylabel('Magnitude'); title('|H(\\omega)|'); grid on; subplot(212); plot(w,angle(H)); xlabel('\\omega(rad/s)'); ylabel('Phase'); title('Phi(\\omega)'); grid on; |H()|0.8Magnitude0.60.40.20012345(rad/s)Phi()67891021Phase0-1-2012345(rad/s)678910由于该因果系统的所有极点不全位于S 平面的左半平面，所以系统是不稳定的。 （3）、已知连续时间系统函数的极点位置分别如下所示（设系统无零点），试用MATLAB绘制6中不同情况下，系统函数的零极点分布图，并绘制相应冲激响应的时域波形，观察并分析系统函数极点位置对冲激响应时域特性的影响。 ○1p=0 MATLAB程序及运行结果如下： H=sym('1/s'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 0]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = 1 121.510.50-6-4-20xPole-Zero Map246Imaginary Axis (seconds-1)10.50-0.5-1-1-0.8-0.6-0.4-0.200.20.40.60.81Real Axis (seconds-1) ○2p=-2 MATLAB程序及运行结果如下： H=sym('1/(s+2)'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 2]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = exp(-2*t) x 104exp(-2 t)420-6-5-4-3-2t-1012Imaginary Axis (seconds-1)Pole-Zero Map10.50-0.5-1-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.20Real Axis (seconds-1) ○3p=2 MATLAB程序及运行结果如下： =sym('1/(s-2)'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 -2]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = exp(2*t) r = 1 x 104exp(2 t)420-2-101t23456Imaginary Axis (seconds-1)Pole-Zero Map10.50-0.5-100.20.40.60.811.21.41.61.82Real Axis (seconds-1) ○4p12j,p22j MATLAB程序及运行结果如下： H=sym('1/(s^2+4)'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 0 4]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = sin(2*t)/2 sin(2 t)/20.50-0.5-6-4-20tPole-Zero Map420-2-4-1-0.8-0.6-0.4-0.200.20.40.60.81246Imaginary Axis (seconds-1)Real Axis (seconds-1) ○5p114j,p214j MATLAB程序及运行结果如下： H=sym('1/(s^2+2*s+17)'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 2 17]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = (sin(4*t)*exp(-t))/4 (sin(4 t) exp(-t))/41050-5-10-6-5-4-3-2-1t01234Imaginary Axis (seconds-1)Pole-Zero Map420-2-4-1.4-1.2-1-0.8-0.6-0.4-0.20Real Axis (seconds-1) ○6p114j,p214j MATLAB程序及运行结果如下： H=sym('1/(s^2-2*s+17)'); h=ilaplace(H); h subplot(211); ezplot(h); b=[1]; a=[1 -2 17]; sys=tf(b,a) subplot(212); pzmap(sys); 结果： h = (sin(4*t)*exp(t))/4 (sin(4 t) exp(t))/4100-10-4-3-2-10t123456Imaginary Axis (seconds-1)Pole-Zero Map420-2-400.20.40.60.811.21.4Real Axis (seconds-1) 由以上六例，可以总结出，在无零点的情况下： 当极点唯一且在原点时，h(t)为常数； 当极点唯一且是负实数时，h(t)为递减的指数函数； 当极点唯一且是正实数时，h(t)为递增的指数函数； 当H（s）有两个互为共轭的极点时，h(t)有一sint因子； 当H（s）有两个互为共轭的极点且他们位于右半平面时，h(t)还有一e因子； 当H（s）有两个互为共轭的极点且他们位于左半平面时，h(t)还有一e因子； tt（4）、已知3个连续时间系统函数 1s2s17s8(2)H(s)2ss17s8(3)H(s)2ss17 (1)H(s)上述三个系统具有相同的极点，只是零点不同，试用MATLAB分别绘制系统的零极点分布图及相应冲激响应的时域波形，观察并分析系统函数零点位置对冲激响应时域特性的影响。 ○1MATLAB程序及运行结果如下： b=[1]; a=[1 2 17]; sys=tf(b,a) subplot(211); pzmap(sys); H=sym('1/(s^2-2*s+17)'); h=ilaplace(H); h subplot(212); ezplot(h); 故h(t)=-0.125*i*exp(-1+4*i)+0.125*i*exp(-1-4*i) h=(-0.125.*i.*exp((-1+4.*i).*t)+0.125.*i.*exp((-1-4.*i).*t)).*heaviside(t); (sin(4 t) exp(t))/4100-10-4-3-2-10t123456Imaginary Axis (seconds-1)Pole-Zero Map420-2-400.20.40.60.811.21.4Real Axis (seconds-1) ○2MATLAB程序及运行结果如下： b=[1 8]; a=[1 2 17]; sys=tf(b,a) subplot(211); pzmap(sys); H=sym('(s+8)/(s^2-2*s+17)'); h=ilaplace(H); h subplot(212); ezplot(h); 故h(t)=((0.5-0.875i)exp((-1+4i)*t)+(0.5+0.875i)exp((-1-4i)*t)).*heaviside(t); h(t)=((0.5-0.875.*i).*exp((-1+4.*i).*t)+(0.5+0.875.*i).*exp((-1-4.*i).*t)).*heaviside(t); Imaginary Axis (seconds-1)Pole-Zero Map420-2-4-8-7-6-5-4-3-2-10Real Axis (seconds-1)exp(t) (cos(4 t) + (9 sin(4 t))/4)500-50-6-4-20t246 ○3MATLAB程序及运行结果如下： b=[1 -8]; a=[1 2 17]; sys=tf(b,a) subplot(211); pzmap(sys); H=sym('(s-8)/(s^2-2*s+17)'); h=ilaplace(H); h subplot(212); ezplot(h); 故h(t)=(0.5+1.125i)*exp((-1+4i)t)+(0.5-1.125i)*exp((-1-4i)*t); h=((0.5+1.125.*i).*exp((-1+4.*i).*t)+(0.5-1.125.*i).*exp((-1-4.*i).*t)).*heaviside(t); Imaginary Axis (seconds-1)Pole-Zero Map420-2-4-2-1012345678Real Axis (seconds-1)exp(t) (cos(4 t) - (7 sin(4 t))/4)500-50-6-4-20t246由以上三例可以看出，当极点不变时，零点分布只影响系统时域响应的幅度和相 位，对时域响应模式没有影响。即不会改变是衰减振荡还是增长振荡。 四、体会与建议 通过此次实验我掌握了用MTALAB实现拉普拉斯变换和拉普拉斯反变换的方法，学习了连续时间系统的复频域分析方法，信号的Laplace变换十分基本而重要，这次的学习使我对Laplace变换有了更深的理解。 实验六 离散时间系统的Z域分析 一、实验目的 ①掌握z变换及其反变换的定义，并掌握MATLAB实现方法。 ②学习和掌握离散时间系统系统函数的定义及z域分析方法。 ③掌握系统零极点的定义，加深理解系统零极点分布于系统特性的关系。 二、实验原理与方法 1、z变换 序列x(n)的z变换定义为 XzZ反变换定义为 nxnzn xnrXzzn1dz 在MATLAB中，可以采用了符号数学工具箱的ztrans函数和iztrans函数计算z变换和z反变换： Z=ztrans(F)求符号表达式F的z变换。 F=iztrans(Z)求符号表达式Z的z反变换。 2、离散时间系统的系统函数 离散时间系统的系统函数H(z)定义为单位抽样响应h(n)的z变换 Hznhnzn 此外，连续时间系统的系统函数还可以由系统输入和输出信号的z变换之比得到 HzYz Xz 由上式描述的离散时间系统的系统函数可以表示为 b0b1z1Hza0a1z1bMzM aNzN 3、离散时间系统的零极点分析 离散时间系统的零点和极点分别指使系统函数分子多项式和分母多项式为零的点，在MATLAB中可以通过函数roots来求系统函数分子多项式和分母多项式的根，从而得到系统的零极点。 此外，还可以利用MATLAB的zplane还是来求解和绘制离散系统的零极点分布图，zplane函数的调用格式为： zplane(b,a) b、a为系统函数的分子、系统函数的几点分母多项式的系数向量（行向量）。 Zplane(z,p) z、p为零极点序列（列向量）。 系统函数是描述系统的重要物理量，研究系统函数的零极点分布不仅可以了解系统单位抽样响应的变化，还可以了解系统的频率特性响应以及判断系统的稳定性： ①系统函数的极点位置决定了系统单位抽样响应h(n)的波形，系统函数零点位置只影响冲激响应的幅度和相位，不影响波形。 ②系统的频率响应取决于系统函数的零极点，根据系统的零极点分布情况，可以通过向量法分析系统的频率响应。 ③因果的离散时间系统稳定的充要条件是H(z)的全部几点都位于单位圆内。 三、实验内容 1、已知因果离散时间系统的系统函数分别为 z22z1H(z)3z0.5z20.005z0.3 （1） MATLAB代码及运行结果如下： b=[1 2 1]; a=[1 0.5 0.005 0.3]; subplot(221); zplane(b,a); [r,p,k]=residue(b,a); r p k subplot(222); impz(b,a,0:10); [H,w]=freqs(b,a); subplot(223); plot(w,abs(H)); xlabel('\\omega(rad/s)'); ylabel('Magnitude');; title('|H(j\\omega)|;') grid on; subplot(224); plot(w,angle(H)); xlabel('\\omega(rad/s)'); ylabel('Phase'); title('\\phi(\\omega)'); grid on; Impulse Response11.510.50-0.5Imaginary Part0-0.5-12-1-0.500.5Real Part|H(j)|;1 Amplitude0.505n (samples)()10642MagnitudePhase05(rad/s)10420-2-405(rad/s)100 r = 0.0098 + 0.0000i 0.4951 - 1.1965i 0.4951 + 1.1965i p = -0.8809 + 0.0000i 0.1905 + 0.5516i 0.1905 - 0.5516i k = [] z3z22H(z)4323z3zz3z1 （2） MATLAB代码及程序运行结果如下： b=[1 -1 0 2]; a=[3 3 -1 3 -1]; subplot(221); zplane(b,a); [r,p,k]=residue(b,a); r p k subplot(222); impz(b,a,0:10); [H,w]=freqs(b,a); subplot(223); plot(w,abs(H)); xlabel('\\omega(rad/s)'); ylabel('Magnitude');; title('|H(j\\omega)|;') grid on; subplot(224); plot(w,angle(H)); xlabel('\\omega(rad/s)'); ylabel('Phase'); title('\\phi(\\omega)'); grid on; Impulse Response140200-20-40Imaginary Part0-0.5-1-10Real Part|H(j)|;1 Amplitude0.505n (samples)()1030-1MagnitudePhase05(rad/s)1021-2-3-405(rad/s)100 r = 0.2263 + 0.0000i -0.2071 + 0.2555i -0.2071 - 0.2555i 0.5213 + 0.0000i p = -1.6462 + 0.0000i 0.1614 + 0.7746i 0.1614 - 0.7746i 0.3234 + 0.0000i k = [] 试采用MATLAB绘出其零极点分布图，求解系统的冲激响应h（n）和频率响应,并判断系统是否稳定。 2、已知离散时间系统系统函数的零点z和极点p分别为： （1）z=0，p=0.25 MATLAB代码及运行结果如下： b=[1 0]; a=[1 0.25]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-101Real PartImpulse Response231 Amplitude0.50-0.50123456n (samples)78910 （2）z=0，p=1 MATLAB代码及运行结果如下： b=[1 0]; a=[1 1]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-101Real PartImpulse Response231 Amplitude0.50-0.5-10123456n (samples)78910 （3）z=0，p=-1.25 MATLAB代码及运行结果如下： b=[1 0]; a=[1 1.25]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-10Real PartImpulse Response12310 Amplitude50-5-100123456n (samples)78910 p0.8e（4）z=0，1j6p0.8e,2j6 MATLAB代码及运行结果如下： b=[1 0]; a=[1 0.4*3^0.5 0.64]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-101Real PartImpulse Response2321 Amplitude0.50-0.5-10123456n (samples)78910 j8（5）z=0，p1epe，1j8 MATLAB代码及运行结果如下： b=[1 0]; a=[1 -2*cos(pi/8) 1]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-101Real PartImpulse Response2324 Amplitude20-2-40123456n (samples)78910 （6）z=0， p11.2ej34， p11.2ej34 MATLAB代码及运行结果如下： b=[1 0]; a=[1 -1.2*2*cos(0.75*pi) 1.44]; subplot(211); zplane(b,a); subplot(212); impz(b,a,0:10); 1Imaginary Part0.50-0.5-1-3-2-101Real PartImpulse Response2324 Amplitude20-2-40123456n (samples)78910 试用MATLAB绘制上述6种不同情况下，系统函数的零极点分布图，并绘制相应单位抽样响应的时域波形，观察分析系统函数极点位置对单位抽样响应时域特性分影响和规律。 3、已知离散时间系统的系统函数分别为： H(Z)（1） z(z2)(z0.8ej6)(z0.8ej6) MATLAB代码及运行结果如下： z=[0 -2]'; p=[0.8*exp(j*pi/6) -0.8*exp(j*pi/6)]'; zplane(z,p); 10.5Imaginary Part0-0.5-1-2-1.5-1-0.5Real Part00.51 H(Z)（2） z(z2)(z0.8ej6)(z0.8ej6) MATLAB代码及运行结果如下： z=[0 2]'; p=[0.8*exp(j*pi/6) -0.8*exp(j*pi/6)]'; zplane(z,p); 10.5Imaginary Part0-0.5-1-1-0.500.5Real Part11.52 上述两个系统具有相同的极点，只是零点不同，试用MATLAB分别绘制上述两个系统的零极点分布图及响应的时域波形，观察分析系统函数零点位置对单位抽样响应时域特性的影响。 四、体会与建议 通过此次实验我掌握了用MTALAB实现Z变化和Z反变换的方法，学习了离散时间系统的复频域分析方法，信号的Z变换十分基本而重要，这次的学习使我对Z变换有了更深的理解。 因篇幅问题不能全部显示，请点此查看更多更全内容