自适应滤波器

1、第三章 自适应数字滤波器,3.1 引 言,自适应数字滤波器和维纳滤波器一样，都是符合某种准则的最佳滤波器。维纳滤波器的参数是固定的，适用于平稳随机信号的最佳滤波，但要设计这种滤波器，必须要求输入信号是平稳的，且具有信号和噪声统计分布规律的先验知识。在实际中， 常常无法知道这些先验知识，且统计特性还会变化，因此实现最佳滤波是困难的,自适应滤波器的特点是：滤波器的参数可以自动地按照某种准则调整到最佳滤波；实现时不需要任何关于信号和噪声的先验统计知识，尤其当输入统计特性变化时，自适应滤波器都能调整自身的参数来满足最佳滤波的需要。 常常将这种输入统计特性未知，调整自身的参数到最佳的过程称为“学习过程”

2、。 将输入信号统计特性变化时，调整自身的参数到最佳的过程称为“跟踪过程”，因此自适应滤波器具有学习和跟踪的性能。 由于自适应滤波器有这些特点，自1967年威德诺(B.Widrow)等人提出自适应滤波器以来，在短短十几年中，自适应滤波器发展很快，已广泛地用于系统模型识别，通信信道的自适应均衡， 雷达与声纳的波束形成，减少或消除心电图中的周期干扰，噪声中信号的检测、跟踪、 增强和线性预测等,3.2 自适应横向滤波器,自适应滤波器的原理框图如图3.2.1所示，图中x(n)称为输入信号，y(n)是输出信号，d(n)称为期望信号，或者称为参考信号、训练信号，e(n)是误差信号。 其中,e(n)=d(n)

3、-y(n,自适应滤波器H(z)的系数根据误差信号，通过一定的自适应算法，不断地进行改变， 使输出y(n)最接近期望信号d(n)。 这里暂时假定d(n)是可以利用的，实际中，d(n)要根据具体情况进行选取， 能够选到一个合适的信号作为期望信号，是设计自适应滤波器的一项有创意的工作。如果真正的d(n)可以获得， 我们将不需要做任何自适应滤波器,图 3.2.1 自适应滤波器原理图,3.2.1 自适应线性组合器和自适应FIR滤波器 1. 自适应滤波器的矩阵表示式 图 3.2.2 表示的是一个有N个权系数的自适应线性组合器， 图中N个权系数w1,w2,wN受误差信号ej的自适应控制。对于固定的权系数，输

4、出yj是输入信号x1j,x2j,xNj的线性组合，因此称它为线性组合器。这里的x1j,x2j,xNj可以理解为是从N个不同的信号源到达的瞬时输入，是一个多输入系统， 也可以是同一个信号源的N个序贯样本，如图 3.2.3 所示。因此它是一个单输入系统， 实际上这种单输入系统就是一个FIR网络结构， 或者说是一个自适应横向滤波器。其输出y(n)用滤波器的单位脉冲相应表示成下式,3.2.1,图 3.2.2 自适应线性组合器,图 3.2.3 自适应FIR滤波器,这里w(n)称为滤波器单位脉冲响应，令：i=m+1,wi=w(i-1), xi=x(n-i+1)，n用j表示，上式可以写成,3.2.2,这里w

5、i也称为滤波器加权系数。用上面公式表示其输出，适合于自适应线性组合器，也适合于FIR滤波器。将上式表示成矩阵形式,3.2.3,式中,误差信号表示为,3.2.4,2. 利用均方误差最小准则求最佳权系数和最小均方误差 误差信号被用来作为权系数的控制信号。下面采用均方误差最小的准则，求最佳权系数。由(3.2.4)式，均方误差为,3.2.5,令,3.2.6,3.2.7,将(3.2.6)、 (3.2.7)式代入(3.2.5)式， 得到,3.2.8,Rdx称为dj与Xj的互相关矩阵，是一个N维列矩阵；Rxx是输入信号的自相关矩阵，特点如下： (1)是对称矩阵，即; (2) 是正定或半正定的，因为对于任意矢

6、量V满足下式,自相关矩阵的主对角线是输入信号的均方值， 交叉项是输入信号的自相关值,3.2.8)式表明，当输入信号和期望信号是平稳随机信号时， 均方误差信号Ee2j是权系数的二次函数，即将(3.2.8)式展开时，公式中的权系数均以它的一次幂或二次幂出现。如果只有一个权系数w1,则Ee2j是w1的口向上的抛物线；如果有两个权系数w1w2,则Eej2是它们的口向上的抛物面；对于两个权系数以上的情况，则属于超抛物面性质。 Eej2在自适应信号处理中是一个重要的函数，经常称它为性能函数。为选择权系数，使性能函数到达它的最小点， 一些有用的自适应方法都是基于梯度法的，我们用 表示Eej2的梯度向量，它是

7、用Eej2对每个权系数求微分而形成的一个列向量， 用公式表示如下,3.2.9,按照(3.2.4)式, 梯度推导如下,3.2.10,还可以用(3.2.8)式对W求导得到,3.2.11,令上式等于0, 得到最佳权矢量W*的表达式,3.2.12,对比第二章维纳滤波器的最佳解，结果是一样的。上式也称为维纳权矢量。当自适应滤波器的权系数满足上式时，均方误差将取最小值。将(3.2.12)式代入(3.2.8)式得到最小均方误差,3.2.13,或者将上式取转置，用下式表示,3.2.14,我们知道，在维纳滤波器中，当滤波器的单位脉冲响应取最佳值时， 其误差信号和输入信号是正交的；这里也有相同的结果， 当权矢量取

