信号与系统

实验教程

（实验报告）

班级:     

姓名:     

目     录

实验一：连续时间信号与系统的时域分析-------------------------------------------------3

一、实验目的及要求---------------------------------------------------------------------------3

二、实验原理-----------------------------------------------------------------------------------3

1、信号的时域表示方法------------------------------------------------------------------4

2、用MATLAB仿真连续时间信号和离散时间信号----------------------------------4

3、LTI系统的时域描述-----------------------------------------------------------------9

三、实验步骤及内容--------------------------------------------------------------------------13

四、实验报告要求-----------------------------------------------------------------------------18

实验二：连续时间信号的频域分析---------------------------------------------------------19

一、实验目的及要求--------------------------------------------------------------------------19

二、实验原理----------------------------------------------------------------------------------19

1、连续时间周期信号的傅里叶级数CTFS---------------------------------------------19

2、连续时间信号的傅里叶变换CTFT--------------------------------------------------20

3、离散时间信号的傅里叶变换DTFT -------------------------------------------------21

4、连续时间周期信号的傅里叶级数CTFS的MATLAB实现------------------------21

5、用MATLAB实现CTFT及其逆变换的计算---------------------------------------25

三、实验步骤及内容----------------------------------------------------------------------27

四、实验报告要求-------------------------------------------------------------------------33

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

一、实验目的及要求--------------------------------------------------------------------------34

二、实验原理----------------------------------------------------------------------------------34

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

2、LTI系统的群延时---------------------------------------------------------------------35

3、用MATLAB计算系统的频率响应--------------------------------------------------36

三、实验步骤及内容----------------------------------------------------------------------37

四、实验报告要求-------------------------------------------------------------------------43

实验一信号与系统的时域分析

一、实验目的

1、熟悉和掌握常用的用于信号与系统时域仿真分析的MATLAB函数；

2、掌握连续时间和离散时间信号的MATLAB产生，掌握用周期延拓的方法将一个非周期信号进行周期信号延拓形成一个周期信号的MATLAB编程；

3、牢固掌握系统的单位冲激响应的概念，掌握LTI系统的卷积表达式及其物理意义，掌握卷积的计算方法、卷积的基本性质；

4、掌握利用MATLAB计算卷积的编程方法，并利用所编写的MATLAB程序验证卷积的常用基本性质；

掌握MATLAB描述LTI系统的常用方法及有关函数，并学会利用MATLAB求解LTI系统响应，绘制相应曲线。

基本要求：掌握用MATLAB描述连续时间信号和离散时间信号的方法，能够编写MATLAB程序，实现各种信号的时域变换和运算，并且以图形的方式再现各种信号的波形。掌握线性时不变连续系统的时域数学模型用MATLAB描述的方法，掌握卷积运算、线性常系数微分方程的求解编程。

二、实验原理

信号（Signal）一般都是随某一个或某几个独立变量的变化而变化的，例如，温度、压力、声音，还有股票市场的日收盘指数等，这些信号都是随时间的变化而变化的，还有一些信号，例如在研究地球结构时，地下某处的密度就是随着海拔高度的变化而变化的。一幅图片中的每一个象素点的位置取决于两个坐标轴，即横轴和纵轴，因此，图像信号具有两个或两个以上的独立变量。

    在《信号与系统》课程中，我们只关注这种只有一个独立变量（Independent variable）的信号，并且把这个独立变量统称为时间变量（Time variable），不管这个独立变量是否是时间变量。

    在自然界中，大多数信号的时间变量都是连续变化的，因此这种信号被称为连续时间信号（Continuous-Time Signals）或模拟信号（Analog Signals），例如前面提到的温度、压力和声音信号就是连续时间信号的例子。但是，还有一些信号的独立时间变量是离散变化的，这种信号称为离散时间信号。前面提到的股票市场的日收盘指数，由于相邻两个交易日的日收盘指数相隔24小时，这意味着日收盘指数的时间变量是不连续的，因此日收盘指数是离散时间信号。

    而系统则用于对信号进行运算或处理，或者从信号中提取有用的信息，或者滤出信号中某些无用的成分，如滤波，从而产生人们所希望的新的信号。系统通常是由若干部件或单元组成的一个整体（Entity）。系统可分为很多不同的类型，例如，根据系统所处理的信号的不同，系统可分为连续时间系统（Continuous-time system）和离散时间系统（Discrete-time system），根据系统所具有的不同性质，系统又可分为因果系统（Causal system）和非因果系统（Noncausal system）、稳定系统（Stable system）和不稳定系统（Unstable system）、线性系统（Linear system）和非线性系统（Nonlinear system）、时变系统（Time-variant system）和时不变系统（Time-invariant system）等等。

    然而，在信号与系统和数字信号处理中，我们所分析的系统只是所谓的线性时不变系统，这种系统同时满足两个重要的基本性质，那就是线性性和时不变性，通常称为线性时不变（LTI）系统。

1. 信号的时域表示方法

1.1将信号表示成独立时间变量的函数

例如              x(t)=sin(ωt) 和 x[n]=n(0.5)nu[n]

分别表示一个连续时间信号和一个离散时间信号。在MATLAB中有许多内部函数，可以直接完成信号的这种表达，例如：

sin()：正弦信号

cos()：余弦信号

exp()：指数信号

1.2用信号的波形图来描述信号

用函数曲线表示一个信号，图1.1就是一个连续时间信号和一个离散时间信号的波形图。

图1.1 连续时间信号与离散时间信号的波形图

1.3将信号用一个数据序列来表示

对于离散时间信号，还可以表示成一个数的序列，例如：

                 x[n]={...., 0.1, 1.1, -1.2, 0, 1.3, ….}

                           ↑n=0

    在《信号与系统》和《数字信号处理》课程中，上述三种信号的描述方法是经常要使用的。

2 用MATLAB仿真连续时间信号和离散时间信号

    在MATLAB中，无论是连续时间信号还是离散时间信号，MATLAB都是用一个数字序列来表示信号，这个数字序列在MATLAB中叫做向量(vector)。通常的情况下，需要与时间变量相对应。

    如前所述，MATLAB有很多内部数学函数可以用来产生这样的数字序列，例如sin()、cos()、exp()等函数可以直接产生一个按照正弦、余弦或指数规律变化的数字序列。

2.1连续时间信号的仿真

程序Program1_1是用MATLAB对一个正弦信号进行仿真的程序，请仔细阅读该程序，并在计算机上运行，观察所得图形。

% Program1_1

% This program is used to generate a sinusoidal signal and draw its plot

clear,                   % Clear all variables

close all,                % Close all figure windows

dt = 0.01;               % Specify the step of time variable

t = -2:dt:2;              % Specify the interval of time

x = sin(2*pi*t);          % Generate the signal

plot(t,x)                % Open a figure window and draw the plot of x(t)

title('Sinusoidal signal x(t)')

xlabel('Time t (sec)')

