MATLAB课程设计任务书（-）
一、名称： MATLAB编程简介 二、目的
熟悉MATLAB编程环境，掌握Help 命令、基本的变量类型、矩阵的基本运算、基本的绘图函数和M-file的建立。 三、内容 （一） 例 题
例1、 Help 命令 help cos help plot help abs help exp help +
例2、 变量和矩阵运算
(1) Matrix——The basic variable type
M=3
M=[1 2 6]
M=[1 2 6; 4 6 7] M13=M(1,3) size(M)
(2) The Colon Operator ( : )
%Creating Array and Vector % v = start: skip: end x1=0:2:10
x2=0:1:10 (or x=0:10) t=-1:0.2:1
?cessing Matrix
A=[1 2 3;4 5 6;7 8 9] A(2:3,1:2) x2(4:8)
(3) Matrix Operations (A±B)
A=[2 3 4; 6 9 8] B=[1 2 3; 5 8 7] C1=A+B C2=A-B C3=A-4
(4) Matrix Operations (A*B A.*B)
% A*B A=[2 3 4; 6 9 8]
1
B=[1 2; 3 5; 8 7] A*B
% A.*B
A=[2 3 4; 6 9 8] B=[1 2 3 ;5 8 7] A.*B
(5) Matrix Operations (B/A ,A\\C, B./A ,A.\\C)
% B/A —— B*inv(A) % A\\C —— inv(A)*C %B./A —— B(i,j)/A(i,j) %A.\\B —— A(i,j)\\B(i,j)
(6) Matrix Operations ( ^ and .^) % ^ Operation
A=[1 2 3; 4 5 6; 7 8 9] b=A^2 % .^ Operation
A=[1 2 3; 4 5 6; 7 8 9] b=A.^2
(7) Matrix Operations ( Ａ′and Ａ. ′)
% Ａ′共轭转置
a=[1+2i 3+4i; 3+2i 5+5i] a′
% Ａ.′非共轭转置 a.′
例3 、 绘图函数plot(x,y) ，stem(k,y) % plot(x,y) x=0:0.01:2; y=sin(2*pi*x); plot(x,y) % stem(k,y) k=0:50;
y=exp(-0.1*k); stem(k,y)
例4、 M file
% y(t)=sin(2t) + sin(5t) -2pi ≤ t ≤ 2pi t =-2*pi:0.02:2*pi; y=sin(2*t) + sin(5*t); plot(t,y)
（二） 练 习 题
1、 基本命令 help plot
help colon help ops
2
help zeros help ones pi*pi-10
sin(pi/4) ans^2 zz=3+4i; conj(zz)
abs(zz) angle(zz) real(zz) imag(zz)
2、Array Indexing
xx=[ones(1,4),[2:2:11],zeros(1,3)] xx(3:7) length(xx)
xx(2:2:length(xx))
xx(3:7)=pi*(1:5)
3、 用以下语句建立M-file
t=-2:0.05:3;
y=sin(2*pi*0.789*t); plot(t,y), grid on
title('TEST PLOT of SINUSOID') xlabel('TIME(sec)')
（-1?t?2)（用M-file实现） 4、 画出以下信号的波形
x1(t)?2cos(2?t?30?) x2(t)?4cos(2?t?60?)
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
3
MATLAB课程设计任务书（二）
一、名称：常见信号的MATLAB表示 二、目的
掌握用MATLAB表示信号与系统中的常见信号， 熟悉MATLAB中一些常用的信号生成函数。 三、内容 （一） 例 题
例1 、 画图表示指数信号y?Aeat, A?1,a??0.4 [MATLAB程序]：
% decaying exponential signal
A=1;a=-0.4; t=0:0.01:10;
ft=A*exp(a*t);
plot(t,ft);grid on;
[运行结果]：
10.90.80.70.60.50.40.30.20.10012345678910例2 、 画图表示正弦信号
y?Asin(?0t??), A?1,?0?2?,???/6
[MATLAB程序]： % sinusoidal signal
A=1;w0=2*pi;phi=pi/6; t=0:0.001:8;
ft=A*sin(w0*t+phi); plot(t,ft);grid on;
[运行结果]：
4
10.80.60.40.20-0.2-0.4-0.6-0.8-1012345678
例3 、 画图表示抽样函数y?sin(?t)/(?t)。（利用函数sinc(t)） [MATLAB程序]： % sample function
t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft);grid on;
[运行结果]：
例4、 画图表示矩形信号，设幅度A=1，宽度为W=2。 [MATLAB程序]：
% rectangular pulse signal t=0:0.001:4; T=1;
ft=rectpuls(t-2*T,2*T);
plot(t,ft);grid on;axis([0 4 -0.5 1.5]);
[运行结果]：
10.80.60.40.20-0.2-0.4-10-8-6-4-202468101.510.50
例5、 画图表示周期形三角波信号，设幅度A=1，宽度为W=4，斜度k=0.5。 [MATLAB程序]：
% periodic triangular pulse signal
t=-5*pi:pi/10:5*pi; x=sawtooth(t,0.5);
plot(t,x);axis([-16 16 -1.5 1.5]);grid on;
[运行结果]：
-0.500.511.522.533.54 5
1.510.50-0.5-1-1.5 例6 、画图表示指数序列y?[MATLAB程序]： % exponential sequence
k=0:10;A=1;a=-0.6; fk=A*a.^k;
stem(k,fk);grid on;
[运行结果]：
10.80.60.40.20-0.2-0.4-0.6-15-10-5051015Aeak, A?1,a??0.6
012345678910例7、 画图表示单位样值序列。 [MATLAB程序]： % unit impulse sequence k=-25:25;
delta=[zeros(1,25),1,zeros(1,25)]; stem(k,delta);
[运行结果]：
10.90.80.70.60.50.40.30.20.10-25-20-15-10-5
0510152025
例8 、利用MATLAB画出f(2t)和
[MATLAB程序]：
%changed triangular pulse signal t=-3:0.001:3;
f(2?2t)的波形。
6
ft1=tripuls(2*t,4,0.5);
subplot(2,1,1);plot(t,ft1);title('f(2t)');grid on; ft2=tripuls((2-2*t),4,0.5);
subplot(2,1,2);plot(t,ft2);title('f(2-2t)');grid on;
[运行结果]：
f(2t)10.80.60.40.20-310.80.60.40.20-3-2-10123-2-10f(2-2t)123
例9 、利用MATLAB计算三角波信号的微分和积分。设三角波信号幅度A=1，宽度为W=4，斜度k=0.5。 [MATLAB程序]： % functri(t) function function yt=functri(t)
yt=tripuls(t,4,0.5);
%differentiation of triangular signal h=0.001;t=-3:h:3;
y1=diff(functri(t))*1/h;
figure(1);plot(t(1:length(t)-1),y1);title('df(t)/dt');grid on;
%integration of triangular signal t=-3:0.1:3;
for x=1:length(t)
y2(x)=quad('functri',-3,t(x)); end
figure(2);plot(t,y2);axis([-3 3 -0.5 2.5]);title('Integral of f(t)');grid on;
[运行结果]：
df(t)/dt0.40.20Integral of f(t)2.521.5-0.2-0.4-0.610.5-0.80-1-1.2-3-2-10123 （二） 练 习 题
1． 利用MATLAB实现下列连续时间信号 （1）f?t??t??t?，取t=0~10;
7
-0.5-3-2-10123
（2）f?t??10sin(100?t), 取t=0~0.2 （3）f?t??Sa??t?cos?20t?，取t=0~5 （4）f?t??4e?0.5tcos?t，取t=0~10
2．用tripuls函数画出信号f(t)??1?t?????t?1????t?1?????3?t?????t?2????t?3???的波形。
3．（1）编写表示信号f(t)?g1?t?0.5???1?t?????t????t?1????g1?t?1.5?的函数；
（2）画出f(t)、f(0.5t)、f(2?0.5t)的波形。 4．利用MATLAB实现下列离散序列 （1）f(k)???k?1? （2）f(k)???k?
（3）f(k)???k?2????k?5? （4）f(k)?k??k?
（5）f(k)?5(0.8)kcos?0.9?k?
5．利用square函数画出下式所示的离散周期方波序列的波形。
n?1 f(k)??n??1f0?k?，其中f0(k)???k??2??k?4????k?10?
6．已知序列f(k)??2?(k?2)???k?1??3??k????k?1????k?2????k?3??2??k?4? （1）用stem函数，画出序列f(k)的波形； （2）画出序列f(k?2)、序列f(k?4)的波形； （3）画出序列f(k)、f(3k)的波形；
3（4）利用fliplr函数实现序列f(?k)，并画出波形。
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
8
MATLAB课程设计任务书（三）
一、名称：连续时间系统时域分析的MATLAB实现 二、目的
掌握应用MATLAB实现对线性时不变连续时间系统的时域分析，即熟悉应用MATLAB 实现微分方程的求解、连续时间信号卷积计算等。 三、内容 （一） 例 题
例1 、已知系统的输入信号为
dy(t)dt22f(t)?10sin2?t，系统的初始状态为零，系统的微分方程为
?2dy(t)dt?100y(t)?f(t)，求y(t)。
[MATLAB程序]：
%solution of differential equation
ts=0; te=5; dt=0.01;
sys=tf([1],[1 2 100]); t=ts:dt:te;
f=10*sin(2*pi*t); 0.2y=lsim(sys,f,t); 0.15plot(t,y); 0.1xlabel('t(sec)'); 0.05ylabel('y(t)'); 0grid on; -0.05 -0.1[运行结果]： -0.15 -0.2 -0.2500.511.522.533.544.5t(sec)
例2 、在上题中，若输入信号为强度为10的冲激信号，求y(t)。 [MATLAB程序]：
%impulse response of LTI system ts=0; te=5;
y(t)5 9
dt=0.01;
sys=tf([10],[1 2 100]); t=ts:dt:te;
y=impulse(sys,t); plot(t,y);
xlabel('t(sec)'); ylabel('y(t)'); grid on;
[运行结果]：
例3、 在上题中，若输入幅度为10的阶跃信号，求y(t)。 [MATLAB程序]：
%step response of LTI system ts=0; te=5; 0.18dt=0.01;
0.16sys=tf([10],[1 2 100]);
0.14t=ts:dt:te;
0.12y=step(sys,t);
0.1plot(t,y);
0.08xlabel('t(sec)');
0.06ylabel('y(t)');
grid on; 0.04 0.02 [运行结果]： 000.511.522.533.544.55 t(sec) 例4 、 已知信号x(t)??(t)??(t?2)，y(t)??(t?3)??(t?5)，利用MATLAB计算x(t)?y(t)，
并画出卷积结果。 [MATLAB程序]：
%convolution of two signal
t=-1:0.01:6; T=1;
ft1=rectpuls(t-1,2*T);
subplot(3,1,1); plot(t,ft1);grid on;axis([-1 6 -0.5 1.5]);ylabel('x(t)'); ft2=rectpuls(t-4,2*T);
subplot(3,1,2); plot(t,ft2);grid on;axis([-1 6 -0.5 1.5]);ylabel('y(t)') z=conv(0.01*ft1,ft2);
t2=(0:length(z)-1)*0.01-2;
subplot(3,1,3); plot(t2,z);grid on;axis([0 10 -0.5 3]);ylabel('x(t)*y(t)')
[运行结果]：
10.80.60.40.2y(t)0-0.2-0.4-0.6-0.800.511.522.5t(sec)33.544.55
y(t)10
1.51x(t)0.50-0.5-11.510123456y(t)0.50-0.5-130123456x(t)*y(t)210012345678910（二） 练 习 题
1． 已知一线性系统的微分方程为 y??(t)?5?y(t?)6y?(t?)?t?()?t (（1） 试求该系统的零状态响应y(t)；
（2） 用lsim求出该系统的零状态响应的数值解，并比较不同的抽样间隔对数值解
的影响。
2．已知一线性系统的微分方程为
dy(t)dt22?4dy(t)dt?25y(t)?e(t)
试求：（1）用impulse函数求该系统的单位冲激响应； （2）用step函数求该系统的单位阶跃响应。
3．下列系统分别为一阶、二阶、三阶BW模拟低通滤波器，用impulse函数分别求出各系统的单位冲激响应，并比较它们的特征。 （1）y?(t)?y(t)?e(t)
（2）y??(t)?2y?(t)?y?t??e(t) （3）y????t??2y??(t)?2y?(t)?y?t??e(t)
4．下列系统分别为BW模拟低通、高通、带通、带阻滤波器，用impulse函数分别求出各系统的单位冲激响应，并比较它们的特征。 （1）y??(t)?2y?(t)?y?t??e(t)
（2）y??(t)?2y?(t)?y?t??e??(t) （3）y??(t)?y?(t)?y?t??e?(t) （4）y??(t)?y?(t)?y?t??e??(t)?e(t)
5．已知某LTI连续系统的单位冲激响应h(t)和各激励信号e(t)的波形如下图所示，试求该系统对各激励信号的零状态响应。
11
h(t) 1 0 2 t 1 e1(t) 1 2 t e2(t) 1 1 2 t e3(t) 0 0 0 2 t
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
MATLAB课程设计任务书（四）
一、名称：连续时间系统频域分析的MATLAB实现 二、目的
掌握应用MATLAB实现对线性时不变连续时间系统的频域分析，掌握应用MATLAB近似计算和绘制信号的频谱、连续时间系统的频率响应（幅频响应和相频响应）。 三、内容
（一） 例 题
例1 、利用MATLAB画出下图所示的周期三角波信号的频谱。
f (t) 1 -0.5 -3 -2 -1 0 0.5 -1 1 2 t
n???4j?sin() n?0经计算，该周期三角波信号的傅立叶级数系数为c(n)??n2?2 2?? 0 n?0
[MATLAB程序]：
%frequency spectrum of periodic triangular signal N=10; n1=-N:-1;
c1=-4*j*sin(n1*pi/2)/pi^2./n1.^2; %计算n=-N到－1时的傅立叶级数系数 c0=0; % 计算n=0时的傅立叶级数系数 n2=1:N;
12
c2=-4*j*sin(n2*pi/2)/pi^2./n2.^2; %计算n=1到N时的傅立叶级数系数 cn=[c1 c0 c2]; n=-N:N;
subplot(2,1,1);stem(n,abs(cn));ylabel('Cn的幅度'); subplot(2,1,2);stem(n,angle(cn));ylabel('Cn的相位'); xlabel('\\omega/\\omega_0');
[运行结果]：
0.50.4Cn的幅度0.30.20.10-1021-8-6-4-20246810Cn的相位0-1-2-10-8-6-4-20246810
?/?0
?1-|t| |t|?1? 0 |t|?1例2 、利用MATLAB采用数值方法近似计算三角波信号f(t)??的频谱。
[MATLAB程序]： %example2 的函数 function y=sf1(t,w)
y=(abs(t)<=1).*(1-abs(t)).*exp(-j*w*t);
% Frequency spectrum of triangular signal using quad8 function
w=linspace(-6*pi,6*pi,512); N=length(w); F=zeros(1,N); for k=1:N
F(k)=quadl(@sf1,-1,1,[],[],w(k)); end
figure(1);
plot(w,real(F));
xlabel('\\omega');ylabel('F(j\\omega)');title('三角波信号近似频谱'); figure(2);
plot(w,real(F)-sinc(w/2/pi).^2); %这里的sinc函数值即为理论计算结果 xlabel('\\omega');ylabel('error');title('计算误差');
[运行结果]：
13
三角波信号近似频谱10.5x 10-9计算误差0.80-0.5-10.6F(j?)error0.4-1.5-2-2.50.20-3-0.2-20-15-10-505101520-3.5-20-15-10-505101520??
例3 、利用MATLAB画出???0.9时F(ej?)?11??e?j?的幅度频谱。
[MATLAB程序]：
% amplitude frequency spectrum of a sequence
b=[1];a1=[1 -0.9];a2=[1 0.9];
w=linspace(0,2*pi,512); %线性均分0到2pi的间隔，共512个点 h1=freqz(b,a1,w); h2=freqz(b,a2,w);
plot(w/pi,abs(h1),w/pi,abs(h2),'k:'); xlabel('\\Omega/\\pi');
legend('\\alpha=0.9','\\alpha=-0.9');
12[运行结果]： ?=0.9?=-0.9 10
8 6
4 2
0 00.20.40.60.811.21.41.61.82?/?
例4 、某连续线性系统的频率响应为：H(j?)?1(j?)?2(j?)?2(j?)?132
利用MATLAB画出该系统的幅频响应|H(j?)|和相频响应?(?)。 [MATLAB程序]：
% frequency response of a system b=[1];a=[1,2,2,1];
w=linspace(0,5,200); H=freqs(b,a,w); subplot(2,1,1); plot(w,abs(H));
set(gca,'xtick',[0 1 2 3 4 5]); set(gca,'ytick',[0 0.4 0.707 1]); xlabel('\\omega(rad/s)'); ylabel('|H(j\\omega)|'); grid on;
subplot(2,1,2); plot(w,angle(H));
set(gca,'xtick',[0 1 2 3 4 5]);
14
xlabel('\\omega(rad/s)'); ylabel('\\phi(rad)'); grid on;
1
0.707[运行结果]：
0.4
（二）练 习 题
|H(j?)|0012345?(rad/s)42?(rad)0-2-4012345
1． 求下图所示周期矩形脉冲信号的傅立叶级数表示式。并用MATLAB求出由前N次谐波合成的信号近似波形。
f(t) …… -2
??(rad/s)1 …… 2
t
2．求方波信号f(t)??n???g??t?nT0?的幅度谱，并画出幅度谱
-0.5 0 0.5
周期矩形脉冲信号
。分别取T0?2?,4?,8?，
讨论周期T0与频谱的关系。
3．利用MATLAB采用数值方法近似计算矩形脉冲信号f(t)??(t?1)??(t?1)的频谱，画出幅度谱。
4．信号f1(t) 和f2(t)如下图所示。
（1） 取t=0:0.05:2.5，画出f(t)?f1(t)?f2(t)cos(50t)的波形。
f1(t) 1 0 1 t 1 0 1 2 t f2(t) （2） 一系统的H(j?)为：
H(j?)?1043422
(j?)?26.131(j?)?3.4142?10(j?)?2.6131?10?j???1034 用freqs画出H(j?)幅频响应和相频响应曲线。
（3）用lsim函数求出信号f(t)和f(t)cos(50t)通过上述系统的响应y1(t)和y2(t)，并根据理论知识解释所得的结果。
四、要求
15
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
MATLAB课程设计任务书（五）
一、名称：连续时间系统复频域分析的MATLAB实现 二、目的
掌握用MATLAB实现连续时间系统S域分析中的零极点求解、拉氏变换及反拉氏变换、由系统函数求冲激响应、频率响应等。 三、内容 （一） 例 题
例1、 用部分分式展开法求F(s)的Laplace反变换：F(s)?s?2s?4s?3s32
[MATLAB程序]：% inverse Laplace transform by partial-fraction expansion
format rat; % 将分数以近似的小整数之比的形式显示，输出格式 num=[1 2]; den=[1 4 3 0];
[r,p]=residue(num,den) % r,p均为列向量
[运行结果]：
r =
-1/6 -1/2 2/3
p =
-3 -1 0
即，由程序已算出
r??[?16,?12,], p??[?3,?1,0]， 23??1/2s?1??1/6s?3则
F(s)?2/3s，
16
反变换为
f(t)?23u(t)?12e?tu(t)?16e?3tu(t)
例2 、 用部分分式展开法求F(s)的Laplace反变换：F(s)?s?2s(s?1)3
[MATLAB程序]：% inverse Laplace transform using conv function and partial-fraction expansion
num=[1 -2];
a=conv([1 0],[1 1]); b=conv([1 1],[1 1]); den=conv(a,b);
[r,p]=residue(num,den) % r,p均为列向量
%第二种解法
% inverse Laplace transform using poly function and partial-fraction expansion num=[1 -2];
den=poly([0,-1,-1,-1]);
[r,p]=residue(num,den) % r,p均为列向量
[运行结果]：
r =
2 2 3 -2
p =
-1 -1 -1 0
即，用程序可算出
r??[2,2,3,?2], p??[?1,?1,?1,0]，
则
F(s)?2s?1?2(s?1)?t2?3(s?1)?t3?2?2s?t，
?2)u(t)
2s3反变换为
f(t)?(2e?2te?1.5te例3、 用部分分式展开法求F(s)的Laplace反变换：F(s)??3s22?5(s?1)(s?s?2)
[MATLAB程序]：% inverse Laplace transform by partial-fraction expansion
num=[2 3 0 5];
den=conv([1 1],[1 1 2]);
[r,p,k]=residue(num,den) % r,p均为列向量 [angle,mag]=cart2pol(real(r),imag(r))
[运行结果]：
r =
-2 + 2024/1785i -2 - 2024/1785i 3
p =
17
-1/2 + 1012/765i -1/2 - 1012/765i -1
k =
2
angle =
1972/751
-1972/751
0
mag =
7895/3434 7895/3434 3
即，用程序可算出
r??[-2.0000 + 1.1339i, -2.0000 - 1.1339i, 3.0000], p??[-0.5000 + 1.3229i, -0.5000 - 1.3229i, -1.0000] , k?2 angle??[2.6258, -2.6258, 0], mag??[2.2991, 2.2991, 3.0000]
?2.2991e?2.6258i所以
F(s)?2?3s?1??t?2?1.1339is?0.5?1.3229i??2?1.1339is?0.5?1.3229i?2?3s?1s?0.5?1.3229i?2.2991e2.6258is?0.5?1.3229i
f(t)?2?(t)?3eu(t)?4.5982e?0.5tcos(1.3229t?2.6258)u(t)
f(t)?e?t例4、 分别利用MATLAB中的laplace和ilaplace函数求：（1）Laplace变换；（2）F(s)?ss22sin(at)u(t)的
?1的Laplace反变换。
[MATLAB程序]：%(1) Laplace transform using laplace function
f=sym('exp(-t)*sin(a*t)'); F=laplace(f)
%(2) Inverse Laplace transform using ilaplace function F=sym('s^2/(s^2+1)'); ft=ilaplace(F)
[运行结果]：
F =
a/((s+1)^2+a^2) ft =
Dirac(t)-sin(t)
即，结果为
例5、 已知系统函数为H(s)?s?1s?2s?22F?a/((s+1)^2+a^2)，
ft?Dirac(t)?sin(t), 即 f(t)??(t)?sin(t)u(t)
，利用MATLAB求出该系统的零极点，并画
出零极点分布图。
[MATLAB程序]：% pole-zero map of H(s) using plot function
b=[1 -1]; a=[1 2 2]; zs=roots(b); ps=roots(a);
plot(real(zs),imag(zs),'o',real(ps),imag(ps),'kx','markersize',12);
18
axis([-2 2 -2 2]); grid on;
legend('零点','极点');
[运行结果]：
21.510.50-0.5-1-1.5-2-2
零点极点-1.5-1-0.500.511.52
例6、 已知系统函数为H(s)?1s?2s?2s?132
，利用MATLAB画出该系统的零极点分
布图；求出该系统的单位冲激响应h(t)的幅频响应，并判断系统是否稳定。
[MATLAB程序]：% impulse response, amplitude frequency response and stability analysis % of LTI H(s)
num=[1]; den=[1 2 2 1]; sys=tf(num,den); poles=roots(den); figure(1); pzmap(sys); t=0:0.02:10;
h=impulse(num,den,t); figure(2); plot(t,h);
xlabel('t(s)');ylabel('h(t)');title('Impulse Response'); [H,w]=freqs(num,den); figure(3);
plot(w,abs(H));
xlabel('ang.freq.\\omega(rad/s)'); ylabel('|H(j\\omega)|');
title('Magnitude Response');
[运行结果]：
19
Impulse ResponsePole-Zero Map10.80.60.40.45零点极点0.40.350.30.25Imaginary Axis0.2h(t)-1.2-1-0.8-0.6-0.4-0.200-0.2-0.4-0.6-0.8-1-1.40.20.150.10.050-0.05012345t(s)678910
10.90.80.70.6|H(j?)|Real Axis
Magnitude Response0.50.40.30.20.100123456ang.freq.?(rad/s)78910（二）练 习 题
1、已知连续时间信号的S域表示式，试用residue求出F(s)的部分分式展开式，并写出f(t)的表达式（f(t)要写成实数）：
(1) F?s??41.6667s?3.7444s?25.7604s?41.666732
(2) F?s??s3?s?5??s2?5s?25?
2、已知某连续时间系统的微分方程为
y???t??4y??t??3y?t??2f??t??f?t?
f?t????t?,y?0???1,y??0???2，试求系统的零输入响应、零状态响应和全响应，并
画出相应的波形。
3、已知系统函数为H?s??s?2s?2s?2s?1321s?2?s?12,试分别画出??0,1,1,2时系统的零极点图。
4如果系统是稳定的，画出系统的幅度响应曲线。 4、已知F?s??，画出该系统的零极点分布图，求出系统的冲激响应、
阶跃响应和频率响应。
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
20
MATLAB课程设计任务书（六）
一、名称：离散时间系统时域分析的MATLAB实现 二、目的
掌握用MATLAB实现离散时间系统的零状态响应、冲激响应、阶跃响应、卷积和的求解等。 三、内容 （一）例 题
例1、 受噪声干扰的信号为
f(k)?s(k)?d(k)，其中s(k)?(2k)0.9k是原始信号，d(k)是噪
MATLAB编程
声。已知M点滑动平均系统的输入输出关系为：y(k)?1MM?1?n?0f(k?n)。试用
实现M点滑动平均系统对噪声干扰的信号去噪。
[提示：噪声信号可用MATLAB中的rand函数产生，将其叠加在原始信号上，即得到受噪声干扰的输入信号f(k)]。
[MATLAB程序]：%signal smoothing by moving average filter
R=51; %length of input signal
%generate(-0.5,0.5)uniformly distributed random numbers d=rand(1,R)-0.5; k=0:R-1;
s=2*k.*(0.9.^k); f=s+d; figure(1);
plot(k,d,'r-.',k,s,'b:',k,f,'k-'); xlabel('Time index k'); legend('d[k]','s[k]','f[k]'); M=5;
b=ones(M,1)/M; a=1;
y=filter(b,a,f); figure(2);
plot(k,s,'b:',k,y,'r-'); xlabel('Time index k');
21
legend('s[k]','y[k]');
[运行结果]：
876558d[k]s[k]f[k]76s[k]y[k]4433210-1210051015202530Time index k35404550051015202530Time index k35404550
例2、用impz函数求下列离散时间系统的单位样值响应h(k)，并与理论值
kkh(k)??(?1)?2(?2), k?0进行比较。y(k)?3y(k?1)?2y(k?2)?f(k) [MATLAB程序]：%impulse response of discrete system
k=0:10; a=[1,3,2]; b=[1];
h=impz(b,a,k); subplot(2,1,1); stem(k,h);
title('单位样值响应的近似值'); grid on;
hk=-(-1).^k+2*(-2).^k; subplot(2,1,2); stem(k,hk); title('单位样值响应的理论值'); grid on;
[运行结果]：
单位样值响应的近似值3000200010000-1000-2000012345678910单位样值响应的理论值3000200010000-1000-2000012345678910
例3、 已知序列x(k)={1,2,3,4;n=0,1,2,3}，y(k)={1,1,1,1,1;n=0,1,2,3,4}，利用MATLAB计算x(k)?y(k)，并画出卷积结果。
[MATLAB程序]：%convolution of two sequences
22
x=[1,2,3,4]; y=[1,1,1,1,1]; z=conv(x,y);
subplot(3,1,1);stem(0:length(x)-1,x);ylabel('x(k)'); subplot(3,1,2);stem(0:length(y)-1,y);ylabel('y(k)');
subplot(3,1,3);stem(0:length(z)-1,z);ylabel('x(k)*y(k)');xlabel('k');
[运行结果]：
4x(k)20100.511.522.53y(k)0.501000.511.522.533.54x(k)*y(k)500123k4567
（二）练 习 题
1．已知f(k)??f?N?1?k?,h?k???h?N?1?k?，用conv函数计算f(k)?h?k?，并总结f(k)?h?k?的类型。（提示：自己构造满足条件的两个序列） 2．利用impz函数，计算系统
y(k)?0.7y?k?1??0.45y?k?2??0.6y?k?3??0.8f?k??0.44f?k?1??0.36f?k?2??0.02f?k?3?
的单位脉冲响应，并画出前31点的图。 3．利用filter 函数，求出系统y(k)?1.85y?k?1??0.85y?k?2??f?k?的单位脉冲响应，并判断系统是否稳定。
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
23
MATLAB课程设计任务书（七）
一、名称：离散时间系统Z域分析的MATLAB实现 二、目的
掌握用MATLAB实现离散时间系统Z域分析中的零极点求解、Z变换及反Z变、由系统函数求冲激响应、频率响应等。 三、实验内容 （一）例 题
例1、 利用MATLAB计算F(z)?1818?3z?1?4z?2?z?3的部分分式展开式。
[MATLAB程序]：% Partial-fraction expansion of F(z)
num=[18];
den=[18 3 -4 -1];
[r,p k]=residuez(num,den)% r,p为列向量
[运行结果]：
r =
0.3600 0.2400 0. 4000
p =
0.5000 -0.3333 -0.3333
k = [ ] 即，由程序已算出
r??[0.3600, 0.2400, 0.4000], p??[0.5000, -0.3333, -0.3333] , k?[]
所以有
F(z)?0.361?0.5z?1?0.241?0.3333z?1?0.4(1?0.3333z?1)2
z1?0.5z?1例2、 已知一离散因果LTI系统的系统函数为：H(z)?该系统的零极点。
说明：先将系统改写为：H(z)?zz32?2z?2?z?2?3?1?0.005z?0.3z?3，求
?2z?12?0.5z?0.005z?0.3，然后用tf2zp函数求系统的零极点。
[MATLAB程序]：% zeros and poles of H(z)
b=[1 2 1];
a=[1 -0.5 -0.005 0.3];
[r,p k]=tf2zp(b,a) % r,p为列向量
24
zplane(b,a);
[运行结果]：
r =
-1 -1
p =
0.5198 + 0.5346i 0.5198 - 0.5346i -0.5396
k =1
求出的零点为：z??1(2阶)， z?-0.5396。 z?0；极点为：z?0.5198 ? 0.5346i，零极点图如下：
10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.512
zz32例3、 已知一离散因果LTI系统的系统函数为：H(z)??2z?12?0.5z?0.005z?0.3。利用
MATLAB画出该系统的零极点分布图，求系统的单位样值响应和幅频响应，并判断系统的稳定性。
[MATLAB程序]：% Impulse response, amplitude frequency response and stability analysis of %LTI H(z)
b=[0 1 2 1];
a=[1 -0.5 -0.005 0.3]; figure(1); zplane(b,a); num=[0 1 2 1];
den=[1 -0.5 -0.005 0.3]; h=impz(num,den) figure(2); stem(h); xlabel('k'); ylabel('h(k)');
title('Impulse Response'); [H,w]=freqz(num,den); figure(3);
plot(w/pi,abs(H));
xlabel('ang.freq.\\Omega(rad/s)'); ylabel('|H(e^j^\\Omega)|');
25
title('Magnitude Response');
[运行结果]：
Impulse Response2.510.80.60.41.52Imaginary Part0.20-0.2-0.4-0.6-0.8-11h(k)0.50-0.5-0.50Real Part0.51-102
8765|H(ej?)|-151015k20253035
Magnitude Response43210
该系统为非稳定系统，∵ 在右半平面存在极点
例4、 分别利用MATLAB中的ztrans和iztrans函数求：（1）f(k)?cos(ak)u(k)的z变换；
（2）F(z)?1(1?z)200.10.20.30.40.50.6ang.freq.?(rad/s)0.70.80.91的反变换。
[MATLAB程序]：% (1)Z transform using ztrans function
f=sym('cos(a*k)'); F=ztrans(f)
% (2)Inverse Z transform using iztrans function F=sym('1/(1+z)^2'); f=iztrans(F)
[运行结果]：
F =
(z-cos(a))*z/(z^2-2*z*cos(a)+1) f =
charfcn[0](n)-(-1)^n+(-1)^n*n
（二）练 习 题
1、用reiduez函数，求出下列各式的部分分式展开式和f(k)
(1) F?z?? (2) F?z??
2z?16z?44z?56z?323z?3z?15z?18z?12432432432432
4z?8.68z?17.98z?26.74z?8.04z?2z?10z?6z?6526
2、已知离散时间系统的差分方程为
2y(k)?y(k?1)?3y(k?2)?2f(k)?f(k?1)
kf(k)?0.5?(k),y(?1)?1,y(?2)?3，试用filter函数求系统的零输入响应、零状态响
应和完全响应
3、利用filter 函数求出系统H?z??断系统是否稳定。
4、用zplane（num,den）函数，画出H?z??2z?16z?44z?56z?323z?3z?15z?18z?1243243211?1.85z?1?0.85z?2的单位脉冲响应，并由h?k?判
的零极点图，
并判断系统的稳定性。
5、已知离散时间LTI系统的单位脉冲响应h?k??ak????k????k?N???,a?0，求系统函数，画出系统的零极点分布图、幅度响应曲线和相位响应曲线。
四、要求
学生对实验练习题编写MATLAB程序并运行，在计算机上输出仿真结果。
27
选做任务（一）
一、名称：滤波器的设计和应用 二、内容：
设计低通、带通、高通三种巴特沃斯滤波器。观察所设计滤波器的幅频、相频曲线。用
同一段音频信号通过不同的滤波器。画出滤波后相应的频谱图，并重放滤波后的声音信号。试听滤波结果，并分析滤波器的作用。 参考结果如图1-1和1-2所示：
图1-1
28
图1-2
三、相关知识：
滤波器是一个以特定方式设计出来的系统，用来改变输入信号的频谱分布。信号通过滤波器后会造成频率分布的变化。对于音频信号来说，频率分布直接影响信号重放时的音色和听感。通过对同一段音频素材不同滤波结果的考察，可以更加直观的理解滤波器的作用。
可以使用巴特沃斯滤波器来完成上述滤波功能。设计巴特沃斯数字低通滤波器可以使用如下语句：
[b,a]=butter(n,wc);
其中[b,a]是描述低通滤波器系统函数的两个行向量：即b表示系统函数分子多项式的系数，a表示分母多项式的系数；n为巴特沃斯滤波器的阶数，滤波器的阶数越大，过渡带越窄；wc为数字滤波器的归一化3dB截止频率（即信号幅度下降至1/2处的频率），其数值在[0,1]之间，取1对应0.5倍的采样频率。
当wc=[w1,w2]时，是设计2n阶的带通滤波器，3dB通带为w1 < w < w2。
[b,a]=butter(n,wc,’ftype’);
用于设计截止频率为wc的巴特沃斯高通或带阻滤波器。当ftype=high时，为n阶高通滤波器；ftype=stop，wc=[w1,w2]时，为2n阶的带阻滤波器。
如果在这个函数输入变量的最后，加一个变量‘s’，比如：[b,a]=butter(n,wc,’s’)。表示设计的是模拟滤波器，详细的使用方法，可以用help butter阅读。
29
除了巴特沃思滤波器（butter），IIR滤波器还有切比雪夫I（cheby1）、切比雪夫II(cheby2)和椭圆滤波器(ellip)，可以使用help学习他们的使用方法。
滤波功能可以使用filter函数实现： c=filter(b,a,y);
其中，[b,a]对应所使用滤波器的系统函数的分子和分母多项式的系数，比如，上面设计巴特沃思滤波器所得到的系数[b,a]。y是需要被滤波的信号。c是滤波后的结果信号。
傅里叶变换可以使用fft函数实现。调用格式是 Y=fft(y,n);
其中，y是需要观测频谱的时域离散信号，n是变换点数，当y的点数小于n时，y后面补n-y个零。Y是变换得到的频谱。离散信号的频谱是周期的，周期为2pi。使用函数fft得到的频谱范围是[0~2pi]，当需要观看[-pi ~ pi]的频谱时，可以使用函数fftshift来转换实现。调用方法是：yy=fftshift(Y);
其中，Y是fft得到的[0~2pi]范围的频谱，yy是转换得到的[-pi ~ pi]的频谱。
音频文件的读取可以使用wavread函数实现，其调用方法为： [y,Fs,bits]=wavread('filename');
其中，返回值y是音频数据，Fs为采样音频数据时使用的采样频率，bits为每个采样样本的比特数。
播放音频数据可以使用sound函数，调用方法为：
sound(y,Fs,bits); 注意采样频率Fs要设置正确，否则音乐会变调。 四、相关函数
[b,a]=butter(n,wc); c=filter(b,a,y); [H,w]=freqz(b,a);
Y=fft(y,n);
yy=fftshift(Y);
[y,Fs,bits]=wavread('filename');
30
任务（二）
一、名称：信号的幅度调制 二、内容
1、仿真幅度调制过程：将一采样信号Sa(t-10)对一个复合信号
s=cos(60*pi*t)+cos(40*pi*t)+cos(20*pi*t)进行幅度调制。后用高通、低通、带通三种不同的滤波器进行滤波，并绘出结果图像。
2、对上述方法得到的调制信号进行解调，恢复出原始的采样信号（解调过程就是：将调制信号和载波再次相乘后经过低通滤波的过程）。
任务1的仿真结果如图2-1和2-2所示（载波信号为周期方波）：
图2-1
31
图2-2
三、相关知识：
幅度调制是指载波幅度随调制信号变化的过程。本实验利用采样信号（Sa(t)）为调制信号，脉冲波（周期方波）信号或正弦信号为载波信号。调制的过程即是两个信号相乘的过程。
脉冲波信号就是周期方波信号，可以使用函数square实现。该函数和正弦、余弦函数很类似。比如：P = square(2*pi*f0*t, duty)
产生基频为f0，即周期为1/f0，占空比为duty的方波信号。
设计巴特沃斯数字滤波器的可以使用函数butter实现；滤波可以使用函数filter实现；
具体使用细节见选做任务（一）。
32