在MATLAB中，对于连续LTI系统的冲激响应和阶跃响应，可分别用控制系统工具箱提供的函数impluse和step来求解。其调用格式为 y=impluse(sys,t) y=step(sys,t)
式中，t表示计算系统响应的抽样点向量，sys是LTI系统模型。 例3-2已知某LTI系统的微分方程为 y’’(t)+ 2y’(t)+100y(t)=10f(t)
求系统的冲激响应和阶跃响应的波形. 解：ts=0;te=5;dt=0.01; sys=tf([10],[1,2,100]);
t=ts:dt:te;
h=impulse(sys,t); figure; plot(t,h);
xlabel('Time(sec)'); ylabel('h(t)');
g=step(sys,t); figure; plot(t,g);
xlabel('Time(sec)');
ylabel('g(t)');
3. 用MATLAB实现连续时间信号的卷积
信号的卷积运算有符号算法和数值算法，此处采用数值计算法，需调用MATLAB 的conv( )函数近似计算信号的卷积积分。连续信号的卷积积分定义是 f(t)?f1(t)?f2(t)?????f1(?)f2(t??)d?
如果对连续信号f1(t)和f2(t)进行等时间间隔?均匀抽样，则f1(t)和f2(t)分别变为离散时间信号f1(m?)和f2(m?)。其中，m为整数。当?足够小时，f1(m?)和f2(m?)既为连续时间信号f1(t)和f2(t)。因此连续时间信号卷积积分可表示为
f(t)?f1(t)?f2(t)??f1(?)f2(t??)d?????lim??0m????f(m?)?f1?
2(t?m?)??采用数值计算时，只求当t?n?时卷积积分f(t)的值f(n?)，其中，n为整数，既
f(n?)?
m?????f(m?)?f1m????2(n??m?)??
???f1(m?)?f2[(n?m)?]9
其中，
m????f(m?)?f1?2[(n?m)?]实际就是离散序列f1(m?)和f2(m?)的卷积和。当?足
够小时，序列f(n?)就是连续信号f(t)的数值近似，既 f(t)?f(n?)??[f1(n)?f2(n)]
上式表明，连续信号f1(t)和f2(t)的卷积，可用各自抽样后的离散时间序列的卷积再乘以抽样间隔?。抽样间隔?越小，误差越小。
例3-3用数值计算法求f1(t)?u(t)?u(t?2)与f2(t)?e?3tu(t)的卷积积分。
解：因为f2(t)?e?3tu(t)是一个持续时间无限长的信号，而计算机数值计算不可能计算真正的无限长信号，所以在进行f2(t)的抽样离散化时，所取的时间范围让f2(t)衰减到足够小就可以了，本例取t?2.5。程序是 dt=0.01; t=-1:dt:2.5;
f1=Heaviside(t)-Heaviside(t-2); f2=exp(-3*t).*Heaviside(t);
f=conv(f1,f2)*dt; n=length(f); tt=(0:n-1)*dt-2; subplot(221), plot(t,f1), grid on;
axis([-1,2.5,-0.2,1.2]); title('f1(t)'); xlabel('t') subplot(222), plot(t,f2), grid on;
axis([-1,2.5,-0.2,1.2]); title('f2(t)'); xlabel('t') subplot(212), plot(tt,f), grid on; title('f(t)=f1(t)*f2(t)'); xlabel('t')
由于f1(t)和f2(t)的时间范围都是从t=-1开始，所以卷积结果的时间范围从 t=-2开始，增量还是取样间隔?，这就是语句tt=(0:n-1)*dt-2的由来。 三、上机实验内容
1． 调试实验原理中所述的相关程序。
2． 已知描述系统的微分方程和激励信号f(t)如下，试用解析法求系统的零状态响应y(t)，并用MATLAB绘出系统零状态响应的时域仿真波形，验证结果是否相同 y’’(t)+ 4y’(t)+4y(t)=f’(t)+3f(t) f(t)= exp(-t)u(t)
3．已知描述系统的微分方程如下，试用MATLAB求系统在0~10秒范围内冲激响应和阶跃响应的数值解，并用绘出系统冲激响应和阶跃响应的时域波形
y’’(t)+3y’(t)+2y(t)=f(t) y’’(t)+ 2y’(t)+2y(t)=f’(t) 4．画出信号卷积积分f1(t)?f2(t)的波形，f1(t)?f2(t)?u(t)?u(t?1)
10
实验三 傅里叶变换、系统的频域分析
（设计型）
一、 实验目的
1、学会用MATLAB实现连续时间信号傅里叶变换 2、学会用MATLAB分析LTI系统的频域特性 3、学会用MATLAB分析LTI系统的输出响应 二、实验原理
1．傅里叶变换的MATLAB求解
MATLAB的symbolic Math Toolbox 提供了直接求解傅里叶变换及逆变换的函数fourier()及ifourier()两者的调用格式如下。 Fourier 变换的调用格式
F=fourier(f)：它是符号函数f的fourier变换默认返回是关于w的函数。
F=fourier(f，v):它返回函数F是关于符号对象v的函数，而不是默认的w，即
F(v)??????f(x)e?jvxdx
Fourier逆变换的调用格式
f=ifourier(F):它是符号函数F的fourier逆变换，默认的独立变量为w，默认返回是
关于x的函数。
f=ifourier(f,u):它的返回函数f是u的函数，而不是默认的x.
注意：在调用函数fourier()及ifourier()之前，要用syms命令对所用到的变量（如t,u,v,w）进行说明,即将这些变量说明成符号变量。 例4-1 求f(t)?e?2t的傅立叶变换
解: 可用MATLAB解决上述问题： syms t
Fw=fourier(exp(-2*abs(t)))
例4-2 求F(jw)?1的逆变换f(t) 1??2解： 可用MATLAB解决上述问题 syms t w
ft=ifourier(1/(1+w^2),t)
2．连续时间信号的频谱图
例4-3 求调制信号f(t)?AG?(t)cos?0t的频谱，式中
A?4,?0?12?,??1??,G?(t)?u(t?)?u(t?) 222解：MATLAB程序如下所示
ft=sym('4*cos(2*pi*6*t)*(Heaviside(t+1/4)-Heaviside(t-1/4))');
11
Fw=simplify(fourier(ft)) subplot(121)
ezplot(ft,[-0.5 0.5]),grid on subplot(122)
ezplot(abs(Fw),[-24*pi 24*pi]),grid
用MATLAB符号算法求傅里叶变换有一定局限，当信号不能用解析式表达时，会提示出错，这时用MATLAB的数值计算也可以求连续信号的傅里叶变换，计算原理是
F(j?)??f(t)e????j?tdt?lim?f(n?)e?j?n??
??0n????当?足够小时，近似计算可满足要求。若信号是时限的，或当时间大于某个给定值时，信号已衰减的很厉害，可以近似地看成时限信号时，n的取值就是有限的，设为N，有
F(k)??f(n?)e?j?kn??n?0N?1,0?k?N,?k?2?k 是频率取样点 N?时间信号取样间隔?应小于奈奎斯特取样时间间隔，若不是带限信号可根据计算精度要求确定一个频率 W0为信号的带宽。
例4-4 用数值计算法求信号f(t)?u(t?1)?u(t?1)的傅里叶变换
解，信号频谱是F(j?)?2Sa(?)，第一个过零点是?，一般将此频率视为信号的带宽，若将精度提高到该值的50倍，既W0=50?，据此确定取样间隔，??R=0.02;t=-2:R:2;
f=Heaviside(t+1)-Heaviside(t-1); W1=2*pi*5;
N=500;k=0:N;W=k*W1/N; F=f*exp(-j*t'*W)*R; F=real(F);
W=[-fliplr(W),W(2:501)]; F=[fliplr(F),F(2:501)]; subplot(2,1,1);plot(t,f); xlabel('t');ylabel('f(t)'); title('f(t)=u(t+1)-u(t-1)'); subplot(2,1,2);plot(W,F); xlabel('w');ylabel('F(w)'); title('f(t)的付氏变换F(w)');
3．用MATLAB分析LTI系统的频率特性
当系统的频率响应H（jw）是jw的有理多项式时，有
1?0.02 2F0jwM)?1?L?b1jw(?)b0B(w)bM(jw)M?bM?1( H(jw)? ?NN?1A(w)aN(jw)?aN?1(jw)?L?a1(jw)?a012