自适应滤波器

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第三章 自适应数字滤波器,第三章 自适应数字滤波器,3.1,引言,3.2,自适应横向滤波器,3.3,自适应格型滤波器,3.4,最小二乘自适应滤波,3.5,自适应滤波的应用,1,3.1,引,言,自适应数字滤波器和维纳滤波器一样，都是符合某种准则的最佳滤波器。维纳滤波器的参数是固定的，适用于平稳随机信号的最佳滤波，但要设计这种滤波器，必须要求输入信号是平稳的，且具有信号和噪声统计分布规律的先验知识。在实际中， 常常无法知道这些先验知识，且统计特性还会变化，因此实现最佳滤波是困难的。,2,自适应滤波器的特点是：滤波器的参数可以自动地按照某种准则调整到最佳滤波；实现时不需要任何关于信号和噪声的先验统计知识，尤其当输入统计特性变化时，自适应滤波器都能调整自身的参数来满足最佳滤波的需要。 常常将这种输入统计特性未知，调整自身的参数到最佳的过程称为“学习过程”。 将输入信号统计特性变化时，调整自身的参数到最佳的过程称为“跟踪过程”，因此自适应滤波器具有学习和跟踪的性能。 由于自适应滤波器有这些特点，自1967年威德诺(B.Widrow)等人提出自适应滤波器以来，在短短十几年中，自适应滤波器发展很快，已广泛地用于系统模型识别，通信信道的自适应均衡， 雷达与声纳的波束形成，减少或消除心电图中的周期干扰，噪声中信号的检测、跟踪、,增强和线性预测等。,3,3.2,自适应横向滤波器,自适应滤波器的原理框图如图3.2.1所示，图中,x,(,n,)称为输入信号，,y,(,n,)是输出信号，,d,(,n,),称为期望信号，或者称为参考信号、训练信号，,e,(,n,),是误差信号。,其中,e,(,n,)=,d,(,n,)-,y,(,n,),自适应滤波器,H,(,z,)的系数根据误差信号，通过一定的自适应算法，不断地进行改变， 使输出,y,(,n,)最接近期望信号,d,(,n,)。 这里暂时假定,d,(,n,)是可以利用的，实际中，,d,(,n,)要根据具体情况进行选取， 能够选到一个合适的信号作为期望信号，是设计自适应滤波器的一项有创意的工作。如果真正的,d,(,n,)可以获得，,我们将不需要做任何自适应滤波器。,4,图,3.2.1,自适应滤波器原理图,5,3.2.1 自适应线性组合器和自适应FIR滤波器,1. 自适应滤波器的矩阵表示式,图 3.2.2 表示的是一个有N个权系数的自适应线性组合器， 图中,N,个权系数,w,1,w,2,w,N,受误差信号,e,j,的自适应控制。对于固定的权系数，输出,y,j,是输入信号,x,1,j,x,2,j,x,Nj,的线性组合，因此称它为线性组合器。这里的,x,1,j,x,2,j,x,Nj,可以理解为是从N个不同的信号源到达的瞬时输入，是一个多输入系统， 也可以是同一个信号源的N个序贯样本，如图 3.2.3 所示。因此它是一个单输入系统， 实际上这种单输入系统就是一个FIR网络结构， 或者说是一个自适应横向滤波器。其输出,y,(,n,),用滤波器的单位脉冲相应表示成下式：,(3.2.1),6,图,3.2.2,自适应线性组合器,7,图,3.2.3,自适应,FIR,滤波器,8,这里,w,(,n,)称为滤波器单位脉冲响应，令：,i,=,m,+1,w,i,=,w,(,i,-1),x,i,=,x,(,n-i,+1),，,n,用,j,表示，上式可以写成,(3.2.2),这里,w,i,也称为滤波器加权系数。用上面公式表示其输出，适合于自适应线性组合器，也适合于,FIR,滤波器。将上式表示成矩阵形式：,(3.2.3),式中,误差信号表示为,(3.2.4),9,2. 利用均方误差最小准则求最佳权系数和最小均方误差,误差信号被用来作为权系数的控制信号。下面采用均方误差最小的准则，求最佳权系数。由,(3.2.4),式，均方误差为,(3.2.5),令,(3.2.6),(3.2.7),10,将,(3.2.6),、,(3.2.7),式代入,(3.2.5),式，,得到,(3.2.8),R,dx,称为,d,j,与,X,j,的互相关矩阵，是一个,N,维列矩阵；,R,xx,是输入信号的自相关矩阵，特点如下：,(1)是对称矩阵，即;,(2),是正定或半正定的，因为对于任意矢量,V,满足下式：,自相关矩阵的主对角线是输入信号的均方值，,交叉项是输入信号的自相关值。,11,(3.2.8)式表明，当输入信号和期望信号是平稳随机信号时， 均方误差信号,E,e,2,j,是权系数的二次函数，即将(3.2.8)式展开时，公式中的权系数均以它的一次幂或二次幂出现。如果只有一个权系数,w,1,则,E,e,2,j,是,w,1,的口向上的抛物线；如果有两个权系数,w,1,w,2,则,E,e,j,2,是它们的口向上的抛物面；对于两个权系数以上的情况，则属于超抛物面性质。,E,e,j,2,在自适应信号处理中是一个重要的函数，经常称它为性能函数。为选择权系数，使性能函数到达它的最小点， 一些有用的自适应方法都是基于梯度法的，我们用 表示,E,e,j,2,的梯度向量，它是用,E,e,j,2,对每个权系数求微分而形成的一个列向量，,用公式表示如下：,12,(3.2.9),按照,(3.2.4),式,梯度推导如下：,(3.2.10),还可以用,(3.2.8),式对,W,求导得到,(3.2.1,1),令上式等于,0,得到最佳权矢量,W,*,的表达式：,(3.2.1,2),13,对比第二章维纳滤波器的最佳解，结果是一样的。上式也称为维纳权矢量。当自适应滤波器的权系数满足上式时，均方误差将取最小值。将(3.2.12)式代入(3.2.8)式得到最小均方误差：,(3.2.13),或者将上式取转置，用下式表示：,(3.2.14),我们知道，在维纳滤波器中，当滤波器的单位脉冲响应取最佳值时， 其误差信号和输入信号是正交的；这里也有相同的结果， 当权矢量取最佳值时，梯度为0，按照,(3.2.10),式：,14,例 3.2.1,一个单输入的二维权矢量自适应滤波器如图 3.2.4所示，图中输入信号与期望信号分别为,这两个信号都是周期性确定性信号， 因为任何正弦函数积的期望值，都可由这个积在一个或多个周期上作时间平均来计算，,可以推导出下面公式,6,:,15,图 3.2.4 两个权的自适应滤波器,16,上式表明性能函数,E,e,j,2,对权函数是二次型的，用(3.2.11)式求梯度向量，得到,求最佳权矢量可以用(3.2.12)式，通过对,R,xx,求逆得到，也可以通过上式，令 ，而求出：,17,用,(3.2.13),式求最小均方误差：,上式说明只要,N,2,不管N取多少，通过对权系数的调整可使均方误差达到0，此时输出信号,y,j,完全等于期望信号,d,j, 例如,N,=2， 按照上面公式，可以求出输入、输出信号以及最佳权系数如下：,18,3.2.2 性能函数表示式及其几何意义,在自适应滤波器的分析研究中，性能函数是一个重要函数， 前面已推导出性能函数用(3.2,.8),式表示，重写如下：,下面我们推导它的其它表示方法以及几何意义。,均方误差是权系数的二次函数，当权系数取最佳值时， 均方误差取最小值，将(3.2.14)式代入(3.2.8)式，可以用最小均方误差表示性能函数，推导如下：,为了表示方便，令,=,E,e,2,j,则,19,将,(3.2.12),式代入上式，得到,(3.2.15),令,V=W-W,*,=,v,1,v,2,v,N,T,(3.2.16),V,称为偏差权向量，它表示权向量对最佳权向量的偏差。这样性能函数可以表示得更简单：,(3.2.17),20,因为,R,xx,是对称的，正定或半正定的，利用它的特征值和特征向量再进一步简化，假设,R,xx,是,N,N,维，它的,N,个特征值为：,1,2,N,，将,R,xx,进行分解，得到,R,xx,=,Q,T,Q,=,Q,T,R,xx,Q,(3.2.18),通过调节使,Q,归一化，即,(3.2.19),(3.2.20),21,式中，,Q,称为正交矩阵或特征矩阵，,q,i,称为特征向量，满足下式：,(3.2.21),(3.2.22),是由特征值组成的对角矩阵，,用下式表示：,(3.2.23),将,(3.2.18),式代入,(3.2.17),式，得到,令,(3.2.24),22,则,(3.2.25),上式将性能函数变成了平方和的形式。再观察(3.2.24)式， 该式将V坐标中的,R,xx,的特征向量变成了,V,坐标中的单位向量。 利用(3.2.24)式将特征向量,q,i,变成,q,i,，再利用,(3.2.20),、,(3.2.21),式，,可得,(3.2.26),23,也就是说，,q,i,为,V,坐标中的第,i,个单位向量，,q,i,亦是,矩阵对应于,i,的特征向量。下面用二维权矢量的情况说明它的几何意义。对于二维权矢量情况，有下面公式：,24,图,3.2.5,二维权矢量性能表面,25,图 3.2.6 等均方误差的椭圆曲线族,26,按照,(3.2.17),式，有,或,当,c,=,min,时，对应椭圆的中心，,V=W-W,*, 则相当于,W,坐标平移到V坐标的原点，即,V,坐标的原点对应W坐标的最佳点,W,*,。这里,v,1,v,2,不是椭圆的主轴。但经过对,R,xx,的分解：,且,V,=,Q,T,V,将性能函数的椭圆族(按照(3.2,.25),式,),变成,27,即,或者,(3.2.27),显然，上式是一个椭圆方程，,v,1,和,v,2,是椭圆族的主轴，如果,1,2,则,v,1,是长轴，,v,2,是短轴。因此(3.2.24)式起坐标旋转的作用，将,v,1,v,2,旋转到主轴上，形成,v,1,v,2,主轴。对于维数,N,2的情况，长轴对应最小特征值，按照上面的椭圆方程长轴正比于；短轴对应于最大特征值，正比于 。另外，,因为,28,得到,(3.2.28),V,中单位矢量就是,V,坐标中的,R,xx,的特征矢量。,29,3.2.3,最陡下降法,1. 最陡下降法的递推公式,将,(3.2.11),式代入,(3.2.29),式，得到,(3.2.30),(3.2.31),在上式两边都减去,W,*,并令,V,j,=,W,j,-,W,*,得到,V,j,+1,=,I,-2,R,xx,V,j,(3.2.32),上式是一个递推公式，由于项不是对角矩阵，计算与分析均复杂。下面仍然采用坐标旋转的方法进行推导。,30,(3.2.33),此时，项已变成对角矩阵，假设起始值是,V,0,，可得到上式的递推解为,(3.2.34),31,再将(3.2.24)式代入，再经过坐标平移，即代入,V,j,=,W,j,-,W,*,式，,最后得到权系数的递推公式：,(3.2.35),上面递推公式中，部分已变成对角矩阵，,这使分析与研究自适应特性变得简单了。,32,2. 收敛条件,由最陡下降法的递推公式不难分析出它的收敛条件，即当迭代次数,j,趋于时，权系数收敛最佳时的条件。按照上式，,显然只有当,(3.2.36),(3.2.37),满足时，才能得到： 。(3.2.37)式即是最陡下降法的收敛条件，式中,max,是,R,xx,的最大特征值。,(3.2.36),式中的,0,表示,0,矢量。,33,3. 过渡过程,过渡过程是指权矢量和性能函数由起始点随迭代次数的增加，进行变化的过程。下面从权矢量和性能函数两方面讨论自适应滤波器的过渡过程。权矢量的过渡过程讨论如下：,按照,(3.2.34),式，权矢量的递推解是,第,i,个权系数递推方程是,(3.2.38),令,(3.2.39),34,将上式代入,(3.2.38),式，得到,(3.2.40),上式说明第,i,个分量,v,i,按指数规律变化，其时常数为,i,=1, 2, 3,N,(3.2.41),因为一般,取得比较小，可以近似为,i,=1, 2, 3,N,(3.2.42),35,因为,所以,再将,(3.2.40),式代入，得到,(3.2.43),36,(3.2.44),式中,(3.2.45),上式说明第,i,个加权系数按照,N,个指数和的规律变化，由初始值收敛到最佳值，其时常数与特征值成反比。下面分析性能函数的过渡过程。按照,(3.2.25),式，性能函数如下式：,(3.2.46),将,(3.2.40),式代入，得到,(3.2.47),37,上式说明性能函数也是按,N,个指数和的规律变化，和加权系数过渡过程不同的是时间常数不同，,它的时常数为,(3.2.48),我们已经知道，性能函数和各个加权系数都是按照N个具有不同时常数的指数和的规律变化的，时常数和特征值成反比，不同的特征值对应的收敛时间是不一样的，但最终的收敛要取决于最慢的指数过程，它的时常数最大，对应最小的特征值，公式如下：,(3.2.49),(3.2.50),38,但为保证收敛，,不能取得太大，受限于最大特征值,max,。 这样，如果特征值比较分散时，即,max,和,min,相差很大时， 使最陡下降法的收敛性能很差。下面分析,值的影响。,值收敛过程影响很大，首先必须选择得足够小，使之满足收敛条件：,但按照(3.2.47)、 (3.2.48)式，它影响收敛速度。 一般希望在保证收敛的条件下，选大一些，使时间常数小一些，收敛的速度快一些。但当,选择得太大时，即使收敛条件满足，也可能形成振动性的过渡特性。 在图 3.2.7 中，图(a)是,较小时的情况；图,(,b,),是,较大时的情况，此时过渡过程已发生振荡。,39,图 3.2.7,值的影响,(,a,) 较小时的情况； （,b,） 较大时的情况,40,3.2.4,最小均方,(LMS),算法,1. LMS算法的权值计算,LMS(Least Mean Square),算法的梯度估计值用一条样本曲线进行计算，公式如下：,(3.2.51),因为,41,所以,(3.2.52),(3.2.53),FIR,滤波器中的第,i,个权系数的计算公式为,(3.2.54),FIR滤波器中的第,i,个权系数的控制电路如图3.2.8所示, LMS自适应滤波器的总框图如图,3.2.9,所示。,42,图,3.2.8 FIR,第,i,个支路的控制电路,43,LMS算法的加权系数按照(3.2.53)式进行控制， 式中加权矢量的改变量是2,e,j,X,j,，梯度的估计值是-2,e,j,X,j,。显然，这是一个随机变量，这说明LMS算法的加权矢量是随机变化的。因此，LMS算法又称为随机梯度法。下面对这种算法的性能进行分析，,主要分析加权矢理和性能函数的平均变化规律以及它们的随机性造成的影响。,按照,(3.2.52),式，,对梯度估计值求统计平均，,得到,(3.2.55),上式说明梯度估计值是无偏估计的，梯度的估计量在理想梯度,j,附近随机变化，权系数也是在理想情况下的权轨迹附近随机变化的。,44,图,3.2.8 LMS自适应滤波器总计算框图,45,2.LMS算法加权矢量的过渡过程,将误差公式,(3.2.4),式代入,(3.2.53),式，得到,(3.2.56),按照,(3.2.53),式，,对加权矢量取统计平均：,(3.2.57),46,类似于最陡下降法的推导，经过坐标平移和旋转，变换到,V,坐标中。其公式推导如下,:,令,V,j,=,W,j,-,W,*,(3.2.58),那么,E,V,j,=,E,W,j,-,W,*,E,V,j,+1,=,E,W,j,+1,-,W,*,(3.2.59),将上面两式代入,(3.2.57),式中，得到,它的递推解是,令,R,xx,=,QQ,T,=,QR,xx,Q,T,(3.2.60),47,得到,(3.2.61),(3.2.62),再将,(3.2.59),、,(3.2.60),和,(3.2.61),式代入上式，得到,E,W,j,=,W,*,+,Q,I,-2,j,Q,-1,(,W,0,-,W,*,),(3.2.63),对比(3.2.35)式，说明LMS算法加权矢量的统计平均值的过渡过程和最陡下降法加权矢量的过渡过程是一样的。换句话说， LMS算法加权矢量是在最陡下降法加权矢量附近随机变化的，,其统计平均值等于最陡下降法加权矢量，那么，其收敛条件同样为,(3.2.64),48,在满足收敛条件的情况下，才有下式：,由于最大的特征值,max,不可能大于,R,的迹(,R,的主对角线元素之和,),，即,因此收敛条件可以表示为,(3.2.65),49,对于横向滤波器，,式中的迹是,NE,x,2,j,，即,N,倍的输入功率，,那么,(3.2.66),实际中，通常,选得很小，选,(3.2.67),同样由,(3.2.62),式，第,i,个分量为,(3.2.68),50,同样引入时常数,i,(3.2.69),(3.2.70),(3.2.71),同样，第,i,个权系数可以表示成,(3.2.72),51,3 . LMS算法性能函数的过渡过程学习过程,由于LMS算法加权矢量的平均值的变化规律与最陡下降法的加权矢量一样，可以推想它的均方误差也会按照最陡下降的均方误差变化规律变化。下面进行推导。,按照,(3.2.4),式，,信号误差为,(3.2.73),52,式中，,e,opt,j,=,d,j,-,X,j,T,W,*,，称为最佳误差信号，它对应于最小均方误差，,即,按照,(3.2.73),式写出均方误差表示式：,假定,X,j,和,V,j,不相关，上式中最后一项为,0,，那么,53,同样，假设加权系数变化很小，,V,j,也变化很小，,E,V,j,V,j,，这样：,类似前面的推导，得到,(3.2.74),(3.2.75),对照最陡下降法性能曲线(3.2.47)式，LMS均方误差变化规律和最陡下降法完全一样，学习曲线同样近似为几个不同时间常数的指数和。,54,4. 稳态误差和失调系数,由上面分析知道，权矢量的平均值可以收敛到它的最佳值， 但权矢量变化过程是随机的，即使其平均值收敛到最佳值，它仍然按照下式：,W,j,+1,=,W,j,+2,e,j,X,j,随机地进行变化，这样使权矢量仍在最佳值附近随机变化， 但均方误差将大于最小均方误差，如图,3.2.10,所示。为此，引入失调系数,M,，,M,定义为,(3.2.76),55,图,3.2.10 LMS算法稳态误差,56,可以推出,5,失调系数为,(3.2.77),或者,M,=,NP,in,(3.2.78),式中，,N,是滤波器的阶数，,P,in,是输入信号功率。上式说明,和输入功率加大都会增加失调系数。 在保证收敛的情况下,加大，会提高收敛速度，也说明为了减小失调系数，,应该适当选择收敛速度，以保证收敛速度和失调系数都满足要求。,57,图 3.2.11 是一个LMS自适应滤波器的计算机结果,5,，阶数,N,=5，其输入是信号加白噪声， 输入信号功率为1，中心频率是0.03,f,s,(,f,s,为采样频率)，噪声功率为0.5, 输入信号自相关函数的特征值为： 5.14、0.853、0.502、0.500、0.500 ，权系数初始值取0，,=0.0065。图中画出了一条样本学习曲线和150条样本学习曲线的平均曲线。该图表明个别学习曲线起伏较大，平均学习曲线起伏很小，计算出的维纳最小均方误差为0.743 96，用LMS算法得到的稳态误差大于该值，按(3.2.77)式计算的失调系数是,4.87%,，按计算机模拟结果测得的失调系数是,5.40%,。,58,图,3.2.11 LMS,算法的学习曲线,59,3.3,自适应格型滤波器,3.3.1 前、后向线性预测误差滤波器,1. 前向线性预测误差滤波器,为了分析简单，假设信号属于实平稳随机信号。前向线性预测误差滤波器直接由信号的线性一步预测导出。在维纳滤波器一章我们已研究了信号的线性一步预测问题，即由,x,(,n,-1),x,(,n-,2), ,x,(,n-p,)预测,x,(,n,)，其估计值,x,(,n,)和预测误差,e,p,(,n,)用下式表示：,(3.3.1),60,由于假设了信号是实的，式中预测误差,e,p,(,n,)和系数,a,p,k,均是实数。(3.3.1)式表明 是由,n,时刻以前的,p,个数据,x,(,n,-1)、,x,(,n,-2),x,(,n-p,)得到的估计，因此称 为前向预测误差。将前向预测误差用 表示，上式重写为,(3.3.2),对上式进行,Z,变换，得到,(3.3.3),令,(3.3.4),61,H,f,(,z,)称为前向预测误差滤波器的系统函数。前向预测误差滤波器的结构图如图 3.3.1所示。,图,3.3.1,前向预测误差滤波器,62,用均方误差最小的准则求前向预测误差滤波器的最佳系数,a,p,k,k,=1, 2,，,p,(3.3.5),将,(3.3.2),式代入上式，得到,k,=1, 2, 3,，,p,(3.3.6),上式表明前向预测误差与用于预测的数据正交，这就是对于前向预测误差的正交原理。按照第二章的推导，前向预测误差滤波器的最佳系数,a,p,k,和信号的自相关函数之间的关系式称为,Yule-Walker,方程式，重写如下：,63,(3.3.7),将上式用矩阵方程表示为,(3.3.8),64,2. 后向线性预测误差滤波器,如果利用,x,(,n,+1),x,(,n,+2), ,x,(,n+p,)数据预测,x,(,n,)，则称为后向预测，其估计值用 表示。这样,(3.3.9),一般前向、后向预测用同一数据进行，即利用,x,(,n,),x,(,n,-1)，,x,(,n,-2)，,x,(,n-p,)进行预测，为此，将上式改为,(3.3.10),这样，前向预测是由,x,(,n-,p),x,(,n-p,+1),x,(,n,-2),x,(,n,-1)预测,x,(,n,)，后向预测是由,x,(,n-p,+1),x,(,n-p,+2),x,(,n,)预测,x,(,n-p,)， 这两种预测数据之间的关系如图 3.3.2 所示。,65,图,3.3.2,前向预测数据之间的关系,66,设后向预测误差用 表示,(,实际表示的是信号在,n-p,时刻的预测误差,),，,这样,(3.3.11),同样，利用最小均方误差的准则，可以得到关于后向预测时的正交原理以及Yule,-Walker,方程，它们分别用下面的,(3.3.12),和,(3.3.13),式表示：,k,=1, 2, 3,p,k,=1, 2, 3,p,(3.3.12),(3.3.13),67,式中， 是后向预测误差的最小误差功率。将(3.3.13)式和(3.3.7)式进行对比，它们极其相似。利用,Toeplitz,矩阵的性质，,可得到以下重要关系：,(3.3.14),(3.3.15),上面两式表明前、后向预测的最小误差功率相等，系数也相等(如果是复数，则是共轭关系,),。由,(3.3.10),、,(3.3.11),、,(3.3.14),式得到,(3.3.16),68,式中，当,k,=0, 1, 2, 3,p,时，,p-k,=,p,p,-1,p,-2, 0,因此也可以写成下式：,由上式画出后向预测误差滤波器的结构图如图,3.3.3,所示。,69,图,3.3.3,后向预测误差滤波器,70,对比图 3.3.1 和图 3.3.3, 或者对比公式(3.3.2)和(3.3.17)， 它们的系数虽然一样，但后向预测误差滤波器的系数排序却是前向预测误差滤波器系数排序的逆转排列。,对,(3.3.16),式进行,Z,变换，得到,(3.3.18),后向预测误差滤波器的系统函数为,(3.3.19),71,将上式与前向预测误差滤波器的系统函数(3.3.4)式对比，得到前、 后向预测误差滤波器的系数函数之间的关系是,为了求解前、后向预测误差滤波器的最佳系数，需要解Yule-Walker方程。可以采用高斯消元法解出,a,p,k,(,k,=1,2, 3,p)以及,2,p,但需要,p,3,量级运算量。利用Yule-Walker方程中的自相关矩阵是一个埃尔米特(Hermitain)和托布列斯(Toeplitz)矩阵的特点，且至少是半正定的，可以有效地减少运算量，这就是下面要推导的,Levinson-Durbin,算法，它的运算量级是,p,2,。,72,3.Levinson-durbin算法, Levinson-Durbin算法首先由一阶AR模型开始，按照(3.3.8)式，一阶AR模型,(,p,=1),的,Yule-Walker,为,由该方程解出：,73,然后增加一阶，即令,p,=2,，按照,(3.3.8),式得到,由上面方程解出：,74,然后令,p,=3, 4,以此类推，,可以得到一般递推公式如下：,(3.3.21),(3.3.22),(3.3.23),(3.3.24),(3.3.25),75,上面(3.3.21)(3.3.25)式就是Levinson-Durbin递推公式， 该式中的,k,p,称为反射系数。在(3.3.24)式中，,2,p,和,2,p,-1,是预测误差的均方值，因此,1-,k,2,p,必须大于等于,0,，这样,k,p,应要求满足下式：,(3.3.26),进而得到 ，即预测误差随递推次数增加而减少。把,k,p,称作反射系数，是类似于传输线的情况，如图3.3.4 所示，第p节的输出功率(即下一级的输入功率,),等于前一级的输出功率减去本级的反射功率，用公式表示如下：,(3.3.27),76,图3.3.4 传输线,77,3.3.2 格型滤波器,1. 由预测误差滤波器导出格型滤波器,将前面已推导的前向预测误差公式,(3.3.2),重写如下,:,再将系数,a,p,k,(,k,=1,2,3,p,)的递推公式(3.3.23)代入上式，并令,k,p,=,a,p,p,，得到,78,将上式与,(3.3.2),式对比，方程式的右边前两项是,p,-1,阶前向预测误差，,即,(3.3.28)方程式的右边最后一项中，因为,k,=1,2,3,p,-1时，,p-k=p,-1,p,-2, ,1，方括号部分可以写成,将上式右边与,(3.3.16),式对比，该部分就是,n,-1,时刻,p,-1,阶的后向预测误差，,即,79,这样由,(3.3.28),式，得到前向预测误差的递推公式，,即,(3.3.29),类似地，得到后向预测误差的递推公式为,(3.3.30),利用(3.3.29)式和(3.3.30)式，组成格型滤波器的第p节的结构图, 如图 3.3.5(a)所示。,80,图,3.3.5,全零点格型滤波器,81,对于,p,=0,的情况，,按照,(3.3.2),式和,(3.3.11),式，,得到,整个预测误差格型滤波器的结构如图 3.3.5(b)所示。由于没有反馈支路，它是一个全零点格型滤波器。 经过变形还可得到其他类型，如全极点格型滤波器、全极点横向滤波器，等等,5,。,82,2. 格型滤波器的性质,(1),各阶后向预测误差相互正交。,用公式表示如下：,设,i,j,按照(3.3.12)式， 与,x,(,n,-,j,+1),x,(,n-j,+2),x,(,n,-,i,),x,(,n-i,+1),x,(,n,)数据正交，但按照(3.3.16)式， 是,x,(,n,-,i,),x,(,n-i,+1),x,(,n,),的线性组合，因此 与 相互正交。,各阶后向预测误差相互正交的结果，使滤波器前后级互相解耦，对于系统最小化问题化为一系列独立的对每一级局部最小化问题。用作自适应滤波时，各级可选用不同的自适应步长， 使收敛速度提高。另外，为提高线性预测性能，需要增加一节或几节，可以只对,新增加的级进行独立的调节，达到输出均方误差最小，无需再调节前面的系数。,83,(2) 平稳随机序列可由自相关函数或反射系数表征。按照Levinson-Durbin递推公式，已知,r,xx,(0),k,1,k,2,k,p,从一阶开始，可以推出全部的预测系数,a,p,1,a,p,2, ,a,p,p,和,2,p,，把得到的这些数据代入Yule-walker方程，可求得信号的自相关函数,r,xx,(0),r,xx,(1),r,xx,(2),r,xx,(,p,)。以上说明平稳随机序列可由自相关函数表征，也可由,r,xx,(0),k,1,k,2,k,p,表征。,(3),前向预测误差滤波器是最小相位滤波器，即它的全部零点在单位圆内。,84,3.,对于复信号的预测误差滤波器和格型滤波器,前向预测误差滤波器的系统函数,H,e,(,z,)以及前向预测误差公式 和实信号情况一样，仍是,(3.3.4),式和,(3.3.2),式，但利用均方误差最小原则求预测系数要用下式求解：,(3.3.32),对于前向预测误差的正交原理，则用下式表示：,(3.3.33),前向预测误差滤波器的预测系数和信号自相关函数之间的Yule-Walker方程仍和,(3.3.8),式一样。,85,后向预测误差和后向预测误差滤波器系统函数分别用下式表示：,(3.3.34),(3.3.35),对于后向预测误差的正交原理为,(3.3.36),对于复信号的,Levinson-Durbin,递推公式为,(3.3.37),86,k,=1, 2, 3,p,-1,(3.3.38),(3.3.39),(3.3.40),(3.3.41),复信号的全零点格型滤波器预测误差递推公式为,(3.3.42),(3.3.43),87,图,3.3.6,复信号预测误差全零点格型滤波器,88,3.3.3,最小均方误差自适应格型滤波器,P,阶格型滤波器由,p,节组成，如果前,m,节的参数,k,i,(,i,=1, 2, 3,m,)为最佳，相应的预测误差功率是最小，而后面的节的参数对前面的最佳参数无影响，因此在m节的基础上再加一节， 则只需根据使第m+1节的预测误差功率最小的原则选择,k,m,+1,即可。 预测误差功率有前向预测误差功率和后向预测误差功率，这里采用使前、后向预测误差功率的和为最小的原则求反射系数。,公式为,(3.3,.44),89,将(3.3.29)、(3.3.30)代入上式， 可以得到,(3.3.45),实际计算时，上式中的统计平均值用时间平均计算，,公式为,(3.3.46),对于复信号情况，公式为,(3.3.47),90,上面两式便是直接利用数据计算反射系数的递推公式。下面讨论公式中的求和限问题，如果输入数据为,x,(,i,),i,=0, 1, 2,n,，,当,p,=1,时，,这里,因此,91,上式中，求和限必须限制在已知的输入数据范围内计算，这样求和限应为,i,=1,2,3, ,n,，计算公式为,当,p,=2,时，,92,按照,(3.3.29),、,(3.3.30),式，得到,将上面两式带入 公式中，可以计算出 ，考虑到输入数据的范围，具体计算公式为,93,再根据 ,按照(3.3.29)、(3.3.30)式计算,e,2,(,i,)、,b,2,(,i,)，按照(,3.3.46),式计算,以此类推。,这样，对于 具体计算公式为,(3.3.48),以上便是直接采用信号数据计算格型滤波器的反射系数以及最小预测误差的方法。但这种算法必须从低阶推起，要求较大的存储时，有较大的计算延迟，使应用受到限制。下面介绍梯度算法，,这种算法可以减少运算量，,且适合非平稳情况。,94,自适应格型滤波器的梯度算法中反射系数的计算，类似于自适应横向滤波器中系数的递推算法，公式为,(3.3.49),式中，,仍然是控制收敛速度和收敛的参数；,kp,表示对方括弧中的部分求梯度，,95,将上式代入,(3.3.49),式中，,得到,式中，,=2,，为步长因子。这部分内容可参考文献,5,。,(3.3.50),96,3.4,最小二乘自适应滤波,本节讨论另外一种以误差的平方和最小作为最佳准则的误差准则最小二乘(Least Square)准则。,定义,(3.4.1),式中,(,n),是误差信号的平方和；,e,j,是,j,时刻的误差信号，,97,d,j,是,j,时刻的期望信号，,X,j,是,j,时刻的输入信号构成的向量， W表示滤波器的权系数构成的向量。通过选择,W,，使,(,n,)取得最小值的滤波称为最小二乘(Least Square,简称LS)滤波，而满足,E,e,2,j,取得最小值的滤波称为最小均方误差(Least Mean Square, 简称LMS)滤波。和LMS滤波相比，LS滤波对非平稳信号的适应性要强许多，这是由于LS滤波总是采用新的准则，在每一个时刻对所有已输入信号而言，重新评估使其误差的平方和最小，因此具有更精确的含义，属于精确分析法。 而LMS滤波是以集合平均为基础的，属于统计分析的方法。,98,3.4.1,最小二乘滤波,1. 最小二乘的基本问题,已知,n,个数据,x,(1),x,(2),x,(,n,)，采用,M,个权的FIR滤波器对数据进行滤波，假设期望信号为,d,(,i,),，如图,3.4.1,所示。,滤波器的输出 是对期望信号,d,(,i,),的估计,其中，,w,k,(,i,),k,=1,2, ,M,，为FIR滤波器在,i,时刻的,M,个系数值(说明滤波器的系数可以变化)，它是一个,M,维的向量，记为,w,M,(,i,)=,w,1,(,i,),w,2,(,i,),w,M,(,i,),T,。同理，输入信号也是一个,M,维的向量，,x,M,(,i,)=,x,1,(,i,),x,2,(i),x,M,(,i,),T,。,n,时刻，估计误差为,99,(3.4.2),误差信号的平方加权和为,上式中，由于数据长度有限，对观测区间以外的数据所做的约定不同，当,i,的取值范围不同时，得到不同的,(,n,),。这里采用前加窗法，约定：,x,(,i,)=0,i,0,得到,(3.4.3),100,图 3.4.1,M,个权的FIR滤波器,101,为了后面叙述方面，引入一些符号。令,M,维向量,w,M,(,n,)和,x,M,(,n,),分别表示,n,时刻的滤波器的权向量和输入信号向量,当,i,=1,2, ,n,时，引入,n,维误差向量,e,(,n,)和期望信号向量,d,(,n,),，以及输入信号构成的,M,n,维矩阵,X,M,(,n,),e,(,n,)=,e,(1),e,(2),e,(,n,),T,d,(,n,)=,d,(1),d,2),d,(,n,),T,X,M,(,n,)=,x,M,(1),x,M,(2),x,M,(,n,),(3.4.5,a,),(3.4.5,b,),(3.4.5,a,),102,为了后面推导方便，引入,n,M,维矩阵,C,，定义,(3.4.5,b,),应用这些符号，期望信号的估计和估计误差可以表示为,(3.4.6),(3.4.7),将,(3.4.7),式代入,(3.4.3),式，得到误差信号能量,103,式中,是加权矩阵，对角线上的元素称为加权因子。为了推导简单起见， 在后面的分析中，取,=I,，则,(3.4.9),104,要使,(,n),取得最小值，满足,(3.4.10),成立的,w,(,n,)就是,w,M,(,n,)的最小二乘估计，记为 。,应用标量求导公式，,计算得到,(3.4.11),将,(3.4.5b),式和,(3.4.7),式代入上式，有,(3.4.12),将,(3.4.11),式展开,(3.4.13),105,引入,M,维向量,p,M,(,n,),以及,M,M,维矩阵,R,M,(,n,),，,(3.4.14),(3.4.15),则,(3.4.13),式可以写为,(3.4.16),可以看出，,R,M,(,n,)类似于输入信号的自相关特性，,p,M,(,n,)类似于输入信号与期望信号的互相关特性。(3.4.16)式与第二章中的维纳-霍夫(Wiener-Hopf)方程相似，不同之处在于维纳,-,霍夫方程中的数学期望符号用求和符号所代替。,106,若矩阵,X,M,(,n,)的秩等于,M,，记做rank,X,M,(,n,)=,M,则,X,M,(,n,),X,M,T,(,n,)非奇异，求解(3.4.16)式，可以得到,w,M,(,n,)的最小二乘估计,w,LS,(,n,),(3.4.17),若rank,X,M,(,n,),M, 则,w,M,(,n,)不能唯一辨识。,在,(3.4.12),式两边同左乘以 ，得,(3.4.18),应用,(3.4.6),式，得到,(3.4.19),107,图,3.4.2,最小二乘估计的几何解释,108,当 存在时，最小二乘的估计值 为,最小二乘估计的误差信号能量,min,为,(3.4.20),综合前面的分析，我们可以把最小二乘问题用模型,z=A,+n,(3.4.21),109,来描述，其中，,z,是观测信号，,n,为噪声信号，,A,可看作一数据矩阵，表征输入与输出之间的关系，,是可调整量。与前面的分析相对应，参见(3.4.7)式，,n,类似于误差信号,e,(,n,)，,z,类似于信号真值,d,(,n,)，,A,类似于信号的估计值 ，且A矩阵与C矩阵相对应， 在数据矩阵A已经确定的情况下，,对,z,的最小二乘估计为,与图3.4.2的信号向量相对应，观测信号z构成的观测值向量与,d,(,n,)相对应，,A,与估计向量 相对应，噪声信号,n,与图中的估计误差,e,(,n,),相对应。,110,令误差信号能量为,J,，并取加权矩阵,=I,，则,(3.4.23),111,2. 最小二乘估计的质量,假设误差向量,n,(,k,)是独立同分布的，具有零均值,方差为,2, 那么，最小二乘估计是无偏估计; 若噪声是高斯噪声， 则最小二乘估计是一致估计。,(1),无偏性：,(3.4.24),(2),一致性：如果噪声是高斯噪声，那么最小二乘参数估计一致收敛，,即,(3.4.25),这里，,W.P.1,表示以概率为,1,的收敛。,112,证明,首先求解最小二乘估计的协方差，令 表示最小二乘估计的误差，,即,那么，当,具有零均值时，最小二乘估计的协方差为,(3.4.26),113,已知噪声为高斯噪声，则,E,nn,T,=,2,I,，得到最小二乘估计的协方差为,式中， 将以概率1收敛于一个正定阵，且,2,是有界的，,因此,最小二乘估计的无偏性已经保证了偏移量为0， 因此， 一致性得证。,前面，所讨论的信号模型都是MA模型，即最小二乘滤波器用FIR滤波器实现，若信号是一ARMA模型，那么怎样来求出信号所对应的最小二乘估计呢? 下面举一例说明。,114,例,3.4.1,已知模型,z,(,k,)+,az,(,k,-1)=,bx,(,k,-1)+,n,(,k,),其中，,n,(,k,)是均值为0，且与,x,(,k,)不相关的噪声信号，设定n(k)是各态遍历的高斯噪声，求,LS,的最小二乘估计。,这是一个近似的ARMA(1, 1)模型，其输出信号为,z,(,k,)， 输入信号是,x,(,k,)。这个模型的近似性表现在不包含当前的输入信号，仅包含前一时刻的输入，并且模型还包括噪声信号,n,(,k,)。 设定n个观测值构成的观测向量表示为,z,n,=,z,(0),z,(1),z,(,n,),T,。,115,和一般的最小二乘估计的模型,z,=,A,+,n,相比，这个问题的关键是如何确定A矩阵。MA模型的输出完全由输入信号所决定， 而ARMA模型的输出还包括以前时刻的输出对当前输出的反馈， 因此，ARMA模型的A矩阵应同时包括输出和输入数据。针对这个例题所选的模型，可以确定,A,矩阵为,其中，,z,(,i,),和,x,(,i,),分别表示,i,时刻的观测值和输入值。,116,解,根据,(3.4.17),式，计算,117,其中，,R,z,表示,z,的自相关函数，假设信号具有遍历性，则 ，同理，得到,R,x,(0),R,xz,(0)，,R,n,(1)等。,考虑,n,(,k,)与,x,(,k,)是统计独立的及,E,n,(,k,)=0， 记=,R,z,(0),R,x,(0)-,R,xz,2,(0),，由已知模型关系,其中，,a,0,和,b,0,为参数真值，得到,118,由此可以看出，只有当,n,(,k,),是白噪声时，,R,n,(1)=0,，才能使得,简单地说，所谓的最小二乘估计，是指在每个时刻对所有已输入的信号而言，误差平方和最小。最小二乘法的原理简单，,编程也易实现，,因此获得了广泛的应用。,119,3.4.2 递推最小二乘法(RLS),递推最小二乘法(Recursive of Least Square，简称RLS)的基本思想就是新的估计值是在老的估计值的基础上修正而成的，为了分析简单，便设,为一向量，且,仅与当前观测值有关，则,新的估计值,(,k,)=,老的估计值,(,k,-1)+,修正项,因此，最小二乘递推算法的关键是得到修正项的表达式。根据最小二乘估计式(3.4.22),用,a,(,i,),表示第,i,步迭代时,A,的取值，,A,k,表示前,k,步,A,的数值构成的向量。,120,定义一个变量,P,：,(3.4.29),其中,那么,(3.4.30),(3.4.31),121,令观测信号向量为,(3.4.32a),(3.4.32b),这里，,z,(,j,)表示,j,时刻观测信号的大小。,把,(3.4.29),式代入,(3.4.22),式，,k,-1,时刻的估计值,(3.4.33),上式两边同时左乘,P,-1,(,k,-1),，得,(3.4.34),122,把,(3.4.33),式和,(3.4.34),式代入,(3.4.22),式，,并应用,(3.4.30),式,(3.4.35),123,为了进一步得到加权矩阵的表达式，可定义,q,(,k,),(3.4.36),则,(3.4.35),式可以写为,(3.4.37),根据矩阵反演公式,(3.4.38),124,这里，,A, C,为一个具体数值，代入,(3.4.31),式，得,(3. 4.39),至此，得到了最小二乘递推算法的一种形式：,(3. 4.40),125,也可以将,(3.4.39),式代入,(3.4.36),式，,整理后得,(3.4.41),因此可以得到递推最小二乘的另一种形式：,(3.4.42),126,3.4.3,线性向量空间,1. 向量空间的基本概念,1) 希尔伯特空间,我们知道，一组数据可以看作是空间中的一个向量，一组向量可以确定一个空间，称为张成一个空间， 如二维空间可以由两个向量,x,=1, 0,T,和,y,=0,1,T,所张成，三维空间可以由三个向量,x,=1, 0, 0,T,y,=0, 1, 0,T,和,z,=0, 0, 1,T,所张成，那么,M,维向量空间，可理解为由,M,个基本向量,u,M,(,n,)所张成(这里的基本向量是指这些向量相互独立且正交)。,127,M,个基本向量,u,M,(,n,)构成一个数据矩阵，记为,U,=,u,1,(,n,),u,2,(,n,), ,u,M,(,n,)，那么矩阵U所张成的空间，记为,U,。因此， 一组数据对应于一个向量，一个(数据)矩阵对应于一个空间。 所谓的线性向量空间,U, 是由若干个向量,u,M,(,n,)线性组合得到的，空间的维数是指用来张成该空间的最少的向量数。,通常，我们所研究的向量空间都是希尔伯特空间，它与内积空间有紧密的联系。我们首先给出内积空间的定义。,128,(1) 内积空间。若,u,(,n,),v,(,n,),w,(,n,)是空间,U,中的任意向量，,且,a,为任一数，任意两个向量的内积为一个数，并满足下面的性质：,则称,U,为内积空间。,129,(2) 希尔伯特空间。若某空间是一个线性内积空间的完整部分，则称该空间为希尔伯特空间，简称希氏空间。所谓的完整，是指不存在这样的向量，可以任意接近该空间但不属于该空间。,自适应滤波这一章的内容所涉及的空间都是希尔伯特空间， 数学中常用的欧几里德空间也是希尔伯特空间,130,2) 子空间,(1) 子空间。所谓的子空间，是指用少于空间维数的向量所张成的空间。如在几何中，,X-Y,平面是立体,X-Y-Z,的一个子空间； 假设由向量,u,1,(,n,)=,u,1,(1),u,1,(2), ,u,1,(,n,),T,和,u,2,(,n,)=,u,2,(1),u,2,(2), ,u,2,(,n,),T,分别张成空间,u,1,和,u,2,，由向量,u,(,n,)=,u,1,T,u,2,T,T,=,u,1,(1),u,1,(2),u,1,(,n,),u,2,(1),u,2,(2),u,2,(,n,),T,张成空间记为,U, 那么,u,1,、,u,2,