系统的复频域分析

实验六 系统的复频域分析

§6.1 MATLAB 函数lsim （用于系统函数） 目的

用lsim 仿真由系统函数表征的因果LTI 系统的输出。 基本题

1．定义系数向量a1和b1用以描述由下面系统函数表征的因果LTI 系统： 2

2

)(1+-=

s s s H 2．定义系数向量a2和b2用以描述由下面系统函数表征的因果LTI 系统： 3

.03

)(2+=

s s H 3．定义系数向量a3和b3用以描述由下面系统函数表征的因果LTI 系统： 8

.02)(3+=

s s

s H 4．利用lsim 和前面部分定义的向量求这些因果LTI 系统对由t=[0:0.1:0.5]，x=cos(t)给出的输入的输出。 以上四题解： t=[0:0.1:0.5]; x=cos(t); b1=[1 -2]; a1=[1 2]; b2=3; a2=[1 0.3] b3=2 a3=[1 0.8] y1=lsim(b1,a1,x,t); subplot(2,2,1) plot(t,y1); y2=lsim(b2,a2,x,t); subplot(2,2,2)

plot(t,y2); y3=lsim(b3,a3,x,t); subplot(2,2,3) plot(t,y3);

§6.2作连续时间的零极点图 目的

这一节要学习如何在一个零极点图上展现有理系统函数的零极点。 基本题

1．下列每个系统函数都对应于稳定的LTI 系统。用roots 求每个系统函数的零极点，如上所示的利用plot 画出零极点图并作适当标注。 (i) 3

25)(2+++=s s s s H (ii) 10

21252)(22++++=s s s s s H

(iii) )

