实验三 连续时间LTI系统的频域分析

一、实验目的

1、掌握系统频率响应特性的概念及其物理意义；

2、掌握系统频率响应特性的计算方法和特性曲线的绘制方法，理解具有不同频率响应特性的滤波器对信号的滤波作用；

3、学习和掌握幅度特性、相位特性以及群延时的物理意义； 4、掌握用MATLAB语言进行系统频响特性分析的方法。

基本要求：掌握LTI连续和离散时间系统的频域数学模型和频域数学模型的MATLAB描述方法，深刻理LTI系统的频率响应特性的物理意义，理解滤波和滤波器的概念，掌握利用MATLAB计算和绘制LTI系统频率响应特性曲线中的编程。

二、实验原理及方法

1 连续时间LTI系统的频率响应

所谓频率特性，也称为频率响应特性，简称频率响应（Frequency response），是指系统在正弦信号激励下的稳态响应随频率变化的情况，包括响应的幅度随频率的变化情况和响应的相位随频率的变化情况两个方面。

x(t)、y(t)分别为系统的时域激励信号和响应信号，h(t)是系统的单位冲激响应，它们三者之间的关系为：y(t)?x(t)*h(t)，由傅里叶变换的时域卷积定理可得到：

Y(j?)?X(j?)H(j?)

或者： H(j?)?3.1

Y(j?) 3.2

X(j?)H(j?)为系统的频域数学模型，它实际上就是系统的单位冲激响应h(t)的傅里叶变换。即

? H(j?)????j?th(t)edt 3.3 ?由于H(j?)实际上是系统单位冲激响应h(t)的傅里叶变换，如果h(t)是收敛的，或者

说是绝对可积（Absolutly integrabel）的话，那么H(j?)一定存在，而且H(j?)通常是复数，因此，也可以表示成复数的不同表达形式。在研究系统的频率响应时，更多的是把它表示成极坐标形式：

H(j?)?H(j?)ej?(?) 3.4

上式中，H(j?)称为幅度频率相应（Magnitude response），反映信号经过系统之后，信号各频率分量的幅度发生变化的情况，?(?)称为相位特性（Phase response），反映信号

经过系统后，信号各频率分量在相位上发生变换的情况。H(j?)和?(?)都是频率?的函数。

对于一个系统，其频率响应为H(j?)，其幅度响应和相位响应分别为H(j?)和?(?)，如果作用于系统的信号为x(t)?ej?0t，则其响应信号为

y(t)?H(j?0)ej?0t3.5

?H(j?0)ej?(?0)ej?0t?H(j?0)ej(?0t??(?0))

若输入信号为正弦信号，即x(t) = sin(?0t)，则系统响应为

y(t)?H(j?0)sin(?0t)?|H(j?0)|sin(?0t??(?0)) 3.6

可见，系统对某一频率分量的影响表现为两个方面，一是信号的幅度要被H(j?)加权，二是信号的相位要被?(?)移相。

由于H(j?)和?(?)都是频率?的函数，所以，系统对不同频率的频率分量造成的幅度和相位上的影响是不同的。

2 LTI系统的群延时

从信号频谱的观点看，信号是由无穷多个不同频率的正弦信号的加权和（Weighted sum）所组成。正如刚才所述，信号经过LTI系统传输与处理时，系统将会对信号中的所有频率分量造成幅度和相位上的不同影响。从相位上来看，系统对各个频率分量造成一定的相位移（Phase shifting），相位移实际上就是延时（Time delay）。群延时（Group delay）的概念能够较好地反映系统对不同频率分量造成的延时。

LTI系统的群延时定义为：

?(?)??3.7

d?(?)d?

群延时的物理意义：群延时描述的是信号中某一频率分量经过线性时不变系统传输处理后产生的响应信号在时间上造成的延时的时间。

如果系统的相位频率响应特性是线性的，则群延时为常数，也就是说，该系统对于所有的频率分量造成的延时时间都是一样的，因而，系统不会对信号产生相位失真（Phase distortion）。反之，若系统的相位频率响应特性不是线性的，则该系统对于不同频率的频率分量造成的延时时间是不同的，因此，当信号经过系统后，必将产生相位失真。

3 用MATLAB计算系统频率响应

在本实验中，表示系统的方法仍然是用系统函数分子和分母多项式系数行向量来表示。实验中用到的MATLAB函数如下：

[H,w] = freqs(b,a)：b,a分别为连续时间LTI系统的微分方程右边的和左边的系数向

量（Coefficients vector），返回的频率响应在各频率点的样点值（复数）存放在H中，系统

默认的样点数目为200点；

Hm = abs(H)：求模数，即进行Hm?H运算，求得系统的幅度频率响应，返回值

