《数字信号处理教程-MATLAB释义与实现》第五章

第五章

变换域中的离散时间系统

15.1z变换傅氏变换有两个缺点，其一，在实际中许多有用的信号,如μ(n)和n*μ(n)...等,它们的离散傅利叶变换不存在。其二，系统对初始条件的暂态响应，或由时变输入引起的系统响应，都无法用离散傅利叶变换方法来计算。为了克服上述的两个缺点，要把离散傅利叶变换方法进行推广，推广后的方法称为z变换。2z变换的定义定义：对于一个序列x(n)，其双边z变换定义为:其中 是复变量，其向径为ρ，幅角为φ，构成一个复平面，称为z平面。序列x(n)的单边z变换则定义为：.3z变换的定义例如，对于序列它的双边z变换是：其单边z变换是:如果序数向量n只取正值，称为右序列，右序列的z变换X(z)中只出现z的负次幂；当采样序数n取负值时，X(z)中将出现z的正次幂。在工程中，人们感兴趣的主要是右序列。4z变换的收敛域z变换的收敛性

1.有限长序列：在除原点外的全z平面上收敛；

2。无限长右序列：在一个半径为r（称为收敛半径）的圆外的全z平面上收敛；

3。无限长左序列：在一个半径为r（也称为收敛半径）的圆内收敛；

4。双向无限序列：右序列和右序列收敛区的（环形）公共区，也可能没有；分别见下图中子图(a),(b),(c),(d)。5z变换的收敛域根据实际情况，只需考虑(a)和(b)两种收敛域。前者是对有限序列的，后者是对右序列的。它们的共同特性，那就是都在z=∞的邻域收敛，

6z反变换单边或双边z变换的反变换定义为。其中，积分路径是在复数平面中处于收敛域中的一条围线。在数字信号处理中，不需要用围线积分来求z反变换。5.2节中，将专门讨论求z反变换的其它方法。

7z反变换的非单值性右序列：和左序列：两个不同的序列具有相同的z变换因此必须规定z变换的收敛域，本书限定研究在∞邻域收敛的右序列。故其反变换唯一。8z变换的重要特性z变换的特性与DTFT和DFT的特性有很多相似之处，其证明都可以类比或从定义直接导出,所以不再重复推证。这里着重讨论几个重要特性的意义和应用。因为讨论都限于右序列，所以也免除了对收敛域的说明。1.线性特性： 设Z[g(n)]=G(z),Z[h(n)]=H(z),α,β为常数，则9z变换的重要特性2。样本的移位：

3。序列卷积：

4。初值定理

10z变换的重要特性5。终值定理：若x(n)是因果序列，且其z变换的极点均在单位圆内部，最多只有一个一阶极点在z=1上，则x(n)在n趋于无穷时的终值等于

在信号处理中有用的z变换特性，主要是上面几个，下面的几个特性可用于解析分析，对工程计算用处不大，供查考。

11z变换的重要特性6。乘以指数序列

7。序列乘以n(z变换微分)8。时域折叠

12z变换的重要特性9。复序列共軛10。序列乘积：

13用z变换计算卷积实例5.1.3

设x1(n)=[2,3,4]，x2(n)=[3,4,5,6]，求它们的z变换及两者的卷积输出y(n)。解：由z变换特性3，可得：X1(z)=2+3z-1+4z-2

和X2(z)=3+4z-1+5z-2+6z-3