常用的图形控制函数

axis([xmin,xmax,ymin,ymax])：图型显示区域控制函数，其中xmin为横轴的显示起点，xmax为横轴的显示终点，ymin为纵轴的显示起点，ymax为纵轴的显示终点。

有时，为了使图形具有可读性，需要在所绘制的图形中，加上一些网格线来反映信号的幅度大小。MATLAB中的grid on/grid off可以实现在你的图形中加网格线。

grid on：在图形中加网格线。

grid off：取消图形中的网格线。

x = input(‘Type in signal x(t) in closed form:’)

在《信号与系统》课程中，单位阶跃信号u(t) 和单位冲激信号δ(t) 是二个非常有用的信号。它们的定义如下

                            1.1(a)

                           1.1(b)

这里分别给出相应的简单的产生单位冲激信号和单位阶跃信号的扩展函数。产生单位冲激信号的扩展函数为：

function y = delta(t)

dt = 0.01;

y = (u(t)-u(t-dt))/dt;

产生单位阶跃信号的扩展函数为：

% Unit step function

function y = u(t)

y = (t＞=0);   % y = 1 for t ＞ 0, else y = 0

    请将这二个MATLAB函数分别以delta 和u为文件名保存在work文件夹中，以后，就可以像教材中的方法使用单位冲激信号δ(t) 和单位阶跃信号u(t)。

2.2离散时间信号的仿真

程序Program1_2用来产生离散时间信号x[n]=sin(0.2πn)。

% Program1_2

% This program is used to generate a discrete-time sinusoidal signal and draw its plot

clear,                  % Clear all variables

close all,               % Close all figure windows

n = -10:10;              % Specify the interval of time

x = sin(0.2*pi*n);         % Generate the signal

stem (n,x)               % Open a figure window and draw the plot of x[n]

title ('Sinusoidal signal x[n]')

xlabel ('Time index n')

    请仔细阅读该程序，比较程序Program1_1和Program1_2中的不同之处，以便自己编程时能够正确使用这种方法方针连续时间信号和离散时间信号。

    程序Program1_3用来仿真下面形式的离散时间信号：

                 x[n]={...., 0.1, 1.1, -1.2, 0, 1.3, ….}

                           ↑n=0

% Program1_3

% This program is used to generate a discrete-time sequence

% and draw its plot

clear,                   % Clear all variables

close all,                % Close all figure windows

n = -5:5;                % Specify the interval of time, the number of points of n is 11.

x = [0, 0, 0, 0, 0.1, 1.1, -1.2, 0, 1.3, 0, 0];   % Generate the signal

stem(n,x,'.')             % Open a figure window and draw the plot of x[n]

grid on,

title ('A discrete-time sequence x[n]')

xlabel ('Time index n')

    由于在程序的stem(n,x,'.') 语句中加有'.'选项，因此绘制的图形中每根棒条线的顶端是一个实心点。

    如果需要在序列的前后补较多的零的话，可以利用函数zeros()，其语法为：

    zeros(1, N)：圆括号中的1和N表示该函数将产生一个一行N列的矩阵，矩阵中的所有元素均为零。利用这个矩阵与序列x[n]进行组合，从而得到一个长度与n相等的向量。

    例如，当     x[n]={ 0.1, 1.1, -1.2, 0, 1.3}   时，为了得到程序Program1_3中的序列，

                           ↑n=0

可以用这个MATLAB语句x = [zeros(1,4) x zeros(1, 2)] 来实现。用这种方法编写的程序如下：

% Program1_4

% This program is used to generate a discrete-time sinusoidal signal

% and draw its plot

clear,                    % Clear all variables

close all,                 % Close all figure windows

n = -5:5;                 % Specify the interval of time

x = [zeros(1,4), 0.1, 1.1, -1.2, 0, 1.3, zeros(1,2)];     % Generate the sequence

stem (n,x,'.')              % Open a figure window and draw the plot of x[n]

grid on,

title ('A discrete-time sequence x[n]')

xlabel ('Time index n')

    离散时间单位阶跃信号u[n]定义为

                                            1.2

离散时间单位阶跃信号u[n]除了也可以直接用前面给出的扩展函数来产生，还可以利用MATLAB内部函数ones(1,N) 来实现。这个函数类似于zeros(1,N)，所不同的是它产生的矩阵的所有元素都为1。

    值得注意的是，利用ones(1,N) 来实现的单位阶跃序列并不是真正的单位阶跃序列，而是一个长度为N单位门(Gate)序列，也就是u[n]-u[n-N]。但是在一个有限的图形窗口中，我们看到的还是一个单位阶跃序列。

    在绘制信号的波形图时，有时我们需要将若干个图形绘制在图一个图形窗口中，这就需要使用MATLAB的图形分割函数subplot()，其用法是在绘图函数stem或plot之前，使用图形分割函数subplot(n1,n2,n3)，其中的参数n1，n2和n3的含义是，该函数将把一个图形窗口分割成n1xn2个子图，即将绘制的图形将绘制在第n3个子图中。

2.3信号的时域变换

2.3.1 信号的时移

    信号的时移可用下面的数学表达式来描述：

    设一个连续时间信号为x(t)，它的时移y(t) 表示为：

                           y(t) = x(t - t0)                           1.3

其中，t0为位移量。若t0为正数，则y(t)等于将x(t)右移t0秒之后的结果。反之，若t0为负数，则y(t)等于将x(t)左移t0秒之后的结果。

    在MATLAB中，时移运算与数学上习惯表达方法完全相同。

    程序Program1_5对给定一个连续时间信号x(t) = e-0.5tu(t)，对它分别左移2秒钟和右移2秒钟得到信号x1(t) = e-0.5(t+2)u(t+2)和x2(t) = e-0.5(t-2)u(t-2)。

% Program1_5

% This program is used to implement the time-shift operation

% on a continuous-time signal and to obtain its time-shifted versions

% and to draw their plots.

clear,close all,

t = -5:0.01:5;

x = exp(-0.5*t).*u(t);          % Generate the original signal x(t)

x1 = exp(-0.5*(t+2)).*u(t+2);   % Shift x(t) to the left by 2 second to get x1(t)

x2 = exp(-0.5*(t-2)).*u(t-2);    % Shift x(t) to the right by 2 second to get x2(t)

subplot(311)

plot(t,x)                    % Plot x(t)

grid on,

title ('Original signal x(t)')

subplot (312)

plot (t,x1)                     % Plot x1(t)

grid on,

title ('Left shifted version of x(t)')

subplot (313)

plot (t,x2)                     % Plot x2(t)

grid on,

title ('Right shifted version of x(t)')

xlabel ('Time t (sec)')

2.3.2 信号的时域反褶

    对一个信号x[n]的反褶运算在数学上表示为

                          y[n] = x[-n]                             1.4

    这种反褶运算，用MATLAB实现起来也是非常简单的。有多种方法可以实现信号的反褶运算。

