Z变换和系统频域特性的MATLAB实现

1、小实验4 Z变换和系统频域特性的MATLAB实现实验目的学习通过Z变换来分析离散系统的频率响应，并用MATLAB实现。加深对系统的零、极点分布概念的理解。2. 实例分析2.1通过Z变换分析求解系统的冲激响应已知用线性常系数差分方程：表示的线性时不变系统其系统函数为： （3-1）上式为两个关于的多项式之比，即为有理分式。同时，式可以表示成部分分式的形式：则可以通过所熟悉的常见序列的Z变换形式求得的Z反变换，从而求得系统的冲激响应函数。MATLAB提供了一个内部函数residuez( )，来计算有理多项式的留数和直接项，residuez( )函数有几种调用方式：（1）R,p,C = residue

2、z(b,a)，在已知以分子行向量b和分母行向量a下，得到列向量R含有留数，列向量p是极点位置，行向量C包含直接项；（2）b,a = residuez(R,p,C)，将部分分式展开式转换到分子行向量b和分母行向量a.MATLAB还提供了一个内部函数impz(b,a,N)，在已知分子行向量b和分母行向量a下，计算N点的单位冲激响应。例3.1 求系统： 的单位冲激响应。解：由上述基本原理和MATLAB提供的函数，将下列指令编辑到 “exe3impz.m” 文件中，可以得到所求系统的冲激响应。% exe3impz.m h(n)求解b=0.001836, 0.007344, 0.011016, 0.00

3、7374, 0.001836;a=1, -3.0544, 3.8291, -2.2925, 0.55075;h,n = impz(b,a,40);stem(n,h,.);ylabel(hn);grid on 运行“exe3impz.m” 文件将产生如图3-1所示的序列。图4-1 系统冲激响应例3.2 求的Z反变换。解：由上述基本原理和MATLAB提供的函数，将下列指令编辑到 “exe3zinver.m”文件中，来求Z反变换。% exe3zinver.m z反变换求解b=0,1,0;a=3,-4,1;R,p,C=residuez(b,a)b,a=residuez(R,p,C)运行结果如下：R =

4、 0.5000 -0.5000p = 1.0000 0.3333C = 0b = -0.0000 0.3333 0a = 1.0000 -1.3333 0.3333因此得到因式分解后的，所以Z反变换的结果为。2.2通过Z变换分析离散系统的频率响应和零、极点分布将式（3-1）表示的的分子、分母进行因式分解，可采用根的形式表示多项式，即其中，为分子多项式的根，称为系统函数的零点，为分母多项式的根，称为系统函数的极点，A为比例 。这样系统函数就表示成了零极点的形式，零极点在Z平面的位置刻画了系统很重要的特性，可以通过系统函数零极点位置估算出系统函数的频率响应，进而判断系统的滤波特性，这是一种非常实用

5、的方法，也称为频率响应的几何确定法。根据频率响应的定义，将代入上式，系统频率响应为：MATLAB提供了一个内部函数zplane(b,a)，在已知系统函数的分子行向量b和分母行向量a下，画出极点和零点。MATLAB还提供了一个内部函数freqz( )来求系统的频率响应，freqz( )函数有几种调用方式：（1）H,w=freqz(b,a,N)，在已知分子行向量b和分母行向量a情况下得到N点的频率向量w和N点的系统复频率响应向量H，频率响应在单位圆的上半圆的N个等分点上求值；（2）H,w=freqz(b,a,N,whole)，在已知分子行向量b和分母行向量a 情况下得到N点的频率向量w和N点的系统

