《数字信号处理》实验指导书(完整)

I《数字信号处理》实验指导书通信教研室安阳工学院二零零九年三月1第1章系统响应及系统稳定性1.1实验目的学会运用MATLAB求解离散时间系统的零状态响应；学会运用MATLAB求解离散时间系统的单位取样响应；学会运用MATLAB求解离散时间系统的卷积和。1.2实验原理及实例分析1.2.1离散时间系统的响应离散时间LTI系统可用线性常系数差分方程来描述，即MjjNiijnxbinya00)()(（1-1）其中，ia（0i，1，…，N）和jb（0j，1，…，M）为实常数。MATLAB中函数filter可对式（13-1）的差分方程在指定时间范围内的输入序列所产生的响应进行求解。函数filter的语句格式为y=filter(b,a,x)其中，x为输入的离散序列；y为输出的离散序列；y的长度与x的长度一样；b与a分别为差分方程右端与左端的系数向量。【实例1-1】已知某LTI系统的差分方程为)1(2)()2(2)1(4)(3nxnxnynyny试用MATLAB命令绘出当激励信号为)()2/1()(nunxn时，该系统的零状态响应。解：MATLAB源程序为a=[3-42];b=[12];n=0:30;2x=(1/2).^n;y=filter(b,a,x);stem(n,y,'fill'),gridonxlabel('n'),title('系统响应y(n)')程序运行结果如图1-1所示。1.2.2离散时间系统的单位取样响应系统的单位取样响应定义为系统在)(n激励下系统的零状态响应，用)(nh表示。MATLAB求解单位取样响应可利用函数filter，并将激励设为前面所定义的impDT函数。例如，求解实例1-1中系统的单位取样响应时，MATLAB源程序为a=[3-42];b=[12];n=0:30;x=impDT(n);h=filter(b,a,x);stem(n,h,'fill'),gridonxlabel('n'),title('系统单位取样响应h(n)')图1-1实例1-1系统的零状态响应3程序运行结果如图1-2所示。MATLAB另一种求单位取样响应的方法是利用控制系统工具箱提供的函数impz来实现。impz函数的常用语句格式为impz(b,a,N)其中，参数N通常为正整数，代表计算单位取样响应的样值个数。【实例1-2】已知某LTI系统的差分方程为)1(2)()2(2)1(4)(3nxnxnynyny利用MATLAB的impz函数绘出该系统的单位取样响应。解：MATLAB源程序为a=[3-42];b=[12];n=0:30;impz(b,a,30),gridontitle('系统单位取样响应h(n)')程序运行结果如图1-3所示，比较图1-2和图1-3，不难发现结果相同。图1-2实例1-1的系统单位取样响应41.2.3离散时间信号的卷积和运算由于系统的零状态响应是激励与系统的单位取样响应的卷积，因此卷积运算在离散时间信号处理领域被广泛应用。离散时间信号的卷积定义为mmnhmxnhnxny)()()(*)()(（1-2）可见，离散时间信号的卷积运算是求和运算，因而常称为“卷积和”。MATLAB求离散时间信号卷积和的命令为conv，其语句格式为y=conv(x,h)其中，x与h表示离散时间信号值的向量；y为卷积结果。用MATLAB进行卷积和运算时，无法实现无限的累加，只能计算时限信号的卷积。例如，利用MALAB的conv命令求两个长为4的矩形序列的卷积和，即)]4()([*)]4()([)(nunununung，其结果应是长为7（4+4-1=7）的三角序列。用向量[1111]表示矩形序列，MATLAB源程序为x1=[1111];x2=[1111];g=conv(x1,x2)g=1234321图1-3系统单位取样响应5如果要绘出图形来，则利用stem命令，即n=1:7;stem(n,g,'fill'),gridon,xlabel('n')程序运行结果如图1-4所示。对于给定函数的卷积和，我们应计算卷积结果的起始点及其长度。两个时限序列的卷积和长度一般等于两个序列长度的和减1。【实例1-3】已知某系统的单位取样响应为88.0nununhn，试用MATLAB求当激励信号为)4()()(nununx时，系统的零状态响应。解：MATLAB中可通过卷积求解零状态响应，即)(*)(nhnx。由题意可知，描述)(nh向量的长度至少为8，描述)(nx向量的长度至少为4，因此为了图形完整美观，我们将)(nh向量和)(nx向量加上一些附加的零值。MATLAB源程序为nx=-1:5;%x(n)向量显示范围(添加了附加的零值)nh=-2:10;%h(n)向量显示范围(添加了附加的零值)x=uDT(nx)-uDT(nx-4);h=0.8.^nh.*(uDT(nh)-uDT(nh-8));y=conv(x,h);ny1=nx(1)+nh(1);%卷积结果起始点%卷积结果长度为两序列长度之和减1,即0到(length(nx)+length(nh)-2)%因此卷积结果的时间范围是将上述长度加上起始点的偏移值ny=ny1+(0:(length(nx)+length(nh)-2));图1-4卷积结果图6subplot(311)stem(nx,x,'fill'),gridonxlabel('n'),title('x(n)')axis([-41603])subplot(312)stem(nh,h','fill'),gridonxlabel('n'),title('h(n)')axis([-41603])subplot(313)stem(ny,y,'fill'),gridonxlabel('n'),title('y(n)=x(n)*h(n)')axis([-41603])程序运行结果如图1-5所示。图1-5利用卷积和法求解系统的零状态响应71.3编程练习1.试用MATLAB命令求解以下离散时间系统的单位取样响应，并判断系统的稳定性。（1）)1()()2()1(4)(3nxnxnynyny（2）)()2(10)1(6)(25nxnynyny2.已知某系统的单位取样响应为10)87(nununhn，试用MATLAB求当激励信号为)5()()(nununx时，系统的零状态响应。附：1.单位取样序列单位取样序列)(n，也称为单位冲激序列，定义为)0()0(01)(nnn（12-1）要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n=0处是取确定的值1。在MATLAB中，冲激序列可以通过编写以下的impDT.m文件来实现，即functiony=impDT(n)y=(n==0);%当参数为0时冲激为1,否则为0调用该函数时n必须为整数或整数向量。【实例2-1】利用MATLAB的impDT函数绘出单位冲激序列的波形图。解：MATLAB源程序为n=-3:3;x=impDT(n);stem(n,x,'fill'),xlabel('n'),gridontitle('单位冲激序列')axis([-33-0.11.1])8程序运行结果如图12-1所示。2.单位阶跃序列单位阶跃序列)(nu定义为)0()0(01)(nnnu（12-2）在MATLAB中，冲激序列可以通过编写uDT.m文件来实现，即functiony=uDT(n)y=n=0;%当参数为非负时输出1调用该函数时n也同样必须为整数或整数向量。【实例2-2】利用MATLAB的uDT函数绘出单位阶跃序列的波形图。解：MATLAB源程序为n=-3:5;x=uDT(n);stem(n,x,'fill'),xlabel('n'),gridontitle('单位阶跃序列')axis([-35-0.11.1])图2-1单位冲激序列9程序运行结果如图12-2所示。图2-2单位阶跃序列10第2章z变换及离散时间LTI系统的z域分析2.1实验目的学会运用MATLAB求离散时间信号的z变换和z反变换；学会运用MATLAB分析离散时间系统的系统函数的零极点；学会运用MATLAB分析系统函数的零极点分布与其时域特性的关系；学会运用MATLAB进行离散时间系统的频率特性分析。2.2实验原理及实例分析2.2.1z正反变换序列nx的z变换定义为nnznxnxzXZ（2-1）其中，符号Z表示取z变换，z是复变量。相应地，单边z变换定义为0nnznxnxzXZ（2-2）MATLAB符号数学工具箱提供了计算离散时间信号单边z变换的函数ztrans和z反变换函数iztrans，其语句格式分别为Z=ztrans(x)x=iztrans(z)上式中的x和Z分别为时域表达式和z域表达式的符号表示，可通过sym函数来定义。【实例2-1】试用ztrans函数求下列函数的z变换。（1）)()cos()(nunanxn；（2）)(])2(2[)(11nunxnn。解：（1）z变换MATLAB源程序为x=sym('a^n*cos(pi*n)');11Z=ztrans(x);simplify(Z)ans=z/(z+a)（2）z变换MATLAB源程序为x=sym('2^(n-1)-(-2)^(n-1)');Z=ztrans(x);simplify(Z)ans=z^2/(z-2)/(z+2)【实例2-2】试用iztrans函数求下列函数的z反变换。（1）65198)(2zzzzX（2）32)2)(1()12112()(zzzzzzX解：（1）z反变换MATLAB源程序为Z=sym('(8*z-19)/(z^2-5*z+6)');x=iztrans(Z);simplify(x)ans=-19/6*charfcn[0](n)+5*3^(n-1)+3*2^(n-1)其中，charfcn[0](n)是)(n函数在MATLAB符号工具箱中的表示，反变换后的函数形式为)()2335()(619)(11nunnxnn。（2）z反变换MATLAB源程序为Z=sym('z*(2*z^2-11*z+12)/(z-1)/(z-2)^3');x=iztrans(Z);simplify(x)ans=-3+3*2^n-1/4*2^n*n-1/4*2^n*n^2其函数形式为)()241241233()(2nunnnxnnn。如果信号的z域表示式)(zX是有理函数，进行z反变换的另一个方法是对)(zX进行部分分式展开，然后求各简单分式的z反变换。设)(zX的有理分式表示为12)()(1)(221122110zAzBzazazazbzbzbbzXnnmm（2-3）MATLAB信号处理工具箱提供了一个对)(zX进行部分分式展开的函数residuez，其语句格式为[R,P,K]=residuez(B,A)其中，B，A分别表示X(z)的分子与分母多项式的系数向量；R为部分分式的系数向量；P为极点向量；K为多项式的系数。若X(z)为有理真分式，则K为零。【实例2-3】试用MATLAB命令对函数321431818)(zzzzX进行部分分式展开，并求出其z反变换。解：MATLAB源程序为B=[18];A=[18,3,-4,-1];[R,P,K]=residuez(B,A)R=0.36000.24000.4000P=0.5000-0.3333-0.3333K=[]从运行结果可知，32pp，表