8、最佳值时，梯度为0，按照(3.2.10)式,例 3.2.1 一个单输入的二维权矢量自适应滤波器如图 3.2.4所示，图中输入信号与期望信号分别为,这两个信号都是周期性确定性信号， 因为任何正弦函数积的期望值，都可由这个积在一个或多个周期上作时间平均来计算， 可以推导出下面公式6,图 3.2.4 两个权的自适应滤波器,上式表明性能函数Eej2对权函数是二次型的，用(3.2.11)式求梯度向量，得到,求最佳权矢量可以用(3.2.12)式，通过对Rxx求逆得到，也可以通过上式，令 ，而求出,用(3.2.13)式求最小均方误差,上式说明只要N2,不管N取多少，通过对权系数的调整可使均方误差达到0，此时

9、输出信号yj完全等于期望信号dj, 例如N=2， 按照上面公式，可以求出输入、输出信号以及最佳权系数如下,3.2.2 性能函数表示式及其几何意义 在自适应滤波器的分析研究中，性能函数是一个重要函数， 前面已推导出性能函数用(3.2.8)式表示，重写如下,下面我们推导它的其它表示方法以及几何意义。 均方误差是权系数的二次函数，当权系数取最佳值时， 均方误差取最小值，将(3.2.14)式代入(3.2.8)式，可以用最小均方误差表示性能函数，推导如下： 为了表示方便，令 =Ee2j, 则,将(3.2.12)式代入上式，得到,3.2.15,令,V=W-W*=v1, v2, , vNT,3.2.16,V

10、称为偏差权向量，它表示权向量对最佳权向量的偏差。这样性能函数可以表示得更简单,3.2.17,因为Rxx是对称的，正定或半正定的，利用它的特征值和特征向量再进一步简化，假设Rxx是NN维，它的N个特征值为： 1,2,N，将Rxx进行分解，得到,Rxx=QTQ,=QTRxxQ,3.2.18,通过调节使Q归一化，即,3.2.19,3.2.20,式中，Q称为正交矩阵或特征矩阵，qi称为特征向量，满足下式,3.2.21,3.2.22,是由特征值组成的对角矩阵， 用下式表示,3.2.23,将(3.2.18)式代入(3.2.17)式，得到,令,3.2.24,则,3.2.25,上式将性能函数变成了平方和的形式

11、。再观察(3.2.24)式， 该式将V坐标中的Rxx的特征向量变成了V坐标中的单位向量。 利用(3.2.24)式将特征向量qi变成qi，再利用(3.2.20)、 (3.2.21)式， 可得,3.2.26,也就是说，qi为V坐标中的第i个单位向量，qi亦是矩阵对应于i的特征向量。下面用二维权矢量的情况说明它的几何意义。对于二维权矢量情况，有下面公式,图 3.2.5 二维权矢量性能表面,图 3.2.6 等均方误差的椭圆曲线族,按照(3.2.17)式，有,或,当c=min时，对应椭圆的中心，V=W-W*, 则相当于W坐标平移到V坐标的原点，即V坐标的原点对应W坐标的最佳点W *。这里, v1v2不是

12、椭圆的主轴。但经过对Rxx的分解,且V=QTV将性能函数的椭圆族(按照(3.2.25)式)变成,即,或者,3.2.27,显然，上式是一个椭圆方程，v1和v2是椭圆族的主轴，如果12,则v1是长轴，v2是短轴。因此(3.2.24)式起坐标旋转的作用，将v1v2旋转到主轴上，形成v1v2主轴。对于维数N2的情况，长轴对应最小特征值，按照上面的椭圆方程长轴正比于；短轴对应于最大特征值，正比于 。另外， 因为,得到,3.2.28,V中单位矢量就是V坐标中的Rxx的特征矢量,3.2.3 最陡下降法,1. 最陡下降法的递推公式 将(3.2.11)式代入(3.2.29)式，得到,3.2.30,3.2.31,

13、在上式两边都减去W *,并令Vj=W j-W*, 得到,Vj+1=I-2RxxVj,3.2.32,上式是一个递推公式，由于项不是对角矩阵，计算与分析均复杂。下面仍然采用坐标旋转的方法进行推导,3.2.33,此时，项已变成对角矩阵，假设起始值是V0，可得到上式的递推解为,3.2.34,再将(3.2.24)式代入，再经过坐标平移，即代入Vj=Wj-W*式， 最后得到权系数的递推公式,3.2.35,上面递推公式中，部分已变成对角矩阵， 这使分析与研究自适应特性变得简单了,2. 收敛条件 由最陡下降法的递推公式不难分析出它的收敛条件，即当迭代次数j趋于时，权系数收敛最佳时的条件。按照上式， 显然只有当

14、,3.2.36,3.2.37,满足时，才能得到： 。(3.2.37)式即是最陡下降法的收敛条件，式中max是Rxx的最大特征值。(3.2.36)式中的0表示0矢量,3. 过渡过程 过渡过程是指权矢量和性能函数由起始点随迭代次数的增加，进行变化的过程。下面从权矢量和性能函数两方面讨论自适应滤波器的过渡过程。权矢量的过渡过程讨论如下： 按照(3.2.34)式，权矢量的递推解是,第i个权系数递推方程是,3.2.38,令,3.2.39,将上式代入(3.2.38)式，得到,3.2.40,上式说明第i个分量v i按指数规律变化，其时常数为,i=1, 2, 3, , N,3.2.41,因为一般取得比较小，可

