信号与系统实验四：信号与系统复频域分析（2 学时）

81 阅读 0 评论 54 点赞

我是靠谱客的博主悦耳枕头，最近开发中收集的这篇文章主要介绍信号与系统实验四：信号与系统复频域分析（2 学时），觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

注明：这是学校专业课实验，利用网站做一个学习和记录。有错的地方希望好心人指点一下。

一、实验目的

1.学会用 MATLAB 进行部分分式展开；

2.学会用 MATLAB 分析 LTI 系统的特性；

3.学会用 MATLAB 进行 Laplace 正、反变换。

4.学会用 MATLAB 画离散系统零极点图；

5.学会用 MATLAB 分析离散系统的频率特性；

二、实验原理

1．用 MATLAB 进行部分分式展开

用 MATLAB 函数 residue 可以得到复杂有理分式 F(s)的部分分式展开式，其调用格式为

$left [ r,p,k right ]= residue(num,den)$

其中，num,den 分别为 F(s)的分子和分母多项式的系数向量，r 为部分分式的系数，p 为极点，k 为 F(s)中整式部分的系数，若 F(s)为有理真分式，则 k 为零。

例 6-1 用部分分式展开法求 F(s)的反变换 $F(s)=frac{s+2}{s^{3}+4s^{2}+3s}$

解:其 MATLAB 程序为

format rat;

num=[1,2];

den=[1,4,3,0];

[r,p]=residue(num,den)

程序中 format rat 是将结果数据以分数形式显示

F(s)可展开为 $F(s)=frac{2/3}{s}+frac{-0.5}{s+1}+frac{-1/6}{s+3}$

所以，F(s)的反变换 $f(t)=left [ frac{2}{3} -frac{1}{2}e^{-t}-frac{1}{6}e^{-3t}right ]u(t)$

2．用 MATLAB 分析 LTI 系统的特性

系统函数 H（s）通常是一个有理分式，其分子和分母均为多项式。计算 H（s）的零极点可以应用 MATLAB 中的 roots 函数，求出分子和分母多项式的根，然后用 plot 命令画图。

在 MATLAB 中还有一种更简便的方法画系统函数 H（s）的零极点分布图，即用 pzmap

函数画图。其调用格式为

pzmap(sys)

sys 表示 LTI 系统的模型，要借助 tf 函数获得，其调用格式为

sys=tf(b,a)

式中，b 和 a 分别为系统函数 H（s）的分子和分母多项式的系数向量。

如果已知系统函数H(s)，求系统的单位冲激响应 $H(t)$ 和频率响应 $H(jomega )$ 可以用以前介绍过的 impulse 和 freqs 函数。

例 6-2 已知系统函数为 $H(s)=frac{1}{s^{3}+2s^{2}+2s+1}$

试画出其零极点分布图，求系统的单位冲激响应 h(t)和频率响应 $H(jomega )$ ，并判断系统

是否稳定。

解：其 MATLAB 程序如下：

num=[1];

den=[1,2,2,1];

sys=tf(num,den);

figure(1);pzmap(sys);

t=0:0.02:10;

h=impulse(num,den,t);

figure(2);plot(t,h)

title('Impulse Response')

[H,w]=freqs(num,den);

figure(3);plot(w,abs(H))

xlabel('omega')

title('Magnitude Response')

3．用 MATLAB 进行 Laplace 正、反变换

MATLAB 的符号数学工具箱提供了计算 Laplace 正、反变换的函数 Laplace 和 ilaplace, 其调用格式为上述两式右端的 f 和 F 分别为时域表示式和 s 域表示式的符号表示，可以应用函数 sym 实现，其调用格式为

S=sym(A)

式中，A 为待分析表示式的字符串，S 为符号数字或变量。

$F=laplace(f)$

$f=ilaplace(f))$

例 6-3 试分别用 Laplace 和 ilaplace 函数求

（1） $f(t)=e^{-t}sin(at)u(t)$ 的 Laplace 变换；

（2） $F(s)=frac{s^{2}}{s^{2}+1}$ 的 Laplace 反变换。

解：

（1）其程序为

%（1）其程序为 
f=sym('exp(-t)*sin(a*t)'); 
F=laplace(f) 
%或者
syms a t 
F=laplace(exp(-t)*sin(a*t)) 
%（2）其程序为
F=sym('s^2/(s^2+1)'); 
ft=ilaplace(F) 
%或 
syms s 
ft= ilaplace(s^2/(s^2+1))

4．离散系统零极点图

