信号与系统及MATLAB实现

《信号与系统及MATLAB实现》实验指导书

前言

长期以来，《信号与系统》课程一直采用单一理论教学方式，同学们依靠做习题来巩固和理解教学内容，虽然手工演算训练了计算能力和思维方法，但是由于本课程数学公式推导较多，概念抽象，常需画各种波形，作题时难免花费很多时间，现在，我们给同学们介绍一种国际上公认的优秀科技应用软件MATLAB，借助它我们可以在电脑上轻松地完成许多习题的演算和波形的绘制。

MATLAB的功能非常强大，我们此处仅用到它的一部分，在后续课程中我们还会用到它，在未来地科学研究和工程设计中有可能继续用它，所以有兴趣的同学，可以对MATLAB 再多了解一些。

MATLAB究竟有那些特点呢？

1．高效的数值计算和符号计算功能，使我们从繁杂的数学运算分析中解脱出来；

2．完备的图形处理功能，实现计算结果和编程的可视化；

3．友好的用户界面及接近数学表达式的自然化语言，易于学习和掌握；

4．功能丰富的应用工具箱，为我们提供了大量方便实用的处理工具；

MATLAB的这些特点，深受大家欢迎，由于个人电脑地普及，目前许多学校已将它做为本科生必须掌握的一种软件。正是基于这些背景，我们编写了这本《信号与系统及MATLAB实现》指导书，内容包括信号的MA TLAB表示、基本运算、系统的时域分析、频域分析、S域分析、状态变量分析等。通过这些练习，同学们在学习《信号与系统》的同时，掌握MATLAB的基本应用，学会应用MATLAB的数值计算和符号计算功能，摆脱烦琐的数学运算，从而更注重于信号与系统的基本分析方法和应用的理解与思考，将课程的重点、难点及部分习题用MATLAB进行形象、直观的可视化计算机模拟与仿真实现，加深对信号与系统的基本原理、方法及应用的理解，为学习后续课程打好基础。另外同学们在进行实验时，最好事先预习一些MATLAB的有关知识，以便更好地完成实验，同时实验中也可利用MATLAB的help命令了解具体语句以及指令的使用方法。

实验一 基本信号在MATLAB 中的表示和运算

一、实验目的

1． 学会用MA TLAB 表示常用连续信号的方法； 2． 学会用MA TLAB 进行信号基本运算的方法； 二、实验原理

1． 连续信号的MATLAB 表示

MATLAB 提供了大量的生成基本信号的函数，例如指数信号、正余弦信号。

表示连续时间信号有两种方法，一是数值法，二是符号法。数值法是定义某一时间范围和取样时间间隔，然后调用该函数计算这些点的函数值，得到两组数值矢量，可用绘图语句画出其波形；符号法是利用MATLAB 的符号运算功能，需定义符号变量和符号函数，运算结果是符号表达的解析式，也可用绘图语句画出其波形图。

例1-1指数信号 指数信号在MATLAB 中用exp 函数表示。

如at Ae t f =)(，调用格式为 ft=A*exp(a*t) 程序是

A=1; a=-0.4;

t=0:0.01:10; %定义时间点

ft=A*exp(a*t); %计算这些点的函数值

plot(t,ft); %画图命令，用直线段连接函数值表示曲线 grid on; %在图上画方格

例1-2 正弦信号 正弦信号在MATLAB 中用 sin 函数表示。

调用格式为 ft=A*sin(w*t+phi) A=1; w=2*pi; phi=pi/6;

t=0:0.01:8; %定义时间点

ft=A*sin(w*t+phi); %计算这些点的函数值 plot(t,ft); %画图命令 grid on; %在图上画方格

例1-3 抽样信号 抽样信号Sa(t)=sin(t)/t 在MA TLAB 中用 sinc 函数表示。

定义为 )/(sin )(πt c t Sa =

t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on;

axis([-10,10,-0.5,1.2]); %定义画图范围，横轴，纵轴 title('抽样信号') %定义图的标题名字

例1-4 三角信号 三角信号在MATLAB 中用 tripuls 函数表示。

调用格式为ft=tripuls(t,width,skew)，产生幅度为1，宽度为width，且以0为中心左右各展开width/2大小，斜度为skew的三角波。width的默认值是1，skew的取值范围是-1~+1之间。一般最大幅度1出现在t=(width/2)*skew的横坐标位置。

t=-3:0.01:3;

ft=tripuls(t,4,0.5);

plot(t,ft); grid on;

axis([-3,3,-0.5,1.5]);

例1-5虚指数信号调用格式是f=exp((j*w)*t)

t=0:0.01:15;

w=pi/4;

X=exp(j*w*t);

Xr=real(X); %取实部

Xi=imag(X); %取虚部

Xa=abs(X); %取模

Xn=angle(X); %取相位

subplot(2,2,1),plot(t,Xr),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]),