方法一，修改绘图函数plot(t,x)和stem(n,x)中的时间变量t和n，即用-t和-n替代原来的t和n，这样绘制出来的图形，看起来就是原信号经时域反褶后的版本。

方法二，直接利用原信号与其反褶信号的数学关系式来实现。这种方法最符合信号反褶运算的实际意义。

方法三，使用MATLAB内部函数fliplr()来实现信号的反褶运算。其用法如下：

    y = fliplr(x)：其中x为原信号x(t)或x[n]，而y则为x的时域反褶。需要说明的是，函数fliplr()对信号作时域反褶，仅仅将信号中各个元素的次序作了一个反转，这种反转处理是独立于时间变量t和n的。因此，如果信号与其时间变量能够用一个数学函数来表达的话，那么建议将时间变量t和n的范围指定在一个正负对称的时间区间即可。

2.3.3 信号的时域尺度变换

    信号x(t)的时域尺度变换在数学描述为

y(t) = x(at)，                      1.5

其中a为任意常数。根据a的不同取值，这种时域尺度变换对信号x(t)具有非常不同的影响。

    当a = 1时，y(t) = x(t)；

    当a = -1时，y(t) = x(-t)，即y(t)可以通过将x(t)反褶运算而得到；

    当a ＞ 1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的压缩而得到；

    当0 ＜ a ＜ 1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的扩展而得到；

    当 -1 ＜ a ＜ 0时，y(t) = x(at)，y(t)是将x(t)在时间轴上的扩展同时翻转而得到；

    当 a ＜ -1时，y(t) = x(at)，y(t)是将x(t)在时间轴上的压缩同时翻转而得到；

    由此可见，信号的时域尺度变换，除了对信号进行时域压缩或扩展外，还可能包括对信号的时域反褶运算。实际上，MATLAB完成式1.5的运算，并不需要特殊的处理，按照数学上的常规方法即能完成。

2.3.4周期信号

    在《信号与系统》课程中，周期信号是一类非常重要的信号。给定一个信号x(t)或x[n]，如果满足

                         x(t) = x(t+kT)                         1.6

                         x[n] = x[n+kN]                        1.7

则该信号叫做周期信号。其中，k为任意整数，T和N为常数，通常称为信号的基本周期或最小周期。

    周期信号可以看作是一个时限的非周期信号经过周期延拓之后形成的。在数字信号处理中，周期延拓这一信号处理方法非常重要。

下面的程序段，就是将一个非周期信号x1(t) = e-2t[u(t)-u(t-2)]经过周期延拓之后而得到一个周期信号：

clear, close all;

t = -4:0.001:4;

T = 2; x = 0;

for k = -2:2;

    x = x+exp(-2*(t-k*T)).*(u(t-k*T)-u(t-(k+1)*T));

end

仔细阅读该程序，可以发现其算法就是：

                  1.8

由于k无法计算到无穷，而是以有限值加以替代，反映到有限宽度图形窗口中得到的效果完全符合要求。

3 LTI系统的时域描述

3.1线性时不变系统

    在分析LTI系统时，有关LTI系统的两个重要的性质是必须首先掌握和理解的。这就是线性性（Linearity）和时不变性（Time-invariance）。所谓线性性就是指系统同时满足齐次性和叠加性。这可以用下面的方法来描述。

    假设系统在输入信号x1(t)作用时的响应信号为y1(t)，在输入信号x2(t)作用时的响应信号为y2(t)，给定两个常数a和b，如果当输入信号为x(t)时系统的响应信号为y(t)，且满足

                   x(t) = x1(t) + x2(t)                         1.9(a)

                   y(t) = y1(t) + y2(t)                         1.9(b)

则该系统具有叠加性（Additivity）。如果满足

                   x(t) = ax1(t)                              1.10(a)

                   y(t) = ay1(t)                              1.10(b)

则该系统具有齐次性（Homogeneity）。一个系统如果是线性系统的话，那么这个系统必须同时具有叠加性和齐次性。

    又假设系统在输入信号x(t)作用时的响应信号为y(t)，对一个给定时间常数t0，如果当输入信号为x(t-t0)时，系统的响应信号为y(t-t0)的话，则该系统具有时不变性。

    同时具有线性性和时不变性的系统，叫做线性时不变系统，简称LTI系统。LTI系统有连续时间LTI系统和离散时间LTI系统之分。连续时间系统的输入和输出信号都必须是连续时间信号，而离散时间系统的输入和输出信号都必须是离散时间信号。

3.2 LTI系统的单位冲激响应和卷积模型

    给定一个连续时间LTI系统，在系统的初始条件为零时，用单位冲激信号δ(t)作用于系统，此时系统的响应信号称为系统的单位冲激响应（Unit impulse response），一般用h(t)来表示。需要强调的是，系统的单位冲激响应是在激励信号为δ(t)时的零状态响应（Zero-state response）。

    离散时间LTI系统的单位冲激响应的定义与连续时间LTI系统的单位冲激响应相同，只是离散时间单位冲激函数δ[n]的定义有所不同。

    系统的单位冲激响应是一个非常重要的概念，对于一个系统，如果我们知道了该系统的单位冲激响应，那么，该系统对任意输入信号的响应信号都可以求得。也就是说，系统的输入信号x(t)、x[n]和输出信号y(t)、y[n]之间的关系可以用一个数学表达式来描述，这个数学表达式为

                                        1.11(a)

                                          1.11(b)

这个表达式就是LTI系统的卷积模型，它是根据系统的线性性和时不变性以及信号可以分解成单位冲激函数经推理得到的。这个表达式实际上告诉了我们一个重要的结论，那就是，任意LTI系统可以完全由它的单位冲激响应h(t)/h[n]来确定。由于系统的单位冲激响应是零状态响应，故按照式1.11求得的系统响应也是零状态响应。式1.11中的积分运算叫做卷积积分，求和运算叫做卷积和，是描述连续时间系统输入输出关系的一个重要表达式。

3.3卷积的计算

    卷积的计算通常可按下面的五个步骤进行（以卷积积分为例）：

1.         该换两个信号波形图中的横坐标，由t改为τ，τ变成函数的自变量；

2.         把其中一个信号反褶，如把h(τ)变成h(-τ)；

3.         把反褶后的信号做移位，移位量是t，这样t是一个参变量。在τ坐标系中，t ＞ 0时图形右移，            t ＜ 0时图形左移。

4.         计算两个信号重叠部分的乘积x(τ)h(t-τ)；

5.         完成相乘后图形的积分。

对于两个时限信号（Time-limited signal），按照上述的五个步骤，作卷积积分运算时，关键是正确确定不同情况下的积分限。只要正确地确定了积分限都能得到正确定积分结果。尽管如此，在时域中计算卷积积分，总体上来说是一项比较困难的工作。

程序convlution_demo用来演示上述作卷积积分运算的五个步骤。本程序较为复杂，不建议读者读懂该程序，只需执行这个程序，观看程序执行过程中有关卷积积分的运算过程，以便于理解这五个步骤。

