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

1、1第二章第二章时域中的离散信号和系统时域中的离散信号和系统22.1 2.1 模拟信号与其采样序列模拟信号与其采样序列 对模拟信号 以采样周期T进行等间隔采样, 得到采样序列x(n) 它的每一个样本点与相应处的模拟信号值相等，而在其它时间位置处的值为零。 这一采样序列与原来的模拟信号是不等价的：不包围面积，量纲不同，物理作用不同。)(txannTxtxnxanTta),()()(3模拟信号与其采样序列的关系模拟信号与其采样序列的关系 要构成一个与模拟信号有等价性的离散序列，应该有 这个序列具有面积的量纲，T很密时，有因此和原始的模拟信号等价。 这个关系式不难从图2.1.1中得到验证)()()(n

2、TxTnxTnxae)(nxe0( )( )limeat nTTx nT x tT宽度 所围面积4模拟信号与其采样序列的关系模拟信号与其采样序列的关系 图 2.1.1 连 续 信 号 与 采 样 信 号 的 面 积 等 价 5模拟信号与其采样序列的关系模拟信号与其采样序列的关系 人们用（2.1.1）式的序列作为离散信号的研究模型，由于抽掉了模拟信号中的主要物理量时间。在数字域处理信号时可只按序号逐次进行运算和存储，这就是抽象为序列模型的好处。 另一方面必须看到模型与原物理过程之间的本质差异。在把序列与模拟信号进行转换和等价比较的时候，必需按等价序列xe(n)的概念（2.1.2）式来考虑。62.

3、2 2.2 基本序列及其运算基本序列及其运算 序列的表示方法完全地表示序列x(n)要用x和n两个向量，x表示序列的值，n表示它的位置或时间次序。 例如x(n)=2,1.2,-1.4,3,1,4,3.1,7（下划双线表示n=0处的采样点）用MATLAB表示时用两个向量： n = -3,-2,-1,0,1,2,3,4; x = 2,1.2,-1.4,3,1,4,3.1,7; 可见n是顺序增加的整数序列，可写成n=ns:nf,序列的长度为lx=nf-ns+1。本例n= -3:4;lx=8。7基本序列及其运算基本序列及其运算 有限长序列和无限长序列 序列位置向量的起点和终点都是有限数称为有限长序列 序

4、列位置向量的起点向- 延伸为无限长序列，简称左序列。 序列位置向量的终点向 延伸也为无限长序列，简称右序列。工程中遇到的序列都属于有限序列，有些情况下可看成右序列，但不会有左序列。复数序列可分解为实部和虚部两个实数序列8 单位脉冲序列(n) x,n = impseq(np,ns,nf) 单位阶跃序列(n) x,n = stepseq(np,ns,nf) 矩形序列 复数指数序列 正余弦序列 x(n)=cos(n)+jsin(n)常用的典型序列常用的典型序列( )( )()NRnnnN()jnxej ne9 function x,n = impseq(np,ns,nf)% ns为序列起始位置，nf

5、为序列终止位置, np为脉冲位置if nsnp | nsnf | npnf error（输入位置参数不满足ns = np = 0单位脉冲序列的函数程序单位脉冲序列的函数程序10基本序列的图形基本序列的图形 11序列的基本运算序列的基本运算 序列的基本运算(a) 序列相乘 (b) 序列相加 (c) 倍乘(d) 右移位(迟延)(e) 左移位(提前)(f) 复制其中(a)属于非线性运算，(e)属于非因果运算图 2.2.2 序 列 运 算 的 基 本 类 型 12 序列相乘相加子程序seqmult,seqadd编法： 将两个序列x1,x2的位置向量n1与n2对准，长度取两者之合；形成一个新的包含n1与

6、n2的公共位置向量n。其起点为两者中小的起点，终点为两者中大的终点， n = min(ns1,ns2):max(nf1,nf2); 将x1,x2补零，在n上延长，构成y1,y2; 各个对应元素相乘或相加； y=y1,*y2; (或y=y1+y2;)序列的序列的相乘和相加相乘和相加运算运算13序列相加的子程序序列相加的子程序seqaddseqadd n = min(min(n1),min(n2):max(max(n1),max(n2); % y(n)的长度 y1= zeros(1,length(n); y2 = y1; % y1,y2初始化 y1(find(n=min(n1)&(n=mi

7、n(n2)&(n=max(n2)=x2; % 把x2补零，扩展为具有y的长度的y2 y = y1+y2; % 序列相加.14例2.2.3：x1=0,1,2,3,4,3,2,1,0; n1=-2:6; x2=2,2,0,0,0,-2,-2; n2=2:8则 ny=-2:8;y1=0,1,2,3,4,3,2,1,0,0,0;y2= 0,0,0,0,2,2,0,0,0,-2,-2; yadd=y1+y2=0,1,2,3,6,5,2,1,0,-2,-2ymult=y1.*y2=0,0,0,0,6,9,0,0,0,0,0序列的序列的相加和相乘相加和相乘运算运算15例例2.2.32.2.3的计算图

8、解的计算图解 图 2.2.3 位置向量不同的两序列的相加和相乘 16序列向量的值不变，只把位置向量右移k拍子程序y,ny=seqshift(x,nx,k)的核心语句： y=x; ny=nx+k 左移时k取负值例：把x=3,7,-2,4; nx=-1:2;k=2代入则 y =3,7,-2,4; ny=1:4序列移位序列移位y(n)=x(n-k)y(n)=x(n-k)17意义为将序列x以n=0点为中心，左右翻转。 例如 x=3,7,-2,4; nx=-1:2则 y=4,-2,7,3; ny=-2:1即x的值及其位置排序都要左右翻转子程序y,ny=seqfold(x,nx)的核心语句：y=flipl

9、r(x);ny=-fliplr(nx)序列的序列的折叠：折叠：y(n)=x(-n)y(n)=x(-n)18序列能量：其定义为它相当于一个序列与自身共軛的相乘 。 序列功率的计算与能量相仿。序列能量的计算序列能量的计算nnxnxnxnxE)()()(219例2.2.2是位置向量相同的序列的运算(a)给出了序列移位、倍乘和相加运算的例(b)给出了序列移位、相乘和相减运算的例(c)给出了随机序列和相加运算的例(d)给出了序列能量运算的例计算结果图示于图2.2.2中。当序列的位置向量不同时，先要将序列放到同样的公共位置上，才能进行运算。例2.2.3给出了这样的算例。基本序列及其综合运算举例基本序列及其

10、综合运算举例20例例2.2.2的计算图解的计算图解 图 2.2.2 序列的简单运算与合成 212.3 2.3 几种重要的序列及其运算几种重要的序列及其运算周期序列：对所有n,满足的 称为周期序列。满足此式最小的N称为基本周期。给定x(n),生成周期序列 (xtide)的MATLAB方法：1.简单复制法：nxtide=0:K*N-1; xtide=x,x,x,2.求余函数(mod)法：nxtide=0:K*N-1; xtide=x(mod(nxtide,N)+1) kNnxnx( )x n( )x n22周期序列及其产生的算例周期序列及其产生的算例例2.3.1 设x=2,3,4,5，求将它延拓6

11、个周期所得的序列。解：x=2,3,4,5;N=length(x); K=6;ny=0:K*N1；y=x(mod(ny,N)1) 图2.3.1 序 列 的 周 期 延 拓 23实序列的对称性实序列的对称性 有限长实数序列，若满足xe(-n)= xe(n) ，称为偶序列，或对称序列。xo(-n)=-xo(n)，称为奇序列，或反对称序列。任何任意实数序列x(n)都可以分解成偶序列和奇序列两个分量: x(n)= xe(n)+ xo(n) (2.3.2) 其中 xe(n)= x(n)+x(-n)/2 (2.3.3a) xo(n)= x(n)-x(-n)/2 (2.3.3b) 24将任意实序列分解为奇偶序

12、列的子程序：function xe,xo,m=evenodd(x,nx)x1,nx1 = seqfold(x,nx);% x1(n)是x (n)的折叠序列,xe,m=seqadd(x/2,nx,x1/2,nx1) % 按(2.3.5)求出序列xexo=seqadd(x/2,nx,-x1/2,nx1) % 按(2.3.5)求出序列xo注意：分解后的位置向量m具有左右对称的特性。这是因为序列经过折叠后，其位置向量与原位置向量正负翻转，再经过相加时求两个位置向量的合，必然形成对称向量。（程序集中的子程序要兼容复数序列，故略有不同。）分解为奇偶序列的子程序分解为奇偶序列的子程序25序列的对称性分解算例

13、序列的对称性分解算例2.3.22.3.2例2.3.2 设x(n)=u(n+2)-u(n-7), 将它分解为偶和奇分量.解：n = -2:10; x = stepseq(-2,-2,10)-stepseq(7,-2,10);x1,n1 = seqfold(x,n);xe,nxe = seqadd(0.5*x,n,0.5*x1,n1);xo,nxo = seqadd(0.5*x,n,-0.5*x1,n1);26 序列的对称性分解算例序列的对称性分解算例2.3.22.3.2 图 2.3.2 实序列分解为奇偶序列 27序列的对称性分解算例序列的对称性分解算例2.3.32.3.3设n1=2:8,x1=s

14、in(n1);n2=-1:3;x2= x2=ones(1,5)(a) 求x1与x2之和。(b) 求x1与x2之积。(c) n3=0:4;x3=6,5,4,3,2; 求它对n=0点折叠后的序列y3。 (d) 利用(c)的结果，构成一个对称序列y4，(e) 利用(c)的结果，构成一个反对称序列y5解: (a) 的核心语句如下，(b)相仿n1=2:8,x1=sin(n1);n2=-1:3;x2= ones(1,5);y1,ny1=seqadd(x1,n1,x2,n2); 28序列的对称性分解算例序列的对称性分解算例2.3.32.3.3(d)和(e)的核心语句：n3=0:4;x3=6,5,4,3,2;

15、 % 将n扩展到负区域ny4=-4:4; %对称序列,把n=0处去一点y4=y3(1:4),x3;%反对称序列,把n=0处置0。 y5=-y3(1:4),0,x3(2:5)29序列的卷积运算序列的卷积运算序列x(n)和h(n)作卷积运算的定义此定义用图解表示时，先把x(k)和h(n-k)按k对齐。这时h(n-k)要解释为对k折叠并移位了n的h(k)。输出y(n)是把x(k)和h(n-k)覆盖区内k相同的样本乘积的求和而得。卷积输出的长度为两个输入长度之和减1： lyly=lx+lh-1=lx+lh-1kknhkxny)()()(30设x = 3,-3,7,0,-1,5,2; nx = -4:2

16、;h = 2,3,0,-5,2,1; nh = -1:4; 求其卷积y(n)。解：图2.3.4。x(k)和h(k)如左上图； x(k)和h(-k)如右上图；设n=-1,得出的x(k)和h(n-k)见左下图，设n=2,得出的x(k)和h(n-k)见右下图。由左下图k= -4: 0五个点上的x和h乘积之和，得到 y(-1)=3*2+ (-3) *(-5)+7*0+0*3+(-1)*2=19由右下图k=-2: 2五个点上的x和h乘积之和，得到 y(2)=7*1+0*2+(-1)*(-5)+5*0+2*3=18对每一个n值，都要这样运算y(n)，很繁。序列的卷积手工运算算例序列的卷积手工运算算例2.3

17、.42.3.431序列的卷积运算算例序列的卷积运算算例 图 2.3.4 卷积序列的折叠和移位 32对每一个m值，求y(m)可用下列语句h 1 ,nh 1 = s eqf o ld ( h , n h ) ; h1(nh1)=h(-nh)h2,nh2= seqshift(h1,nh1,m); h2(nh2)=h1(-nh1-m)y1,ny1= seqmult(x,nx,h2,nh2); y1=h2.*xy=sum(y1)用MATLAB函数conv可计算所有m值的卷积y=conv(x,h)不过它没有给出y的位置向量n.序列卷积用序列卷积用MATLABMATLAB运算算例运算算例33自编带有位置向量

18、的卷积子程序convwthn.mfunction y,ny = convwthn (x,nx,h,nh)nys = nx(1)+nh(1); nyf = nx(end) + nh(end); % MATLAB允许用end表示最后一个下标 y = conv(x,h); ny = nys:nyf;将例2.3.4中的 x,nx,h,nh代入，执行y,ny = convwthn (x,nx,h,nh)得到y= 6 3 5 6 19 -31 30 18 -27 -1 9 2ny=-5 -4 -3 -2 -1 0 1 2 3 4 5 6序列卷积用序列卷积用MATLABMATLAB运算算例运算算例34序列的

19、相关序列的相关 序列互相关的定义：它提供了两个序列之间相关程度的度量。它的定义与卷积相仿，只是不必将一个序列折叠。相关输出长度也是两个输入长度之和减一，lrxy=lx+ly-1自相关的定义：）（8 . 3 . 2)()()(nxyknynxkrnxxknxnxkr)()()(35序列的相关序列的相关MATLAB提供的计算序列相关性的函数xcorr.mrxy=xcorr(x,y) 它也不给出rxy的位置向量。当要知道rxy的位置向量nrxy时，可把相关计算变换为卷积计算，然后调用convwthn函数：根据定义：所以可先把x(n)折叠，求出x1(n)=x(-n)x1,nx1=seqfold(x,n

20、x)，再求相关序列rxy,nrxyrxy,nrxy = convwthn(y,ny,x1,nx1) )()()()()(nxnyknynxkrnxy36序列的相关例序列的相关例2.3.62.3.6例2.3.6 令x(n) = 3, 5, -7, 2, -1, -3, 2为一原型序列（可看作雷达发射信号）,令y(n)为x(n)加入噪声干扰并移位后的序列（可看作雷达回波信号） :y(n)=x(n-2)+w(n) 计算两者的相关性，可以判断被噪声淹没的回波的延迟，从而确定目标的距离。程序hc236执行的结果见图2.3.5。得知相关函数最大值在m=2处，即延迟量为2。37序列的相关例序列的相关例2.3

21、.62.3.6 图2.3.5 用 相 关 计 算 求 延 迟 382.4 2.4 时域表述的离散系统时域表述的离散系统 设x1(n)和x2(n)分别为系统的输入，y1(n)和y2(n)分别为其输出，T.为系统算子。y1(n)=Tx1(n), y2(n)=Tx2(n)则对线性系统一定满足：Tx1(n)+x2(n)=y1(n)+y2(n)；Tax1(n)=ay1(n)若系统 y(n)=Tx(n)，对于任何整数n0，都有：y(n-n0)=Tx(n-n0),则该系统称为时(移)不变离散系统。(LTILinear Time (Shift) - Invariant)39时域离散系统时域离散系统 由集总参数

22、部件构成的线性时（移）不变离散系统(LTILLinear Time-Invariant-Lumped) 可以用线性差分方程表示：式中 ,x(n) 和 y(n) 分别是系统的输入序列和输出序列,均为常数向量。对于物理可实现的系统，NM。 1, )()(000ainxbknyaMiiNkk40差分方程的解法差分方程的解法递推算法递推算法只把当前输出y(n)放在左端，将差分方程改写为：这个式子说明：除当前输入x(n)外，只要知道n时刻前的N个y的值 y(n-N),y(n-1) 及M个x的值 x(n-M),x(n-1) ,就可以求出当前的y(n)。由此可见，如果在起始点n=0处，知道了以前N个y和M个

23、x（统称初始条件)，即可依次递推计算所有的输出y。NkkMiiknyainxbny10)()()(41差分方程的解法经典算法把方程的解写成通解和特解两部分 y(n)= yH(n)+ yP(n) (2.4.8)通解为其中zk，k=1,2,.,N是特征代数方程的N个根。只有这些根都在单位圆以内，则由(2.4.6)描述的系统是稳定的。 特解部分由方程右端的输入函数确定。（从略）)9 . 4 . 2()(1NknkkHzcny00Nkkka z42递推算法举例例 2.4.4 设系统用差分方程 y(n)-y(n-1)+ 0.9y(n-2)=x(n)描述,输入序列为x=n,（1）设初始条件为 y(-2)=

24、0,y(-1)=0,求输出序列 y(n)。解：递推过程如下： y(0)= y(-1)- 0.9y(-2)+x(0)=1y(1)= y(0)- 0.9y(-1)+x(1)=1y(2)= y(1)- 0.9y(0)+x(2)=0.143用MATLAB递推算法差分方程的递推在MATLAB中用filter子程序来实现。调用的最简单形式为:y = filter(b,a,x) 其中 b = b0, b1, ., bM; a = a0, a1, ., aN 为差分方程(2.4.6)的系数数组,x 是输入序列。在上例中，语句为：b=1;a=1,-1,0.9; x=1,zeros(1,7);y1=filter(

25、b,a,x) 就可解出y。44 用时移算子用时移算子z -1表示离散系统表示离散系统定义z为时移算子，则左移运算y(n)=x(n+1)可以写成yn=zxn，简写为y(n) =x(n+1) =z x(n) 如果x=2,4,6,8，nx=-1,0,1,2, 则y与x相同，但ny=-2,-1,0,1。 n=0时的y要取 n=1时的x，这意味着要预知未来，物理上是无法实现的。定义z的逆运算子z -1，它的作用是延迟，即右移运算y1(n)= z -1x(n) =x(n-1) 在上例中就意味着y1=x,ny1=0,1,2,3,即y(0)=x(-1)。 n=0时的y要取 n=-1时的x，这是可以做到的。 4

26、5z -1是延迟一拍，因而z -k就是把信号延迟k拍。所以有x(n-k)=z-kx(n)， 于是差分方程(2.4.6)就可以表为从形式上推理，可以得到 y(n)=H(z)x(n) 这是离散系统的传递函数形式，但把z 1多项式放在分母上的运算没有定义，所以现在这个式子可以看，不能算。 用时移算子用时移算子z -1表示离散系统表示离散系统1, )()(000anxzbnyzaMiiiNkkkNNNNMNMNzazazaazbzbzbbzH1111011110)(462.5 离散系统的脉冲响应脉冲响应：系统的初始状态为零，输入为单位脉冲，其输出y(n) 定义为系统的脉冲响应h(n)。用公式表示为 h

27、(n)=T(n)设系统的输入为x(n)，则它的输出为 最后因此，系统的输出是其输入与脉冲响应的卷积。)()()()()(nhnxmnhmxnym )()()()()(mmmnTmxmnmxTny47用MATLAB计算离散脉冲响应(1)输入脉冲序列,用filter函数求输出 delta=1,zeros(1,L-1); h=filter(b,a,delta)(2) 用信号处理工具箱专门函数impzh,t=impz(b,a,L) 例如求例2.4.4中的脉冲响应，可以用下列语句：b=1,0,0;a=1,-1,0.9; h,t=impz(b,a,L)48用卷积法解析求系统输出例2.5.1，设系统脉冲响应

28、为h(n)=(0,9)nu(n)，输入x(n)=u(n)-u(n-10) ，求输出。解：h(n)和x(n)的波形见图2.5.1。将h(n)折叠，向左平移，它与x(n)的相互覆盖将出现三种情况：情况1， n0 ，h(-n)位于x(n)的左方，互不覆盖；y=0情况2，0n9 ,h(-n)逐步移入x(n)的区域，覆盖区逐渐加大，直到与x(n)全部覆盖；情况3，n9，10个x(n)脉冲始终全部与h(-n)覆盖。这三段的积分算式不同，因此形成的曲线也不同。49用卷积法解析求系统输出图上表示的是情况3，x(n)完全和h(-n)覆盖。此时卷积的计算为把相互覆盖的部分(子图3)相乘并求和，得到子图4中最右边的

29、样本。从子图4可以看出卷积的结果由三段线段构成：第一段为零；第二段上升；第三段下降直到零（未画全）。 图2.5.1 卷积的图示 50用向量矩阵乘法进行卷积计算 卷积计算可写成如下的MATLAB乘积形式y=x*H 其中(L=Nx+Nh-1)H的各列依次下移表示图2.5.1中h(k-n)的逐次右移。)(,),2(),1 (,)(,),2(),1 (NxxxxxLyyyy)(00)2()22(0) 1() 12() 11 (NxLhLhhLhhhH51H矩阵有Nx行和L=Nx+Nh-1列，各列依次向右下方移动。构成Toeplitz矩阵。它由第一列向量的转置C和第一行向量R完全确定：因此可以按照这个规

30、律来构成卷积矩阵。用下列语句：C=h(1),zeros(1,Nx-1);R=h,zeros(1,Nx-1)； H=toeplitz(C,R 用向量矩阵乘法进行卷积计算 . R(2) 1)( C(2) C(3). R(3) R(2) 1)( C(2). R(4) R(3) R(2) 1)(),(CCCRCToeplitzH52若：h = 2,3,0,-5,2,1; 则 H=toeplitz(h(1),zeros(1,Nx-1),h,zeros(Nx-1)得到用向量矩阵乘法进行卷积计算 2 3 0 -52 1 0 0 0 0 0 0 0 2 3 0 -52 1 0 0 0 0 0 0 0 2 3 0 -52 1 0 0 0 0 H 0 0 0 2 3 0 -5 2 1 0 0 0 0 0 0 0 2 3 0 -5 2 1 0 0 0 0 0 0 0 2 3 0 -5 2 1 0 0 0 0 0 0 0 2 3 0 -5