15、以近似为,i=1, 2, 3, , N,3.2.42,因为,所以,再将(3.2.40)式代入，得到,3.2.43,3.2.44,式中,3.2.45,上式说明第i个加权系数按照N个指数和的规律变化，由初始值收敛到最佳值，其时常数与特征值成反比。下面分析性能函数的过渡过程。按照(3.2.25)式，性能函数如下式,3.2.46,将(3.2.40)式代入，得到,3.2.47,上式说明性能函数也是按N个指数和的规律变化，和加权系数过渡过程不同的是时间常数不同， 它的时常数为,3.2.48,我们已经知道，性能函数和各个加权系数都是按照N个具有不同时常数的指数和的规律变化的，时常数和特征值成反比，不同的特征

16、值对应的收敛时间是不一样的，但最终的收敛要取决于最慢的指数过程，它的时常数最大，对应最小的特征值，公式如下,3.2.49,3.2.50,但为保证收敛，不能取得太大，受限于最大特征值max。 这样，如果特征值比较分散时，即max和min相差很大时， 使最陡下降法的收敛性能很差。下面分析值的影响。 值收敛过程影响很大，首先必须选择得足够小，使之满足收敛条件,但按照(3.2.47)、 (3.2.48)式，它影响收敛速度。 一般希望在保证收敛的条件下，选大一些，使时间常数小一些，收敛的速度快一些。但当选择得太大时，即使收敛条件满足，也可能形成振动性的过渡特性。 在图 3.2.7 中，图(a)是较小时的

17、情况；图(b)是较大时的情况，此时过渡过程已发生振荡,3.2.4 最小均方(LMS)算法,1. LMS算法的权值计算 LMS(Least Mean Square)算法的梯度估计值用一条样本曲线进行计算，公式如下,3.2.51,因为,所以,3.2.52,3.2.53,FIR滤波器中的第i个权系数的计算公式为,3.2.54,FIR滤波器中的第i个权系数的控制电路如图3.2.8所示, LMS自适应滤波器的总框图如图 3.2.9 所示,图 3.2.8 FIR第i个支路的控制电路,LMS算法的加权系数按照(3.2.53)式进行控制， 式中加权矢量的改变量是2ejXj，梯度的估计值是-2ejXj。显然，这

18、是一个随机变量，这说明LMS算法的加权矢量是随机变化的。因此，LMS算法又称为随机梯度法。下面对这种算法的性能进行分析， 主要分析加权矢理和性能函数的平均变化规律以及它们的随机性造成的影响。 按照(3.2.52)式， 对梯度估计值求统计平均， 得到,3.2.55,上式说明梯度估计值是无偏估计的，梯度的估计量在理想梯度j附近随机变化，权系数也是在理想情况下的权轨迹附近随机变化的,图 3.2.8 LMS自适应滤波器总计算框图,2.LMS算法加权矢量的过渡过程 将误差公式(3.2.4)式代入(3.2.53)式，得到,3.2.56,按照(3.2.53)式， 对加权矢量取统计平均,3.2.57,类似于最