借助MATLAB的内部函数conv()可以很容易地完成两个信号的卷积积分运算。其语法为：y = conv(x,h)。其中x和h分别是两个作卷积运算的信号，y为卷积结果。

为了正确地运用这个函数计算卷积，这里有必要对conv(x,h)做一个详细说明。conv(x,h)函数实际上是完成两个多项式的乘法运算。例如，两个多项式p1和p2分别为：

  和 

这两个多项式在MATLAB中是用它们的系数构成一个行向量来表示的，如果用x来表示多项式p1，h表示多项式p2，则x和h分别为

                   x = [1  2  3  4]

                   h = [4  3  2  1]

在MATLAB命令窗口依次键入

＞＞ x = [1  2  3  4];

＞＞ h = [4  3  2  1];

＞＞ y=conv(x,h)

在屏幕上得到显示结果：

y =  4    11    20    30    20    11     4

这表明，多项式p1和p2的乘积为：

    正如前所述，用MATLAB处理连续时间信号时，独立时间变量t的变化步长应该是很小的，假定用符号dt表示时间变化步长，那么，用函数conv()作两个信号的卷积积分时，应该在这个函数之前乘以时间步长方能得到正确的结果。也就是说，正确的语句形式应为：y = dt*conv(x,h)。

对于定义在不同时间段的两个时限信号x(t)，t0 ≤ t ≤ t1，和h(t)，t2 ≤ t ≤ t3。 如果用y(t)来表示它们的卷积结果，则y(t)的持续时间范围要比x(t)或h(t)要长，其时间范围为t0+t2 ≤ t ≤ t1+t3。这个特点很重要，利用这个特点，在处理信号在时间上的位置时，可以很容易地将信号的函数值与时间轴的位置和长度关系保持一致性。

根据给定的两个连续时间信号x(t) = t[u(t)-u(t-1)]和h(t) = u(t)-u(t-1)，编写程序，完成这两个信号的卷积运算，并绘制它们的波形图。范例程序如下：

% Program1_6

% This program computes the convolution of two continuou-time signals

clear;close all;

t0 = -2;    t1 = 4;   dt = 0.01;

t = t0:dt:t1;

x = u(t)-u(t-1);

h = t.*(u(t)-u(t-1));

y = dt*conv(x,h);         % Compute the convolution of x(t) and h(t)

subplot(221)

plot(t,x), grid on, title('Signal x(t)'), axis([t0,t1,-0.2,1.2])

subplot(222)

plot(t,h), grid on, title('Signal h(t)'), axis([t0,t1,-0.2,1.2])

subplot(212)

t = 2*t0:dt:2*t1;          % Again specify the time range to be suitable to the

                       % convolution of x and h.

plot(t,y), grid on, title('The convolution of x(t) and h(t)'), axis([2*t0,2*t1,-0.1,0.6]),

xlabel('Time t sec')

在有些时候，做卷积和运算的两个序列中，可能有一个序列或者两个序列都非常长，甚至是无限长，MATLAB处理这样的序列时，总是把它看作是一个有限长序列，具体长度由编程者确定。实际上，在信号与系统分析中所遇到的无限长序列，通常都是满足绝对可和或绝对可积条件的信号。因此，对信号采取这种截短处理尽管存在误差，但是通过选择合理的信号长度，这种误差是能够减小到可以接受的程度的。若这样的一个无限长序列可以用一个数学表达式表示的话，那么，它的长度可以由编程者通过指定时间变量n的范围来确定。

例如，对于一个单边实指数序列x[n] = 0.5nu[n]，通过指定n的范围为0 ≤n ≤ 100，则对应的x[n]的长度为101点，虽然指定更宽的n的范围，x[n]将与实际情况更相符合，但是，注意到，当n大于某一数时，x[n]之值已经非常接近于0了。对于序列x[n] = 0.5nu[n]，当n = 7时，x[7] = 0.0078，这已经是非常小了。所以，对于这个单边实指数序列，指定更长的n的范围是没有必要的。当然，不同的无限长序列具有不同的特殊性，在指定n的范围时，只要能够反映序列的主要特征就可以了。

3.4  用线性常系数微分方程描述LTI系统

线性常系数微分方程或差分方程是描述LTI系统的另一个时域模型。一个连续时间LTI系统，它的输入信号x(t)输出信号y(t)关系可以用下面的微分方程来表达

                                        1.12

式1.12中，max (N, M)定义为系统的阶。式1.12描述了LTI系统输入信号和输出信号的一种隐性关系（Implicit relationship）。为了求得系统响应信号的显式表达式（Explicit expression），必须对微分方程和差分方程求解。

    在MATLAB中，一个LTI系统也可以用系统微分方程的系数来描述，例如，一个LTI连续时间系统的微分方程为

                           

MATLAB则用两个系数向量num = [1]和den = [1 3 2]来描述该系统，其中num和den分别表示系统微分方程右边和左边的系数，按照微分运算的降阶排列。

MATLAB的内部函数impulse()，step()，initial()，lsim() 可以用来计算并绘制连续时间LTI系统的单位冲激响应，单位阶跃响应，零输入响应和任意信号作用于系统的零状态响应。这些函数的用法描述如下：

h= impulse(num, den, T) 和  impulse(num, den, T)

s = step(num, den, T)    和  step(num, den, T)

y = lsim(num, den, x, t)  和  lsim (num, den, x, t)

函数impulse()，step()用来计算由num和den表示的LTI系统的单位冲激响应和单位阶跃响应，响应的时间范围为0~T，其中den和num分别为系统微分方程左右两边的系数向量，T为指定的响应的终点时间。h和s的点数默认值为101点。由此可以计算时间步长为dt = T/(101-1)。不带返回值的函数如impulse(num, den, T)和step(num, den, T)将直接在屏幕上绘制系统的单位冲激响和单位阶跃响应曲线。带返回值的函数如y = lsim(num, den, x, t)和y = lsim(num, den, x, t)，用来计算由num和den表示的LTI系统在输入信号x作用下的零状态响应。其中t为指定的时间变化范围，x为输入信号，它们的长度应该是相同的。如带返回参数y，则将计算的响应信号保存在y中，若不带返回参数y，则直接在屏幕上绘制输入信号x和响应信号y的波形图。

例如，编写程序，计算并绘制由下面的微分方程表示的系统的单位冲激响应h(t)，单位阶跃响应s(t)。

                         

MATLAB范例程序如下：

% Program1_7

% This program is used to compute the impulse response h(t) and the step response s(t) of a

% continuous-time LTI system

clear, close all;

num = input('Type in the right coefficient vector of differential equation：');

den = input('Type in the left coefficient vector of differential equation：');

t = 0:0.01:8;

x = input('Type in the expression of the input signal x(t)：');

subplot(221),  impulse(num,den,8);

subplot(222),  step(num,den,8)