6、复频率响应向量H，频率响应用环绕整个单位圆的N个等分点求值；（3）H=freqz(b,a,w)，在已知分子行向量b和分母行向量a情况下得到在频率向量w上的频率响应向量H。例3.3 画出因果系统的零极点图和频率响应曲线，并求出脉冲响应。解：由上述的基本原理和MATLAB提供的函数，将下列指令编辑到 “exe3freqz.m”文件中，即可求得系统的零极点图和频率响应。% exe3freqz.m 系统求解b=1,0;a=1,-0.9; subplot(2,2,1);zplane(b,a);H,w=freqz(b,a,100);magH=abs(H);phaH=angle(H); subplot(2,

7、2,3);plot(w/pi,magH);grid onxlabel(Frequency in Pi uints);ylabel(Magnitude);title(Magnitude response); subplot(2,2,4);plot(w/pi,phaH/pi);xlabel(Frequency in Pi uints);ylabel(Phase);title(Phase response); h,n=impz(b,a,40);subplot(2,2,2);stem(n,h,.); ylabel(hn);grid on运行“exe3freqz.m” 文件将产生如图3-2所示的序列。图

8、4-2 系统零极点及频率响应例3.4定性画出的零极点图和频率响应曲线。解：由上述基本原理和MATLAB提供的函数，将下列指令编辑到 “exe3freqz2.m” 文件中，即可求得系统的零极点图和频率响应。% exe3freqz2.m 系统求解N=8;b=1,0,0,0,0,0,0,0,-1;a=1,0,0,0,0,0,0,0,0; subplot(2,1,1);zplane(b,a);title(零点-极点图); H,w=freqz(b,a,100);magH=abs(H);phaH=angle(H);subplot(2,2,3);plot(w/pi,magH);grid onxlabel(F

9、requency in Pi uints);ylabel(Magnitude);title(Magnitude response); subplot(2,2,4);plot(w/pi,phaH/pi);xlabel(Frequency in Pi uints);ylabel(Phase);title(Phase response);运行“exe3freqz2.m” 文件将产生如图3-3所示的序列。图4-3 零极点及频率响应图2.3 由求解零极点及由零极点得到的表达式例3.5 已知系统函数有5个零点，一个在，另外4个是及其共轭与镜像，求的系数并进行验证。解：MATLAB提供了一个内部函数tf2z

10、p( )，来求的零极点和增益；提供了一个内部函数zp2tf( )，在零极点已知时求和的系数；提供了一个内部函数roots( )，来求一个多项式的根；提供了一个内部函数poly( )，由给定的根求出相应多项式的系数。因此，将下列指令编辑到 “exe2zphz.m” 文件中，可求解和验证。% exe3zphz.m 零极点与系统函数的转换B=poly(0.5*exp(j*pi/3),0.5*exp(-j*pi/3),2*exp(j*pi/3),2*exp(-j*pi/3),0.8);L=length(B);A=zeros(1,L);A(1)=1;Z,P,K=tf2zp(B,A);sort(Z)b,a

11、=zp2tf(Z,P,K)Z1=roots(B)运行“exe3zphz.m” 文件结果如下，可自行比较验证。ans = 0.2500 - 0.4330i 0.2500 + 0.4330i 0.8000 1.0000 - 1.7321i 1.0000 + 1.7321ib = 1.0000 -3.3000 7.2500 -6.7000 3.0000 -0.8000a = 1 0 0 0 0 0Z1 = 1.0000 + 1.7321i 1.0000 - 1.7321i 0.8000 0.2500 + 0.4330i 0.2500 - 0.4330i 3. 实验内容3.1求系统的零极点和频率响应。

12、b=0.0528, 0.0797, 0.1295, 0.1295, 0.797, 0.0528;a=1, -1.8107, 2.494, -1.8801, 0.9537,-0.2336;subplot(2,2,1);zplane(b,a);H,w=freqz(b,a,1000);magH=abs(H);phaH=angle(H); subplot(2,2,3);plot(w/pi,magH);grid onxlabel(Frequency in Pi uints);ylabel(Magnitude);title(Magnitude response); subplot(2,2,4);plot(

13、w/pi,phaH/pi);xlabel(Frequency in Pi uints);ylabel(Phase);title(Phase response); h,n=impz(b,a,40);subplot(2,2,2);stem(n,h,.); ylabel(hn);grid on3.2 求出矩形序列的零极点和频率响应。3.3 已知离散系统差分方程为：求其Z变换，画出零极点图，并判断系统的稳定性。b=1,-0.5,0;a=1,0.75,0.125;subplot(2,2,1);zplane(b,a);H,w=freqz(b,a,100);magH=abs(H);phaH=angle(H)