存于Hm之中。

real(H)：求H的实部； imag(H)：求H的虚部；

phi = atan(-imag(H)./(real(H)+eps))：求相位频率相应特性，atan()用来计算反正切值；

或者

phi = angle(H)：求相位频率相应特性；

tao = grpdelay(num,den,w)：计算系统的相位频率响应所对应的群延时。 计算频率响应的函数freqs()的另一种形式是：

H = freqs(b,a,w)：在指定的频率范围内计算系统的频率响应特性。在使用这种形式

的freqs/freqz函数时，要在前面先指定频率变量w的范围。

例如在语句H = freqs(b,a,w)之前加上语句：w = 0:2*pi/256:2*pi。

下面举例说明如何利用上述函数计算并绘制系统频率响应特性曲线的编程方法。 假设给定一个连续时间LTI系统，下面的微分方程描述其输入输出之间的关系

d2y(t)dy(t)?3?2?x(t)

dt2dt编写的MATLAB范例程序，绘制系统的幅度响应特性、相位响应特性、频率响应的实

部和频率响应的虚部。程序如下：

% Program3_1

% This Program is used to compute and draw the plots of the frequency response % of a continuous-time system

b = [1]; % The coefficient vector of the right side of the differential equation a = [1 3 2]; % The coefficient vector of the left side of the differential equation [H,w] = freqs(b,a); % Compute the frequency response H Hm = abs(H); % Compute the magnitude response Hm phai = angle(H); % Compute the phase response phai

Hr = real(H); % Compute the real part of the frequency response

Hi = imag(H); % Compute the imaginary part of the frequency response subplot(221)

plot(w,Hm), grid on, title('Magnitude response'), xlabel('Frequency in rad/sec') subplot(223)

plot(w,phai), grid on, title('Phase response'), xlabel('Frequency in rad/sec') subplot(222)

plot(w,Hr), grid on, title('Real part of frequency response'), xlabel('Frequency in rad/sec') subplot(224)

plot(w,Hi), grid on, title('Imaginary part of frequency response'), xlabel('Frequency in rad/sec')

三、实验内容及步骤

实验前，必须首先阅读本实验原理，了解所给的MATLAB相关函数，读懂所给出的全部

范例程序。实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。并结合范例程序所完成的工作，进一步分析程序中各个语句的作用，从而真正理解这些程序。

实验前，一定要针对下面的实验项目做好相应的实验准备工作，包括事先编写好相应的实验程序等事项。

给定三个连续时间LTI系统，它们的微分方程分别为

d2y(t)dy(t)dx(t)?1?25y(t)?系统1： Eq.3.1 2dtdtdt系统2： 系统3：

dy(t)dx(t)?y(t)??x(t) Eq.3.2 dtdtd6y(t)d5y(t)d4y(t)d3y(t)d2y(t)dy(t)?10?48?148?306?401?262y(t)?262x(t)65432dtdtdtdtdtdt Eq.3.3

Q3-1 修改程序Program3_1，并以Q3_1存盘，使之能够能够接受键盘方式输入的微分方

程系数向量。并利用该程序计算并绘制由微分方程Eq.3.1、Eq.3.2和Eq.3.3描述的系统的幅

度响应特性、相位响应特性、频率响应的实部和频率响应的虚部曲线图。 抄写程序Q3_1如下： clear, close all;

a = input('微分方程左边的系数：'); b = input('微分方程右边的系数：'); [H,w] = freqs(b,a); Hm = abs(H); phai = angle(H); Hr = real(H); Hi = imag(H); subplot(221) plot(w,Hm), grid on,

title('Magnitude response'), xlabel('Frequency in rad/sec') subplot(223) plot(w,phai), grid on,

title('Phase response'),

xlabel('Frequency in rad/sec') subplot(222) plot(w,Hr), grid on,

title('Real part of frequency response'),

xlabel('Frequency in rad/sec') subplot(224) plot(w,Hi), grid on,

title('Imaginary part of frequency response'), xlabel('Frequency in rad/sec')

执行程序Q3_1，绘制的系统1的频率响应特性曲线如下：

Magnitude response11Real part of frequency response0.50.5005Frequency in rad/secPhase response10210-1-205Frequency in rad/sec10510Frequency in rad/secImaginary part of frequency response0.5000-0.505Frequency in rad/sec10

从系统1的幅度频率响应曲线看，系统1是低通、高通、全通、带通还是带阻滤波器？ 答：带通滤波器。

执行程序Q3_1，绘制的系统2的频率响应特性曲线如下：

Magnitude response110.501-0.5-1Real part of frequency response1105Frequency in rad/secPhase response104321005Frequency in rad/sec10510Frequency in rad/secImaginary part of frequency response100.5005Frequency in rad/sec10从系统2的幅度频率响应曲线看，系统2低通、高通、全通、带通还是带阻滤波器？ 答：低通滤波器。