subplot(223),  plot(t,x);

subplot(224),  y=lsim(num,den,x,t);plot(t,y)

三、实验内容及步骤

实验前，必须首先阅读本实验原理，读懂所给出的全部范例程序。实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。并结合范例程序应该完成的工作，进一步分析程序中各个语句的作用，从而真正理解这些程序。

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

Q1-1：修改程序Program1_1，将dt改为0.2，再执行该程序，保存图形，看看所得图形的效果如何？

dt = 0.01时的信号波形                            

dt = 0.2时的信号波形

此处粘贴图形                                  此处粘贴图形

这两幅图形有什么区别，哪一幅图形看起来与实际信号波形更像？

答：

Q1-2：修改程序Program1_1，并以Q1_2为文件名存盘，产生实指数信号x(t)=e-0.5t。 要求在图形中加上网格线，并使用函数axis()控制图形的时间范围在0~2秒之间。然后执行该程序，保存所的图形。

修改Program1_1后得到的程序Q1_2如下：             

信号x(t)=e-0.5t的波形图

Q1-3：修改程序Program1_1，并以Q1_3为文件名存盘，使之能够仿真从键盘上任意输入的一个连续时间信号，并利用该程序仿真信号x(t)=e-2t。

修改Program1_1后得到的程序Q1_3如下：                    

信号x(t)=e-2t的波形图

                                                               此处粘贴图形

Q1-4：将实验原理中所给的单位冲激信号和单位阶跃信号的函数文件在MATLAB文件编辑器中编写好，并分别以文件名delta和u存入work文件夹中以便于使用。

抄写函数文件delta如下：                       抄写函数文件u如下：

Q1-5：修改程序Program1_4，并以Q1_5为文件名存盘，利用axis()函数，将图形窗口的横坐标范围改为-2≤n≤5，纵坐标范围改为-1.5≤ x ≤1.5。

修改Program1_4后得到的程序Q1_5如下：                         

信号的波形图

                                                               此处粘贴图形

Q1-6：仿照前面的示例程序的编写方法，编写一个MATLAB程序，以Q1_6为文件名存盘，使之能够在同一个图形窗口中的两个子图中分别绘制信号x[n]=0.5|n| 和x(t)=cos(2πt)[u(t)-u(t-3)]。要求选择的时间窗能够表现出信号的主要部分（或特征）。

编写的程序Q1_6如下：

信号x[n]=0.5|n| 的波形图和信号x(t)=cos(2πt)[u(t)-u(t-3)]的波形图

Q1-7：根据示例程序的编程方法，编写一个MATLAB程序，以Q1_7为文件名存盘，由给定信号

                                  x(t) = e-0.5tu(t)

求信号y(t) = x(1.5t+3)，并绘制出x(t) 和y(t)的图形。

编写的程序Q1_7如下：

    信号x(t)的波形图           信号y(t) = x(1.5t+3) 的波形图

此处粘贴图形                    此处粘贴图形

Q1-8：给定一个离散时间信号x[n] = u[n] – u[n-8]，仿照示例程序Program1_5，编写程序Q1_8，产生x[n]的左移序列x1[n] = x[n+6]和右移序列x2[n] = x[n-6]，并在同一个图形窗口的三个子图中分别绘制这三个序列的图形。

编写的程序Q1_8如下：

信号波形图

此处粘贴图形

Q1-9：给定两个离散时间序列

                  x[n] = 0.5n{u[n]-u[n-8]}

                  h[n] = u[n]-u[n-8]

编写程序Q1_10，计算它们的卷积，并分别绘制x[n]、h[n]和它们的卷积y[n]的图形。

编写的程序Q1_10如下：

信号x[n]、h[n]和y[n]的波形图

此处粘贴图形

四、实验报告要求

1、按要求完整书写你所编写的全部MATLAB程序

2、详细记录实验过程中的有关信号波形图（存于自带的U盘中），图形要有明确的标题。全部的MATLAB图形应该用打印机打印，然后贴在本实验报告中的相应位置，禁止复印件。

3、实事求是地回答相关问题，严禁抄袭。

本实验完成时间：    年     月    日

实验二  连续时间信号的频域分析

一、实验目的

1、掌握连续时间周期信号的傅里叶级数的物理意义和分析方法；

2、观察截短傅里叶级数而产生的“Gibbs现象”，了解其特点以及产生的原因；

3、掌握连续时间傅里叶变换的分析方法及其物理意义；

4、掌握各种典型的连续时间非周期信号的频谱特征以及傅里叶变换的主要性质；

5、学习掌握利用MATLAB语言编写计算CTFS、CTFT和DTFT的仿真程序，并能利用这些程序对一些典型信号进行频谱分析，验证CTFT、DTFT的若干重要性质。

基本要求：掌握并深刻理傅里叶变换的物理意义，掌握信号的傅里叶变换的计算方法，掌握利用MATLAB编程完成相关的傅里叶变换的计算。

二、实验原理及方法

1、连续时间周期信号的傅里叶级数CTFS分析

任何一个周期为T1的正弦周期信号，只要满足狄利克利条件，就可以展开成傅里叶级数。

其中三角傅里叶级数为：

                2.1

或：                              2.2

其中，称为信号的基本频率（Fundamental frequency），分别是信号的直流分量、余弦分量幅度和正弦分量幅度，为合并同频率项之后各正弦谐波分量的幅度和初相位，它们都是频率的函数，绘制出它们与之间的图像，称为信号的频谱图（简称“频谱”），－图像为幅度谱，－图像为相位谱。

三角形式傅里叶级数表明，如果一个周期信号x(t)，满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related）的正弦信号所组成，其中每一个不同频率的正弦信号称为正弦谐波分量(Sinusoid component)，其幅度（amplitude）为。也可以反过来理解三角傅里叶级数：用无限多个正弦谐波分量可以合成一个任意的非正弦周期信号。

指数形式的傅里叶级数为：

                            2.3

其中，为指数形式的傅里叶级数的系数，按如下公式计算：

                                            2.4

指数形式的傅里叶级数告诉我们，如果一个周期信号x(t)，满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related）的周期复指数信号所组成，其中每一个不同频率的周期复指数信号称为基本频率分量，其复幅度（complex amplitude）为。这里“复幅度（complex amplitude）”指的是通常是复数。

上面的傅里叶级数的合成式说明，我们可以用无穷多个不同频率的周期复指数信号来合成任意一个周期信号。然而，用计算机（或任何其它设备）合成一个周期信号，显然不可能做到用无限多个谐波来合成，只能取这些有限个谐波分量来近似合成。

假设谐波项数为N，则上面的和成式为：

                            2.5

显然，N越大，所选项数越多，有限项级数合成的结果越逼近原信号x(t)。本实验可以比较直观地了解傅里叶级数的物理意义，并观察到级数中各频率分量对波形的影响包括“Gibbs”现象：即信号在不连续点附近存在一个幅度大约为9%的过冲，且所选谐波次数越多，过冲点越向不连续点靠近。这一现象在观察周期矩形波信号和周期锯齿波信号时可以看得很清楚。

