数字信号处理实验指导书

前 言数字信号处理是一门理论和工程实践密切结合的课程。为了加深对教学内容的理解，应在学习理论的同时，加强上机实验，深入理解和消化基本理论，锻炼初学者独立解决问题的能力。本课程实验要求学生运用MATLAB编程完成一些数字信号处理的基本功能。 MATLAB是一高效的工程计算语言，它将计算、可视化和编程等功能集于一个易于使用的环境。在MATLAB环境中描述问题计编制求解问题的程序时，用户可以按照符合人们科学思维的方式和数学表达习惯的语言形式来书写程序。MATLAB广泛应用于工业，电子，医疗和建筑等众多领域。其典型应用主要包括以下几个方面：数学计算；算法开发；数据采集；系统建模和仿真；数据分析和可视化科学和工程绘图；应用软件开发（包括用户界面）。 1 / 100； 实验1 用MATLAB产生时域离散信号一、.实验目的：1、了解常用时域离散信号及其特点2、掌握用MATLAB产生时域离散信号的方法二、.实验原理：1、时域离散信号的概念在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x(n)表示，其幅度可以在某一范围内连续取值。由于信号处理设备或装置（如计算机、专用的信号处理芯片等）均以有限位的二进制数来表示信号的幅度，因此，信号的幅度也必须离散化。我们把时间和幅度均取离散值的信号称为时域离散信号或数字信号。在MATLAB语言中，时域离散信号可以通过编写程序直接产生。2、常用时域离散信号的生成1) 单位抽样序列单位抽样序列的表示式为 或 以下三段程序分别用不同的方法来产生单位抽样序列。例1-1 用MATLAB的关系运算式来产生单位抽样序列。n1= -5;n2=5;n0=0;n=n1:n2;x=n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x); xlabel(时间(n);ylabel(幅度x(n);title(单位脉冲序列); 运行结果如图1-1所示：图1-1例1-2 用zeros函数和抽样点直接赋值来产生单位抽样序列。n1=-5;n2=5;k=0;n=n1:n2;nt=length(n);（取参数的长度）nk=abs(k-n1)+1;取绝对值x=zeros(1,nt);（行、列）（zeros生成零矩阵）x(nk)=1;绘图部分的程序及作图结果与例1-1相同。例1-3 生成移位的单位脉冲序列。 n1=-5;n2=5;n0=2;n=n1:n2;x=(n-n0)=0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(时间(n);ylabel(幅度x(n);title(单位脉冲序列); 运行结果如图1-2所示图1-22) 单位阶跃序列 单位阶跃序列表示式为 或 以下三段程序分别用不同的方法来产生单位阶跃序列。例1-4用MATLAB的关系运算式来产生单位阶跃序列。n1=-2;n2=8;n0=0;n=n1:n2;x=n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(时间(n);ylabel(幅度x(n);title(单位阶跃序列);box运行结果如图1-3所示：图1-3例1-5 用zeros和ones函数来产生单位阶跃序列。n1=-2;n2=8;k=0;n=n1:n2;nt=length(n);nk=abs(k-n1)+1;x=zeros(1,nk-1),ones(1,nt-nk+1);绘图部分的程序及作图结果与例1-4相同。.例1-6生成移位的单位阶跃序列。n1=-10;n2=10;n0=4;n=n1:n2;x=(n-n0)=0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(时间(n);ylabel(幅度x(n);title(移位的单位阶跃序列);box程序运行结果如图1-4所示：图1-43) 实指数序列实指数序列的表示式为 x(n)=an 其中a为实数例1-7 编写产生a=1/2和a=2的实指数连续信号和离散序列的程序n1=-10;n2=10;a1=0.5;a2=2;na1=n1:0;x1=a1.na1;na2=0:n2;x2=a2.na2;subplot(2,2,1);plot(na1,x1);(图形窗口)title(实指数信号（a1）);subplot(2,2,3);stem(na1,x1,filled);title(实指数序列（a1）);subplot(2,2,4);stem(na2,x2,filled);title(实指数序列（a1）);box程序运行结果如图1-5所示：图1-54）复指数序列复指数序列的表示式为 x(n)=e(+j)n当=0时，x(n)为实指数序列;当=0时，x(n)为虚指数序列，即ejn=cos(n)+jsin(n)其实部为余弦序列，虚部为正弦序列。例1-8 编写程序产生=-0.1，=0.6的复指数连续信号与离散序列。n1=30;a=-0.1;w=0.6;n=0:n1;x=exp(a+j*w)*n);subplot(2,2,1);plot(n,real(x);title(复指数信号的实部);subplot(2,2,3);stem(n,real(x),filled);title(复指数序列的实部);subplot(2,2,2);plot(n,imag(x);title(复指数信号的虚部);subplot(2,2,4);stem(n,imag(x),filled);title(复指数序列的虚部);box程序运行结果如图1-6所示图1-65）正（余）弦序列正（余）弦序列的表示式为 x(n)=Umsin(0n+)例1-9 已知一时域周期性正弦信号的频率为1Hz，振幅值为1V。编写程序在图形窗口上显示两个周期的信号波形，并对该信号的一个周期进行32点采样获得离散信号。f=1;Um=1;nt=2;N=32;T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=Um*sin(2*f*pi*tn);subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(t);subplot(2,1,2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(n);box程序运行结果如图1-7所示图1-76）矩形波序列MATLAB提供有专门函数square用于产生矩形波。其调用格式如下：x=square(t) ：类似于sin(t),产生周期为2，幅值为1的方波。x=square(t，duty)：产生指定周期的矩形波，其中duty用于指定占空比。将square的参数t换成n，且n取整数，则可以获得矩形序列。例1-10 一个周期性矩形信号频率为5kHz，信号幅度在02V之间，占空比为0.25，。编写程序生成该信号，要求在图形窗口上显示2个周期的信号波形；对信号的一个周期进行16点采样获得离散信号。f=5000;nt=2;N=16;T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=square(2*f*pi*tn,25)+1;subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(t);subplot(2,1,2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(n);box程序运行结果如图1-8所示;三、.实验内容：1、阅读并上机验证实验原理部分的例题程序，理解每一条语句的含义。改变例题中的有关参数（如信号的频率、周期、幅度、显示时间的取值范围、采样点数等），观察对信号波形的影响。2、编写程序，产生以下离散序列：（1）f(n)=(n) (-3n4)（2）f(n)=u(n) (-5n5)（3）f(n)= e(0.1+j1.6)n (0n16)（4）f(n)=3sin(n/4) (0n=0;gn=filter(b,a,x2,xi);subplot(1,2,1);stem(n,hn,k);title(系统的单位序列响应);ylabel(h(n);xlabel(n);axis(0,N,1.1*min(hn),1.1*max(hn);subplot(1,2,2);stem(n,gn,k);title(系统的单位阶跃响应);ylabel(g(n);xlabel(n);axis(0,N,1.1*min(gn),1.1*max(gn);程序运行结果如图2-2所示：图2-24、用MATLAB实现线性卷积1）用函数conv进行卷积运算：求解两个序列的卷积和，关键在于如何确定卷积结果的时宽区间。MATLAB提供的求卷积函数conv默认两个序列的序号均从n=0开始，卷积结果y对应的序列的序号也从n=0开始。例2-3 已知两个序列f1=0.8n （0n20），f2=u(n) （0n10），求两个序列的卷积和。n1=0:20;f1=0.8.n1;subplot(2,2,1);stem(n1,f1,filled);title(f1(n);n2=0:10;N2=length(n2);f2=ones(1,N2);subplot(2,2,2);stem(n2,f2,filled);title(f2(n);y=conv(f1,f2);subplot(2,1,2);stem(y,filled);程序运行结果如图2-3所示：图2-32）非零起始序列的卷积运算：当两个序列不是从0开始时，必须对conv函数稍加扩展。由卷积原理可知，若待卷积的两个序列序号分别为x(n)；nx=nxs:nxf，h(n)；nh=nhs:nhf，则卷积和y(n)的序号起点和终点分别为：nys=nxs+nhs，nyf=nxf+nhf。据此可定义通用卷积函数convu：functiony,ny=convu(h,nh,x,nx)nys=nh(1)+nx(1);nyf=nh(end)+nx(end);y=conv(h,x);ny=nys:nyf;例2-4 已知序列f1=0.5n (0n10)，f2=u(n+2) (-2n10)，求两个序列的卷积和。程序清单如下：n1=0:10;f1=0.5*n1;n2=-2:10;nt=length(n2);f2=ones(1,nt);y,ny=convu(f1,n1,f2,n2);subplot(2,2,1);stem(n1,f1);subplot(2,2,2);stem(n2,f2);subplot(2,1,2);stem(ny,y);程序运行结果如图2-4所示：图2-43）卷积积分的动态过程演示：为了更深入地理解序列卷积的原理，下面提供一段演示卷积和的动态过程的程序。例2-5 动态演示例2-3两个序列f1=0.8n （0n20），f2=u(n) （0nlf1 nf2=0;nf1=lf2-lf1;elseif lf2lf1 nf1=0;nf2=lf1-lf2;else nf2=0;lf1=0;end;lt=m;u=zeros(1,lt),f2,zeros(1,nf2),zeros(1,lt);t1=(-lt+1:2*lt);f1=zeros(1,2*lt),f1,zeros(1,nf1);hf1=fliplr(f1);N=length(hf1);y=zeros(1,3*lt);for k=0:2*lt p=zeros(1,k),hf1(1:N-k); y1=u.*p; yk=sum(y1); y(k+lt+1)=yk; subplot(4,1,1);stem(t1,u); subplot(4,1,2);stem(t1,p); subplot(4,1,3);stem(t1,y1); subplot(4,1,4);stem(k,yk); axis(-20,50,0,5);hold on pause(2);end程序运行结果如图2-5所示：图2-55、离散LSI系统时域响应的求解：MATLAB提供了多种方法求解离散LSI系统的响应：1）用conv函数进行卷积积分，求任意输入的系统零状态响应；例2-6 已知描述某因果系统的差分方程为6y(n)+2y(n-2)=x(n)+3x(n-1)+3x(n-2)+x(n-3)，满足初始条件y(-1)=0，x(-1)=0。在该系统的输入端加一个矩形脉冲序列，其占空比为0.25，一个周期取16个采样点，求该系统的响应。程序清单如下：N=16;n=0:N-1;x=ones(1,N/4),zeros(1,3*N/4);subplot(2,2,1);stem(n,x);a=1,0,1/3,0,;b=1/6,1/2,1/2,1/6;hn=impz(b,a,n);subplot(2,2,2);stem(n,hn);y=conv(x,hn);subplot(2,1,2);stem(y);程序运行结果如图2-6所示：图2-62）用dlsim函数求任意输入的系统零状态响应；例2-7 已知某IIR数字低通滤波器的系统函数为输入两个正弦叠加的序列，求该系统的响应。程序清单如下：nx=0:8*pi;x=sin(nx/2)+sin(10*nx)/3;subplot(3,1,1);stem(nx,x);a=1,-0.34319,0.60439,-0.20407;b=0.1321,0.3963,0.3963,0.1321;nh=0:9;h=impz(b,a,nh);subplot(3,1,2);stem(nh,h);y=dlsim(b,a,x);subplot(3,1,3);stem(y);程序运行结果如图2-7所示图2-73）用filtic和filter函数求任意输入的系统完全响应。例2-8 已知描述某系统的差分方程为y(n)-1.5y(n-1)+0.5 y(n-2)=x(n) n0，满足初始条件y(-1)=4，y(-2)=10，求系统输入为x(n)=(0.25)nu(n)时的零输入、零状态及全响应。解 为了更深入地理解filtic和filter函数的用法，先用经典法求得系统完全响应表达式：，并编写程序绘出其图形，以便与用MATLAB函数求解的结果进行对比。程序清单如下：a=1,-1.5,0.5;b=1;N=20;n=0:N-1;x=0.25.n;x0=zeros(1,N);y01=4,10;xi=filtic(b,a,y01);y0=filter(b,a,x0,xi);xi0=filtic(b,a,0);y1=filter(b,a,x,xi0);y=filter(b,a,x,xi);y2=(1/3)*(1/4).n+(1/2).n+(2/3).*ones(1,N);subplot(2,3,1);stem(n,x);title(输入信号x(n);subplot(2,3,2);stem(n,y0);title(系统的零输入响应);subplot(2,3,3);stem(n,y1);title(系统的零状态响应);subplot(2,2,3);stem(n,y);title(用filter求得的完全响应);subplot(2,2,4);stem(n,y2);title(经典法求得的完全响应);程序运行结果如图2-8所示：图2-8三、实验内容：1、输入并运行例题程序，理解每一条语句的含义。2、已知描述某离散LSI系统的差分方程为2y(n)-3y(n-1)+y(n-2)=x(n-1)，分别用impz和dstep函数、filtic和filter函数两种方法求解系统的单位序列响应和单位阶跃响应。3、编写程序描绘下列序列的卷积波形：（1）f1(n=u(n),f2(n)=u(n-2), (0n10)（2）x(n)=sin(n/2),h(n)=(0.5)n (-3n4)4、已知某离散LSI系统的单位序列响应为h(n)=3(n-3)+0.5(n-4)+0.2(n-5)+0.7(n-6)-0.8(n-7)求输入为x(n)=e-0.5nu(n)时的系统响应。5、已知描述某离散LSI系统的差分方程为y(n)=0.7y(n-1)+2x(n)-x(n-2)，求输入为x(n)=u(n-3)时的系统响应。四、实验预习：1、认真阅读实验原理部分，明确实验目的，复习有关离散LSI系统的理论知识。2、读懂实验原理部分的例题程序，熟悉与本实验有关的MATLAB函数。3、根据实验内容预先编写实验程序，并思考本实验提出的有关MATLAB函数在调用时应注意哪些问题。五、实验报告：1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、列出本实验提出的有关MATLAB函数在调用时应注意的问题。实验3 离散LSI系统的频域分析一、实验目的1、加深对离散系统变换域分析z变换的理解，掌握使用MATLAB进行z变换和逆z变换的常用函数的用法。2、了解离散系统的零极点与系统因果性和稳定性的关系，熟悉使用MATLAB进行离散系统的零极点分析的常用函数的用法。3、加深对离散系统的频率响应特性基本概念的理解，掌握使用MATLAB进行离散系统幅频响应和相频响应特性分析的常用方法。二、实验原理1、z变换和逆z变换（1）用ztrans函数求无限长序列的z变换。该函数只给出z变换的表达式，而没有给出收敛域。另外，由于这一函数还不尽完善，有的序列的z变换还不能求出，逆z变换也存在同样的问题。例7-1 求以下各序列的z变换 x1(n)=an x2(n)=n x3(n)=n(n-1)/2 x4(n)=ejon x5(n)=1/n(n-1)程序清单如下：syms w0 n z a;x1=n*an;X1=ztrans(x1)x2=sin(w0*n);X2=ztrans(x2)x3=exp(-a*n)*sin(w0*n);X3=ztrans(x3)程序运行结果如下：X1 =z/a/(z/a-1) X2 =z/(z-1)2X3 =1/2*z*(z+1)/(z-1)3-1/2*z/(z-1)2X4 =z/exp(i*w0)/(z/exp(i*w0)-1)X5 =z/(z-1)-ztrans(1/n,n,z)（2）用iztrans函数求无限长序列的逆z变换。例3-2 求下列函数的逆z变换。 程序清单如下：syms n z a;X1=z/(z-a);x1=iztrans(X1)X2= z/(z-a)2;x2=iztrans(X2)X3=z/z-exp(j*w0);x3=iztrans(X3)X4=(1-z-3)/(1-z-1);x4=iztrans(X4)程序运行结果如下：x1 =1x2 =an*nx3 =1/2*n2-1/2*nx4 =iztrans(1-z(-n)/(1-1/z),z,n)2、离散系统的零极点分析（系统极点位置对系统响应的影响）例3-3 研究z右半平面的实数极点对系统的影响。已知系统的零极点增益模型分别为：求这些系统的零极点分布图以及系统的单位序列响应，判断系统的稳定性。程序清单如下：z1=0;p1=0.85;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(极点在单位圆内);subplot(3,2,2);impz(b1,a1,20);z2=0;p2=1;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(极点在单位圆上);subplot(3,2,4);impz(b2,a2,20);z3=0;p3=1.5;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(极点在单位圆外);subplot(3,2,6);impz(b3,a3,20);程序运行结果如图3-1所示。由图可见，这三个系统的极点均为实数且处于z平面的右半平面。由图可知，当极点位于单位圆内，系统的单位序列响应随着频率的增大而收敛；当极点位于单位圆上，系统的单位序列响应为等幅振荡；当极点位于单位圆外，系统的单位序列响应随着频率的增大而发散。由此可知系统1、2为稳定系统。图3-1例3-4 研究z左半平面的实数极点对系统的影响。已知系统的零极点增益模型分别为：求这些系统的零极点分布图以及系统的单位序列响应，判断系统的稳定性。程序清单如下：z1=0;p1=-0.85;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(极点在单位圆内);subplot(3,2,2);impz(b1,a1,20);z2=0;p2=-1;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(极点在单位圆上);subplot(3,2,4);impz(b2,a2,20);z3=0;p3=-1.5;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(极点在单位圆外);subplot(3,2,6);impz(b3,a3,20);程序运行结果如图3-2所示。由图可见，这三个系统的极点均为实数且处于z平面的左半平面。由图可知，当极点位于单位圆内，系统的单位序列响应随着频率的增大而收敛；当极点位于单位圆上，系统的单位序列响应为等幅振荡；当极点位于单位圆外，系统的单位序列响应随着频率的增大而发散。由此可知系统1、2为稳定系统。图3-2例3-5 研究z右半平面的复数极点对系统响应的影响已知系统的零极点增益模型分别为：求这些系统的零极点分布图以及系统的单位序列响应，判断系统的稳定性。程序清单如下：z1=0.3,0;p1=0.5+0.7j,0.5-0.7j;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(极点在单位圆内);subplot(3,2,2);impz(b1,a1,20);z2=0.3,0;p2=0.6+0.8j,0.6-0.8j;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(极点在单位圆上);subplot(3,2,4);impz(b2,a2,20);z3=0.3,0;p3=1+j,1-j;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(极点在单位圆外);subplot(3,2,6);impz(b3,a3,20);程序运行结果如图3-3所示。由图可见，这三个系统的极点均为复数且处于z平面的右半平面。由图可知，当极点位于单位圆内，系统的单位序列响应随着频率的增大而收敛；当极点位于单位圆上，系统的单位序列响应为等幅振荡；当极点位于单位圆外，系统的单位序列响应随着频率的增大而发散。由此可知系统1、2为稳定系统。图3-3由以上三例可得结论：系统只有在极点处于单位圆内才是稳定的。例3-6 已知某离散时间系统的系统函数为求该系统的零极点及零极点分布图，并判断系统的因果稳定性。程序清单如下：b=0.2,0.1,0.3,0.1,0.2;a=1,-1.1,1.5,-0.7,0.3;rz=roots(b)rp=roots(a)subplot(2,1,1);zplane(b,a);title(系统的零极点分布图);subplot(2,1,2);impz(b,a,20);title(系统的单位序列响应);xlabel(n);ylabel(h(n);程序运行结果如下：rz = -0.5000 + 0.8660i -0.5000 - 0.8660i 0.2500 + 0.9682i 0.2500 - 0.9682irp = 0.2367 + 0.8915i 0.2367 - 0.8915i 0.3133 + 0.5045i 0.3133 - 0.5045i图3-4由零极点分布图可见，该系统的所有极点均在单位圆内，因此该系统是一个因果稳定系统。3、离散系统的频率响应（1）离散系统的频率响应的基本概念已知稳定系统传递函数的零极点增益模型为则系统的频响函数为其中，系统的幅频特性为系统的相频特性为由以上各式可见，系统函数与频率响应有着密切的联系。适当地控制系统函数的零极点分布，可以改变离散系统的频响特性： 在原点（z=0）处的零点或极点至单位圆的距离始终保持不变，其值|ej|=1，所以，对幅度响应不起作用； 单位圆附近的零点对系统幅度响应的谷值位置及深度有明显影响； 单位圆内且靠近单位圆附近的极点对系统幅度的峰值位置及大小有明显的影响。（2）系统的频响特性分析例3-7 已知某离散时间系统的系统函数为求该系统在0频率范围内的相对幅频响应与相频响应。程序清单如下：b=0.1321,0,-0.3963,0,0.3963,0,-0.1321;a=1,0,0.34319,0,0.60439,0,0.20407;freqz(b,a);程序运行结果如图3-5所示。该系统是一个IIR数字带通滤波器。其中幅频特性采用归一化的相对幅度值，以分贝（dB）为单位。图3-5例3-8 已知某离散时间系统的系统函数为求该系统在0频率范围内的绝对幅频响应与相频响应。程序清单如下：b=0.2,0.1,0.3,0.1,0.2;a=1,-1.1,1.5,-0.7,0.3;n=(0:500)*pi/500;h,w=freqz(b,a,n);subplot(2,1,1);plot(n/pi,abs(h);grid;axis(0,1,1.1*min(abs(h),1.1*max(abs(h);xlabel(omega/pi);ylabel(幅度);subplot(2,1,2);plot(n/pi,angle(h);grid;axis(0,1,1.1*min(angle(h),1.1*max(angle(h);xlabel(omega/pi);ylabel(相位);程序运行结果如图3-6所示。该系统为一低通滤波器。图3-6例3-9 已知某离散时间系统的系统函数为求该系统在0频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及零极点分布图。程序清单如下：b=0.1,-0.4,0.4,-0.1;a=1,0.3,0.55,0.2;n=(0:500)*pi/500;h,w=freqz(b,a,n);db=20*log10(abs(h);subplot(2,2,1);plot(w/pi,abs(h);grid;axis(0,1,1.1*min(abs(h),1.1*max(abs(h);title(幅频特性（V）);xlabel(omega/pi);ylabel(幅度(V);subplot(2,2,2);plot(w/pi,angle(h);grid;axis(0,1,1.1*min(angle(h),1.1*max(angle(h);xlabel(omega/pi);ylabel(相位);title(相频特性);subplot(2,2,3);plot(w/pi,db);gridaxis(0,1,-100,5);title(幅频特性（dB）);subplot(2,2,4);zplane(b,a);title(零极点分布);程序运行结果如图3-7所示：图3-7（3） 一个求解频率响应的实用函数。在实际使用freqz进行离散系统频响特性分析时。通常需要求解幅频响应、相频响应、群时延，幅频响应又分为绝对幅频和相对幅频两种表示方法。下面定义函数freqz_m，利用该函数，可方便求出上述各项。freqz_m函数定义如下：functiondb,mag,pha,grd，w=freqz_m(b,a);H,w=freqz(b,a,1000,whole);H=(H(1:501);w=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);例3-10 已知某离散时间系统的系统函数为求该系统在0频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。程序清单如下：b=0.1321,0,0.3963,0,0.3963,0,0.1321;a=1,0,-0.34319,0,0.60439,0,-0.20407;db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);gridaxis(0,1,1.1*min(mag),1.1*max(mag);title(幅频特性（V）);xlabel(omega/pi);ylabel(幅度(V);subplot(2,2,2);plot(w/pi,pha);grid;axis(0,1,1.1*min(pha),1.1*max(pha);xlabel(omega/pi);ylabel(相位);title(相频特性);subplot(2,2,3);plot(w/pi,db);gridaxis(0,1,-100,5);title(幅频特性（dB）);subplot(2,2,4);plot(w/pi,grd);gridaxis(0,1,0,10)title(群时延);程序运行结果如图3-8所示：图3-8三、实验内容：1、输入并运行例题程序，理解每条语句的含义。2、求以下各序列的z变换：3、求下列函数的逆z变换4、求一下系统函数所描述的离散系统的零极点分布图，并判断系统的稳定性（1）（2）5、已知某离散时间系统的系统函数为求该系统在0频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。四、实验预习：1、认真阅读实验原理部分，明确实验目的，复习有关离散LSI系统频率响应的理论知识。2、读懂实验原理部分的例题程序，熟悉与本实验有关的MATLAB函数。3、根据实验内容预先编写实验程序，并思考本实验提出的有关MATLAB函数在调用时应注意哪些问题。4、预习思考题： 系统函数零极点的位置与系统单位序列响应有何关系？ 离散系统的零极点对系统幅频响应有何影响？五、实验报告：1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、列出本实验提出的有关MATLAB函数在调用时应注意的问题。 3、给出预习思考题答案。实验4 DFS、DFT与FFT一、实验目的：1、加深对周期序列DFS、有限长序列DFT和FFT的基本概念及理论的理解。2、掌握用MATLAB语言求解DFS、DFT、FFT的以及相应反变换的方法。3、观察周期序列的重复周期数对序列频谱特性的影响。 二、实验原理：1、周期序列的离散傅里叶级数（DFS）（1）DFS的基本概念：离散时间序列x(n)满足x(n)= x(n+rN)，称为离散周期序列，用表示。其中，N为信号的周期，x(n)称为离散周期序列的主值。周期序列可以用离散傅里叶级数（DFS）表示：其中，是周期序列DFS第k次谐波分量的系数，也称为周期序列的频谱，可表示为以上两式也是周期序列的一对傅里叶级数变换对。令，以上DFS变换对又可以写成：与连续周期信号的傅里叶级数相比，周期序列的离散傅里叶级数有以下特点： 连续周期信号的傅里叶级数由无穷多个与基波频率成整数倍的谐波分量叠加而成，而周期为N的周期序列的傅里叶级数仅有N个独立的谐波分量。 周期序列的频谱也是一个以N为周期的周期序列。（2）周期序列的DFS和IDFS例4-1 已知一个周期性矩形序列的脉冲宽度占整个周期的1/4，一个周期的采样点数为16点，编程显示3个周期的序列波形，并： 用傅里叶级数求信号的幅度谱和相位谱。 求傅里叶级数逆变换的图形，并与原序列进行比较。程序清单如下：N=16;xn=ones(1,N/4),zeros(1,3*N/4);xn=xn,xn,xn;n=0:3*N-1;k=0:3*N-1;Xk=xn*exp(-j*2*pi/N).(n*k);x=(Xk*exp(j*2*pi/N).(n*k)/N;subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,3*N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFS|X(k)|);axis(-1,3*N,1.1*min(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,3*N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,3*N,1.1*min(angle(Xk),1.1*max(angle(Xk);程序运行结果如图4-1所示：图4-1由离散傅里叶级数逆变换图形可见，与原序列相比，幅度扩大了32倍。这是因为周期序列为原主值序列周期的3倍，做逆变换时未做处理。可将逆变换程序改为：x=(Xk*exp(j*2*pi/N).(n*k)/3*3*N;由上例可见，周期序列的DFS和IDFS是依据变换公式编程的，无论信号序列如何变化，求解的公式总是一样的。因此，可将其编写成通用子程序： 离散傅里叶级数变换通用子程序dfs.mfunctionXk=dfs(xn,N)n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;Xk=xn*WN.nk;离散傅里叶级数逆变换通用子程序idfs.mfunctionxn=idfs(Xk,N)n=0:N-1;k=0:N-1;WN=exp(j*2*pi/N);nk=n*k;xn=(Xk*WN.nk)/N;例4-2 利用上述两个子程序，重做例4-1程序清单如下：N=16;xn=ones(1,N/4),zeros(1,3*N/4);Xk=dfs(xn,N);x=idfs(Xk,N);subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,3*N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFS|X(k)|);axis(-1,3*N,1.1*min(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,3*N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,3*N,1.1*min(angle(Xk),1.1*max(angle(Xk);程序运行结果如图4-2所示。由于子程序仅适用于对主值区间进行变换，周期次数无法传递给子程序，因此程序执行结果仅显示一个周期的变换情况。图4-2（2）周期重复次数对序列频谱的影响理论上讲，周期序列不满足绝对可积条件，因此不能用傅里叶级数来表示。实际处理时可先取K个周期进行处理，然后令K趋于无穷大，分析其极限情况。根据这一分析思路，可以观察序列由非周期到周期变化时，频谱由连续谱逐渐向离散谱过渡的过程。例4-3 已知一矩形序列的脉冲宽度占整个周期的1/2，一个周期的采样点数为10点，用DFS求序列的重复周期数分别为1、4、7、10时的幅度谱。程序清单如下：xn=ones(1,5),zeros(1,5);Nx=length(xn);Nw=1000;dw=2*pi/Nw;k=floor(-Nw/2+0.5):(Nw/2+0.5);for r=0:3 K=3*r+1; nx=0:(K*Nx-1); x=xn(mod(nx,Nx)+1); Xk=x*(exp(-j*dw*nx*k)/K;subplot(4,2,2*r+1);stem(nx,x);axis(0,K*Nx-1,0,1.1);ylabel(x(n);subplot(4,2,2*r+2);plot(k*dw,abs(Xk);axis(-4,4,0,1.1*max(abs(Xk);ylabel(X(k);end程序运行结果如图4-3所示。由图可见，序列的重复周期数越多，频谱约是向几个频点集中。当序列的周期数趋于无穷大时，频谱转化为离散谱。图4-32、离散傅里叶变换（DFT）（1）DFT与IDFT在实际中常常使用有限长序列。如果有限长序列为x(n)，则该序列的离散傅里叶变换对可表示为从离散傅里叶变换定义式可以看出，有限长序列在时域上是离散的，在频域上也是离散的，式中，即仅在单位圆上N个等间距的点上取值，这为使用计算机进行处理带来了方便。由有限长序列的傅里叶变换和逆变换定义可知，DFT和DFS的变换公式非常相似，因此，在程序编写上也基本一致。例4-4 已知x(n)=0,1,2,3,4,5,6,7，求其DFT和IDFT。要求： 画出序列傅里叶变换对应的|X(k)|和argX(k)图形。 画出x(n)图形，并与IDFTX(k)图形进行比较。程序清单如下：xn=0,1,2,3,4,5,6,7;N=length(xn);n=0:N-1;k=0:N-1;Xk=xn*exp(-j*2*pi/N).(n*k);x=(Xk*exp(j*2*pi/N).(n*k)/N;subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFT|X(k)|);axis(-1,N,1.1*min(x),1.1*max(x);subplo