19、陡下降法的推导，经过坐标平移和旋转，变换到V坐标中。其公式推导如下: 令,Vj=Wj-W,3.2.58,那么,EVj=EWj-W,EVj+1=EWj+1-W,3.2.59,将上面两式代入(3.2.57)式中，得到,它的递推解是,令,Rxx=QQ T, =QRxxQT,3.2.60,得到,3.2.61,3.2.62,再将(3.2.59)、(3.2.60)和(3.2.61)式代入上式，得到,EWj=W*+QI-2j Q-1(W0-W*,3.2.63,对比(3.2.35)式，说明LMS算法加权矢量的统计平均值的过渡过程和最陡下降法加权矢量的过渡过程是一样的。换句话说， LMS算法加权矢量是在最陡下降

20、法加权矢量附近随机变化的， 其统计平均值等于最陡下降法加权矢量，那么，其收敛条件同样为,3.2.64,在满足收敛条件的情况下，才有下式,由于最大的特征值max不可能大于R的迹(R的主对角线元素之和)，即,因此收敛条件可以表示为,3.2.65,对于横向滤波器， 式中的迹是NEx2j，即N倍的输入功率， 那么,3.2.66,实际中，通常选得很小，选,3.2.67,同样由(3.2.62)式，第i个分量为,3.2.68,同样引入时常数i,3.2.69,3.2.70,3.2.71,同样，第i个权系数可以表示成,3.2.72,3 . LMS算法性能函数的过渡过程学习过程 由于LMS算法加权矢量的平均值的变

21、化规律与最陡下降法的加权矢量一样，可以推想它的均方误差也会按照最陡下降的均方误差变化规律变化。下面进行推导。 按照(3.2.4)式， 信号误差为,3.2.73,式中，eoptj=dj-XjTW*，称为最佳误差信号，它对应于最小均方误差， 即,按照(3.2.73)式写出均方误差表示式,假定Xj和Vj不相关，上式中最后一项为0，那么,同样，假设加权系数变化很小，Vj也变化很小，EVjVj，这样,类似前面的推导，得到,3.2.74,3.2.75,对照最陡下降法性能曲线(3.2.47)式，LMS均方误差变化规律和最陡下降法完全一样，学习曲线同样近似为几个不同时间常数的指数和,4. 稳态误差和失调系数

22、由上面分析知道，权矢量的平均值可以收敛到它的最佳值， 但权矢量变化过程是随机的，即使其平均值收敛到最佳值，它仍然按照下式,Wj+1=Wj+2ejXj,随机地进行变化，这样使权矢量仍在最佳值附近随机变化， 但均方误差将大于最小均方误差，如图 3.2.10 所示。为此，引入失调系数M，M定义为,3.2.76,图 3.2.10 LMS算法稳态误差,可以推出5失调系数为,3.2.77,或者 M=NPin,3.2.78,式中，N是滤波器的阶数，Pin是输入信号功率。上式说明和输入功率加大都会增加失调系数。 在保证收敛的情况下加大，会提高收敛速度，也说明为了减小失调系数， 应该适当选择收敛速度，以保证收敛

23、速度和失调系数都满足要求,图 3.2.11 是一个LMS自适应滤波器的计算机结果5，阶数N=5，其输入是信号加白噪声， 输入信号功率为1，中心频率是0.03fs(fs为采样频率)，噪声功率为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

24、.40,图 3.2.11 LMS算法的学习曲线,3.3 自适应格型滤波器,3.3.1,由于假设了信号是实的，式中预测误差ep(n)和系数ap,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,Hf(z)称为前向预测误差滤波器的系统函数。前向预测误差滤波器的结构图如图 3.3.1所示,图 3.3.1 前向预测误差滤波器,用均方误差最小的准则求前向预测误差滤波器的最佳系数ap,k,k=1, 2, ，p,3.3.5,将

25、(3.3.2)式代入上式，得到,k=1, 2, 3,，p,3.3.6,上式表明前向预测误差与用于预测的数据正交，这就是对于前向预测误差的正交原理。按照第二章的推导，前向预测误差滤波器的最佳系数ap,k和信号的自相关函数之间的关系式称为Yule-Walker方程式，重写如下,3.3.7,将上式用矩阵方程表示为,3.3.8,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.

26、3.10,图 3.3.2 前向预测数据之间的关系,设后向预测误差用 表示(实际表示的是信号在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,式中， 是后向预测误差的最小误差功率。将(3.3.13)式和(3.3.7)式进行对比，它们极其相似。利用Toeplitz矩阵的性质， 可得到以下重要关系,3.3.14,3.3.15,上面两式表明前、后向预测的最小误差功率相

27、等，系数也相等(如果是复数，则是共轭关系)。由(3.3.10)、 (3.3.11)、 (3.3.14)式得到,3.3.16,式中，当k=0, 1, 2, 3, , p时， p-k=p, p-1, p-2, , 0, 因此也可以写成下式,由上式画出后向预测误差滤波器的结构图如图 3.3.3 所示,图 3.3.3 后向预测误差滤波器,对比图 3.3.1 和图 3.3.3, 或者对比公式(3.3.2)和(3.3.17)， 它们的系数虽然一样，但后向预测误差滤波器的系数排序却是前向预测误差滤波器系数排序的逆转排列。 对(3.3.16)式进行Z变换，得到,3.3.18,后向预测误差滤波器的系统函数为,3

28、.3.19,将上式与前向预测误差滤波器的系统函数(3.3.4)式对比，得到前、 后向预测误差滤波器的系数函数之间的关系是,为了求解前、后向预测误差滤波器的最佳系数，需要解Yule-Walker方程。可以采用高斯消元法解出ap,k(k=1,2, 3,p)以及2p,但需要p3量级运算量。利用Yule-Walker方程中的自相关矩阵是一个埃尔米特(Hermitain)和托布列斯(Toeplitz)矩阵的特点，且至少是半正定的，可以有效地减少运算量，这就是下面要推导的Levinson-Durbin算法，它的运算量级是p2,3.Levinson-durbin算法 Levinson-Durbin算法首先由

29、一阶AR模型开始，按照(3.3.8)式，一阶AR模型(p=1)的Yule-Walker为,由该方程解出,然后增加一阶，即令p=2，按照(3.3.8)式得到,由上面方程解出,然后令p=3, 4, , 以此类推， 可以得到一般递推公式如下,上面(3.3.21)(3.3.25)式就是Levinson-Durbin递推公式， 该式中的kp称为反射系数。在(3.3.24)式中，2p和2p-1是预测误差的均方值，因此1-k2p必须大于等于0，这样kp应要求满足下式,3.3.26,进而得到 ，即预测误差随递推次数增加而减少。把kp称作反射系数，是类似于传输线的情况，如图3.3.4 所示，第p节的输出功率(即

30、下一级的输入功率)等于前一级的输出功率减去本级的反射功率，用公式表示如下,3.3.27,图3.3.4 传输线,3.3.2 格型滤波器 1. 由预测误差滤波器导出格型滤波器 将前面已推导的前向预测误差公式(3.3.2)重写如下,再将系数ap,k(k=1,2,3,p)的递推公式(3.3.23)代入上式，并令kp=ap,p，得到,将上式与(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阶的后向预

31、测误差， 即,这样由(3.3.28)式，得到前向预测误差的递推公式， 即,3.3.29,类似地，得到后向预测误差的递推公式为,3.3.30,利用(3.3.29)式和(3.3.30)式，组成格型滤波器的第p节的结构图, 如图 3.3.5(a)所示,图 3.3.5 全零点格型滤波器,对于p=0的情况， 按照(3.3.2)式和(3.3.11)式， 得到,整个预测误差格型滤波器的结构如图 3.3.5(b)所示。由于没有反馈支路，它是一个全零点格型滤波器。 经过变形还可得到其他类型，如全极点格型滤波器、全极点横向滤波器，等等5,2. 格型滤波器的性质 (1) 各阶后向预测误差相互正交。 用公式表示如下,

32、设ij,按照(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)的线性组合，因此 与 相互正交。 各阶后向预测误差相互正交的结果，使滤波器前后级互相解耦，对于系统最小化问题化为一系列独立的对每一级局部最小化问题。用作自适应滤波时，各级可选用不同的自适应步长， 使收敛速度提高。另外，为提高线性预测性能，需要增加一节或几节，可以只对新增加的级进行独立的调节，达到输出均方误差最小，无需再调节前面的系数,2) 平稳随机序列可由自相关函数或反射系数表征。按照Levins

33、on-Durbin递推公式，已知rxx(0),k1,k2,kp,从一阶开始，可以推出全部的预测系数ap,1,ap,2, ,ap,p和2p，把得到的这些数据代入Yule-walker方程，可求得信号的自相关函数rxx(0), rxx(1), rxx(2),rxx(p)。以上说明平稳随机序列可由自相关函数表征，也可由rxx(0),k1,k2,kp表征。 (3) 前向预测误差滤波器是最小相位滤波器，即它的全部零点在单位圆内,3.3.32,对于前向预测误差的正交原理，则用下式表示,3.3.33,前向预测误差滤波器的预测系数和信号自相关函数之间的Yule-Walker方程仍和(3.3.8)式一样,后向预

34、测误差和后向预测误差滤波器系统函数分别用下式表示,3.3.34,3.3.35,对于后向预测误差的正交原理为,3.3.36,对于复信号的Levinson-Durbin递推公式为,3.3.37,k=1, 2, 3, , p-1,3.3.38,3.3.39,3.3.40,3.3.41,复信号的全零点格型滤波器预测误差递推公式为,3.3.42,3.3.43,图 3.3.6 复信号预测误差全零点格型滤波器,3.3.3 最小均方误差自适应格型滤波器,P阶格型滤波器由p节组成，如果前m节的参数ki(i=1, 2, 3,m)为最佳，相应的预测误差功率是最小，而后面的节的参数对前面的最佳参数无影响，因此在m节的

35、基础上再加一节， 则只需根据使第m+1节的预测误差功率最小的原则选择km+1即可。 预测误差功率有前向预测误差功率和后向预测误差功率，这里采用使前、后向预测误差功率的和为最小的原则求反射系数。 公式为,3.3.44,将(3.3.29)、(3.3.30)代入上式， 可以得到,3.3.45,实际计算时，上式中的统计平均值用时间平均计算， 公式为,3.3.46,对于复信号情况，公式为,3.3.47,上面两式便是直接利用数据计算反射系数的递推公式。下面讨论公式中的求和限问题，如果输入数据为x(i), i=0, 1, 2, , n， 当p=1时,这里,因此,上式中，求和限必须限制在已知的输入数据范围内计

36、算，这样求和限应为i=1,2,3, , n，计算公式为,当p=2时,按照(3.3.29)、(3.3.30)式，得到,将上面两式带入 公式中，可以计算出 ，考虑到输入数据的范围，具体计算公式为,再根据 ,按照(3.3.29)、(3.3.30)式计算e2(i)、b2(i)，按照(3.3.46)式计算 , 以此类推。 这样，对于 具体计算公式为,3.3.48,以上便是直接采用信号数据计算格型滤波器的反射系数以及最小预测误差的方法。但这种算法必须从低阶推起，要求较大的存储时，有较大的计算延迟，使应用受到限制。下面介绍梯度算法， 这种算法可以减少运算量， 且适合非平稳情况,自适应格型滤波器的梯度算法中反

37、射系数的计算，类似于自适应横向滤波器中系数的递推算法，公式为,3.3.49,式中，仍然是控制收敛速度和收敛的参数；kp表示对方括弧中的部分求梯度,将上式代入(3.3.49)式中， 得到,式中，=2，为步长因子。这部分内容可参考文献5,3.3.50,3.4 最小二乘自适应滤波,本节讨论另外一种以误差的平方和最小作为最佳准则的误差准则最小二乘(Least Square)准则。 定义,3.4.1,式中, (n)是误差信号的平方和；ej是j时刻的误差信号,dj是j时刻的期望信号，Xj是j时刻的输入信号构成的向量， W表示滤波器的权系数构成的向量。通过选择W，使(n)取得最小值的滤波称为最小二乘(Lea

38、st Square,简称LS)滤波，而满足Ee2j取得最小值的滤波称为最小均方误差(Least Mean Square, 简称LMS)滤波。和LMS滤波相比，LS滤波对非平稳信号的适应性要强许多，这是由于LS滤波总是采用新的准则，在每一个时刻对所有已输入信号而言，重新评估使其误差的平方和最小，因此具有更精确的含义，属于精确分析法。 而LMS滤波是以集合平均为基础的，属于统计分析的方法,3.4.1 最小二乘滤波,1. 最小二乘的基本问题 已知n个数据x(1),x(2),x(n)，采用M个权的FIR滤波器对数据进行滤波，假设期望信号为d(i)，如图3.4.1所示,滤波器的输出 是对期望信号d(i)

39、的估计,其中，wk(i),k=1,2, ,M，为FIR滤波器在i时刻的M个系数值(说明滤波器的系数可以变化)，它是一个M维的向量，记为wM(i)=w1(i),w2(i),wM(i)T。同理，输入信号也是一个M维的向量，xM(i)=x1(i), x2(i), , xM(i)T。n时刻，估计误差为,3.4.2,误差信号的平方加权和为,上式中，由于数据长度有限，对观测区间以外的数据所做的约定不同，当i的取值范围不同时，得到不同的(n)。这里采用前加窗法，约定,x(i)=0 i0,得到,3.4.3,图 3.4.1 M个权的FIR滤波器,为了后面叙述方面，引入一些符号。令M维向量wM(n)和xM(n)分

40、别表示n时刻的滤波器的权向量和输入信号向量,当i=1,2, ,n时，引入n维误差向量e(n)和期望信号向量d (n)，以及输入信号构成的Mn维矩阵XM(n,e(n)=e(1), e(2), , e(n)T,d(n)=d(1), d2), , d(n)T,XM(n)=xM(1),xM(2),xM(n,3.4.5a,3.4.5b,3.4.5a,为了后面推导方便，引入nM维矩阵C，定义,3.4.5b,应用这些符号，期望信号的估计和估计误差可以表示为,3.4.6,3.4.7,将(3.4.7)式代入(3.4.3)式，得到误差信号能量,式中,是加权矩阵，对角线上的元素称为加权因子。为了推导简单起见， 在后

41、面的分析中，取=I，则,3.4.9,要使(n)取得最小值，满足,3.4.10,成立的w(n)就是wM(n)的最小二乘估计，记为 。 应用标量求导公式， 计算得到,3.4.11,将(3.4.5b)式和(3.4.7)式代入上式，有,3.4.12,将(3.4.11)式展开,3.4.13,引入M维向量pM(n)以及MM维矩阵RM(n,3.4.14,3.4.15,则(3.4.13)式可以写为,3.4.16,可以看出，RM(n)类似于输入信号的自相关特性，pM(n)类似于输入信号与期望信号的互相关特性。(3.4.16)式与第二章中的维纳-霍夫(Wiener-Hopf)方程相似，不同之处在于维纳-霍夫方程中

42、的数学期望符号用求和符号所代替,若矩阵XM(n)的秩等于M，记做rankXM(n)=M,则XM(n)XMT (n)非奇异，求解(3.4.16)式，可以得到wM(n)的最小二乘估计wLS(n,3.4.17,若rankXM(n)M, 则wM(n)不能唯一辨识。 在(3.4.12)式两边同左乘以 ，得,3.4.18,应用(3.4.6)式，得到,3.4.19,图 3.4.2 最小二乘估计的几何解释,当 存在时，最小二乘的估计值 为,最小二乘估计的误差信号能量min为,3.4.20,综合前面的分析，我们可以把最小二乘问题用模型,z=A+n,3.4.21,来描述，其中，z是观测信号，n为噪声信号，A可看作

43、一数据矩阵，表征输入与输出之间的关系，是可调整量。与前面的分析相对应，参见(3.4.7)式，n类似于误差信号e(n)，z类似于信号真值d(n)，A类似于信号的估计值 ，且A矩阵与C矩阵相对应， 在数据矩阵A已经确定的情况下，对z的最小二乘估计为,与图3.4.2的信号向量相对应，观测信号z构成的观测值向量与d(n)相对应，A与估计向量 相对应，噪声信号n与图中的估计误差e(n)相对应,令误差信号能量为J，并取加权矩阵=I，则,3.4.23,2. 最小二乘估计的质量 假设误差向量n(k)是独立同分布的，具有零均值,方差为2, 那么，最小二乘估计是无偏估计; 若噪声是高斯噪声， 则最小二乘估计是一致

44、估计。 (1) 无偏性,3.4.24,2) 一致性：如果噪声是高斯噪声，那么最小二乘参数估计一致收敛， 即,3.4.25,这里，W.P.1表示以概率为1的收敛,证明 首先求解最小二乘估计的协方差，令 表示最小二乘估计的误差， 即,那么，当具有零均值时，最小二乘估计的协方差为,3.4.26,已知噪声为高斯噪声，则EnnT=2I，得到最小二乘估计的协方差为,式中， 将以概率1收敛于一个正定阵，且2是有界的， 因此,最小二乘估计的无偏性已经保证了偏移量为0， 因此， 一致性得证。 前面，所讨论的信号模型都是MA模型，即最小二乘滤波器用FIR滤波器实现，若信号是一ARMA模型，那么怎样来求出信号所对应

45、的最小二乘估计呢? 下面举一例说明,例 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个观测值构成的观测向量表示为zn=z(0), z(1), , z(n)T,和一般的最小二乘估计的模型z=A+n相比，这个问题的关键是如何确定A矩阵。MA模型的输出完全由输入信号所决定，

46、 而ARMA模型的输出还包括以前时刻的输出对当前输出的反馈， 因此，ARMA模型的A矩阵应同时包括输出和输入数据。针对这个例题所选的模型，可以确定A矩阵为,其中，z(i)和x(i)分别表示i时刻的观测值和输入值,解 根据(3.4.17)式，计算,其中，Rz表示z的自相关函数，假设信号具有遍历性，则 ，同理，得到Rx(0),Rxz(0)，Rn(1)等。 考虑n(k)与x(k)是统计独立的及En(k)=0， 记=Rz(0)Rx(0)-Rxz2 (0)，由已知模型关系,其中，a0和b0为参数真值，得到,由此可以看出，只有当n(k)是白噪声时，Rn(1)=0，才能使得,简单地说，所谓的最小二乘估计，是

47、指在每个时刻对所有已输入的信号而言，误差平方和最小。最小二乘法的原理简单， 编程也易实现， 因此获得了广泛的应用,3.4.2 递推最小二乘法(RLS) 递推最小二乘法(Recursive of Least Square，简称RLS)的基本思想就是新的估计值是在老的估计值的基础上修正而成的，为了分析简单，便设为一向量，且仅与当前观测值有关，则,因此，最小二乘递推算法的关键是得到修正项的表达式。根据最小二乘估计式(3.4.22),用a(i)表示第i步迭代时A的取值，Ak表示前k步A的数值构成的向量,定义一个变量P,3.4.29,其中,那么,3.4.30,3.4.31,令观测信号向量为,3.4.32

48、a,3.4.32b,这里，z(j)表示j时刻观测信号的大小。 把(3.4.29)式代入(3.4.22)式， k-1时刻的估计值,3.4.33,上式两边同时左乘P-1(k-1)，得,3.4.34,把(3.4.33)式和(3.4.34)式代入(3.4.22)式， 并应用(3.4.30)式,3.4.35,为了进一步得到加权矩阵的表达式，可定义q(k,3.4.36,则(3.4.35)式可以写为,3.4.37,根据矩阵反演公式,3.4.38,这里，A, C为一个具体数值，代入(3.4.31)式，得,3. 4.39,至此，得到了最小二乘递推算法的一种形式,3. 4.40,也可以将(3.4.39)式代入(3

49、.4.36)式， 整理后得,3.4.41,因此可以得到递推最小二乘的另一种形式,3.4.42,3.4.3 线性向量空间,1. 向量空间的基本概念 1) 希尔伯特空间 我们知道，一组数据可以看作是空间中的一个向量，一组向量可以确定一个空间，称为张成一个空间， 如二维空间可以由两个向量x=1, 0T和y=0,1T所张成，三维空间可以由三个向量x=1, 0, 0T,y=0, 1, 0T和z=0, 0, 1T所张成，那么,M维向量空间，可理解为由M个基本向量uM(n)所张成(这里的基本向量是指这些向量相互独立且正交,M个基本向量uM(n)构成一个数据矩阵，记为U=u1(n), u2(n), ,uM(n

50、)，那么矩阵U所张成的空间，记为U。因此， 一组数据对应于一个向量，一个(数据)矩阵对应于一个空间。 所谓的线性向量空间U, 是由若干个向量uM(n)线性组合得到的，空间的维数是指用来张成该空间的最少的向量数。 通常，我们所研究的向量空间都是希尔伯特空间，它与内积空间有紧密的联系。我们首先给出内积空间的定义,1) 内积空间。若u(n),v(n),w(n)是空间U中的任意向量， 且a为任一数，任意两个向量的内积为一个数，并满足下面的性质,则称U为内积空间,2) 希尔伯特空间。若某空间是一个线性内积空间的完整部分，则称该空间为希尔伯特空间，简称希氏空间。所谓的完整，是指不存在这样的向量，可以任意接

51、近该空间但不属于该空间。 自适应滤波这一章的内容所涉及的空间都是希尔伯特空间， 数学中常用的欧几里德空间也是希尔伯特空间,2) 子空间 (1) 子空间。所谓的子空间，是指用少于空间维数的向量所张成的空间。如在几何中，X-Y平面是立体X-Y-Z的一个子空间； 假设由向量u1(n)=u1(1),u1(2), , u1(n)T和u2(n)=u2(1), u2(2), , u2(n)T分别张成空间u1和u2，由向量u(n)=u1T,u2T T=u1(1),u1(2),u1(n),u2(1),u2(2),u2(n)T张成空间记为U, 那么u1、u2都是U的子空间。 所谓相互正交的两个子空间，是指在这两个

52、子空间上的任意向量u(n), v(n)的内积都等于0，即,3.4.43,2) 两子空间之和。由空间U1的向量u1(n)和子空间U2的向量u2(n)的所有线性组合所张成的空间称为这两个子空间之和，记为U1U2, 简记为U1U2,2. 投影与投影矩阵 这里所说的投影的含义与几何中投影的直观含义相同， M+1维空间中的某一点与原点构成的向量在M维子空间的投影， 可以用该向量与一加权矩阵相乘来表示，这个加权矩阵称作投影矩阵或投影算子。当然，投影矩阵也可以是一个向量。具体的定义如下： 假设在n时刻得到的数据向量为x(n)=x(1),x(2), , x(n-1),x(n)T,U表示一个数据矩阵或数据向量，

53、记为U，由U张成的空间记为U, 定义,3.4.44,称PU为在U上的投影矩阵或投影算子。若U是一个Mn维的矩阵， 则PU是一个MM的方阵,投影定理：给定希尔伯特空间中的子空间U和向量x， 在U中必有唯一向量PUx,使得在U中的任意向量y， 均满足,3.4.45,PUx称为x在U中的投影向量。该式可以理解为向量x-PUx与子空间U正交，或者由该向量张成的一维子空间x-PUx与子空间U正交。(3.4.45)式启示我们， 在进行空间分解时，可以采用正交分解的方法，使得张成空间所需的向量数目最少。 一个M维线性向量空间U可以分解为一组相互正交的子空间,定理：若U1, U2为希尔伯特空间H中的相互正交的

54、子空间，那么，对H中的任一向量x，均满足,3.4.46,即当子空间相互正交时，向量在子空间和上的投影等于在各子空间上投影之和。 若已知一个向量在某子空间U上的投影矩阵为PU,则该向量在U的正交子空间上的投影矩阵称为正交投影矩阵，记为,3.4.47,图 3.4.3 所示是投影矩阵和正交投影矩阵的含义，从向量x在向量u上的投影和正交投影中可以看出，投影矩阵和正交投影矩阵的实质是一加权向量。对于维数较高的子空间，投影矩阵和正交投影矩阵实质上是一加权矩阵,图 3.4.3 投影矩阵和正交投影矩阵,容易证明，投影矩阵和正交投影矩阵具有下面的性质： (1) 幂等性，即,3.4.48,2) 反身性，也称对称性