2)(102(1252)(2

2+++++=

s s s s s s H

若干不同的信号能有相同的拉普拉斯变换有理表达式，但有不同的收敛域。例如，具有单位冲激响应为)()(),()(t u e t h t u e t h at ac at c --==--的因果和反因果LTI 系统就有相同的分子和分母多项式的有理系统函数为 a s a

s s H c ->+=

)R e ( 1

)( a s a

s s H ac -<+=

)Re( 1

)( 然而，它们有不同的系统函数，因为 它们有不同的收敛域。 解：代码如下 b1=[1 5]; a1=[1 2 3]; zs1=roots(b1) ps1=roots(a1) subplot(3,2,1)

plot(real(zs1),imag(zs1),'o'); grid

axis([-10 10 -10 10]) subplot(3,2,2)

plot(real(ps1),imag(ps1),'x'); grid

axis([-10 10 -10 10])

b2=[2 5 12]; a1=[1 2 10]; zs2=roots(b2) ps2=roots(a2) subplot(3,2,3)

plot(real(zs2),imag(zs2),'o'); grid

axis([-10 10 -10 10]) subplot(3,2,4)

plot(real(ps2),imag(ps2),'x'); grid

axis([-10 10 -10 10])

b3=[1 5];

a3=[1 2 3];

zs3=roots(b3)

ps3=roots(a3)

subplot(3,2,5)

plot(real(zs3),imag(zs3),'o'); grid

axis([-10 10 -10 10]) subplot(3,2,6)

plot(real(ps3),imag(ps3),'x'); grid

axis([-10 10 -10 10])

2．对1中每个有理表达式，确定它们的收敛域。 解：代码如下

i.Re(s)>-1;ii. Re(s)为任意值 i.Re(s)>-2 3．对输入和输出满足下面微分方程：

)(5)

(2)()(3)(2

2t x dt t dx dt t x d t y dt t dy ++=- 的因果LTI 系统，求系统函数的零点和极点，并完成一幅适当标注的零极点图。 解：代码如下 b=[1 -3]; a=[1 2 5]; zs=roots(b) ps=roots(a) subplot(2,1,1)

plot(real(zs),imag(zs),'o'); grid

axis([-10 10 -10 10]) subplot(2,1,2)

plot(real(ps),imag(ps),'x'); grid

axis([-10 10 -10 10])

中等题

下面的练习需要用到函数plotpz，其具体实现如下。函数plotpz是用于对分子和分母多项式的系数存入向量a和b中的LTI系统画出零极点图用的。这个函数除了画图以外，还产生零极点的值。

function [ps,zs]=plotpz(b,a,ROC)

ps=roots(a);

zs=roots(b);

ps=ps(:);

zs=zs(:);

Maxl=max(abs(imag([ps;zs;j])));

MaxR=max(abs(imag([ps;zs;i])));

plot(1.5*[-MaxR MaxR],[0 0],'w');

hold on

text(1.5*MaxR,0,'Re');

plot([0 0],1.5*[-Maxl Maxl],'w');

text(0,1.5*Maxl,'Im');

plot(real(zs),imag(zs),'bo');

plot(real(ps),imag(ps),'bx');

if (nargin('plotpz')>2)%nargin是用来判断输入变量个数的函数

if(any(real(ps)

ipole=max(real(ps(real(ps)

plot([ipole ipole],1.5*[-Maxl Maxl],'r--');

end

if(any(real(ps)>ROC))

rpole=min(real(ps(real(ps)>ROC)));

plot([rpole rpole],1.5*[-Maxl Maxl],'r--');

end

text(ROC,-1.25*Maxl,'ROC');

end

axis('equal');

grid;

hold on

任选宗量ROC用来在图中指出收敛域。在该系统的收敛域内用选取ROC是一个点，plotpz将适当地给系统的收敛域标注。例如，试试执行下列命令

>> b=[1 -1];

>> a=[1 3 1];

>> [ps,zs]=plotpz(b,a,1);

>> [ps,zs]=plotpz(b,a,-2);

b=[1 -1];

a=[1 3 1];

[ps,zs]=plotpz(b,a,1);

[ps,zs]=plotpz(b,a,-2);

4．说明plotpz如何在ROC内画出一个单一的点确定一个有理变换的收敛域。解：通过any(real(ps)ROC)))函数找出边界，然后用plot画出收敛域

§6.3 MATLAB 函数freqz 目的

学习用函数freqz 绘制LTI 系统的频率响应。 基本题

1．利用freqz 画出微分方程)(5)()()(4)(3

2

22

2t x dt t x d t y dt t dy dt t y d +=++描述的一个因果、

稳定的LTI 系统的频率响应的幅值和相位。 解：代码如下： a=[3 4 1]; b=[1 0 5]; freqz(b,a);

§6.4 系统的时域和频域特性 目的

考虑几个由线性常系数微分方程描述的稳定的LTI 系统，对这些系统要求计算它们的单位冲激响应和频率响应。尽管用频率响应或单位脉冲响应都足以完全表征一个LTI 系统，但是将会明白有时候既从时域，又从频域来考虑系统特性是很有利的。 基本题

考虑由下面微分方程给出的一类因果LTI 系统

)()()

(00t x a t y a dt

t dy =+ 其中00>a 以保证稳定性。定义系统I 是满足上式，30=a 的系统，系统II 是3

10=a 的系统。

1．用解析法导出对应于上式的稳定LTI 系统的频率响应，同时确定这个频率响应的幅值和相位。 解：代码如下：

a1=[1 3];

b1=3; freqs(b1,a1) a2=[1 1/3]; b2=1/3; freqs(b2,a2)

2．定义w=linspace(0,10)，利用freqz计算系统I和系统II在w频率上的频率响应，在单一的图上画出这两个频率响应的幅值。这两个幅值图与解析表达式中频率响应的幅值一致吗？

解：代码如下：

w=linspace(0,10);

a=[1 3]

b=3;

freqz(b,a,w)

w=linspace(0,10);

a=[1 1/3];

b=1/3;

freqz(b,a,w)

结果分析：这两个幅值图与解析表达式中频率响应的幅值一致

3．用函数impulse计算系统I和系统II在向量t=linspace(0,5)所定义的时间样本点上的单位冲激响应。

解：

t=linspace(0,5);

b=3;

a=[1 3];

impulse(b,a,t);

t=linspace(0,5); b=1/3;

a=[1 1/3]; impulse(b,a,t);

4．单位冲激响应随时间衰减的速率与频率响应幅值随频率下降的速率之间是什么关系？CTFT中的何种性质说明这一关系？

解：单位冲激响应随时间衰减的速率与频率响应幅值随频率下降的速率之间是一对傅立叶变换关系

中等题

5．巴特沃兹滤波器就属于这样一类连续时间频率选择滤波器，它们能用线性常系数微分方程实现。为了确定一个二阶3

ω的巴特沃兹滤波器的系数，键入

=

c

>> wc=3;

>> [b2,a2]=butter(2,wc,'s');

宗量's'表明所要求的是一个连续时间滤波器而不是离散时间滤波器；如有缺损，则由butter产生的是一个离散时间滤波器。向量b2和a2包含的是二阶微分方程的系数，其格式与 2.3节所用相同。用freqz计算并画出相应LTI系统在w=linspace(0,10)频率点上频率响应的幅值。在一张图上要包括2中已计算出的系统I的频率响应的幅值。

解：代码如下：

wc=3;

[b2,a2]=butter(2,wc,'s');

w=linspace(0,10);

H2=freqs(b2,a2,w)

a1=[1 3];

b1=3;

H1=freqs(b1,a1,w)

plot(w,abs(H2),'r');

hold on

plot(w,abs(H1));

6．利用频率响应幅值，哪个系统更逼近于3=

ω的理想低通滤波器，是系统I还

c

是二阶巴特沃兹滤波器？将每一系统的频率响应相位与该理想低通滤波器的频率响应相位比较，情况怎样？

解：二阶巴特沃兹滤波器更逼近于3

ω的理想低通滤波器；将每一系统的频率

=

c

响应相位与该理想低通滤波器的频率响应相位比较，相位基本保持一个固定值不变

深入题

7．用impulse计算二阶巴特沃兹滤波器在时间样本t=linspace(0,5)上的单位冲激响应，在一张图上画出系统I和二阶巴特沃兹滤波器的单位冲激响应。

解：代码如下：

a1=[1 3];

b1=3;

t=linspace(0,5)

h1=impulse(b1,a1,t);

wc=3;

[b2,a2]=butter(2,wc,'s');

h2=impulse(b2,a2,t);

plot(t,h1,'r')

hold on;

plot(t,h2);

§6.5 二阶系统的极点位置

目的

考查二阶系统阻尼比和无阻尼自然频率对系统零极点分布以及频率响应的影响。

基本题

考查

ω保持为1，ζ取不同值时的频率响应和极点位置。

n

1．定义)(1s H ～)(4s H 分别是1=n ω，而1,41,0=ζ和2时，由上式确定的系统函数。定义向量a1～a4分别是)(1s H ～)(4s H 分母多项式的系数向量。求出并画出每个系统极点的位置。 解：代码如下： a1=[1 0 1]; a2=[1 1/2 1]; a3=[1 2 1]; a4=[1 4 1]; ps1=roots(a1); ps2=roots(a2); ps3=roots(a3); ps4=roots(a4); subplot(4,1,1);

plot(real(a1),imag(a1),'x'); subplot(4,1,2);

plot(real(a2),imag(a2),'x'); subplot(4,1,3);

plot(real(a3),imag(a3),'x'); subplot(4,1,4);

plot(real(a4),imag(a4),'x');

2．定义omega=[-5:0.1:5]是这些频率，在这些频率上要计算这4个系统的频率响应。用freqz 函数计算并画出在1中所定义的这4个系统的)(ωj H 。对于1<ζ和1≥ζ的频率响应定性来看是怎样的不同？你能解释系统的极点位置为何会导致这种差别？另外，从几何上你能论证为什么对所有这4种系统都有相同的0)(=ωωj H 值吗？

解：代码如下： a1=[1 0 1]; a2=[1 1/2 1]; a3=[1 2 1]; a4=[1 4 1]; b=1;

omega=[-5:0.1:5]; h1=abs(freqs(b,a1,omega)); h2=abs(freqs(b,a2,omega)); h3=abs(freqs(b,a3,omega)); h4=abs(freqs(b,a4,omega)); subplot(4,1,1); plot(h1'); subplot(4,1,2); plot(h2); subplot(4,1,3); plot(h3); subplot(4,1,4); plot(h4);

结果分析： 中等题

下面将要画出 2)(2

22n

n n

ac s s s H ωζωω++=

当改变ζ

和n ω时极点位置的轨迹，并

且看看这些参数变化是如何影响系统的频率响应的。

3．首先保持1=n ω，在100≤≤ζ范围内改变ζ。定义zetarange=[0 logspace(-1,1,99)]以得到在100≤≤ζ内100个按对数分开的点。定义azata 是一个3×100的矩阵，其中每一列是当ζ为zetarange 中对应列的值时，)(s H 分母多项式的系数。定义zetapoles 是2×100的矩阵，其中每一列就是azeta 对应列的根。在单独的图上，画出zetspoles 每一行实部对于虚部的图，并叙述画图的轨迹。在图上指出下列各点：1,1,0=ζ和2。为了在图中得到具有等长度轴的平方纵横比，可以键入 >> axis('equal'); >> axis([-4 0 -2 2]);

定性叙述一下，当ζ从0改变到1，然后又从1改变到10时，估计频率响应会如何变化。 解：代码如下：

wn=1;

N=100;

zetarange=[0 logspace(-1,1,99)];

azata=ones(3,100);

azata(2,:)=zetarange;

zetapoles=zeros(2,100);

for n=1:N;

zetapoles(:,n)=roots(azata(:,n));

end;

subplot(2,1,1)

plot(real(zetapoles(1,:)),imag(zetapoles(1,:))); axis([-4 0 -2 2]);

subplot(2,1,2);

plot(real(zetapoles(2,:)),imag(zetapoles(2,:))); axis([-4 0 -2 2]);