为求x1(n)和x2(n)的卷积，先把X1(z)和X2(z)相乘Y(z)=(2+3z-1+4z-2)(3+4z-1+5z-2+6z-3）

=6+17z-1+34z-2+43z-3+38z-4+24z-5

它的反变换就是卷积输出序列：

y(n)＝6173443382414z变换的计算实例

其实多项式相乘和卷积计算相仿，在MATLAB中用的是同一个函数conv。在本例中要求z-3的系数，可以把第二组系数反过来排列,与第一组系数对齐如下。

234 6543 12+15+16+0=43

把对应项相乘（空项看作零）并逐项相加，得到上面的结果43。求z的其他幂次的系数时只需把第二组系数向左或向右移位即可，所以其计算和卷积过程相同。15z变换的计算实例5.1.4求长度为N的方波函数x(n)=RN(n)的z变换解：RN(n)是一个因果的有限长序列。当N趋向于无穷时， 。这就是阶跃函数的z变换。这类常用序列的z变换示于表5.1.1中

16z变换的计算实例5.1.5

用z变换性质和z变换表求下面序列的z变换。解：这是一个解析推导题。依次运用z变换的移位特性、微分特性和查表，得到175.2z反变换和差分方程的解

由定义（5.1.13）可知，求z反变换需要在一复围线上求积分：工程上计算z反变换并不需要求积分。解决这个问题有三种方法。（1）极点留数法；（2）部分分式法；（3）幂级数法（长除法）。MATLAB为这几种方法提供了简便易用的函数。

18z反变换和差分方程的解

用留数定理求z反变换：如果在围线c内的极点用zk表示，根据留数定理：等式右端求和号内是被积函数 在极点zk处的留数，所以z反变换是围线c内所有极点的留数之和。

19z反变换和差分方程的解

各个极点的留数求法如下如果是单极点，则

如果是N重极点，则手工计算留数也是相当麻烦的。MATLAB提供了计算多项式留数的函数。20z反变换和差分方程的解

例5.2.1设 求其z反变换。解：因为X(z)具有位于z=a处单极点，故解得：21z反变换和差分方程的解

例5.2.2设 求其z反变换

解：因为X(z)在z=0处有三重极点，故按这个式子，当n=1,n=2时，x都是0。只在n＝3处，有zn-3=00＝1，故有x(3)=3×2×1＝6。22z反变换和差分方程的解

可以看到，留数定理没有给出n=0处的x(0)。这是因为n=0时，z=0处并不是极点，所以用留数定理不能解决求x(0)的问题。正确的方法应该是用初值定理求x(0)，在本例中可得到：按照定义，用等式左右系数相比较的方法，一眼就可以看出x(0)=0,x(1)=0,x(2)=0,x(3)=6，可以校验留数定理的正确性。23用部分分式法求z反变换用部分分式法求z反变换过程总结如下：给定如果M≥N，可得式子右边第一项是真有理分式部分，第二项是直接多项式（无穷项）部分。

24用部分分式法求z反变换将真有理分式部分X1(z)进行部分分式展开，得此处pk是X（z）的第k个极点，Rk是该极点pk处的留数。假设它们都是单极点，则其单项的反变换由表查出为25用部分分式法求z反变换因此，(5.2.7)式的反变换x(n)可求出为：

其中留数Rk由下式给出：由此z反变换的计算就化成了一个代数问题，便于用MATLAB求解。26用部分分式法求z反变换MATLAB中的极点留数计算函数residuez，其基本调用格式为：

[r,p,C]=residuez(b,a)其中b和a为分子和分母的系数向量，p为分母的根向量，也就是X(z)的极点向量；r为对应于根向量中各个根的留数向量，C为直接项多项式系数向量，仅在M≥N时存在。Z反变换为27部分分式法求z反变换例5.2.3计算下式的反变换解:先用函数poly求出分母多项式的系数，再用函数residuez求X(z)的极点和留数。

b=1;a=poly([0.9,0.9,-0.7]); [r,p,C]=residuez(b,a)得到 r=[0.2461;0.5625;0.1914] p=[0.9000;0.9000;-0.7000] C=[]说明X(z)可分解为如下的部分分式：28部分分式法求z反变换例5.2.3其中第二项的形式是为了利用时域移位性(5.1.14),由此可得到反变换式：要计算出前8点的值，可以用MATLAB列出程序：n=0:7;x=（r(1)+r(2))*p(1).^n+r(2)*n.*p(2).^n+r(3)*p(3).^n29部分分式法求z反变换例5.2.3把这两条语句和前面的两条语句联在一起执行，得到：x=[1.001.101.661.752.052.112.252.25]如果把X(z)看作系统函数，则求X(z)的反变换就是求它的脉冲响应h(n)。函数名impz.m。其调用方式为h=impz(b,a,L)键入h=impz(1,poly([0.9,0.9,-0.7]),8)可以得到与上面的x相同的结果，

30部分分式法求z反变换例5.2.4例5.2.4计算下式的反变换要求求出的序列是因果的，并且不含复数。解：从表5.1.1中看到，当z变换的分母由实系数组成时，复数极点将成对共軛出现，结果本应不含虚数。为消除由于计算误差引入的微小虚部，程序末行用了real

函数b=[1,0.5];a=[1,-1.2,0.7];[r,p,k]=residuez(b,a) %求极点留数x=real(r(1)*p(1).^n+r(2)*p(2).^n)%取实部31幂级数法（长除法）求z反变换

