信号与系统实验三 四.doc

实验三、 连续信号的采样与恢复
一、实验目的
(1) 加深理解采样对信号的时域和频域特性的影响；
(2) 加深对采样定理的理解和掌握，以及对信号恢复的必要性；
(3) 掌握对连续信号在时域的采样与重构的方法。
二、实验原理
(1) 信号的采样
信号的采样原理图如下图所示，其数学模型表示为：
=
其中的f(t)为原始信号，δTs(t)为理想的开关信号（冲激采样信号）δTs(t) =，fs(t)为采样后得到的信号称为采样信号。由此可见，采样信号在时域的表示为无穷多冲激函数的线性组合，其权值为原始信号在对应采样时刻的定义值。

令原始信号f(t)的复里叶变换为F(ejw)=FT(f(t))，则采样信号fs(t) 的复里叶变换Fs(ejw)=FT(fs(t))=。由此可见，采样信号fs(t)的频谱就是将原始信号f(t)的频谱在频率轴上以ws为周期进行周期延拓后的结果（幅度为原频谱的1/Ts）。如果原始信号为有限带宽的信号（即当|w|>|wm|时，有F(ejw)=0），则有：如果取样频率ws≥2wm时，频谱不发生混叠；否则会出现频谱混叠。
(2) 信号的重构
设信号f(t)被采样后形成的采样信号为fs(t)，信号的重构是指由fs(t)经过内插处理后，恢复出原来的信号f(t)的过程。因此又称为信号恢复。
由前面的介绍可知，在采样频率ws≥2wm的条件下，采样信号的频谱Fs(ejw)是以ws为周期的谱线。选择一个理想低通滤波器，使其频率特性H(jw)满足：
H(jw)=
式中的wc称为滤波器的截止频率，满足wm≤wc≤ws/2。将采样信号通过选择的理想低通滤波器，输出信号的频谱将与原信号的频谱相同。根据信号的时域表示与频率表示的一一对应关系可得，经过理想滤波器还原得到的信号即为原信号本身。信号重构的原理图见下图。

通过以上分析，得到如下的时域采样定理：一个带宽为wm的带限信号f(t)，可唯一地由它的均匀取样信号fs(nTs)确定，其中，取样间隔Ts<π/wm, 该取样间隔又称为奈奎斯特（Nyquist）间隔。
根据时域卷积定理，求出信号重构的数学表达式为：
=
=
=
式中的抽样函数Sa(wct)起着内插函数的作用，信号的恢复可以视为将抽样函数进行不同时刻移位后加权求和的结果，其加权的权值为采样信号在相应时刻的定义值。利用MATLAB中的抽样函数Sinc(t)=sin(πt)/(πt)来表示Sa(t)，有Sa(t)=Sinc(t/π)，于是，信号重构的内插公式也可表示为：
=
(3) 模拟低通滤波器的设计
在任何滤波器的设计中，第一步是确定滤波器阶数N及适当的截止频率Ωc。对于巴特沃斯滤波器，可使用MATLAB命令buttord来确定这些参数，设计滤波器的函数为butter。其调用形式为
[N,wn]=buttord(wp,ws,rp,rs,’s’)
[b,a]=butter[N,wn,’s’]
其中，wp,ws,rp,rs为设计滤波器的技术指标：通带截止频率、阻带截止频率、通带最大衰减和阻带最小衰减；’s’表示设计滤波器的类型为模拟滤波器；N,wn为设计得到的滤波器的阶数和3dB截止频率，b,a为滤波器系统函数的分子和分母多项式的系数矢量，假定系统函数的有理分式表示为：