title('实部');

subplot(2,2,3),plot(t,Xi),axis([0,15,-(max(Xa)+0.5),max(Xa)+0.5]),

title('虚部');

subplot(2,2,2), plot(t,Xa),axis([0,15,0,max(Xa)+1]),title('模');

subplot(2,2,4),plot(t,Xn),axis([0,15,-(max(Xn)+1),max(Xn)+1]),title('相角');

%subplot(m,n,i) 命令是建立m行n列画图窗口，并指定画图位置i

例1-6复指数信号调用格式是f=exp((a+j*b)*t)

t=0:0.01:3;

a=-1;b=10;

f=exp((a+j*b)*t);

subplot(2,2,1),plot(t,real(f)),title('实部')

subplot(2,2,3),plot(t,imag(f)),title('虚部')

subplot(2,2,2),plot(t,abs(f)),title('模')

subplot(2,2,4),plot(t,angle(f)),title('相角')

例1-7 矩形脉冲信号矩形脉冲信号可用rectpuls函数产生，

调用格式为y=rectpuls(t,width),幅度是1，宽度是width，以t=0为对称中心。

t=-2:0.01:2;

width=1;

ft=2*rectpuls(t,width);

plot(t,ft)

grid on;

例1-8 单位阶跃信号单位阶跃信号u(t)用“t>=0”产生，调用格式为ft=(t>=0) t=-1:0.01:5;

ft=(t>=0);

plot(t,ft); grid on;

axis([-1,5,-0.5,1.5]);

例1-9 正弦信号符号算法

syms t %定义符号变量t

y=sin(pi/4*t) %符号函数表达式

ezplot(y,[-16,16]) %符号函数画图命令

或者

f=sym('sin(pi/4*t)') %定义符号函数表达式

ezplot(f,[-16,16])

例1-10单位阶跃信号MA TTLAB符号数学函数Heaviside表示阶跃信号，但要画图需在工作目录创建Heaviside的M文件

function f=Heaviside(t)

f=(t>0);

保存，文件名是Heaviside ，调用该函数即可画图，例

t=-1:0.01:3;

f=heaviside(t);

plot(t,f)

axis([-1,3,-0.2,1.2])

或者

y=sym('Heaviside(t)');

ezplot(y,[-1,5]);grid on

2．信号基本运算的MATLAB实现

信号基本运算是乘法、加法、尺度、反转、平移、微分、积分，实现方法有数值法和符号法

例1-11 以f(t)为三角信号为例，求f(2t) , f(2-2t)

t=-3:0.001:3;

ft=tripuls(t,4,0.5);

subplot(3,1,1);

plot(t,ft); grid on;

title ('f(t)');

ft1= tripuls(2*t,4,0.5);

subplot(3,1,2);

plot(t,ft1); grid on;

title ('f(2t)');

ft2= tripuls(2-2*t,4,0.5);

subplot(3,1,3);

plot(t,ft2); grid on;

title ('f(2-2t)');

例1-12 已知f1(t)=sinwt , f2(t)=sin8wt , w=2pi , 求f1(t)+f2(t)和f1(t)f2(t) 的波形图w=2*pi;