由于右序列得到的z变换在z=∞的邻域解析。所以，可以把有理分式在z=∞的邻域展开为z–1的幂级数。幂级数展开可以用多项式相除的方法实现。例5.2.5用长除法求下列z变换多项式的反变换32幂级数法求z反变换

正确的长除法：系数按z-1的升幂排列，意味着在z=∞的邻域展开幂级数，用于右序列。33幂级数法求z反变换错误的方法：系数按z的升幂排列，意味着在z=0的邻域展开幂级数，只能用于左序列。34幂级数法求z反变换用MATLAB实现长除法：多项式除法是乘法的逆运算，它的函数名为deconv.m。其调用方法为：

[q,r]=deconv(b,a)

其中b为分子系数向量，a为分母系数向量，q为商的系数向量，r为余数的系数向量，多项式都按z-1的升幂排列。长除的目的是求q,商q的长度为M-N+1.如果长度不够,可以給分子向量b补零，使其长度M增加到Mbz。若希望商q的长度为Nq,补零的数目应为

Mbz-M=Nq+N-M-135幂级数法求z反变换长除求反变换的程序的核心语句为：b=?;a=?;Nq(输出向量的长度)=?M=length(b);N=length(a); x=deconv([b,zeros(1，Nq+N-M-1)],a)例5.2.6用长除法算例5.2.3,求出8点输出。解：

键入x=deconv([3,zeros(1,8+4-1-1)],[1,-1.1,-0.45,0.567])得到x=[1.001.101.661.752.052.112.252.25]36幂级数法求z反变换例5.2.7用长除法算例5.2.4,求出6点输出。解：写出MATLAB程序hc527如下b=[1,0.5];a=[1,-1.2,0.7];N=length(a);M=length(b);Nq=6;x=deconv([b,zeros(1,Nq+N-M-1)],a)运行结果仍为

x=[1.001.701.33990.4180-0.4364-0.8163]所以，求数值结果时，长除法比部分分式法方便得多。

37用z变换解差分方程

z变换的一个特点是可以求初始条件引起的响应，这时假定输入为零，所以也称为为零输入响应；把初始状态为零而加了输入造成的响应称为零状态响应。这两种响应的叠合就是方程的全解。为了推导初始条件引起的输出，采用单边z变换来分析n=0附近的问题。单边z变换时移特性如下38用z变换解差分方程上面的结果可用于求解具有非零初始条件和非零输入x(n)的如下差分方程其初始条件:y(i),i=-1,...,-N和

x(i),i=-1,...,-M。解的步骤：先对方程作单边z变换，把方程中的Z+[y(n-k)]和Z+[x(n-i)]都按前式展开；然后整理出Y+(z)的解；最后作z反变换求y(n)。

39用z变换解差分方程例5.2.8求解其中初始条件为：y(-1)=4和y(-2)=10。解：对差分方程的两边同时进行单边z变换，得到代入初始条件并整理，得40用z变换解差分方程进行部分分式展开，得到作z反变换后，得到解为这就是该系统在给定初始条件和输入信号x(n)下的全响应。即差分方程的全解。41用z变换解差分方程差分方程的全响应。可以用几种形式分解：(1)。通解和特解；从微分方程理论出发的分解。(2).暂态响应和稳态响应；从输出信号持续时间上进行分解。(3).零输入响应和零状态响应：从输出的诱发原因上分解。书上给出了本例的结果按三种情况进行分解的表达式。其总的响应是相同的。本书通常只求总的响应。42用z变换解差分方程函数filter也用于求全响应。当输入和初始条件给出时，调用格式如下

y=filter(b,a,x,xic) 其中，xic是等效的初始条件输入序列。它用filtic函数来计算。其调用格式为

xic=filtic(b,a,Y,X)b和a是滤波器的分子分母系数数组，Y和X是初始条件数组，分别由y(n)和x(n)的初始条件

Y=[y(-1),y(-2),...,y(-N)] X=[x(-1),x(-2),...,x(-M)]确定43用z变换解差分方程例5.2.8要计算例5.2.8的全响应，可用下列程序a=[1,-1.5,0.5];b=1;n=[0:7];x=(1/4).^n;Y=[4,10];xic=filtic(b,a,Y)%计算xicy1=filter(b,a,x,xic)y2=(1/3)*(1/4).^n+(1/2).^n+(2/3)*ones(1,8)％与解析式结果对比