2、连续时间信号傅里叶变换----CTFT

    傅里叶变换在信号分析中具有非常重要的意义，它主要是用来进行信号的频谱分析的。傅里叶变换和其逆变换定义如下：

                                            2.6

                           2.7

    连续时间傅里叶变换主要用来描述连续时间非周期信号的频谱。按照教材中的说法，任意非周期信号，如果满足狄里克利条件，那么，它可以被看作是由无穷多个不同频率（这些频率都是非常的接近）的周期复指数信号ejwt的线性组合构成的，每个频率所对应的周期复指数信号ejwt称为频率分量（frequency component），其相对幅度为对应频率的|X(jw)|之值，其相位为对应频率的X(jw)的相位。

X(jw)通常为关于的复函数，可以按照复数的极坐标表示方法表示为：

X(jw)=| X(jw)|ejÐ X(jw)

其中，| X(jw)|称为x(t)的幅度谱，而ÐX(jw)则称为x(t)的相位谱。

    给定一个连续时间非周期信号x(t)，它的频谱也是连续且非周期的。对于连续时间周期信号，也可以用傅里变换来表示其频谱，其特点是，连续时间周期信号的傅里叶变换时有冲激序列构成的，是离散的——这是连续时间周期信号的傅里叶变换的基本特征。

3、离散时间序列的傅里叶变换---DTFT

给定一个非周期离散时间序列x[n]，它的傅里叶变换定义为

                                           2.8

其反变换定义为                          2.9

    式2.8称为离散时间傅里叶正变换，式2.9称为离散时间傅里叶反变换。由式2.9可以看出，一个非周期离散时间序列，总是可以被看作是由无穷多个不同频率的加权的基本频率分量ejωn组合而成的。对序列中的频率为ω的频率分量来说，其权为X(ejω)，通常是复数，也可以将它表示为

                  X(ejω) = | X(ejω)|ejÐX(ejω)

|X(ejω)|称为序列的幅度谱，而ÐX(ejω)称为序列的相位谱，它们都是频率ω的周期函数。

4、连续周期信号的傅里叶级数CTFS的MATLAB实现

4.1 傅里叶级数的MATLAB计算

设周期信号x(t)的基本周期为T1，且满足狄里克利条件，则其傅里叶级数的系数可由式2.4计算得到。式2.4重写如下：

基本频率为：           

对周期信号进行分析时，我们往往只需对其在一个周期内进行分析即可，通常选择主周期（Principle period）。假定x1(t)是x(t)中的主周期，则

计算机不能计算无穷多个系数，所以我们假设需要计算的谐波次数为N，则总的系数个数为2N+1个。在确定了时间范围和时间变化的步长即T1和dt之后，对某一个系数，上述系数的积分公式可以近似为：

          

对于全部需要的2N+1个系数，上面的计算可以按照矩阵运算实现。MATLAB实现系数计算的程序如下：

dt = 0.01;

T = 2;  t = -T/2:dt:T/2;  w0 = 2*pi/T;

x1 = input(‘Type in the periodic signal x(t) over one period x1(t)=’);

N = input(‘Type in the number N=’);

k = -N:N;   L = 2*N+1;

ak = x1*exp(-j*k*w0*t’)*dt/T;

需要强调的是，时间变量的变化步长dt的大小对傅里叶级数系数的计算精度的影响非常大，dt越小，精度越高，但是，计算机计算所花的时间越长。

例题2-1：给定一个周期为T1 = 2s的连续时间周期方波信号，如图所示，其一个周期内的数学表达式为：

                

解：首先，我们根据前面所给出的公式，计算该信号的傅里叶级数的系数。

                      

因为：w0 = 2π/T1 = π，代入上式得到：

在MATLAB命令窗口，依次键入：

＞＞ k = -10:10;

＞＞ ak = ((-j).^k).* (sin((k+eps)*pi/2)./((k+eps)*pi))  % The expression of ak

ak =

  Columns 1 through 4

  -0.0000                  0 + 0.0354i    -0.0000                  0 + 0.0455i

  Columns 5 through 8

  -0.0000                  0 + 0.0637i    -0.0000                  0 + 0.1061i

  Columns 9 through 12

  -0.0000                  0 + 0.3183i     0.5000                  0 - 0.3183i

  Columns 13 through 16

  -0.0000                  0 - 0.1061i     -0.0000                  0 - 0.0637i

  Columns 17 through 20

  -0.0000                  0 - 0.0455i     -0.0000                  0 - 0.0354i

  Column 21

  -0.0000         

从MATLAB命令窗口，我们得到了该周期信号从到共21个系数。

紧接着再键入以下命令：

＞＞ subplot(221)

＞＞ stem(k,abs(ak),'k.')

＞＞ title('The Fourier series coefficients')

＞＞ xlabel('Frequency index k')

就得到一幅如右图所示的描述与k之间的关系的图形。

以上是我们通过手工计算得到的这个周期信号的傅里叶级数表达式及其频谱图，下面给出完成傅里叶级数系数计算的相应MATLAB范例程序。

% Program2_1

% This program is used to evaluate the Fourier series coefficients ak of a periodic square wave

clear, close all

T = 2;  dt = 0.00001;  t = -2:dt:2;

x1 = u(t) - u(t-1-dt);  x = 0;

for m = -1:1                    % Periodically extend x1(t) to form a periodic signal

    x = x + u(t-m*T) - u(t-1-m*T-dt);

end

w0 = 2*pi/T;

N = 10;                       % The number of the harmonic components

L = 2*N+1;

for k = -N: N;                  % Evaluate the Fourier series coefficients ak

    ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt;

end

phi = anglel(ak);                % Evaluate the phase of ak

执行程序Program2_1后，就完成了信号的傅里叶级数的系数的计算，在命令窗口键入

＞＞ ak

命令窗口就可以显示傅里叶级数的21个系数：

ak =

  Columns 1 through 4

   0.0000 + 0.0000i   0.0000 + 0.0354i   0.0000 - 0.0000i   0.0000 + 0.0455i

  Columns 5 through 8

   0.0000 - 0.0000i   0.0000 + 0.0637i   0.0000 - 0.0000i   0.0000 + 0.1061i

  Columns 9 through 12

   0.0000 - 0.0000i   0.0000 + 0.3183i   0.5000             0.0000 - 0.3183i

  Columns 13 through 16

   0.0000 + 0.0000i   0.0000 - 0.1061i   0.0000 + 0.0000i   0.0000 - 0.0637i

  Columns 17 through 20

   0.0000 + 0.0000i   0.0000 - 0.0455i   0.0000 + 0.0000i   0.0000 - 0.0354i

  Column 21

   0.0000 - 0.0000i

将这里的ak之值同前面手工计算得到的ak比较，可见两者是完全相同的。