离散系统可以用下述差分方程描述： $sum_{i=0}^{N}a_{i}y(k-i)=sum_{m=0}^{M}b_{m}f(k-m)$

Z变换后可得系统函数： $H(z)=frac{Y(z)}{F(z)}=frac{b_{0}+b_{1}z^{-1}+...+b_{M}z^{-M}}{a_{0}+a_{1}z^{-1}+...+a_{N}z^{-N}}$

用 MATLAB 提供的 root 函数可分别求零点和极点，调用格式是

p=[a0,a1…an],q=[b0,b1…bm,0,0…0], 补 0 使二者维数一样。画零极点图的方法有多种，可以用 MATLAB 函数[z,p,k]=tf2zp(b,a)和 zplane(q,p)，也可用 plot 命令自编一函数 ljdt.m,话图时调用。

%注意：建立函数时以lidt.m为文件名，切函数有参数不足报错不用细究
function ljdt(A,B)

% The function to draw the pole-zero diagram for discrete system

p=roots(A)； %求系统极点

q=roots(B)； %求系统零点

p=p'； %将极点列向量转置为行向量

q=q'；%将零点列向量转置为行向量

x=max(abs([p q 1]))； %确定纵坐标范围

x=x+0.1；

y=x； %确定横坐标范围

clf

hold on

axis([-x x -y y]) %确定坐标轴显示范围

w=0:pi/300:2*pi；

t=exp(i*w)；

plot(t) %画单位园

axis('square')

plot([-x x],[0 0])  %画横坐标轴

plot([0 0],[-y y]) %画纵坐标轴

text(0.1,x,'jIm[z]')

text(y,1/10,'Re[z]')

plot(real(p),imag(p),'x') %画极点

plot(real(q),imag(q),'o') %画零点

title('pole-zero diagram for discrete system') %标注标题

hold off

end

例 6-4 求系统函数零极点图 $H(z)=frac{z+1}{3z^{5}-z^{4}+1}$

%在建立好上述函数的情况下，另起脚本，运行代码
a=[3 -1 0 0 0 1];

b=[1 1];

ljdt(a,b)

p=roots(a)

q=roots(b)

pa=abs(p)

5．离散系统的频率特性

离散系统的频率特性可由系统函数求出，既令 $z=e^{jomega }$ ,MATLAB 函数 freqz 可计算频率特性，调用格式是：

[H，W]=freqz(b,a,n)，

b 和 a 是系统函数分子分母系数，n是 $0-pi$ 范围 n 个等份点，默认值 512，H 是频率响应函数值，W 是相应频率点； [H，W]=freqz(b,a,n,’whole’), n 是 $0-2pi$ 范围 n 个等份点freqz(b,a,n)，直接画频率响应幅频和相频曲线；

例 6-5 系统函数 $H(z)frac{z-0.5}{z}$

运行如下语句，可得 10 个频率点的计算结果

A=[1 0];

B=[1 -0.5];

[H,W]=freqz(B,A,10)

继续运行如下语句，可将 400 个频率点的计算结果用 plot 语句画幅频和相频曲线

B=[1 -0.5];

A =[1 0];

[H,w]=freqz(B,A,400,'whole');

Hf=abs(H);

Hx=angle(H);

clf

figure(1)

plot(w,Hf)

title('离散系统幅频特性曲线')

figure(2)

plot(w,Hx)

title('离散系统相频特性曲线')

%还可用 freqz 语句直接画图，注意区别

A=[1 0];

B=[1 -0.5];

[H,W]=freqz(B,A,10)

继续运行如下语句，可将 400 个频率点的计算结果用 plot 语句画幅频和相频曲线

B=[1 -0.5];

A =[1 0];

[H,w]=freqz(B,A,400,'whole');

Hf=abs(H);

Hx=angle(H);

clf

figure(1)

plot(w,Hf)

title('离散系统幅频特性曲线')

figure(2)

plot(w,Hx)

title('离散系统相频特性曲线')

%还可用 freqz 语句直接画图，注意区别

%A=[1 0];

%B=[1 -0.5];

%freqz(B,A,400)

例 6-6 用几何矢量法，自编程序画频率响应

原理：频率响应 $H(e^{jomega })=frac{prod_{j-1}^{M}e^{jomega }-q^{j}}{prod_{i=1}^{N}e^{jomega }-p^{i}}$

编程流程：定义 Z 平面单位圆上 k 个频率等分点；求出系统函数所有零点和极点到

这些等分点的距离；求出系统函数所有零点和极点到这些等分点的矢量的相角；求出单位圆上各频率等分点的

