⑵ 连续时间系统冲激响应和阶跃响应的求解
在MATLAB中, 求解系统冲激响应可应用控制系统工具箱提供的函数impulse，求解阶跃响应可利用函数step。其调用形式为：
y?impulse(sys,t)y?step(sys,t)
式中，t表示计算系统响应的抽样点向量，sys是LTI系统模型。
【例14】 在例13所述力学系统中，若外力f(t)是强度为10的冲激信号，求物体的位移y(t)。
解：由已知条件，系统的输入信号为f(t)?10?(t)，系统的微分方程为：
d2y(t)dy(t)?2?100y(t)?10?(t) dt2dt计算物体位移y(t)的MATLAB源程序如下：
% program exa_14.m, Impulse response of LTI system
ts=0;te=5;dt=0.01;
sys=tf([10],[1 2 100]); t=ts:dt:te;
y=impulse(sys,t); plot(t,y)
xlabel('Time(sec)'); ylabel('y(t)') 程序运行结果如图所示：
0.210.80.60.4
y(t)0-0.2-0.4-0.6-0.800.511.522.5Time(sec)33.544.5515
⑶ 离散时间系统零状态响应的求解
大多LTI离散时间系统都可用如下线性常系数差分方程描述：
?ay[k?i]??bii?0j?0NMjf[k?j]
其中，f[k],y[k]分别表示系统的输入和输出，N是差分方程阶数。已知差分方程的N个初始状态和输入f[k]，就可以通过编程由下式迭代计算出系统的输出：
y[k]???(ai/a0)y[k?i]??(bj/a0)f[k?j]
i?1j?0NM在零初始状态下，MATLAB信号处理工具箱提供了一个filter函数,可以计算由差分方程描述的系统的响应。其调用形式为：
y?filter(b,a,f)
式中，b?[b0,b1,b2,,bM],a?[a0,a1,a2,,aN] 分别是差分方程左、右端的系数
向量，f表示输入序列，y表示输出序列。注意，输出序列的长度与输入序列的长度相同。
【例15】 受噪声干扰的信号为f[k]?s[k]?d[k]，其中s[k]?(2k)0.9k是原始信号，d[k]是噪声。已知M点滑动平均（moving average）系统的输入输出关系为：
1y[k]?M?f[k-n]
n?0M-1试利用MATLAB编程实现用M点滑动平均系统对受噪声干扰的信号去噪。 解：系统的输入信号f[k]含有有用信号s[k]和噪声信号d[k]。噪声信号d[k]可以用rand函数产生，将其叠加在有用信号s[k]上，即得到受噪声干扰的输入信号
f[k]。对信号f[k]去噪的源程序如下（取M?5）：
% program exa_15.m, 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;
16
k=0:R-1;
s=2*k.*(0.9.^k); f=s+d;
figure(1);
plot(k,d,'r-.',k,s,'b--',k,f,'g-'); 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'); legend('s[k]','y[k]'); 程序运行结果如图所示：
876543210-1d[k]s[k]f[k]051015202530Time index k35404550876543210s[k]y[k]051015202530Time index k3540455017
⑷ 离散时间系统单位脉冲响应的求解
在MATLAB中, 求解离散时间系统单位脉冲响应,可应用信号处理工具箱提供的函数impz，其调用形式为：
h?impz(b,a,k)
式中，b?[b0,b1,b2,,bM],a?[a0,a1,a2,,aN] 分别是差分方程左、右端的系数
向量，k表示输出序列的取值范围，h就是系统的单位脉冲响应。 【例16】 用impz函数求离散时间系统：
y[k]?3y[k?1]?2y[k?2]?f[k]
的单位脉冲响应h[k]，并与理论值 h[k]??(?1)k?2(?2)k,k?0 进行比较。 解：源程序如下：
% program exa_16.m, Impulse response of discrete system k=0:10
a=[1 3 2]; b=[1];
h=impz(b,a,k); subplot(211) stem(k,h)
title('单位脉冲响应的近似值') hk=-(-1).^k+2*(-2).^k; subplot(212) stem(k,hk)
axis(0 10 -2000 3000]) title('单位脉冲响应的理论值') 程序运行结果如图所示：
18
单位脉冲响应的近似值3000200010000-1000-2000012345678910单位脉冲响应的理论值400020000-2000012345678910⑸ 离散卷积的计算
卷积是用来计算系统零状态响应的有力工具。MATLAB信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv,其调用形式为：
c?conv(a,b)
式中，a、b为待卷积两序列的向量表示，c是卷积结果。向量c的长度为向量a、
b的长度之和减1，即 length(c)=length(a)+length(b)-1。
k?【例17】 已知序列x[k]?[1,2,3,4;0,1,2,3],yk?[][1,1,1,1,1;k?,0,1,4]，
利用MATLAB编程计算x[k]?y[k]并画出卷积结果。 解：源程序如下：
% program exa_17.m, Linear Convolution
x=[1,2,3,4]; y=[1,1,1,1,1]; z=conv(x,y) N=length(z); stem(0:N-1,z) 程序运行结果：
19