执行程序Q3_1，绘制的系统3的频率响应特性曲线如下：

Magnitude response110.50.50-0.50-1Real part of frequency response05Frequency in rad/secPhase response10420-2-405Frequency in rad/sec10510Frequency in rad/secImaginary part of frequency response100.50-0.5-105Frequency in rad/sec10

从系统3的幅度频率响应曲线看，系统3是低通、高通、全通、带通还是带阻滤波器？ 答：带阻滤波器。

这三个系统的幅度频率响应、相位频率相应、频率响应的实部以及频率响应的虚部分别具有何种对称关系？请根据傅里叶变换的性质说明为什么会具有这些对称关系？ 答：

Q3-2 编写程序Q3_2，使之能够能够接受键盘方式输入的输入信号x(t)的数学表达式，系

统微分方程的系数向量，计算输入信号的幅度频谱，系统的幅度频率响应，系统输出信号y(t)的幅度频谱，系统的单位冲激响应h(t)，并按照下面的图Q3-2的布局，绘制出各个信号的时域和频域图形。

图Q3-2

你编写的程序Q3_2抄写如下： clear, close all; t = 0:0.01:40; T =0.01; dw =0.1;

w=-4*pi:dw:4*pi;

a = input('微分方程左边的系数：'); b = input('微分方程右边的系数：'); x = input('表达式的输入信号x（t）：');

subplot(323); impulse(b,a,40); axis([0 40 -0.2 1]);

grid on,

title('系统单位冲击响应h(t)')

subplot(321), plot(t,x)

title('输入信号x(t):'); xlabel('t/s');

subplot(325), y=lsim(b,a,x,t); plot(t,y)

title('输出信号y(t)'); xlabel('t/s');

X=x*exp(-j*t'*w)*T; X1=abs(X); subplot(322); plot(w,X1),

axis([-4*pi 4*pi 0 20]); grid on,

title('输入信号x(t)的幅度频谱') xlabel('频率 弧度/秒')

Y=y'*exp(-j*t'*w)*T; Y1=abs(Y); subplot(326) plot(w,Y1),

axis([-4*pi 4*pi 0 20]); grid on,

title('输出信号y(t)的幅度频谱'); xlabel('频率 弧度/秒')

[H,w] = freqs(b,a); Hm = abs(H);

phai = angle(H); subplot(324) plot(w,Hm), grid on,

title('系统的幅度频率响应:'), xlabel('频率 弧度/秒')

执行程序Q3_2，输入信号x(t) = sin(t) + sin(8t)，输入由Eq.3.3描述的系统。得到的图形如下：

此处粘帖执行程序Q3_2所得到的图形

输入信号x(t):210-1-2010203040t/s系统单位冲击响应h(t)1edu0.5itlpAm0010203040Time (sec)输出信号y(t)10.50-0.5-1010203040t/s

输入信号x(t)的幅度频谱20151050-10-50510频率 弧度/秒系统的幅度频率响应:10.500510频率 弧度/秒输出信号y(t)的幅度频谱20151050-10-50510频率 弧度/秒

请手工绘制出信号x(t) = sin(t) + sin(8t) 的幅度频谱图如下：

输入信号x(t)的幅度频谱350300250200150100500-15-10-50频率 弧度/秒51015你手工绘制的信号x(t) = sin(t) + sin(8t) 的幅度频谱图与执行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度频谱图是否相同？如不同，是何原因造成的？ 答：不同。

执行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度频谱图实际上是另外一个信号x1(t)的幅度

ejw?e?jw?ej8w?e?j8w频谱，这个信号的时域数学表达式为 x1(t) =

2j请利用傅里叶变换的相关性质计算并绘制信号x1(t)的幅度频谱图。 计算过程：

手工绘制的x1(t) 的幅度频谱图如下：

结合所学的有关滤波的知识，根据上面所得到的信号的时域和频域图形，请从时域和频域两个方面解释滤波的概念。

答：滤波在时域的理解即是输入信号与单位冲激响应的卷积，也即是说，当前时刻的滤波结果，是对当前及之前若干个输入信号加上不同系数的权值，然后再累加的结果。权系数的不同，代表了不同的滤波特性。

滤波在频域的理解即是首先对原信号f(x, y)进行傅里叶变换，将信号从时域变换到频域，获得频域信号F(u, v)，然后用滤波器H(u, v)和F(u, v)相乘，改变原信号的频谱成分，对改变后的频域信号进行傅里叶逆变换，将信号从频域变换回时域，获得滤波后的图像。

四、实验体会

本实验完成时间： 年 月 日