y1和y2的结果是相同的44用z变换解差分方程例5.2.9

设 x(n)=cos(n

π/3)u(n),y(-1)=-2;y(-2)=-3;x(-1)=1;x(-2)=1,求解差分方程要求列出解的函数表达式，然后用数值验证。解：先解析求解，然后用MATLAB计算。对差分方程作单边z变换，代入初始条件，得到：将 代入并化简，得到Y+(z)。

45用z变换解差分方程例5.2.9化简和部分分式展开的程序hc529：b=[1,1,1]/3;a=[1,-0.95,0.9025];Y=[-2,-3];X=[1,1];xic=filtic(b,a,Y,X)bxplus=[1,-0.5];axplus=[1,-1,1];％X(z)系数ayplus=conv(a,axplus)％Yplus(z)的分母byplus=conv(b,bxplus)+conv(xic,axplus)％分子

[R,p,C]=residuez(byplus,ayplus)Mp=abs(p),Ap=angle(p)/pi％极坐标形式

46用z变换解差分方程例5.2.9求出R,p,C,Mp和Ap后，把Y+(z)展开为部分分式：并根据表5.1.1，写出其反变换：前两项对应稳态响应，后两项是暂态响应。

47用z变换解差分方程例5.2.9若只要数值结果，求y(n)前8个点的程序为。n=[0:7];x=cos(pi*n/3);y=filter(b,a,x,xic)%下面为对它的复数表示式进行验算A=real(2*R(1));B=imag(2*R(2));C=real(2*R(3));D=imag(2*R(4));y1=A*cos(pi*n/3)+B*sin(pi*n/3)+((0.95).^n).*(C*cos(pi*n/3)+D*sin(pi*n/3))执行这个程序，可以证明y和y1是相同的。复数计算程序太复杂，实际上不必采用。485.3z域中对系统的描述在z域中，离散系统的系统函数定义为它的脉冲响应h(n)的z变换因为h(n)通常都是右序列，所以求和的下限可以取为零。根据系统的输出序列y(n)等于其输入序列x(n)与系统脉冲响应序列h(n)的卷积，可以得出下列重要关系式：

Y(z)=H(z)X(z)49z域中对系统的描述证明如下：对卷积式两边取z变换，可得：改变求和次序，引入下标变量n=n-k，并利用h(n)的因果性，得到

50z域中对系统的描述另一种定义方法。根据LTI系统的差分方程：

对两边取z变换，利用z变换的特性，可以得到把系统函数H(z)定义为输出序列的z变换Y(z)与输入序列的z变换X(z)之比，即得

51z域中对系统的描述这是系统函数的负幂形式，其中把a0归一化为1。在信号处理中，通常都采用这种形式。z是定义在一个复数平面上的变量 ，而H(z)通常是z的多项式有理分式，所以它也是一个复数。可以用其幅度和相位表示。

如果变量z的取值限于z平面的单位圆圆周上，ρ=1,故

。则系统函数就成为频率响应：

52z域中对系统的描述所以，频率响应是z变换表示的系统函数的特殊情况。或者说傅立叶变换是z变换的特例。为了方便，常常统称为系统函数,只是用括号中的自变量取z或取ejω（也可以直接用ω）加以区分。这样，我们就有了四种表示离散系统特性的方法：差分方程；脉冲响应；频率响应；系统函数；它们之间的转换关系列成表5.2.2。

53几种描述方法的变换关系差分方程脉冲响应h(n)频率响应

H(jω)系统函数H(z)差分方程递推，

DTFT后，解Y/X

Z变换后，解Y/X脉冲响应h(n)FIR情况:b(n)=h(n)

H(jω)=DTFT[h(n)]Z变换频率响应

H(jω)分子、分母系数对应方程系数

IDTFT[H(jω)]系统函数H(z)H(z)用z-1表示，用时延代替z-1展开Z反变换54z域中对系统的描述根据上述算式，系统输出、系统函数和系统输入的z变换Y(z),H(z),X(z)之间满足简单的乘法关系，故可以用下列框图(a)表示。

在实际中，人们往往不用Y(z)和X(z)而直接采用了序列y(n)和x(n)自身，画成框图(b)，应该注意,当H(z)含有分母部分时,它只有形式的意义。