三、程序示例
示例1：选取门信号f(t)= g2(t)为被采样信号。利用MATLAB实现对信号f(t)的采样，显示原信号与采样信号的时域和频域波形。
因为门信号并非严格意义上的有限带宽信号，但是，由于其频率f>1/τ的分量所具有的能量占有很少的比重，所以一般定义fm=1/τ为门信号的截止频率。其中的τ为门信号在时域的宽度。在本例中选取fm=0.5，临界采样频率为fs=1，过采样频率为fs>1（为了保证精度，可以将其值提高到该值的50倍），欠采样频率为fs<1。
% 显示原信号及其Fourier变换示例
R=0.01;%采样周期
t=-4:R:4;
f=rectpuls(t,2)
w1=2*pi*10; % 显示[-20*pi 20*pi]范围内的频谱
N=1000; % 计算出2*1000+1个频率点的值
k=0:N;
wk=k*w1/N;
F=f*exp(-j*t'*wk)*R; % 利用数值计算连续信号的Fourier变换
F=abs(F); % 计算频谱的幅度
wk=[-fliplr(wk),wk(2:1001)];
F=[fliplr(F),F(2:1001)]; % 计算对应负频率的频谱
figure;
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(wk,F);xlabel('w'); ylabel('F(jw)');
title('f(t)的Fourier变换');
程序运行后的结果见下图。

%显示采样信号及其Fourier变换示例
R=0.25; % 可视为过采样
t=-4:R:4;
f=rectpuls(t,2);
w1=2*pi*10;
N=1000;
k=0:N;
wk=k*w1/N;
F=f*exp(-j*t'*wk); % 利用数值计算采样信号的Fourier变换
F=abs(F);
wk=[-fliplr(wk),wk(2:1001)]; % 将正频率扩展到对称的负频率
F=[fliplr(F),F(2:1001)]; %将正频率的频谱扩展到对称的负频率的频谱
figure;
subplot(2,1,1)
stem(t/R,f); % 采样信号的离散时间显示
xlabel('n'); ylabel('f(n)');
title('f(n)');
subplot(2,1,2)
plot(wk,F); % 显示采样信号的连续的幅度谱
xlabel('w'); ylabel('F(jw)');
title('f(n)的Fourier变换');
程序运行后的结果如下图。

示例2：利用MATLAB实现对示例1中采样信号的重构，并显示重构信号的波形。
%采样信号的重构及其波形显示示例程序
Ts=0.25; % 采样周期，可修改
t=-4:Ts:4;
f=rectpuls(t,2); % 给定的采样信号
ws=2*pi/Ts;
wc=ws/2;
Dt=0.01;
t1=-4:Dt:4; % 定义信号重构对应的时刻，可修改
fa=Ts*wc/pi*(f*sinc(wc/pi*(ones(length(t),1)*t1-t'*ones(1,length(t1)))));
%信号重构
figure
plot(t1,fa);
xlabel('t'); ylabel('fa(t)');
title('f(t)的重构信号');
t=-4:0.01:4;
err=fa-rectpuls(t,2);
figure; plot(t,err);
sum(abs(err).^2)/length(err); % 计算重构信号的均方误差
程序运行后的结果见下图。

示例3：通过频率滤波的方法，利用MATLAB实现对示例1中采样信号的重构，并显示重构信号的波形。
%采用频率滤波的方法实现对采样信号的重构
Ts=0.25; % 采样周期
t=-4:Ts:4;
f=rectpuls(t,2);
w1=2*pi*10;
N=1000;
k=0:N;
wk=k*w1/N;
F=f*exp(-j*t'*wk); % 利用数值计算连续信号的Fourier变换
wk=[-fliplr(wk),wk(2:1001)];
F=[fliplr(F),F(2:1001)];
Tw=w1/N; % 频率采样间隔
w=-2*pi*10:Tw:2*pi*10;
H=Ts*rectpuls(w,2.*pi/Ts); % 理想低通滤波器频率特性
%下面两行为可修改程序
%[b,a]=butter(M, Wc,'s'); % 确定系统函数的系数矢量，M, Wc为设定滤波器的阶数和3dB截止频率
%H=Ts*freqs(b,a,w);
Fa=F.*H; % 采样信号通过滤波器后的频谱
Dt=0.01;
t1=-4:Dt:4;
fa=Tw/(2*pi)*(Fa*exp(j*wk'*t1));
% 利用数值计算连续信号的Fourier逆变换
figure
subplot(2,1,1) ; plot(t1,fa);xlabel('t');ylabel('fa(t)');title('f(t)的重构信号');
err=fa-rectpuls(t1,2);
subplot(2,1,2);plot(t1,err); xlabel('t');ylabel('err(t)');title('f(t)的重构误差信号');
sum(abs(err).^2)/length(err)
程序运行后的结果见下图。

四、实验内容与步骤
(1) 修改示例中的门信号宽度、采样周期等参数，重新运行程序，观察得到的采样信号时域和频域特性，以及重构信号与误差信号的变化。
(2) 将原始信号分别修改为抽样函数Sa(t)、正弦信号sin(20*pi*t)+cos(40*pi*t)、指数信号e-2tu(t)时，在不同采样频率的条件下，观察对应采样信号的时域和频域特性，以及重构信号与误差信号的变化。
(3) 利用频域滤波的方法（将采样信号通过一个（butterworth）低通滤波器），修改实验中的部分程序，完成对采样信号的重构？
五、实验报告要求
整理并给出“实验内容与步骤”（1）、（2）、（3）中的程序代码与产生的图形，并回答下面的问题：
(1) 根据实验内容与步骤(1)，说明信号在时域宽度的变化对其频率特性的影响，总结信号在时域的宽度与在频域的宽度的关系；
(2) 根据实验内容与步骤(2)和(3)，运用采样定理的知识，说明采样周期的变化对重构信号质量的影响；
根据实验内容与步骤(3)，比较与步骤(2)重构信号的波形，看看重构信号相对于原信号在时域是否有延时？为什么？如何设计一段程序修正信号的延时，使得重构信号与原始信号基本对齐？

实验四、 离散时间系统的时域、Z域分析
一、实验目的
(1) 加深对线性时不变离散系统中零状态响应概念的理解，掌握其求解方法；
(2) 深刻理解卷积和运算，掌握离散序列求卷积和的计算方法；
(3) 掌握求给定离散系统的单位序列响应和单位阶跃序列响应的方法；
(4) 加深理解和掌握序列信号的求Z变换和逆Z变换的方法；
(5) 加深理解和掌握离散系统的系统函数零点、极点分布与系统时域特性、系统稳定性的关系。
二、实验原理
(1)
线性时不变 (LTI) 离散时间系统用常系数线性差分方程进行描述：

其中，f[k]和y[k]分别表示系统的输入和输出，N=max(n,m)是差分方程的阶数。
在已知差分方程的初始状态以及输入的条件下，可以通过编程由下式迭代算出系统的输出：

系统的零状态响应就是在系统初始状态为零条件下微分方程的解。在零初始状态下，MATLAB控制系统工具箱提供了一个filter函数，可以计算差分方程描述的系统的响应，其调用形式为：y=filter(b,a,f)
其中，、分别是系统差分方程左、右端的系数向量，f表示输入向量，y表示输出向量。注意，输出序列的长度与输入序列的长度相同。
(2)
离散系统的冲激响应、阶跃响应分别是输入信号为和所对应的零状态响应。MATLAB控制系统工具箱专门提供了两个函数求解离散系统的冲激响应和阶跃响应。
冲激响应：h = impz(b, a, K)，其中的h表示系统的单位序列响应，、分别是系统差分方程左、右端的系数向量，K表示输出序列的时间范围。
阶跃响应：g = stepz(b, a, N)，其中的g表示系统的单位阶跃序列响应，b和a的含义与上相同，N表示输出序列的长度。
(3)
卷积是信号与系统中一个最基本、也是最重要的概念之一。在时域中，对于LTI 连续时间系统，其零状态响应等于输入信号与系统冲激响应的卷积；而利用卷积定理，这种关系又对应频域中的乘积。
MATLAB信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv。设向量a、b代表待卷积的两个序列，则c = conv(a, b)，就是a与b卷积后得到的新序列。
我们知道两个序列卷积以后，一般而言所得新序列的时间范围、序列长度都会发生变化。例如设f1(n)长度为5，－3≤n≤1；f2(n)长度为7，2≤n≤8；则卷积后得到的新序列长度为5+7-1=11，且有在－1≤n≤9时，新序列的值不为零。
(4)
如有序列f[k](k=0,±1,±2,…,)，z为复变量，则函数

称为序列f[k]的双边Z变换。如果上式的求和只在k的非负值域进行，则称为序列的单边Z变换。
MATLAB的符号数学工具箱提供了计算Z正变换的函数ztrans和计算逆Z变换的函数iztrans，其调用形式为：
F=ztrans(f)
f=iztrans(F)
上面两式中，右端得f和F分别为时域表示式和Z域表示式的符号表示，可利用函数sym来实现，其调用形式为S=sym(A)。式中的A为待分析的表示式的字符串，S为符号化的数字或变量。
(5)
线性时不变离散系统可用其Z域的系统函数H(Z)表示，其通常具有如下有理分式的形式

为了能从系统的Z域表示方便地得到其时域表示式，可将H(z)展开为部分分式和的形式，在对其求逆Z变换。MATLAB的信号处理工具箱提供了对H(z)进行部分分式展开的函数residuez，其调用形式为：[r, p, k]=residuez(B, A)
式中的B和A分别为H(z)的分子多项式和分母多项式的系数向量，r为部分分式的向量，p为极点向量，k为多项式的系数向量。由此借助于residuez函数将上述有理函数H(z)分解为：

进一步通过上面介绍的求逆Z变换的方法求出系统的单位序列响应。
(6)
通过系统函数的表达式，可以方便地求出系统函数的零点和极点。系统函数的零点和极点的位置对于系统的时域特性和频域特性有重要影响。位于Z平面的单位圆上和单位圆外的极点将使得系统不稳定。系统函数的零点将使得系统的幅频响应在该频率点附近出现极小值，而其对应的极点将使得系统的幅频响应在该频率点附近出现极大值。
在MATLAB中可以借助函数tf2zp来直接得到系统函数的零点和极点的值，并通过函数zplane来显示其零点和极点的分布。利用MATLAB中的impz函数和freqz函数可以求得系统的单位序列响应和频率响应。假定系统函数H(z)的有理分式形式为：

tf2zp函数的调用形式如下：

式中的b和a分别表示H(z)中(z的正幂次表示)的分子多项式和分母多项式的系数向量，该函数的作用是将H(z)转换为用零点、极点和增益常数组成的表示式，即：

zplane函数的调用形式如下：
zplane(B,A)
式中的B和A分别表示H(z)中(z的负幂次表示)的分子多项式和分母多项式的系数向量，该函数的作用是在Z平面画出单位圆以及系统的零点和极点。
freqz函数的调用形式如下：
[H,w]=freqz(B,A)
式中的B和A分别表示H(z)中(z的负幂次表示)的分子多项式和分母多项式的系数向量，H表示频率响应矢量，w为频率矢量。
三、程序示例
示例1：已知系统的差分方程为（M点滑动平均系统），输入，其中，是有用信号，是噪声信号。求系统的零状态响应输出。
%信号的滑动平均滤波程序示例
R=100;%输入信号的长度
d=rand(1,R)-0.5;%产生(-0.5,0.5)均匀分布的随机噪声
k=0:R-1;
s=2*k.*(0.9.^k);
f=s+d;
figure;
subplot(1, 2, 1)
plot(k,d,'r-',k,s,'b:',k,f,'k-');
xlabel('time index k');legend('d[k]','s[k]','f[k]');
m=5; b=ones(m,1)/m; a=1;
y=filter(b,a,f);
subplot(1, 2, 2)
plot(k,s,'b:',k,y,'r-');
xlabel('time index k');legend('s[k]','y[k]');
程序运行结果如下图所示。左图中的三条曲线分别为噪声信号d[k]、有用信号s[k]和受干扰的输入信号f[k]，右图中的三条曲线分别为原始的有用信号s[k]和去噪信号y[k]。由此可见，该系统实现了对受噪声干扰信号的去噪处理。
示例2：利用函数impz和stepz求下列离散系统的单位序列响应h[k]和单位阶跃响应g[k]，并与其理论值比较：

%计算单位响应和阶跃响应示例；
k=0:10;
a=[1 3 2];
b=[1];
h=impz(b,a,k);
g=stepz(b,a,length(k));
figure;
subplot(2,2,1)
stem(k,h);
grid on;
title('单位响应的近似值');
subplot(2,2,3)
stem(k,g);
title('单位阶跃响应的近似值');
grid on;
hk=-(-1).^k+2*(-2).^k;
gk=1/6-(-1).^k/2+4*(-2).^k/3;
subplot(2,2,2)
stem(k,hk);
grid on;
title('单位响应的理论值');
subplot(2,2,4)
stem(k,gk);
title('单位阶跃响应的理论值');
grid on;
程序运行结果如下图所示。图(a)为调用函数的结果显示，图(b)为理论计算值输出显示。

(a) (b)
示例3：已知序列x[k]={1,2,3,4,5; k=-1,0,1,2,3}, h[k]={1,1,1,1,1; k=0,1,2,3,4}，利用conv函数计算两个序列卷积后的新序列并显示结果。
%利用函数conv(a,b)计算两序列的卷积和；
k1=-1:3;
x=[1 2 3 4 5];
k2=0:4;
h=[1 1 1 1 1];
figure(1);
subplot(1,2,1)
stem(k1,x);
grid on;
xlabel('输入序列x[k]');
subplot(1,2,2)
stem(k2,h);
grid on;
xlabel('单位序列响应h[k]');
y=conv(x,h);
k=k1(1)+k2(1):k1(length(k1))+k2(length(k2));
figure;
stem(k,y);
grid on;
xlabel('输出响应y[k]');
程序运行后的结果如下图所示。图(a)为被卷积的两序列x[k]和h[k]，图(b)为卷积后的结果显示。

(a) (b)
示例4：已知序列=ak*u[k]，序列f2[k]的Z域函数为F2(z)=z/(z-1/2)2。
求：(1) 序列f1[k]的Z变换；(2) F2(z)的逆Z变换。
%计算序列的Z变换和逆Z变换示例；
f1=sym('a^k');
F1=ztrans(f1);
F1
F2=sym('z/(z-1/2)^2');
f2=iztrans(F2);
f2
程序运行后的结果为：
F1=
z/a/(z/a-1)
f2 =
2*(1/2)^n*n
注：即有F1(z)=z/(z-a)和f2[k]=2*(1/2)k*k*u[k]
示例5：已知因果系统的系统函数为。利用MATLAB：(1)计算H(z)的部分分式展开形式；(2) 求系统的单位序列响应并显示波形。
%(1) Z域函数的部分分式展开示例程序
B=[1];
A=[1 -0.75 0.125];
[r,p,k]=residuez(B,A);
运行结果为：
r = 2, -1
p = 0.5000, 0.2500
k = []
由运行结果说明：该系统有两个极点p(1)=0.5和p(2)=0.25；展开的多项式系数全为零。
%(2) 通过逆Z变换求系统单位响应示例
B=[1];
A=[1 -0.75 0.125];
n=0:5;
F1=sym('2/(1-z^(-1)/2.)');
f1=iztrans(F1)
F2=sym('-1/(1-z^(-1)/4.)');
f2=iztrans(F2)
h=f1+f2
hn=impz(B,A,n);
figure;
stem(n,hn);
运行结果为：
h = 2.*.5000^n-1.*.2500^n

示例6：已知一离散因果系统的系统函数为：

利用MATLAB：(1) 求出系统函数的零点和极点，并在Z平面显示它们的分布；(2)求出系统的单位序列响应并显示；(3)求出系统的频率响应，并画出幅频响应和相频响应特性曲线。
%求离散系统零、极点并显示的示例程序
b=[1 2 1];
a=[1 -0.5 -0.005 0.3];
[z,p,k] = tf2zp(b,a)
B=[0 1 2 1];
A=[1 -0.5 -0.005 0.3];
figure;
zplane(B,A);
程序运行后，求得零点和极点的值为：
z = -1, -1
p = 0.5198 + 0.5346i, 0.5198 - 0.5346i, -0.5396
k = 1
以及画出零点和极点在Z平面的分布图如下图所示。图中圆圈表示零点的位置，旁边的数字2表示有一个处在改点的2阶零点，叉叉表示极点位置。

%求解系统单位序列响应及频率响应示例程序
B=[0 1 2 1];
A=[1 -0.5 -0.005 0.3];
k=0:40;
h=impz(B,A,k);
figure(1);
stem(k,h);
xlabel('k');
ylabel('h[k]');
title('Impulse response');
[H,w]=freqz(B,A);
figure(2);
subplot(2,1,1)
plot(w/pi,abs(H));
xlabel('ang.freq.\Omega(rad/s)');
ylabel('|H(e^j^\Omega)|');
title('Magnitude response');
subplot(2,1,2)
plot(w/pi,angle(H));
xlabel('ang.freq.\Omega(rad/s)');
ylabel('Angle');
title('Angle response');
程序运行后的结果分别显示单位响应如下图(a)，以及画出系统幅频响应和相频响应如下图(b)。

(a)

(b)
四、实验内容与步骤
(1) 已知系统的差分方程为，输入为。计算系统的零状态响应、单位序列响应和阶跃响应，并画出相应的图形（选取k=0:10）。
(2) 已知系统的单位序列响应为=，输入信号为。利用MATLAB计算：
a.
b.
画出、、和的波形。
(3) 已知因果离散系统的系统函数为。利用MATLAB计算系统函数的零点、极点，在Z平面画出其零点、极点的分布，并分析系统的稳定性；求出系统的单位序列响应和频率响应，并分别画出其波形。
五、实验报告要求
整理并给出“实验内容与步骤”（1）、（2）、（3）中的程序代码与产生的图形；并回答下面的问题。
(1) 在“实验内容与步骤”(1)，其单位序列响应取值不为零的长度应该等于多少？如果将方程左边的系数0.7和0.1都置为0，其单位序列响应的长度将会如何改变？
(2) 在“实验内容与步骤”(2)，比较输出和的波形，可得到怎样的结论？
(3) 在“实验内容与步骤”(3)，由系统函数的零点、极点可否判定该系统是否稳定？为什么？系统的单位序列响应具有怎样的特征，由此应该判定系统是否稳定？根据系统的幅频特性可判定该系统是什么系统？