14、;极点在单位圆内，所以系统是稳定的4. 思考题 试从理论上计算分析例子3.4的零极点分布和幅频响应，并把结果同MATLAB程序运行结果对比。附录资料：MATLAB Cell函数使用技巧谈谈MATLAB中cell函数如果p为一个数，那么h(1)=p,是没有问题的。如果p为一个向量，那么h(1,:)=p是没有问题的。如果p是一个矩阵的话，上面的两种赋值方法都是会有错误的。那么要如何处理呢？这时就用到了cell数据类型了。cell的每个单元都可以存储任何数据，比如传递函数等。当然，存储矩阵更是没有问题的了。但是用cell数据类型之前，要先初始化。a=cell(n,m)那么就把a初始化为一个n行m列的

15、空cell类型数据。如何赋值呢？a1,1=rand(5)那么a的1行1列的单元中存储的就是一个随机的55的方阵了。那么要用第一个单元中的方阵中的某个值呢？可以如下引用：a1,1(2,3)就可以了，引用cell单元时要用,再引用矩阵的某个数据就要用()了。cell单元中的每个单元都是独立的，可以分别存储不同大小的矩阵或不同类型的数据。下面举个例子：a=cell(2,2);%预分配a1,1=cellclass;a1,2=1 2 2;a2,1=a,b,c;a2,2=9 5 6;a1,1ans =cellclassa1,2ans = 1 2 2a2,:ans =abcans = 9 5 6 b=a1,

16、1b =cellclass元胞数组：元胞数组是MATLAB的一种特殊数据类型，可以将元胞数组看做一种无所不包的通用矩阵，或者叫做广义矩阵。组成元胞数组的元素可以是任何一种数据类型的常数或者常量，每一个元素也可以具有不同的尺寸和内存占用空间，每一个元素的内容也可以完全不同，所以元胞数组的元素叫做元胞（cell）。和一般的数值矩阵一样，元胞数组的内存空间也是动态分配的。（1）元胞数组的创建 a=matlab,20;ones(2,3),1:10a = matlab 20 2x3 double 1x10 double b=matlab,20;ones(2,3),1:10b = matlab 20 2x

17、3 double 1x10 double c=10c = 10c(1,2)=2c = 10 2c(2,2)=5c = 10 2 5isequal(a,b)ans = 1whosName Size Bytes Class Attributesa 2x2 388 cell ans 1x1 1 logical b 2x2 388 cell c 2x2 208 cell 用cell函数创建元胞数组，创建的数组为空元胞。cell函数创建空元胞数组的主要目的是为数组预先分配连续的存储空间，节约内存占用，提高执行效率。 a=cell(1)a = b=cell(1,2)b = c=cell(3,3)c = d

18、=cell(2,2,2)d(:,:,1) = d(:,:,2) = whosName Size Bytes Class Attributesa 1x1 4 cell ans 1x1 1 logical b 1x2 8 cell c 3x3 36 cell d 2x2x2 32 cell （2）元胞数组的数据获得从元胞数组中读取数据，可保存为一个标准的数组或一个新的单元数组，或取出数组进行计算。元胞数组中数据的访问，可通过元胞内容的下标进行，用元胞数组名加大括号。大括号中数值表示元胞的下标。如a1，2表示元胞数组中第一行第二列的元胞。 a=20,matlab;ones(2,3),1:3a = 2