55、， 即,3.4.49,3) 正交性，即,3.4.50,3. 数据向量的扩充及投影矩阵的更新 若数据矩阵U张成的空间是U，它的投影矩阵和正交投影矩阵分别为PU、 ，当一个新的向量u增加到U中时产生一新矩阵，新矩阵记为(U,u)，由(U,u)所张成的空间为U,u，该空间的投影矩阵和正交投影矩阵分别记为PU,u和 ,u，那么由于向量u的增加带来的投影矩阵的变化可以表示为,PU,u=PU+校正项,PU,u相当于对新数据矩阵U, u进行分解，其中PU是原来的投影矩阵，校正项为新的投影矩阵与原投影矩阵的差或称为误差项。对PU,u进行正交分解，可以得到,3.4.51,说明w向量所张成的空间W与U正交。一般情

56、况下,u与U不一定正交，即uw。因此需要进行的分解就是把向量u向空间U和与U垂直的空间W上进行投影，换言之，w向量就是向量u在与U垂直空间上的投影,图 3.4.4 子空间U,u及其等效子空间U,w,采用正交分解方法，保证了U,u和U,w包含的信息是一样的(如图 3.4.4 所示)。那么空间W上的投影矩阵Pw可以计算， 即,3.4.53,利用投影矩阵的反身性， 上式可以化简为,分别用PU左乘和右乘(3.4.53)式，可以证明Pw与PU满足正交性,把前面得到的投影矩阵的变化应用于实际问题中， 可以看到，当由数据矩阵U变为(U, u)，投影矩阵和正交投影矩阵的更新公式为,3.4.54,若某一向量y(

57、n)在U上的投影存在，记为PUy，当数据矩阵发生变化时， 同一向量y(n)在新空间U,u的投影分别是,3.4.55,研究自适应滤波器的递推算法，当然要研究输出信号的变化，某时刻滤波器的估计误差及误差信号能量都是标量。鉴于标量可以用一广义内积来表示，因此，当空间发生变化时， 标量的更新用投影向量与某一向量内积的变化来表示， 这里假设向量z不变，另一个向量由PUy变化为PU,uy，根据内积运算的性质， 可以得到标量的更新公式,3.4.56,4. 用向量空间法描述最小二乘问题 首先介绍怎样由数据向量构造数据矩阵。已知一组数据x(1),x(2), ,x(n)T构成一个向量x(n)，通过对x(n)平移，

