1、实验三、 连续信号的采样与恢复一、实验目的(1) 加深理解采样对信号的时域和频域特性的影响；(2) 加深对采样定理的理解和掌握，以及对信号恢复的必要性；(3) 掌握对连续信号在时域的采样与重构的方法。二、实验原理(1) 信号的采样信号的采样原理图如下图所示，其数学模型表示为：=)()(ttffTssnssnTtf)()其中的 f(t)为原始信号， Ts(t)为理想的开关信号（冲激采样信号） Ts(t) =，f s(t)为采样后得到的信号称为采样信号。由此可见，采样信号在时域的表nTt)(示为无穷多冲激函数的线性组合，其权值为原始信号在对应采样时刻的定义值。 )(tf)(tfs)(tTs令原始信
2、号 f(t)的复里叶变换为 F(ejw)=FT(f(t)，则采样信号 fs(t) 的复里叶变换 Fs(ejw)=FT(fs(t)= 。由此可见，采样信号 fs(t)的频谱就是将原始信号 f(t)的频谱nsswjFT)(1在频率轴上以 ws 为周期进行周期延拓后的结果（幅度为原频谱的 1/Ts） 。如果原始信号为有限带宽的信号（即当|w |wm|时，有 F(ejw)=0） ，则有：如果取样频率 ws2w m 时，频谱不发生混叠；否则会出现频谱混叠。(2) 信号的重构设信号 f(t)被采样后形成的采样信号为 fs(t)，信号的重构是指由 fs(t)经过内插处理后，恢复出原来的信号 f(t)的过程。
3、因此又称为信号恢复。由前面的介绍可知，在采样频率 ws2w m 的条件下，采样信号的频谱 Fs(ejw)是以 ws为周期的谱线。选择一个理想低通滤波器，使其频率特性 H(jw)满足：H(jw)=csT|,0式中的 wc 称为滤波器的截止频率，满足 wmw cw s/2。将采样信号通过选择的理想低通滤波器，输出信号的频谱将与原信号的频谱相同。根据信号的时域表示与频率表示的一一对应关系可得，经过理想滤波器还原得到的信号即为原信号本身。信号重构的原理图见下图。低通滤波器)(jFs )(jF通过以上分析，得到如下的时域采样定理：一个带宽为 wm 的带限信号 f(t)，可唯一地由它的均匀取样信号 fs(
4、nTs)确定，其中，取样间隔 Ts1/ 的分量所具有的能量占有很少的比重，所以一般定义 fm=1/ 为门信号的截止频率。其中的 为门信号在时域的宽度。在本例中选取 fm=0.5，临界采样频率为 fs=1，过采样频率为 fs1（为了保证精度，可以将其值提高到该值的 50 倍） ，欠采样频率为 fs1。% 显示原信号及其 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*
5、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 变换);程序运行后的结果见下图。-4 -3 -2 -1 0 1 2 3 4
6、00.20.40.60.81tf(t)f(t)=u(t+1)-u(t-1)-80 -60 -40 -20 0 20 40 60 8000.511.52wF(jw)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
7、: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 变换);程序运行后的结果如下图。-20 -15 -10 -5 0 5 10 15 2000.20.40.60.81nf(n)f(n)-80 -60 -40 -20 0 20 40 60 8002468wF(jw)
8、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);%信号重构figureplot(t1,fa);xlabel(t); ylabel(fa(t)
9、;title(f(t)的重构信号 );t=-4:0.01:4;err=fa-rectpuls(t,2);figure; plot(t,err);sum(abs(err).2)/length(err); % 计算重构信号的均方误差程序运行后的结果见下图。-4 -3 -2 -1 0 1 2 3 4-0.500.511.5tfa(t)f(t)似似似似似示例 3：通过频率滤波的方法，利用 MATLAB 实现对示例 1 中采样信号的重构，并显示重构信号的波形。%采用频率滤波的方法实现对采样信号的重构Ts=0.25; % 采样周期t=-4:Ts:4;f=rectpuls(t,2);w1=2*pi*10;N
10、=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; % 采样信号通过滤
11、波器后的频谱Dt=0.01;t1=-4:Dt:4;fa=Tw/(2*pi)*(Fa*exp(j*wk*t1);% 利用数值计算连续信号的 Fourier 逆变换figuresubplot(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)程序运行后的结果见下图。-4 -3 -2 -1
12、0 1 2 3 4-0.500.511.5tfa(t)f(t)似似似似似-4 -3 -2 -1 0 1 2 3 4-0.6-0.4-0.200.20.40.6terr(t)f(t)似似似似似似似四、实验内容与步骤(1) 修改示例中的门信号宽度、采样周期等参数，重新运行程序，观察得到的采样信号时域和频域特性，以及重构信号与误差信号的变化。(2) 将原始信号分别修改为抽样函数 Sa(t)、正弦信号 sin(20*pi*t)+cos(40*pi*t)、指数信号e-2tu(t)时，在不同采样频率的条件下，观察对应采样信号的时域和频域特性，以及重构信号与误差信号的变化。(3) 利用频域滤波的方法（将采样
13、信号通过一个（butterworth）低通滤波器） ，修改实验中的部分程序，完成对采样信号的重构？五、实验报告要求整理并给出“实验内容与步骤” （1） 、 （2） 、 （3）中的程序代码与产生的图形，并回答下面的问题：(1) 根据实验内容与步骤(1)，说明信号在时域宽度的变化对其频率特性的影响，总结信号在时域的宽度与在频域的宽度的关系；(2) 根据实验内容与步骤(2)和(3) ，运用采样定理的知识，说明采样周期的变化对重构信号质量的影响；根据实验内容与步骤(3)，比较与步骤(2) 重构信号的波形，看看重构信号相对于原信号在时域是否有延时？为什么？如何设计一段程序修正信号的延时，使得重构信号与原
14、始信号基本对齐？实验四、 离散时间系统的时域、Z 域分析一、实验目的(1) 加深对线性时不变离散系统中零状态响应概念的理解，掌握其求解方法；(2) 深刻理解卷积和运算，掌握离散序列求卷积和的计算方法；(3) 掌握求给定离散系统的单位序列响应和单位阶跃序列响应的方法；(4) 加深理解和掌握序列信号的求 Z 变换和逆 Z 变换的方法；(5) 加深理解和掌握离散系统的系统函数零点、极点分布与系统时域特性、系统稳定性的关系。二、实验原理(1)线性时不变 (LTI) 离散时间系统用常系数线性差分方程进行描述：mjni jkfbikya00 其中，fk和 yk分别表示系统的输入和输出，N=max(n,m)
15、是差分方程的阶数。在已知差分方程的初始状态以及输入的条件下，可以通过编程由下式迭代算出系统的输出： mjjnii jkfabikyaky010 )/()/(系统的零状态响应就是在系统初始状态为零条件下微分方程的解。在零初始状态下，MATLAB 控制系统工具箱提供了一个 filter 函数，可以计算差分方程描述的系统的响应，其调用形式为：y =filter(b,a,f)其中， 、 分别是系统差分方程左、右端的系数向量，,.210mbb,.210naaf 表示输入向量，y 表示输出向量。注意，输出序列的长度与输入序列的长度相同。(2)离散系统的冲激响应、阶跃响应分别是输入信号为 和 所对应的零状态
16、响应。)(kuMATLAB 控制系统工具箱专门提供了两个函数求解离散系统的冲激响应和阶跃响应。冲激响应：h = impz(b, a, K)，其中的 h 表示系统的单位序列响应，、 分别是系统差分方程左、右端的系数向量，K 表示,.210mbb.,210na输出序列的时间范围。阶跃响应：g = stepz(b, a, N)，其中的 g 表示系统的单位阶跃序列响应，b 和 a 的含义与上相同，N 表示输出序列的长度。(3)卷积是信号与系统中一个最基本、也是最重要的概念之一。在时域中，对于 LTI 连续时间系统，其零状态响应等于输入信号与系统冲激响应的卷积；而利用卷积定理，这种关系又对应频域中的乘积
17、。MATLAB 信号处理工具箱提供了一个计算两个离散序列卷积和的函数 conv。设向量a、b 代表待卷积的两个序列，则 c = conv(a, b)，就是 a 与 b 卷积后得到的新序列。我们知道两个序列卷积以后，一般而言所得新序列的时间范围、序列长度都会发生变化。例如设 f1(n)长度为 5，3n1；f 2(n)长度为 7，2n8；则卷积后得到的新序列长度为 5+7-1=11，且有在 1n9 时，新序列的值不为零。(4)如有序列 fk(k=0,1,2,)，z 为复变量，则函数 kkzfZF)(称为序列 fk的双边 Z 变换。如果上式的求和只在 k 的非负值域进行，则称为序列的单边Z 变换。M
18、ATLAB 的符号数学工具箱提供了计算 Z 正变换的函数 ztrans 和计算逆 Z 变换的函数 iztrans，其调用形式为：F=ztrans(f)f=iztrans(F)上面两式中，右端得 f 和 F 分别为时域表示式和 Z 域表示式的符号表示，可利用函数 sym来实现，其调用形式为 S=sym(A)。式中的 A 为待分析的表示式的字符串， S 为符号化的数字或变量。(5)线性时不变离散系统可用其 Z 域的系统函数 H(Z)表示，其通常具有如下有理分式的形式 )()(210 zABazabbzHnm为了能从系统的 Z 域表示方便地得到其时域表示式，可将 H(z)展开为部分分式和的形式，在对
19、其求逆 Z 变换。MATLAB 的信号处理工具箱提供了对 H(z)进行部分分式展开的函数residuez，其调用形式为：r, p, k =residuez(B, A)式中的 B 和 A 分别为 H(z)的分子多项式和分母多项式的系数向量，r 为部分分式的向量，p 为极点向量，k 为多项式的系数向量。由此借助于 residuez 函数将上述有理函数 H(z)分解为：)(11 1()2()()()( nmzkkznprzprzAB 进一步通过上面介绍的求逆 Z 变换的方法求出系统的单位序列响应。(6)通过系统函数的表达式，可以方便地求出系统函数的零点和极点。系统函数的零点和极点的位置对于系统的时域
20、特性和频域特性有重要影响。位于 Z 平面的单位圆上和单位圆外的极点将使得系统不稳定。系统函数的零点将使得系统的幅频响应在该频率点附近出现极小值，而其对应的极点将使得系统的幅频响应在该频率点附近出现极大值。在 MATLAB 中可以借助函数 tf2zp 来直接得到系统函数的零点和极点的值，并通过函数 zplane 来显示其零点和极点的分布。利用 MATLAB 中的 impz 函数和 freqz 函数可以求得系统的单位序列响应和频率响应。假定系统函数 H(z)的有理分式形式为： )1()2()1(1nazambbzHnnmtf2zp 函数的调用形式如下： ),(,bptfkz式中的 b 和 a 分别
21、表示 H(z)中(z 的正幂次表示)的分子多项式和分母多项式的系数向量，该函数的作用是将 H(z)转换为用零点、极点和增益常数组成的表示式，即： )()2()1()( npzzpmkzHzplane 函数的调用形式如下：zplane(B,A)式中的 B 和 A 分别表示 H(z)中(z 的负幂次表示)的分子多项式和分母多项式的系数向量，该函数的作用是在 Z 平面画出单位圆以及系统的零点和极点。freqz 函数的调用形式如下：H,w=freqz(B,A)式中的 B 和 A 分别表示 H(z)中(z 的负幂次表示)的分子多项式和分母多项式的系数向量，H 表示频率响应矢量，w 为频率矢量。三、程序示
22、例示例 1：已知系统的差分方程为 （M 点滑动平均系统） ，输入10mjjkfky，其中， 是有用信号， 是噪声信号。求系统的零状kdskfks9.*)2( kd态响应输出。%信号的滑动平均滤波程序示例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(dk,sk,fk);m=5; b=ones(m,1)/m; a=1;y=filt
23、er(b,a,f);subplot(1, 2, 2)plot(k,s,b:,k,y,r-);xlabel(time index k);legend(sk,yk);程序运行结果如下图所示。左图中的三条曲线分别为噪声信号 dk、有用信号 sk和受干扰的输入信号 fk，右图中的三条曲线分别为原始的有用信号 sk和去噪信号 yk。由此可见，该系统实现了对受噪声干扰信号的去噪处理。示例 2：利用函数 impz 和 stepz 求下列离散系统的单位序列响应 hk和单位阶跃响应 gk，并与其理论值比较： 213 kfyky0 20 40 60 80 100-1012345678time index kdks
24、kfk0 20 40 60 80 100-1012345678time index kskyk%计算单位响应和阶跃响应示例；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 o
25、n;title(单位响应的理论值);subplot(2,2,4)stem(k,gk);title(单位阶跃响应的理论值);grid on;程序运行结果如下图所示。图(a)为调用函数的结果显示，图(b)为理论计算值输出显示。0 2 4 6 8 10-2000-10000100020003000 似似似似似似似似0 2 4 6 8 10-1000-500050010001500 似似似似似似似似似似0 2 4 6 8 10-2000-10000100020003000 似似似似似似似似0 2 4 6 8 10-1000-500050010001500 似似似似似似似似似似(a) (b)示例 3：已
26、知序列 xk=1,2,3,4,5; k=-1,0,1,2,3, hk=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(输入序列 xk);subplot(1,2,2)stem(k2,h);grid on;xlabel(单位序列响应 hk);y=conv(x,h);k=k1(1)+k2(1):k1(length(k1)
27、+k2(length(k2);figure;stem(k,y);grid on;xlabel(输出响应 yk);程序运行后的结果如下图所示。图(a)为被卷积的两序列 xk和 hk，图(b)为卷积后的结果显示。-1 0 1 2 300.511.522.533.544.55似似似似xk 0 1 2 3 400.10.20.30.40.50.60.70.80.91似似似似似似hk -1 0 1 2 3 4 5 6 7051015似似似似yk(a) (b)示例 4：已知序列=a k*uk，序列 f2k的 Z 域函数为 F2(z)=z/(z-1/2)2。求：(1) 序列 f1k的 Z 变换；(2) F2
28、(z)的逆 Z 变换。%计算序列的 Z 变换和逆 Z 变换示例；f1=sym(ak);F1=ztrans(f1);F1F2=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)和 f2k=2*(1/2)k*k*uk示例 5：已知因果系统的系统函数为 。利用 MATLAB：(1)计算 H(z)3/1(2/)(zzH的部分分式展开形式；(2) 求系统的单位序列响应并显示波形。%(1) Z 域函数的部分分式展开示例程序B=1;A=1 -0.75 0.125;r,p,k=res
29、iduez(B,A);运行结果为：r = 2, -1p = 0.5000, 0.2500k = 由运行结果说明：该系统有两个极点 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+f2hn=impz(B,A,n);figure;stem(n,hn);运行结果为：h = 2.*.5000n-1.*.2500n0 0.5 1 1.5
30、2 2.5 3 3.5 4 4.5 500.10.20.30.40.50.60.70.80.91示例 6：已知一离散因果系统的系统函数为： 3.0501)(23zzH利用 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);程序运行后，求得零点和极点的
31、值为：z = -1, -1p = 0.5198 + 0.5346i, 0.5198 - 0.5346i, -0.5396k = 1 以及画出零点和极点在 Z 平面的分布图如下图所示。图中圆圈表示零点的位置，旁边的数字 2 表示有一个处在改点的 2 阶零点，叉叉表示极点位置。-1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.812Real PartImaginaryPart%求解系统单位序列响应及频率响应示例程序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);xl
32、abel(k);ylabel(hk);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(ejOmega)|);title(Magnitude response);subplot(2,1,2)plot(w/pi,angle(H);xlabel(ang.freq.Omega(rad/s);ylabel(Angle);title(Angle response);程序运行后的结果分别显示单位响应如下图(a)，以及画出
33、系统幅频响应和相频响应如下图(b)。0 5 10 15 20 25 30 35 40-1-0.500.511.522.5khkImpulse response(a)0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 102468ang.freq.(rad/s)|H(ej)|Magnitude response0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-4-2024ang.freq.(rad/s)AngleAngle response(b)四、实验内容与步骤(1) 已知系统的差分方程为 ，输入为1271.07. kffkyky。计算系统
34、的零状态响应 、单位序列响应 和阶跃响应 ，)4.0(kukf hkg并画出相应的图形（选取 k=0:10） 。(2) 已知系统的单位序列响应为 = ，输入信号为kh5ku。利用 MATLAB 计算：)5()5.0ukkfa. 1fhyb. 22kfhky画出 、 、 和 的波形。1yk(3) 已知因果离散系统的系统函数为 。利用 MATLAB 计算系统函25.04)(2zzH数的零点、极点，在 Z 平面画出其零点、极点的分布，并分析系统的稳定性；求出系统的单位序列响应和频率响应，并分别画出其波形。五、实验报告要求整理并给出“实验内容与步骤” （1） 、 （2） 、 （3）中的程序代码与产生的图形；并回答下面的问题。(1) 在“实验内容与步骤”(1)，其单位序列响应取值不为零的长度应该等于多少？如果将方程左边的系数 0.7 和 0.1 都置为 0，其单位序列响应的长度将会如何改变？ (2) 在“实验内容与步骤”(2)，比较输出 和 的波形，可得到怎样的结论？1ky2(3) 在“实验内容与步骤”(3)，由系统函数的零点、极点可否判定该系统是否稳定？为什么？系统的单位序列响应具有怎样的特征，由此应该判定系统是否稳定？根据系统的幅频特性可判定该系统是什么系统？