19、0 matlab 2x3 double 1x3 doublestr=a(1,2)str = matlabclass(str)ans =cellstr=a1,2str =matlabclass(str)ans =char（）和有着本质的区别，大括号用于表示元胞的内容，小括号表示指定的元胞。a = 20 matlab 2x3 double 1x3 doublea2,1(2,2)ans = 0.9134a2,1(2,3)ans = 0.0975a1,2(2)ans =a使用元胞的下标，可将一个元胞数组的子集赋值给另一个变量，创建新的元胞数组。 a=1,2,3;4,5,6;7,8,9a = 1 2 3

20、 4 5 6 7 8 9 b=a(2:3,2:3)b = 5 6 8 9 c=a(1:3,2:3)c = 2 3 5 6 8 9本例使用元胞下标的方式创建了新的元胞数组b和c，通过结果看出b和c就是元胞数组a的一部分。（3）元胞数组的删除和重塑要删除单元数组中的行或列，可以用冒号表示单元数组中的行或列，然后对其赋一个空矩阵即可。a=20,matlab;ones(2,3),1:3a = 20 matlab 2x3 double 1x3 doublea(1,:)=a = 2x3 double 1x3 double a=20,matlab;ones(2,3),1:3;a1=a = matlab 2x

21、3 double 1x3 doublea(1)=a = 2x3 double matlab 1x3 doublea(2)=a = 2x3 double 1x3 doublea(1,2)=? A null assignment can have only one non-colon index.a(1)=a = 1x3 double元宝数组和其他数组一样，也可以通过reshape函数改变形状，改变后的元胞数组与原元胞数组的元素个数相同，不能通过改变形状来添加或删除元胞数组中的元素。 a=cell(4,4)a = size(a)ans = 4 4 b=reshape(a,2,8)b = size(

22、b)ans = 2 8（5）元胞数组中的操作函数 cell：创建空的元胞数组cellfun：为元胞数组的每个元胞执行指定的函数celldisp：显示所有元胞的内容cellplot：利用图形方式显示元胞数组 cell2mat：将元胞数组转变成为普通的矩阵 mat2cell：将数值矩阵转变成为元胞数组 num2cell：将数值数组转变成为元胞数组 deal：将输入参数赋值给输出 cell2struct：将元胞数组转变成为结构 struct2cell：将结构转变为元胞数组iscell：判断输入是否为元胞数组 a=20,matlab,3-7i;ones(2,3),1:3,0a = 20 matlab

23、3.0000 - 7.0000i 2x3 double 1x3 double 0 b=cellfun(isreal,a)b = 1 1 0 1 1 1 c=cellfun(length,a)c = 1 6 1 3 3 1 d=cellfun(isclass,a,double)d = 1 0 1 1 1 1（函数的应用）cellfun函数的主要功能是对元胞数组的元素（元胞）分别指定不同的函数，不过，能够在cellfun函数中使用的函数ushuliang是有限的。能在cellfun中使用的函数：isempty：若元胞元素为空，则返回逻辑真islogical：若元胞元素为逻辑类型，则返回逻辑真isr

24、eal：若元胞元素为实数，则返回逻辑真 length：元胞元素的长度ndims：元胞元素的维数prodofsize：元胞元素包含的元素个数（7）元胞数组的嵌套元胞数组的元胞中包含其他的元胞数，称为嵌套元胞数组，没有嵌套结构的元胞则称为页元胞。使用嵌套的大括号或cell函数，或直接用赋值表达式，都可以创建嵌套单元数组，另外还可以访问嵌套元胞数组的子数组、元胞或元胞的元素。 a=cell(1,2)a = a(1,2)=cell(2,2)a = 2x2 cella(1,1)=magic(3);a1,2(1,1)=1 2 3;4 5 6;7 8 9;a1,2(2,1)=2-i;4+7i;a1,2(2,2)=cell(1,2);a1,22,2(2)=5;cellplot(a)（8）元胞数组与数值数组间的转化应用循环，可以将元胞数