t=0:0.01:3;

f1=sin(w*t); f2=sin(8*w*t); subplot(211)

plot(t,f1+1,':',t,f1-1,':',t,f1+f2) grid on,title('f1(t)+f2(t))') subplot(212)

plot(t,f1,':',t,-f1,':',t,f1.*f2) grid on,title('f1(t)*f2(t)')

符号算法也可实现上述运算，以信号的微积分运算为例说明符号算法应用 微分的调用格式为 diff(function,’variable’,n) 积分的调用格式为 int(function,’variable’,a,b)

式中function 表示要微分或积分的函数，variable 表示运算变量，n 表示求导阶数，默认值是求一阶导数，a 是积分下限，b 是积分上限，a b 默认是求不定积分。 例1-13 求一阶导数的例题，已知)sin(21ax y =，x x x y ln sin 2= clear

syms a x y1 y2 %定义符号变量a ， x ，y1， y2 y1=sin(a*x^2); %符号函数y1 y2=x*sin(x)*log(x); %符号函数y2

dy1=diff(y1,’x’) %无分号直接显示结果 dy2=diff(y2) %无分号直接显示结果

例1-14 求积分的例题，dx x ax x )2(5

+-?，?+102)

1(dx x xe x

clear

syms a x y3 y4

y3=x^5-a*x^2+sqrt(x)/2; y4=(x*exp(x))/(1+x)^2; iy3=int(y3,'x') iy4=int(y4,0,1) 三、上机实验内容 1． 验证实验原理中程序 2． 画出信号波形

（1）)()2()(2t u e t f t

--= （2）)]2()()[cos 1()(--+=t u t u t t f π

3．信号)()2()(2t u e

t f t

--=，求)2(t f 、)2(t f -波形

实验二 离散信号与系统的时域分析

一、实验目的

1．学会用MA TLAB 表示常用离散信号的方法； 2．学会用MA TLAB 实现离散信号卷积的方法； 3． 学会用MA TLAB 求解离散系统的单位响应； 4． 学会用MA TLAB 求解离散系统的零状态响应； 二、实验原理

1．离散信号的MA TLAB 表示

表示离散时间信号f(k)需要两个行向量，一个是表示序号k=[ ]，一个是表示相应函数值f=[ ]，画图命令是stem 。 例2-1正弦序列信号 正弦序列信号可直接调用MATLAB 函数cos ，例)cos(?ω+k ，当

ωπ/2是整数或分数时，才是周期信号。画)8/cos(?π+k ，)2cos(k 波形程序是：

k=0:40;

subplot(2,1,1)

stem(k,cos(k*pi/8),'filled') title('cos(k*pi/8)') subplot(2,1,2)

stem(k,cos(2*k),'filled') title('cos(2*k)')

例2-2 单位序列信号 ?

??≠==000

1)(k k k δ

本题先建立一个画单位序列)(0k k +δ的M 函数文件，画图时调用。M 文件建立方法：file / new / m-file 在文件编辑窗输入程序，保存文件名用函数名。

function dwxulie(k1,k2,k0) % k1 , k2 是画图时间范围，k0是脉冲位置 k=k1:k2; n=length(k); f=zeros(1,n);

f(1,-k0-k1+1)=1; stem(k,f,'filled') axis([k1,k2,0,1.5]) title('单位序列δ(k)') 保存文件名dwxulie.m

画图时在命令窗口调用，例：dwxulie(-5,5,0) 例2-3 单位阶跃序列信号

??

?<≥=0

00

1)(k k k ε

本题也可先建立一个画单位阶跃序列)(0k k +ε的M 函数文件，画图时调用。 function jyxulie(k1,k2,k0) k=k1:-k0-1; kk=-k0:k2; n=length(k); nn=length(kk) u=zeros(1,n);

uu=ones(1,nn); stem(kk,uu,'filled') hold on

stem(k,u,'filled') hold off

title('单位阶跃序列') axis([k1 k2 0 1.5]) 保存文件名jyxulie.m

画图时在命令窗口调用，例：jyxulie(-3,8,0)

例2-4 实指数序列信号 k ca k f =)( ，c 、 a 是实数。 建立一个画实指数序列的M 函数文件，画图时调用。 function dszsu(c,a,k1,k2) %c ：指数序列的幅度 %a ：指数序列的底数

%k1：绘制序列的起始序号 %k2：绘制序列的终止序号 k=k1:k2; x=c*(a.^k); stem(k,x,'filled') hold on

plot([k1,k2],[0,0]) hold off

调用该函数画信号：)()4

5()(1k k f k

ε= ，)()4

3()(2k k f k

ε-=波形。 dszsu(1,5/4,0,40) dszsu(1,-3/4,0,40)

2 .离散信号的卷积和

两个有限长序列f1，f2卷积可调用MATLAB 函数conv ，调用格式是f=conv(f1,f2), f 是卷积结果，但不显示时间序号，可自编一个函数dconv 给出f 和k ，并画图。 function [f,k]=dconv(f1,f2,k1,k2) %The function of compute f=f1*f2

% f: 卷积和序列f(k)对应的非零样值向量

% k ： 序列f(k)的对应序号向量 % f1: 序列f1(k)非零样值向量 % f2: 序列f2(k)的非零样值向量 % k1: 序列f1(k)的对应序号向量 % k2: 序列f2(k)的对应序号向量 f=conv(f1,f2) %计算序列f1与f2的卷积和f k0=k1(1)+k2(1); %计算序列f 非零样值的起点位置 k3=length(f1)+length(f2)-2; %计算卷积和f 的非零样值的宽度

k=k0:k0+k3 %确定卷积和f 非零样值的序号向量 subplot(2,2,1) stem(k1,f1) %在子图1绘序列f1(k)时域波形图 title('f1(k)') xlabel('k') ylabel('f1(k)') subplot(2,2,2) stem(k2,f2) %在图2绘序列f2(k)时波形图 title('f1(k)') xlabel('k') ylabel('f2(k)') subplot(2,2,3) stem(k,f); %在子图3绘序列f(k)的波形图 title('f(k)f1(k)与f2(k)的卷积和f(k)') xlabel('k') ylabel('f(k)')

h=get(gca,'position'); h(3)=2.5*h(3);

set(gca,'position',h) %将第三个子图的横坐标范围扩为原来的2.5倍 例2-5求卷积和，

)

1()(2)1()(1-+++=k k k k f δδδ)2()1()()1()2()(2-+-+++++=k k k k k k f δδδδδ

f1=[1 2 1];

k1=[-1 0 1]; f2=ones(1,5); k2=-2:2;

[f, k]=dconv(f1,f2,k1,k2)

由运行结果知，f 的长度等于f1和f2长度之和减一， f 的起点是f1和f2的起点之和，f 的终点是f1和f2的终点之和。 3． 离散系统的单位响应

MATLAB 提供画系统单位响应函数impz ，调用格式是 impz(b,a) 式中b 和a 是表示离散系统的行向量；

impz(b,a,n) 式中b 和a 是表示离散系统的行向量，时间范围是0~n ；

impz(b,a,n1,n2) 时间范围是n1~n2 ；y=impz(b,a,n1,n2) 由y 给出数值序列；

例2-6已知 )())2(9.0)1()(k f k y k y k y =-+-- 求单位响应。

a=[1,-1,0.9]; b=[1]; impz(b,a) impz(b,a,60) impz(b,a,-10:40)

4． 离散系统的零状态响应

MATLAB 提供求离散系统零状态响应数值解函数filter ，调用格式为filter(b,a,x),式中b 和a 是表示离散系统的向量，x 是输入序列非零样值点行向量，输出向量序号同x 一样。 例2-7 已知 )1()())2(5.0)1(25.0)(-+=-+--k f k f k y k y k y ， )()2