55z域中信号通过系统求输出例5.3.1设系统函数为输入信号x=[2,3,4,5]，用z变换计算输出y(n)。解：先求输入信号的z变换X(z)，所以最后对Y(z)求z反变换得到y(n)=Z-1[Y(z)]。56z域中信号通过系统例5.3.1程序hc531

x=[2,4,3,5]; %输入序列及初始序号nfx=length(x)-1; %计算序列终止序号b=-3;a=[2,-2.2,0.5];%分子分母系数,B=conv(-3,x);A=a;%输入与分子系数相乘[r,p,k]=residuez(B,A) %求r,p及k 得到r,p,k后，就可以写出输出函数的表达式57z域中信号通过系统例5.3.1如果只要求给出输出的数据和图形，则可以那本题就可以简单地用工具箱函数filter来解，程序为

x=[2,4,3,5,zeros(1,10)];b=-3;a=[2,-2.2,.5];y=filter(b,a,x);stem([0:ny],x)stem([0:ny],y)得到的波形如右图。58正幂系统函数及零极增益形式把分子分母同乘以(此处设N≥M，若N≤M，则同乘以），得到系统函数的正幂形式：

其分子分母系数向量可能变化， ，前或后会加几个零。把分子分母各分解因式后，得到

59正幂系统函数及零极增益形式zl为系统的零点；pk为系统的极点；b0则相当于增益。这种系统函数称为零极增益形式。LT1系统可在z域中用零极点图的形式来描述。这在设计简单的滤波器时很重要，只要正确地配置零极点就可达到一定的设计要求。分别对正幂有理函数H(z)的分子分母多项式使用MATLAB的roots函数，就可求得其零、极点。roots的逆向函数是poly，用它可由多项式的根求得其系数。用zplane(b,a)函数,可由给定的分子行向量b和分母行向量a绘制成系统的零极点图。

60z域中对系统的描述例5.3.2设系统的差分方程为：y(n)+0.3y(n-1)=2x(n)+4x(n-2)-2x(n-4)求此系统的系统函数，分别求出它的负幂形式、正幂形式和零极增益形式。解：对两边作z变换，得到所以其系统函数的负幂形式为它的系数向量为61z域中对系统的描述若化成正幂形式，分子分母同乘以z4，得到因此，它的系数向量为可以看到两种形式的系数向量是不同的。因为正幂形式的系数向量从z的最高幂一直排到零次幂（即常数项），低次项有一个空缺时要在系数向量中放一个0。而负幂形式的系数向量从z的零次幂（即常数项）一直排到最低幂（负的最高幂），最低幂以后就不管了，因此尾部不会出现0。

62z域中对系统的描述现对例5.3.2的系统函数的分子分母分别求根：z=roots([2,0,4,0,-2]),p=roots([1,0.3,0,0,0])得到z=[1.5538i;-1.5538i;0.6436;-0.6436] p=[0；0；0；-0.3]故其系统函数的零极增益形式为此系统函数分子分母中z的最高次均为4，有四个零点和四个极点，增益项为2。

63z域中对系统的描述为了得到系统零极点在z平面上的分布。键入a=[1,0.3,0,0,0];b=[2,0,4,0,-2];zplane(b,a) ％画零极点分布图图中。用‘x’表示系统的四个极点，其中原点为三重极点；用‘o’表示它的四个零点，两个在实轴上，两个在虚轴上形成共軛。。

64z域中对系统的描述本例系统函数的分子分母多项式的最高正幂次数相同，所以极点和零点的数目相等。若分母的最高正幂次数高于分子，零极点图上看到的极点数目会比零点多。实际上此系统的H(∞)=0，意味着在z=∞处还有零点。把这些无穷远处的零点也算上，系统的零极点数目总是相同的。对右序列，H(z)在z=∞处收敛，所以在该处只会有零点，不可能有极点。

注意zplane函数的调用方法规定：a和b必须取正幂系数，且必须用行向量输入。如果输入的是零极点向量z和p，那末它们必须取列向量形式。65z域中对系统的描述利用MATLAB的三维图形功能，可以画出系统函数H(z)的对数幅特性在z平面上的分布，如右图，从图上也可以大体看出零点和极点的位置。66z域中对系统的描述上述图形是用下列MATLAB语句得到的a=[1,0.3,0,0,0];b=[1,0,2,0,-1];[X,Y]=meshgrid([-2:0.05:2]+0.001,[-2:0.05:2]+0.001);H=polyval(b,X+j*Y)./polyval(a,X+j*Y);mesh(Y,X,20*log10(abs(H)))程序第二行语句是为了设置平面自变量矩阵X和Y,每个方向各81个点，自变量总数为81×81=6561个点,在meshgrid函数中加的0.001是为了避开真正H(z)=∞。第三行语句是为了求H(z),因为函数freqz只能求单位圆上的H(z),不能求全z平面的H(z)。所以只能从基本的多项式求值函数polyval出发来求。67零极点分布与系统的因果性零极点分布与系统的因果性：因果系统的充分必要条件：系统的脉冲响应必须是右序列。前面指出，右序列在变换域的特征是它的z变换在无穷远处收敛。因此因果系统的极点不可能在无穷远处，只能在z平面上一个有界区域内。如果系统函数用的是多项式分式形式，则因果性的要求是N≥M，即在正幂形式时分母上z的最高次数（在负幂形式时为z-1的最低次数）大于分子上z或z-1的对应次数。68零极点分布与系统的稳定性不同的极点位置所对应的脉冲响应如下图。

其规则是：69零极点分布与系统的稳定性在单位圆内的单极点和重极点：当n→∞时，脉冲响应趋向于零。在单位圆外的单极点和重极点：当n→∞时，脉冲响应趋向于无穷大。在单位圆上的单极点：当n→∞时，脉冲响应趋向于常数和等幅振荡。

在单位圆上的重极点：当n→∞时，脉冲响应趋向于无穷大。所以LTI系统稳定的充要条件是：系统的极点位于单位圆内部。705.4z平面上的谱分析零极点与频率特性的关系

如果系统是稳定的，则单位圆处于收敛域内，就可以算出单位圆圆周上的H(z)。。在单位圆周上， ，所以其频率特性可用零极增益方法写成：71零极点分布与频谱分析因子 可解释为成为z复平面中由零点 指向单位圆周上处的向量，因子 可看成为复平面中由极点指向单位圆周上 处的向量，如图示。因此幅频响应函数可看成所有零点到单位圆周上给定点ejω的矢量长度的积，除以所有极点到该点矢量长度的积，再乘以b0。也可类比地看出相频响应函数的几何意义。72零极点分布与频谱分析可以看出系统频率响应的几何意义。系统幅特性应该等于两虚线的乘积除以两个实线长度的乘积。

73零极点分布与频谱分析例5.4.1设系统有一对位于0.5±0.7j的共軛极点，另有位于0.3和1.5的两个实数零点，增益为2。则其系统函数为用zplane([0.3;1.5],[0.5+0.7j,0.5-0.7j]语句，可以画出它的零极点如图5.4.1所示。所以这个系统的脉冲响应是右序列，是因果的。而且它的极点都在单位圆内，所以是稳定的。

74零极点分布与频谱分析如果极点靠单位圆很近时，幅频率特性在靠近该极点附近会出现大的谐振峰，这是由于极点到它的距离很近，分母迅速减小的缘故。从零点的分布还可以导出最小相位系统的概念。若零点在单位圆内，当ω沿反时针方向由0→

π时，由零点产生的向量β的相角ψ=0→π；若零点在单位圆外左方，相角ψ=0→某θ→0；若零点在单位圆外右方，则ψ=π→某θ→π；所以如果系统所有的零点都在单位圆内，当ω沿反时针方向由0→

π时，它们引起的正相移最大，或负相移最小，称为最小相位系统。

75z平面上的谱分析非单位圆周上的频谱分析

在很多实际应用中，并非整个单位圆上的频谱都很有意义。如果信号的极点离单位圆较远，可以使采样点轨迹沿一条接近这些极点的弧线或圆周进行，则采样结果将会在极点对应的频率上出现明显的尖峰，这样就能较准确地测定出极点频率。因此往往有必要计算半径为r<1的圆上的频谱,或者计算螺旋线上的频谱。

76z平面上的谱分析对均匀分布在以原点为圆心的任何圆上的N点( ,k=0,1,…,N-1)频谱采样，其结果为上式说明，要计算x（n）在半径为r的圆上的N点等间隔频谱分量，可以先对x(n)乘以r-n，再计算N点DFT（FFT）即可。77z平面上的谱分析例5.4.2

设序列信号 取64点,求其频谱特性，及在z平面上r=0.9的中心圆上的频谱。解：此序列的z变换为 ，即在z=0.91处有一个单实极点。所以r=0.9的圆在该点附近H(z)会有很大的幅特性。程序hc542:n=0:63;p=0.91;x=p.^n;%给定序列参数H=freqz(x,1,64); %求序列的频率响应Hx1=x.*(0.9).^-n; %将序列乘r-n得到x1，

H1=freqz(x1,1,64);%求r圆上的频谱H178z平面上的谱分析程序运行的结果见右图。注意子图(d)上的峰值比子图(c)的峰值大了好多倍。79z平面上螺旋线的谱分析Chirp-z变换

要分析z平面上M点频谱采样值。设 式中A和W为复数，用极坐标形式表示为式中，A0和W0为实数。当k=0时，有由此可见，(5.4.5)式中，A决定谱分析起始点z0的位置，W0的值决定分析路径的螺旋趋势，φ0则表示两个相邻分析点之间的夹角。

80z平面上螺旋线的谱分析例5.4.3设序列信号为：取64点，求其频谱及在z平面上r=0.9的中心圆上的频谱，并求它在一根通过1+0i和0.4500+0.7794i两点的z平面螺旋线上的频谱。解：此序列的z反变换为 可以求出它有两个共軛极点0.4500±0.7794i。在调用czt函数时，主要的难点是确定参数W，使所取的螺旋线通过所关心的极点附近。

81z平面上的谱分析今设定点数M=64已知，螺旋线在上半平面，相角的增量＝π/M。今设第一点到关心点之间的相角为π/3，对应于Mx=(π/3)/(π/M)步，两点的模值比为R=1/0.9，W0的取法应为 ,程序的核心语句为：M=64,n=0:M-1;dw=pi/M;%

定点数和频率分辨率x=0.9.^n.*(cos(pi*n/3)+2*0.2887*sin(pi*n/3));A=1;W0=0.9^(-3/M),W=W0*exp(-j*dw);%参数A,Wz=A*W.^(-(0:M-1));%螺旋线轨迹方程

G=czt(x,M,W,A);%

调用czt函数

82z平面上的谱分析

图(a)为序列x(n)图(b)为它的频谱，图(c)为它在r<1的圆上的频谱，图(d)为所取的螺旋线轨迹， 图(e)和(f)为它在螺旋线轨迹上的频谱，两图的坐标比例尺不同。835.5理想滤波器理想的低通数字滤波器应该具有如下频率特性。其曲线见右图。从幅特性看，它可以让|ω|≤ωc的频谱分量100%地通过，而让|ω|>ωc的频谱分量100%地衰减。从相特性看，它使通过的信号产生了一个相位迟延-ω×τ,τ是一个正的无量纲常数，所以也常用拍数n0表示。

84理想滤波器线性相位的重要性迟延也可以看成是一种失真。无延迟的理想滤波器的脉冲响应h(n)是一个双向序列，因而是不可实现的非因果系统。这点时延通常是可以容忍的，工程中重要的往往是不要有波形失真。如果让所有的频谱分量都迟延同样的时间，那么整个波形将原封不动地保持它的形状，仅仅迟延了一定的时间而已。因为相位＝频率×时延如果时延相同，那么由于频率不同造成的相位迟延就与频率成线性关系。考虑到因果性，时延不可能是正的，所以相位也必须是负的。85理想滤波器

理想滤波器的脉冲响应，用傅立叶反变换求得令ωc=1，及n0＝6，h(n)可以用下列语句画出

wc=1;n0=6;n=-[10:30]+1e-5; h=sin(n-n0)./(n-n0)/pi;stem(n,h)%绘图得到的脉冲响应是一个双向序列。考虑到因果性，只能截取出n≥0的部分,而考虑到线性相位的要求，截取的脉冲序列必须对称。86理想滤波器图(a)中n0=6，截取的长度只能达到13。图(b)n0=15，截取出n≥0的部分达到31，滤波器的幅频率特性更接近理想滤波器，付出的代价则是时延增长到15拍。

87理想滤波器

高通滤波器、带通滤波器和带阻滤波器的理想频率特性见图。也可以求出它们的脉冲响应。它们频率特性和脉冲响应都可以由低通滤波器特性叠合而成。88理想滤波器线性相位条件对脉冲响应对称性的要求：先任意取两项来分析。令可以把这个式子看作复平面上两个向量的相加。这两个向量的模分别为h1和h2，而其相角分别为θ1=-ωn1和θ2=-ωn2,设合成向量的模为h3，相角为θ3。提出的问题是：如果θ1和θ2都加k倍,要使θ3也加k倍，则h1和h2应满足什么关系？这是一个简单的三角关系如下图。

89理想滤波器如果h1=h2，则三角形OAC是一个等腰三角形，其两个底角相等，由此可得：

θ3=θ1+θ1=θ1+(θ1-θ2)/2=(θ1+θ2)/2=-ω(n1+n2,)/2=-τω其中τ=(n1+n2)/2=n0

可见此时θ3与ω呈线性关系。它的迟延量为原来两向量迟延量的平均数。

90理想滤波器当N个向量合成，要保证所有向量的相角增长与合成向量的增长满足线性关系，那就不仅要它们的大小成对出现，而且要求各对合成的向量都在同一方向。即合成的向量都取同样的n0=(n1+n2)/2。结论是，合成的向量H(jω)的相角-ωn0中的n0应该就是序列的对称中心点，每一对向量都应该关于它对称配置。若序列位置从n=0到n=N-1，对称条件应为h(0)=h(N-1),h(1)=h(N-2),…，对称中心在n0=(N-1)/2处。当N是奇数时，n0为整数；当N是偶数时，n0含有小数0.5。91理想滤波器在反对称的情况下，h1=-h2,上述几何关系仍正确，所以相应的论证也成立，它相当于把h2的相角加了一个π，从而θ3=-ω(n1+n2,)/2–π/2=–π/2-τω所以反对称的序列构成的滤波器，其相位中除了线性项-τω外，还有一个等于–π/2的常数项。同样它也可以根据N为奇数或偶数，分成两类。其对称中心也在n0=(N-1)/2处，当N是奇数时，n0为整数；当N是偶数时，n0含有小数0.5。

92理想滤波器归纳起来，能够保证线性相位的滤波器，其脉冲响应只有四种类型：类型I。N奇数，脉冲响应对称；n0为整数。类型II。N偶数，脉冲响应对称；n0含0.5。类型III。N奇数，脉冲响应反对称;n0为整数。类型IV。N偶数，脉冲响应反对称；n0含0.5。这四种类型的滤波器用途各有不同，设计时也各有特点，本书将在第七章FIR滤波器设计中深入讨论。93理想滤波器

右图为三种截断条件下的滤波器幅频相频特性：上：n0=6,N=13中：n0=15,N=31下：n0=6,N=14可见：N愈大幅特性愈理想；不对称相特性愈非线性。945.6低阶数字滤波器一阶数字滤波器它的传递函数具有如下形式它在z=b处有一个零点，而在z=a处有一个极点。为了保证稳定性，极点应该位于单位圆之内，因而|a|<1。零点可以在单位圆内，也可以在单位圆外。本节讨论最小相位系统，则零点应在单位圆内，|b|<1

。

95一阶数字滤波器取单位圆上的一点 ，从极点和零点向该点引向量α和β。则有因而系统的频率特性为：先设极点在原点，a=0，于是p为常数1，若零点为b=-1时，从ω=0增加到π时，q将从2变到0，

H(ω)的模也将同样变化，它是一个滤波作用很弱的低通滤波器，

96一阶数字滤波器 若零点放在b=1处，它就成为一个高通滤波器。 总之离零点远的频率，其幅度就大，这是一条规则。 对于极点，就反过来，离极点近的频率，其幅度就大。97一阶数字滤波器设极点在a=0.8处，零点在原点b=0，这时q=1为常数，所以零点对幅特性没有影响。系统的频率特性将与p成反比。可以想象，ω＝0的频点离极点最近，所以该处的幅频特性应该最大。当从ω＝0增加到π时，p将从0.2变到1.8，H(ω)的模也将从5变到0.55。它的幅频特性见图5.6.1(c)中的实线。

98一阶数字滤波器又设零点固定在b=-1处，而让a在0与1之间变动，研究此滤波器的通带如何变化。此时滤波器系统函数为：它满足H(ej0)=0。在取a=0.7,0.8,0.9三种数值时，它的频率特性见图5.6.1(d)。可以计算出它的通带ωp近似公式为

ωp=1-a

这也就是极点离开单位圆的最小距离。99一阶数字滤波器的设计例5.6.1设有一连续信号

，要求设计一阶低通数字滤波器，滤