再次特别提示：程序中，时间变量的变化步长dt的大小对傅里叶级数系数的计算精度的影响非常大，dt越小，精度越高，本程序中的dt之所以选择0.00001就是为了提高计算精度。但是，计算机所花的计算时间越长。

在程序Program2_1中添加相应的计算| ak |和绘图语句，就可以绘制出信号的幅度谱和相位谱的谱线图。

4.2 周期信号的合成以及Gibbs现象

从傅里叶级数的合成式（Synthesis equation）

可以看出，用无穷多个不同频率和不同振幅的周期复指数信号可以合成一个周期信号。然而，我们无法用计算机实现对无穷多个周期复指数信号的合成。但是，用有限项来合成却是可行的，在实际应用中，多半也就是这么做的。然而，这样做的一个必然结果，就是引入了误差。

如果一个周期信号在一个周期有内断点存在，那么，引入的误差将除了产生纹波之外，还将在断点处产生幅度大约为9%的过冲（Overshot），这种现象被称为吉伯斯现象（Gibbs phenomenon）。

为了能够观察到合成信号与原信号的不同以及Gibbs现象，我们可以利用前面已经计算出的傅里叶级数的系数，计算出截短的傅里叶级数：

这个计算可用L = 2N+1次循环来完成：

其中r作为循环次数，x2在循环之前应先清零。完成这一计算的MATLAB程序为：

x2 = 0;  L = 2*N+1;

for r = 1:L;

      x2 = x2+ak(r)*exp(j*(r-1-N)*w0*t);

end;

完成了所有的计算之后，就可以用绘图函数：plot()和stem()将计算结果包括x1, x2, abs(ak)和angle(ak)以图形的形式给出，便于我们观察。

观察吉伯斯现象的最好的周期信号就是图2-1所示的周期方波信号，这种信号在一个周期内有两个断点，用有限项级数合成这个信号时，吉伯斯现象的特征非常明显，便于观察。

例题2-2：修改程序Program2_1，使之能够用有限项级数合成例题2-1所给的周期方波信号，并绘制出原始周期信号、合成的周期信号、信号的幅度谱和相位谱。

为此，只要将前述的for循环程序段和绘图程序段添加到程序Program2_1中即可，范例程序如下：

% Program2_2

% This program is used to compute the Fourier series coefficients ak of a periodic square wave

clear,close all

T = 2;  dt = 0.00001;  t = -2:dt:2;

x1 = u(t)-u(t-1-dt);  x = 0;

for m = -1:1

    x = x + u(t-m*T) - u(t-1-m*T-dt);   % Periodically extend x1(t) to form a periodic signal

end

w0 = 2*pi/T;

N = input('Type in the number of the harmonic components N = :');

L = 2*N+1;

for k = -N:1:N;

    ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt;

end

phi = angle(ak);

y=0;

for q = 1:L;     % Synthesiz the periodic signal y(t) from the finite Fourier series

    y = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T);

end;

subplot(221),

plot(t,x),   title('The original signal x(t)'),   axis([-2,2,-0.2,1.2]),

subplot(223),

plot(t,y),   title('The synthesis signal y(t)'),   axis([-2,2,-0.2,1.2]),   xlabel('Time t'),

subplot(222)

k=-N:N;   stem(k,abs(ak),'k.'),   title('The amplitude |ak| of x(t)'),   axis([-N,N,-0.1,0.6])

subplot(224)

stem(k,phi,'r.'),   title('The phase phi(k) of x(t)'),   axis([-N,N,-2,2]),   xlabel('Index k')

在用这个程序观察吉伯斯现象时，可以反复执行该程序，每次执行时，输入不同之N值，比较所的图形的区别，由此可以观察到吉伯斯现象的特征。

5 用MATLAB实现CTFT及其逆变换的计算

5.1 用MATLAB实现CTFT的计算

MATLAB进行傅里叶变换有两种方法，一种利用符号运算的方法计算，另一种是数值计算，本实验要求采用数值计算的方法来进行傅里叶变换的计算。严格来说，用数值计算的方法计算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号（Time limited signal），也就是当时间|t|大于某个给定时间时其值衰减为零或接近于零，这个条件与前面提到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。计算机只能处理有限大小和有限数量的数。

采用数值计算算法的理论依据是：

                 

若信号为时限信号，当时间间隔T取得足够小时，上式可演变为：

上式用MATLAB表示为：

                   X=x*exp(j*t’*w)*T

其中X为信号x(t)的傅里叶变换，w为频率Ω，T为时间步长。

相应的MATLAB程序：

T = 0.01; dw = 0.1;       %时间和频率变化的步长

t = -10:T:10;

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

X(jw)可以按照下面的矩阵运算来进行：

X=x*exp(-j*t’*w)*T;  %傅里叶变换

X1=abs(X);          %计算幅度谱

phai=angle(X);       %计算相位谱

为了使计算结果能够直观地表现出来，还需要用绘图函数将时间信号x(t)，信号的幅度谱|X(jw)|和相位谱Ð X(jw)分别以图形的方式表现出来，并对图形加以适当的标注。

5.2 用MATLAB实现傅里叶逆变换

连续时间傅里叶逆变换可用式2.7进行计算。式2.7重写如下：

从定义式可看出，其计算方法与傅里叶变换是一样的，因此可以采用同样的矩阵运算的方法来计算，即

                         x(t)=X(jw)*exp(jw’*t)*dw

具体的MATLAB函数如下：

t = -5:0.01;5;    % 指定信号的时间范围，此范围应根据信号的持续时间确定。

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

X = input(‘Type in the expression of X(jw)’);

x = X* exp(jw’*t)*dw;

然后用绘图函数就可以绘制出逆变换得到的时域信号波形图。

三、实验内容和要求

实验前，必须首先阅读本实验原理，读懂所给出的全部范例程序。实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。并结合范例程序应该完成的工作，进一步分析程序中各个语句的作用，从而真正理解这些程序。

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

Q2-1 编写程序Q2_1，绘制下面的信号的波形图：

   

其中，w0 = 0.5π，要求将一个图形窗口分割成四个子图，分别绘制cos(w0t)、cos(3w0t)、cos(5w0t) 和x(t) 的波形图，给图形加title，网格线和x坐标标签，并且程序能够接受从键盘输入的和式中的项数。

抄写程序Q2_1如下：

执行程序Q2_1所得到的图形如下：

Q2-2 给程序Program2_1增加适当的语句，并以Q2_2存盘，使之能够计算例题2-1中的周期方波信号的傅里叶级数的系数，并绘制出信号的幅度谱和相位谱的谱线图。

通过增加适当的语句修改Program2_1而成的程序Q2_2抄写如下：

执行程序Q2_2得到的图形

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

Q2-3 反复执行程序Program2_2，每次执行该程序时，输入不同的N值，并观察所合成的周期方波信号。通过观察，你了解的吉伯斯现象的特点是：