58、 并在前端添0，可以构成一系列向量z-jx(n)(0jn)。这些向量的组合可构成数据矩阵，由这些数据矩阵分别张成相应的向量空间,其中，数据矩阵X1,M(n)和X0,M-1(n)分别为,3.4.57,这里，X1,M(n)是一个nM维的矩阵,这里，X0,M-1(n)是一个nM维的矩阵。根据X0,M-1(n)的定义， 可以得到X0,M-1(n-1)是一个(n-1)M的矩阵，其表达式为,比较上式和(3.4.57)式， 可以看出,3.4.58,将数据矩阵X1,M(n)和X0,M-1(n)张成的向量空间，分别记为X1,M(n)、X0,M-1(n)。 仍用前加窗的方法讨论最小二乘问题。参照图 3.4.1，

59、已知数据x(1), x(2), ,x(n)T，依次通过M个权的FIR滤波器， 第i个误差信号为,1in,将n个误差信号全部写出来，得到,写成向量的形式,根据(3.4.17)式得到最小二乘意义下， 的最佳解及其最佳估计,根据投影矩阵的定义式(3.4.44)，子空间X0,M-1(n)投影矩阵P0,M-1(n)为,得到,5. 抽取参量和角参量 1) 抽取向量(n) 抽取向量(n)用于提取一个向量最新时刻的分量，也称为单位时间向量，或现时向量，为了书写方便，简写为(n,n)=0, 0， ，1T,3.4.59,可以看出，(n)表征的是当前数据向量的方向，因此可以把数据空间分解为过去与当前两个数据子空间。

60、假设一个n维数据向量x(n)=x(1), x(2), ,x(n)T，它的现时分量的值可以表示为,x(n)=(n), x(n,3.4.60,n)到(n)的投影矩阵P(n)和正交投影矩阵 分别为,3.4.61,3.4.62,设定U是一个nM维矩阵，取u=(n)，考虑当U增加了这样一个列向量后，投影矩阵和正交投影矩阵发生的变化。这里， (U,(n)是n (M+1)维的矩阵，因此PU(n)仍是一个nn维的矩阵，将u的取值代入(3.4.54)式，得到,3.4.63,3.4.64,2) 角参量U(n) 我们先给出角参量的定义， 然后解释其名称的由来。角参量定义为,根据投影矩阵的反身性,3.4.65,取U(