$left | H(e^{jomega}) right |$ 和 $varphi (omega )$

画指定范围内的幅频与相频。若要画零极点图，可调用 ljdt.m 函数。

%与例4相同，建立一个函数以dplxy.m 命名，
function dplxy(k,r,A,B)

%The function to draw the frequency response of discrete system

p=roots(A); %求极点

q=roots(B); %求零点

figure(1)

ljdt(A,B) %画零极点图

w=0:l*pi/k:r*pi;

y=exp(i*w);%定义单位圆上的 k 个频率等分点

N=length(p);%求极点个数

M=length(q);%求零点个数

yp=ones(N,1)*y;%定义行数为极点个数的单位圆向量

yq=ones(M,1)*y;%定义行数为零点个数的单位圆向量

vp=yp-p*ones(1,r*k+1);%定义极点到单位圆上各点的向量

vq=yq-q*ones(1,r*k+1);%定义零点到单位圆上各点的向量

Ai=abs(vp);%求出极点到单位圆上各点的向量的模

Bj=abs(vq);%求出零点到单位圆上各点的向量的模

Ci=angle(vp);%求出极点到单位圆上各点的向量的相角

Dj=angle(vq);%求出零点到单位圆上各点的向量的相角

fai=sum(Dj,1)-sum(Ci,1); %求系统相频响应

H=prod(Bj,1)./prod(Ai,1); %求系统幅频响应

figure(2)

plot(w,H); %绘制幅频特性曲线

title('离散系统幅频特性曲线')

xlabel('角频率')

ylabel('幅度')

figure(3)

plot(w,fai)

title('离散系统的相频特性曲线')

xlabel('角频率')

ylabel('相位')

end

已知系统函数 $H(z)=frac{5/4(1-z^{-1})}{1-1/4z^{-1}}$ ，画频率响应和零极点图。

%同例4，在上述dplxy.m函数建好情况下，另建脚本运行
A=[1 -1/4];

B=[5/4 -5/4];

dplxy(500,2,A,B) %绘制系统 2π 频率范围内 500 个频率点的幅频和相频特性曲线

三。上机实验

1.验证实验原理中的相关程序。

2 ．求信号 $f(t)=te^{-3t}u(t)$ 的拉普拉斯变换

%上机实验 f=t*exp(-3*t)u(t)的拉普拉斯变化
f=str2sym('t*exp(-3*t)'); 
F=laplace(f) 
%或
syms  t 
F=laplace(t*exp(-3*t))

3 ．求函数 $F(s)=frac{s^{3}+5s^{2}+9s+7}{s^{2}+3s+2}$ 的反变换

%求F（s）=s^3+s^2+9*s+7/(s^2+3*s+2)的反变化
%反变换:
F=str2sym('(s^3+5*s^2+9*s+7)/(s^2+3*s+2)'); 
ft=ilaplace(F)

4 ．已知连续系统的系统函数如下，试用 MATLAB 绘制系统的零极点图，并根据零极点图

判断系统的稳定性

$H(s)=frac{s^{2}+s+2}{3s^{3}+5s^{2}+4s-6}$

% 已知连续系统的系统函数如下，试用 MATLAB 绘制系统的零极点图，并根据零极点图，判断系统的稳定性
%H(s)=s^2+s+2/(3*s^3+5*s^2+4*s-6)

num=[1,1,2];%用向 量表示 分子系数
den=[3,5,4,-6];%用向量表示分母系数
sys=tf(num,den);%系统传递函数
set(gcf,'color','w');
pzmap(sys);%绘制其零极点图
p=pole(sys);
z=zero(sys);%零极点图： X为极点， O为零点
%由零极点图看出：系统不稳定，因为有极点在右半平面

5 ．系统函数是 $1+5z^{-1}+5z^{-2}+z^{-3}$ 求频率响应

第一种方法：

（总的求他的频率响应图，有问题请指正）

num=[1,5,5,1]; 
den=[1];
sys=tf(num,den); 
figure(1);pzmap(sys); 
[H,w]=freqs(num,den); 
figure(2);plot(w,abs(H)) 
xlabel('omega') 
title('频率响应')

第二种方法：

（分别求频率响应的幅频和相频）

num=[1,5,5,1]; 
den=[1,0,0,0];
[H,w]=freqz(num,den,400,'whole'); 
Hf=abs(H);
Hx=angle(H);
clf
figure(1) 
plot(w,Hf)
title('幅频特性') 
figure(2)
plot(w,Hx)
title('相频特性')