1、周期信号的傅里叶级数与GIBBS现象

给定如下两个周期信号：

Q2-4 分别手工计算x1(t) 和x2(t) 的傅里叶级数的系数。

信号x1(t) 在其主周期内的数学表达式为：

计算x1(t) 的傅里叶级数的系数的计算过程如下：

通过计算得到的x1(t)的傅里叶级数的系数的数学表达式是：

信号x2(t) 在其主周期内的数学表达式为：

计算x2(t) 的傅里叶级数的系数的计算过程如下：

通过计算得到的x1(t)的傅里叶级数的系数的数学表达式是：

用MATLAB帮助你计算出你手工计算的傅里叶级数的系数ak从-10到10共21个系数。

从命令窗口上抄写x1(t)的21个系数如下：

从命令窗口上抄写x2(t)的21个系数如下：

2. 连续时间非周期信号的傅里叶变换

给定两个时限信号:

        

Q2-5 利用单位阶跃信号u(t)，将x1(t) 表示成一个数学闭式表达式，并手工绘制x1(t) 和x2(t) 的时域波形图。

信号x1(t) 的闭式数学表达式为：x1(t) = ：

手工绘制的x1(t)的时域波形图                    手工绘制的x2(t)的时域波形图

Q2-6 编写MATLAB程序Q2_10，能够接受从键盘输入的时域信号表达式，计算并绘制出信号的时域波形、幅度谱。

程序Q2_10抄写如下

执行程序Q2_10，输入信号x1(t)的数学表达式，得到的信号时域波形、幅度谱和相位谱如下：

执行程序Q2_10，输入信号x2(t)的数学表达式，得到的信号时域波形、幅度谱和相位谱如下：

Q2-7 修改程序Q2_10，并以程序Q2_11为文件名存盘，要求能够接受从键盘输入的时域信号表达式，计算其傅里叶变换，并分别绘制其傅里叶变换的实部、虚部、幅度频谱和相位频谱的图形。

编写的程序Q2_11如下：

选定适当的信号，该信号的时域表达式为：

执行你编写好的MATLAB程序Q2_11，输入你选定的信号的数学表达式，绘制出的该信号的傅里叶变换的图形如下：

Q2-8 修改程序Q2_11，并以Q2_12存盘，要求程序能接受从键盘输入信号的时域表达式，计算并绘制信号的时域波形、信号的幅度频谱和相位频谱图。

编写的程序Q2_12如下：

四、实验报告要求

1、按要求完整书写你所编写的全部MATLAB程序

2、详细记录实验过程中的有关信号波形图（存于自带的U盘中），图形要有明确的标题。全部的MATLAB图形应该用打印机打印，然后贴在本实验报告中的相应位置，禁止复印件。

3、实事求是地回答相关问题，严禁抄袭。

本实验完成时间：    年     月    日

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

一、实验目的

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

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

3、学习和掌握幅度特性、相位特性以及群延时的物理意义；

4、掌握用MATLAB语言进行系统频响特性分析的方法。

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

二、实验原理及方法

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

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

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

                           3.1

或者：                                           3.2

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

                                 3.3

由于H(jw)实际上是系统单位冲激响应h(t)的傅里叶变换，如果h(t)是收敛的，或者说是绝对可积（Absolutly integrabel）的话，那么H(jw)一定存在，而且H(jw)通常是复数，因此，也可以表示成复数的不同表达形式。在研究系统的频率响应时，更多的是把它表示成极坐标形式：

                                             3.4

上式中，称为幅度频率相应（Magnitude response），反映信号经过系统之后，信号各频率分量的幅度发生变化的情况，称为相位特性（Phase response），反映信号经过系统后，信号各频率分量在相位上发生变换的情况。和都是频率w的函数。

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

           3.5

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

                            3.6

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

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

2  LTI系统的群延时

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

LTI系统的群延时定义为：

                                          3.7

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

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

3 用MATLAB计算系统频率响应

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

[H,w] = freqs(b,a)：b,a分别为连续时间LTI系统的微分方程右边的和左边的系数向量（Coefficients vector），返回的频率响应在各频率点的样点值（复数）存放在H中，系统默认的样点数目为200点；

Hm = abs(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系统，下面的微分方程描述其输入输出之间的关系

                         

编写的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系统，它们的微分方程分别为

系统1：                     Eq.3.1

系统2：                            Eq.3.2

系统3：

                                                    Eq.3.3

Q3-1 修改程序Program3_1，并以Q3_1存盘，使之能够能够接受键盘方式输入的微分方程系数向量。并利用该程序计算并绘制由微分方程Eq.3.1、Eq.3.2和Eq.3.3描述的系统的幅度响应特性、相位响应特性、频率响应的实部和频率响应的虚部曲线图。

抄写程序Q3_1如下：

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

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

答：

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

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

答：

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

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

答：

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

答：

Q3-2 编写程序Q3_2，使之能够能够接受键盘方式输入的输入信号x(t)的数学表达式，系统微分方程的系数向量，计算输入信号的幅度频谱，系统的幅度频率响应，系统输出信号y(t)的幅度频谱，系统的单位冲激响应h(t)，并按照下面的图Q3-2的布局，绘制出各个信号的时域和频域图形。

图Q3-2

你编写的程序Q3_2抄写如下：

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

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

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

你手工绘制的信号x(t) = sin(t) + sin(8t) 的幅度频谱图与执行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度频谱图是否相同？如不同，是何原因造成的？

答：

执行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度频谱图实际上是另外一个信号x1(t)的幅度频谱，这个信号的时域数学表达式为   x1(t) =                            

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

    计算过程：

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

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

答：

Q3-3 编写程序Q3_3，能够接受从键盘输入的系统微分方程系数向量，并分别绘制所给三个系统的群延时曲线图。

抄写程序Q3_3如下：

系统Eq.3.1的群延时曲线图               系统Eq.3.3的群延时曲线图

根据上面的群延时曲线图可以看出，对系统Eq.3.1，当频率为5弧度/秒时，群延时为     秒，当频率为10弧度/秒时，群延时为     秒，如何解释这两个群延时时间？

根据上面的群延时曲线图，说明这两个系统是否会造成对信号的相位失真？为什么？

从系统Eq.3.3的群延时曲线图中可以看出，当信号的频率为1弧度/秒时，系统Eq.3.3对这一频率的信号的延时是     秒。所以，执行程序Q3_2时，当作用于系统Eq.3.3的输入信号为x(t) = sin(t) + sin(8t)时，其输出信号y(t)的数学表达式为：

四、实验报告要求

1、按要求完整书写你所编写的全部MATLAB程序

2、详细记录实验过程中的有关信号波形图（存于自带的U盘中），图形要有明确的标题。全部的MATLAB图形应该用打印机打印，然后贴在本实验报告中的相应位置，禁止复印件。

3、实事求是地回答相关问题，严禁抄袭。

本实验完成时间：    年     月    日