1

()(k k f k

ε= 求零状态响应， 范围0~20。 a=[1 -0.25 0.5]; b=[1 1]; t=0:20; x=(1/2).^t; y=filter(b,a,x) subplot(2,1,1) stem(t,x)

title('输入序列') subplot(2,1,2) stem(t,y)

title('响应序列') 三、上机实验内容 1．验证实验原理中程序

2．已知)2(2)1(3)()2()1(2)(2-+-+=-+--k f k f k f k y k y k y ，画单位响应波形。 3．已知)()2(25.0)1()(k f k y k y k y =-+-+，输入)()(k t f ε=，画输出波形，范围0~15。

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

一、实验目的

1．学会用MA TLAB 求解连续系统的零状态响应； 2. 学会用MATLAB 求解冲激响应及阶跃响应； 3．学会用MA TLAB 实现连续信号卷积的方法； 二、实验原理

1．连续时间系统零状态响应的数值计算

我们知道，LTI 连续系统可用如下所示的线性常系数微分方程来描述，

()

()0

()()N

M

i j i j i j a y

t b f t ===∑∑

在MA TLAB 中，控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim 。其调用格式

y=lsim(sys,f,t)

式中，t 表示计算系统响应的抽样点向量，f 是系统输入信号向量，sys 是LTI 系统模型，用来表示微分方程，差分方程或状态方程。其调用格式

sys=tf(b,a)

式中，b 和a 分别是微分方程的右端和左端系数向量。例如，对于以下方程:

''''''''''''32103210()()()()()()()()a y t a y t a y t a y t b f t b f t b f t b f t +++=+++

可用32103210[,,,];[,,,];a a a a a b b b b b == (,)sys tf b a = 获得其LTI 模型。

注意，如果微分方程的左端或右端表达式中有缺项，则其向量a 或b 中的对应元素应为零，不能省略不写，否则出错。

例3-1 已知某LTI 系统的微分方程为 y’’(t)+ 2y’(t)+100y(t)=f(t)

其中，'

(0)(0)0,()10sin(2)y y f t t π===，求系统的输出y(t).

解：显然，这是一个求系统零状态响应的问题。其MATLAB 计算程序如下： ts=0;te=5;dt=0.01; sys=tf([1],[1,2,100]); t=ts:dt:te;

f=10*sin(2*pi*t); y=lsim(sys,f,t); plot(t,y);

xlabel('Time(sec)'); ylabel('y(t)');

2．连续时间系统冲激响应和阶跃响应的求解 在MA TLAB 中，对于连续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 实现连续时间信号的卷积

信号的卷积运算有符号算法和数值算法，此处采用数值计算法，需调用MA TLAB 的conv( )函数近似计算信号的卷积积分。连续信号的卷积积分定义是 τττd t f f t f t f t f ?

∞

∞

--=

*=)()()()()(2121

如果对连续信号)(1t f 和)(2t f 进行等时间间隔?均匀抽样，则)(1t f 和)(2t f 分别变为离散时间信号)(1?m f 和)(2?m f 。其中，m 为整数。当?足够小时，)(1?m f 和)(2?m f 既为连续时间信号)(1t f 和)(2t f 。因此连续时间信号卷积积分可表示为

∑?∞

-∞

=→?∞

∞

-?

??-??=-=*=m m t f

m f d t f f t f t f t f )()()()()()()(2

1

02121lim

τ

ττ

采用数值计算时，只求当?=n t 时卷积积分)(t f 的值)(?n f ，其中，n 为整数，既

∑∑∞

-∞

=∞

-∞

=?-???=?

??-???=

?m m m n f m f m n f

m f n f ]

)[()()()()(212

1

其中，

∑∞

-∞

=?-??m m n f

m f ])[()(2

1

实际就是离散序列)(1?m f 和)(2?m f 的卷积和。当?足

够小时，序列)(?n f 就是连续信号)(t f 的数值近似，既 )]()([)()(21n f n f n f t f *?=?≈

上式表明，连续信号)(1t f 和)(2t f 的卷积，可用各自抽样后的离散时间序列的卷积再乘以抽样间隔?。抽样间隔?越小，误差越小。

例3-3用数值计算法求)2()()(1--=t u t u t f 与)()(32t u e t f t -=的卷积积分。

解：因为)()(32t u e t f t -=是一个持续时间无限长的信号，而计算机数值计算不可能计算真正的无限长信号，所以在进行)(2t f 的抽样离散化时，所取的时间范围让)(2t f 衰减到足够小就可以了，本例取5.2=t 。程序是 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')

由于)(1t f 和)(2t f 的时间范围都是从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．画出信号卷积积分)()(21t f t f *的波形，)1()()()(21--==t u t u t f t f

实验四 傅里叶变换、系统的频域分析

一、 实验目的

1、学会用MA TLAB 实现连续时间信号傅里叶变换

2、学会用MA TLAB 分析LTI 系统的频域特性

3、学会用MA TLAB 分析LTI 系统的输出响应 二、实验原理

1．傅里叶变换的MATLAB 求解

MA TLAB 的symbolic Math Toolbox 提供了直接求解傅里叶变换及逆变换的函数fourier()及ifourier()两者的调用格式如下。 Fourier 变换的调用格式

F=fourier(f)：它是符号函数f 的fourier 变换默认返回是关于w 的函数。

F=fourier(f ，v):它返回函数F 是关于符号对象v 的函数，而不是默认的w ，即

()()jvx F v f x e dx +∞

--∞

=?

Fourier 逆变换的调用格式

f=ifourier(F):它是符号函数F 的fourier 逆变换，默认的独立变量为w ，默认返回是

关于x 的函数。

f=ifourier(f,u):它的返回函数f 是u 的函数，而不是默认的x.

注意：在调用函数fourier()及ifourier()之前，要用syms 命令对所用到的变量（如t,u,v,w ）进行说明,即将这些变量说明成符号变量。 例4-1 求2()t

f t e

-=的傅立叶变换

解: 可用MA TLAB 解决上述问题： syms t

Fw=fourier(exp(-2*abs(t)))

例4-2 求2

1

()1F jw ω

=

+的逆变换f(t) 解： 可用MA TLAB 解决上述问题 syms t w

ft=ifourier(1/(1+w^2),t)

2．连续时间信号的频谱图

例4-3 求调制信号t t AG t f 0cos )()(ωτ=的频谱，式中

)2

()2()(,21,12,40τττπωτ--+==

==t u t u t G A 解：MATLAB 程序如下所示

ft=sym('4*cos(2*pi*6*t)*(Heaviside(t+1/4)-Heaviside(t-1/4))'); Fw=simplify(fourier(ft))

subplot(121)

ezplot(ft,[-0.5 0.5]),grid on subplot(122)

ezplot(abs(Fw),[-24*pi 24*pi]),grid

用MA TLAB 符号算法求傅里叶变换有一定局限，当信号不能用解析式表达时，会提示出错，这时用MATLAB 的数值计算也可以求连续信号的傅里叶变换，计算原理是

∑?∞

-∞

=-→-∞

∞

-==n n j t

j e n f dt e

t f j F ττωτωτω)(lim )()(0

当τ足够小时，近似计算可满足要求。若信号是时限的，或当时间大于某个给定值时，信号已衰减的很厉害，可以近似地看成时限信号时，n 的取值就是有限的，设为N ，有

k N N

k e n f k F k N n n j k τ

π

ωτ

ττω2,0,)()(1

=

≤≤=∑-=- 是频率取样点 时间信号取样间隔τ应小于奈奎斯特取样时间间隔，若不是带限信号可根据计算精度要求确定一个频率 W0为信号的带宽。

例4-4 用数值计算法求信号)1()1()(--+=t u t u t f 的傅里叶变换

解，信号频谱是)(2)(ωωSa j F =，第一个过零点是π，一般将此频率视为信号的带宽，若将精度提高到该值的50倍，既W0=50π，据此确定取样间隔，02.00

21

=

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．用MA TLAB 分析LTI 系统的频率特性

当系统的频率响应H （jw ）是jw 的有理多项式时，有

11101

110

()(

)()()()()()()()M M M M N N N N b jw b jw b jw b B w H jw A w a jw a jw a jw a ----++++==++++L L

MATLAB 信号处理工具箱提供的freqs 函数可直接计算系统的频率响应的数值解。其调用格式如下

H=freqs(b,a,w)

其中，a 和b 分别是H(jw)的分母和分子多项式的系数向量，w 为形如w1:p:w2的向量，定义系统频率响应的频率范围，w1为频率起始值，w2为频率终止值,p 为频率取样间隔。H 返回w 所定义的频率点上，系统频率响应的样值。

例如，运行如下命令，计算0~2pi 频率范围内以间隔0.5取样的系统频率响应的样值 a=[1 2 1]; b=[0 1];

h=freqs(b,a,0:0.5:2*pi)

例 4-5 三阶归一化的butterworth 低通滤波器的频率响应为

32

1

()()2()2()1

H jw jw jw jw =

+++ 试画出该系统的幅度响应()H jw 和相位响应()?ω。

解 其MA TLAB 程序及响应的波形如下 w=0:0.025:5; b=[1];a=[1,2,2,1]; H=freqs(b,a,w); subplot(2,1,1);

plot(w,abs(H));grid; xlabel('\omega(rad/s)'); ylabel('|H(j\omega)|'); title('H(jw)的幅频特性'); subplot(2,1,2);

plot(w,angle (H));grid; xlabel('\omega(rad/s)'); ylabel('\phi(\omega)'); title('H(jw)的相频特性');

4．用MA TLAB 分析LTI 系统的输出响应

例 4-6已知一RC 电路如图所示 系统的输入电压为f(t)，输出信号为电阻两端的电压y(t).当RC=0.04，f(t)=cos5t+cos100t, t -∞<<+∞ 试求该系统的响应y(t)

解 由图可知 ，该电路为一个微分电路，其频率响应为 ()11R jw

H jw R jwC jw RC

=

=

++ 由此可求出余弦信号0cos t ω通过LTI 系统的响应为

000()()c o s (())

y t H j w t ω?ω=+ 计算该系统响应的MATLAB 程序及响应波形如下

RC=0.04;

t=linspace(-2,2,1024); w1=5;w2=100;

H1=j*w1/(j*w1+1/RC); H2=j*w2/(j*w2+1/RC); f=cos(5*t)+cos(100*t);

y=abs(H1)*cos(w1*t+angle(H1))+ abs(H2)*cos(w2*t+angle(H2)); subplot(2,1,1); plot(t,f); ylabel('f(t)'); xlabel('Time(s)'); subplot(2,1,2); plot(t,y); ylabel('y(t)'); xlabel('Time(s)'); 三、 上机实验内容

1.验证实验原理中所述的相关程序；

2.试用MA TLAB 求单边指数数信号()()at

f t e u t -=的傅立叶变换，并画出其波形；

3.设2

1

()0.08()0.41

H jw jw jw =

++，试用MATLAB 画出该系统的幅频特性()H jw 和相频特性()?ω，并分析系统具有什么滤波特性。

-

+

-

+ f(t)

y(t) R

C

实验五 信号抽样与恢复

一、实验目的

学会用MA TLAB 实现连续信号的采样和重建 二、实验原理 1．抽样定理

若)(t f 是带限信号，带宽为m ω, )(t f 经采样后的频谱)(ωs F 就是将)(t f 的频谱

)(ωF 在频率轴上以采样频率s ω为间隔进行周期延拓。因此，当s ω≥m ω时，不会发生频

率混叠；而当 s ω

2．信号重建

经采样后得到信号)(t f s 经理想低通)(t h 则可得到重建信号)(t f ，即：

)(t f =)(t f s *)(t h

其中：)(t f s =)

(t f ∑∞∞

--)(s

nT t δ=∑∞

∞

--)()(s

s

nT t nT f δ

)()(t Sa T t h c c

s ωπ

ω=

所以：

)(t f =)(t f s *)(t h =∑∞

∞

--)()(s s nT t nT f δ*)(t Sa T c c

s

ωπ

ω =π

ω

c

s T ∑∞

∞

--)]([)(s

c

s

nT t Sa nT f ω

上式表明，连续信号可以展开成抽样函数的无穷级数。

利用MATLAB 中的t t t c ππ)

sin()(sin =

来表示)(t Sa ，有 )(sin )(π

t c t Sa =，所以可以得到在MA TLAB 中信号由)(s nT f 重建)(t f 的表达式如下：

)(t f =π

ω

c

s T ∑

∞

∞

--)]([

sin )(s c

s nT t c nT f π

ω 我们选取信号)(t f =)(t Sa 作为被采样信号，当采样频率s ω=2m ω时，称为临界采样。我们取理想低通的截止频率c ω=m ω。下面程序实现对信号)(t f =)(t Sa 的采样及由该采样

信号恢复重建)(t Sa ：

例5-1 Sa(t)的临界采样及信号重构；

wm=1; %信号带宽

wc=wm; %滤波器截止频率 Ts=pi/wm; %采样间隔

ws=2*pi/Ts; %采样角频率 n=-100:100; %时域采样电数 nTs=n*Ts %时域采样点 f=sinc(nTs/pi);

Dt=0.005;t=-15:Dt:15;

fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); %信号重构 t1=-15:0.5:15; f1=sinc(t1/pi); subplot(211); stem(t1,f1); xlabel('kTs'); ylabel('f(kTs)');

title('sa(t)=sinc(t/pi)的临界采样信号'); subplot(212); plot(t,fa) xlabel('t'); ylabel('fa(t)');

title('由sa(t)=sinc(t/pi)的临界采样信号重构sa(t)'); grid;

例5-2 Sa(t)的过采样及信号重构和绝对误差分析

程序和例4-1类似，将采样间隔改成Ts=0.7*pi/wm , 滤波器截止频率该成wc=1.1*wm ， 添加一个误差函数 wm=1;

wc=1.1*wm; Ts=0.7*pi/wm; ws=2*pi/Ts; n=-100:100; nTs=n*Ts

f=sinc(nTs/pi);

Dt=0.005;t=-15:Dt:15;

fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t)))); error=abs(fa-sinc(t/pi)); %重构信号与原信号误差 t1=-15:0.5:15; f1=sinc(t1/pi); subplot(311); stem(t1,f1); xlabel('kTs');

ylabel('f(kTs)');

title('sa(t)=sinc(t/pi)的采样信号');

subplot(312);

plot(t,fa)

xlabel('t');

ylabel('fa(t)');

title('由sa(t)=sinc(t/pi)的过采样信号重构sa(t)');

grid;

subplot(313);

plot(t,error);

xlabel('t');

ylabel('error(t)');

title('过采样信号与原信号的误差error(t)');

例5-3 Sa(t)的欠采样及信号重构和绝对误差分析

程序和例4-2类似，将采样间隔改成Ts=1.5*pi/wm , 滤波器截止频率该成wc=wm=1

三、上机实验内容

1．验证实验原理中所述的相关程序；

2．设f(t)=0.5*(1+cost)*(u(t+pi)-u(t-pi)) ，由于不是严格的频带有限信号，但其频谱大部分集中在[0，2]之间，带宽wm可根据一定的精度要求做一些近似。试根据以下两种情况用MATLAB实现由f(t)的抽样信号fs(t)重建f(t) 并求两者误差，分析两种情况下的结果。

(1) wm=2 , wc=1.2wm , Ts=1;

(2) wm=2 , wc=2 , Ts=2.5