CN1084978C - 具有用于最佳适应的卡尔曼滤波器的回声消除方法 - Google Patents

Abstract

卡尔曼滤波器被应用来获得一种快速收敛，低复杂度的回声(y(t))消除器(18)。一种平均的对角线Riccati等式允许以少量的标量Riccati等式来预先计算或在线计算最佳时变自适应增益(t)。

Description

具有用于最佳适应的卡尔曼滤波器的回声消除方法

                        技术领域

该发明涉及一种在电信系统中自适应滤波器类型的回声消除方法。

                        发明背景

在无线电通信系统中，回声控制是一个重要的问题。在公共交换电话网(PSTN)中，由于阻抗不匹配，回声常常在四线到两线的连接处产生。更精确地，在用户和交换机之间的本地回路使用两线连接，而交换机之间的连接使用四线。在两线和四线连接之间的耦合由一种称之为混合电路的电气设备形成。理想地，混合电路应该将所有流入信号能量从远端(远处的说话者)直接传送到近的终端。由于是全双工的连接，当远端说话者在讲话时，近端说话者也必须可能讲话。不幸地是，这意味着由于在混合电路中阻抗不匹配，一些来自于远端的信号能量会被反射回远端。这些混合电路的数目是很多的，因此手动调整它们是不可行的。这种情况是对称的，即，当近端说话时，在远端会产生同样的效果。

在解决回声消除问题的尝试中，曾建议使用自适应滤波器，尤其是相应于NLMS算法修改的滤波器。这种算法由于它的相对低的复杂性而吸引人。然而，该算法的一个缺点是它的慢收敛性。

                        发明概要

因此，提供一个具有改进的收敛时间和低复杂度的回声消除方法是该发明的一个目标。根据本发明在电信系统的混合电路中消除由一个输入信号X(t)产生的回声信号的方法，包括步骤：在一个精细脉冲响应卡尔曼滤波器中对所述输入信号进行滤波，所述卡尔曼滤波器有一个时间变化卡尔曼增益矢量K(t)，所述矢量确定所述滤波器的时间变化滤波器抽头，所述增益矢量K(t)和以下矢量成正比： $\left[\begin{array}{c}{P}_1\left(t\right)x(t-1)\\ P_1\left(t\right)x(t-2)\\ \underset{.}{\overset{.}{.}}\\ P_n\left(t\right)x(t-n)\end{array}\right]$ 其中n表示滤波器抽头数，t表示离散时间，并且P_i(t)，i＝1…n，表示矩阵P(t)的第i个对角线元素，所述矩阵P(t)是所有t的对角，独立于所述输入信号，并且每个所述对角元素满足以下形式的相应标量等式： ${\stackrel{-}{P}}_i(t+1)={\stackrel{-}{P}}_i\left(t\right)+S(t-i)r_i-\frac{S(t-i){\stackrel{-}{P}}_i^2\left(t\right)}{\mathrm{\α}+S(t-i){\stackrel{-}{P}}_i\left(t\right)},i=1...n$ 其中r_I和α是恒定参数，而s(t-i)是单位步长函数；并且从所述回声信号中减去所述滤波器信号。

附图简要描述

该发明以及进一步的目标和它的优越之处可以通过参考下面与附图一起进行的描述而更好地理解。其中：

图1是说明回声产生过程的方框图。

图2是说明用于回声消除的一种自适应滤波器方法的方框图。及

图3是分别相应于以前技术运行的回声消除滤波器的收敛时间和相应于该发明运行的回声消除滤波器收敛时间的比较。

                优选实施方案的描述

图1表示说明回声产生过程的方框图。说话者A在两线连接上讲话。在混合电路10，这一连接被变换成四线连接，语音被送到四线传送路径上的混合电路12，并具有一个由方块14表示的延迟。在混合电路12的另一边的两线连接上，语音到达听者B。同样地，来自B的语音将通过一个四线返回路径到达A，并具有一个由方块16表示的延迟。然而，由于在混合电路10，12中的阻抗误匹配，在返回路径上，A的回声也被送回给A。

因此，流入信号能量的一部分通常被反射回发送方。如果回路延迟较短，即小于10-15ms，说话者不能区别反射信号和来自他/她的终端旁边的音调。然而，当回路延迟大于10-15ms，说话者可以辨别来自回路另一端的回声。感知上，回声会听起来令人讨厌。如果(回声)大，它可以使对话实际上不可进行。注意到，在无线电通信系统中，回声也会由于其它机构产生。例如，在蜂窝电话的免提设备中，回声由扬声器和紧邻的麦克风之间的声耦合产生。

混合电路是一种动态电子电路，即，它必须从数学上被描述为一种动态系统，这意味着来自混合电路的反射信号也依赖于以前的流入信号。这对于设计象回声消除器这样的回声控制设备具有重要的意义。

图2说明了适于采用该发明方法的回声消除器的基本原理。这种思想是简单构造一个电路，该电路在当它被来自远端的信号激励时，模拟混合电路的特性。因为混合电路可以由一个线性动态系统来准确描述，回声消除器被构造成为一个滤波器。出于实际的原因，采用了一种时间离散有限脉冲响应滤波器(FIR)。在电话系统中，采样率是8kHz。此外，使得滤波器脉冲响应具有64ms的最长持续时间是很普通的。这对应于滤波器中的512抽头。

实验表明，混合电路具有变化很大的脉冲响应。因此，有必要使得回声消除器的FIR滤波器自适应。这意味着通过一些测量图2中残留回声

的机制来调整脉冲响应参数(它等于FIR滤波器的抽头值)。调整机制努力调整FIR滤波器的阶数，使得残留回声的方差最小。在下一部分中，描述了一些普通的自适应算法。

            自适应滤波，NLMS和卡尔曼滤波

在[2]中详细描述了自适应滤波算法。这里给出一个简略的概括。通常通过一个适当的损耗函数的最小化来引入自适应。因此，复习一下最普通的算法即，最陡梯度降落搜索方法和更复杂的二阶牛顿法是有意义的。出于该目的，考虑最小化标量函数V(x)的问题，其中x是一个矢量。

最陡降落方法基于一种思想，这种思想建立一种迭代，其中对所搜寻的最小值的下一次估计等于前一次估计加上在V(x)负梯度方向上的一个小的修正。(朝着最小值下降)。这里给出迭代： $\left(1\right)---x_{n+1}=x_n-\mathrm{\μ}{\frac{\∂V^T\left(x\right)}{\∂x}}_{|x=x_n}.$

这里，μ是个常数，n表示迭代次数。用户需提供一个适当的开始值。

牛顿法可以开始以估计值x_n附近的V(x)的泰勒级数展开式来理解。因此有 $\left(2\right)--V\left(x\right)\≈V\left(x_n\right)+{\frac{\∂V\left(x\right)}{\∂x}}_{|x=x_n}(x-x_n)+\frac{1}{2}{(x-x_n)}^T{\frac{{\∂}^{2}V\left(x\right)}{{\∂x}^2}}_{|x=x_n}(x-x_n).$

其中，T表示转置，在最小值情形下，V(x)/x应该等于零。因此，利用近似式(2)来判定下一个估计值使得这一点保持正确。对(2)式微分并使结果等于零，得到(利用²V(x)/x²的对称性) $\left(3\right)--0={\frac{\∂V^T\left(x\right)}{\∂x}}_{|x=x_{n+1}}\≈{\frac{\∂V^T\left(x\right)}{\∂x}}_{|x=x_n}+{\frac{{\∂}^{2}V\left(x\right)}{\∂x^2}}_{|x=x_n}{(x-x_n)}_{|x=x_{n+1}}$ $={\frac{\∂V^T\left(x\right)}{\∂x}}_{|x=x_n}+{\frac{{\∂}^{2}V\left(x\right)}{\∂x^2}}_{|x=x_n}(x_{n+1}-x_n).$

从迭代式中解出x_n+1 $\left(4\right)--x_{n+1}=x_n-{\left({\frac{{\∂}^{2}V\left(x\right)}{\∂x^2}}_{|x=x_n}\right)}^{-1}{\frac{\∂V^T\left(x\right)}{\∂x}}_{|x=x_n}$

这就是牛顿迭代。

注意到，上述V(x)是一个标量，V(x)/x是一个与x维数相同的矢量。而²V(x)/x²是一个方阵，它的行数和列数等于x的维数((²V(x)/x²|x＝x_n)^-1可以被看作为常数标量μ对依赖于时间的矩阵形式的一个产生式)。因此，众所周知的是

在每一次迭代中，牛顿迭代要比最陡降落法复杂的多。

然而，牛顿迭代却有一个重要的有利之处：

至少在接近最小值时，牛顿法比最陡降落法收敛的快的多。

下一步是应用上述方法来构造自适应FIR滤波算法。出于该目的，必须构造一个损耗函数。用y(t)来表示混合电路的输出，用y^(t，h)来表示自适应FIR滤波器的输出。这里，t表示离散时间(此后t代替了上述的迭代索引n)。h表示FIR滤波器的输出依赖于集中在参数矢量h中的未知的抽头值h_i(将被自适应算法优化)。残留回声由下式给出： $\left(5\right)--\mathrm{\ϵ}(t,h)=y\left(t\right)-\hat{y}(t,h).$ 为了继续，需要确定FIR滤波器，它由下式给出。 $\left(6\right)--\hat{y}(t,h)=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{h}_{i}x(t-i)=C\left(t\right)h$

(7)               C(t)＝(x(t)…x(t-n))

(8)               h＝(h_o…h_n)^T.

现在可以构造一个损耗函数，根据该函数可以构造一种最陡降落类型的算法。该算法被称作为最小均方算法(LMS)。并普遍用于很多应用。该损耗函数为 $\left(9\right)--V\left(h\right)=\frac{1}{2}{\mathrm{\ϵ}}^2(t,h)$

其中，目标是通过最小化误差ε(t，h)来估计h。

由于最小化应该采用算法(1)。可以得到下述迭代 $\left(10\right)--\hat{h}(t+1)=\hat{h}\left(t\right)-\mathrm{\μ}{\frac{\∂V^T\left(h\right)}{\∂h}}_{|h=\hat{h}\left(t\right)}$

并且，因此有必要计算V(h)的梯度。可以利用(5)-(8)来完成这一点。考虑梯度矢量中的第j个元素 $\left(1\right)--\frac{\∂V\left(h\right)}{\∂h_j}=\frac{1}{2}\frac{\∂{\mathrm{\ϵ}}^2(t,h)}{\∂h_j}=\frac{\∂\mathrm{\ϵ}(t,h)}{\∂h_j}\mathrm{\ϵ}(t,h)$ $=\frac{\∂}{\∂h_j}(y\left(n\right)-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{h}_{i}x(t-i))\mathrm{\ϵ}(t,h)=-x(t-j)\mathrm{\ϵ}(t,h).$

因此，根据下式更改第j阶抽头值。 $\left(12\right)--{\hat{h}}_j(t+1)={\hat{h}}_j\left(t\right)+\mathrm{\μx}(t-j)\mathrm{\ϵ}{(t,h)}_{|h=\hat{h}\left(t\right)}.$

所有抽头的梯度等式可以被叠加在每一个之上并公式化为下述矢量：见(5)-(8) $\left(13\right)--{\frac{\∂V^T\left(h\right)}{\∂h}}_{|h=\hat{h}\left(t\right)}=-C^T\left(t\right)\mathrm{\ϵ}{(t,h)}_{|h=\hat{h}\left(t\right)}.$

这给出下面定义LMS算法的矢量迭代： $\hat{h}(t+1)=\hat{h}\left(t\right)+\mathrm{\μ}{C}^T\left(t\right)\stackrel{.}{\mathrm{\ϵ}}{(t,h)}_{|h=\hat{h}\left(t\right)}$ $\left(14\right)--\hat{y}{(t+1,h)}_{|h=\hat{h}\left(t\right)}=C(t+1)\hat{h}\left(t\right)$ $\mathrm{\ϵ}{(t+1,h)}_{|h=\hat{h}\left(t\right)}=y(t+1)-\hat{y}{(t+1,h)}_{|h=\hat{h}\left(t\right)}.$

应该注意到，在LMS算法中的更改步长是与C(t)和ε(t，h)|h＝h^(t)的大小成比例的。尤其是，稳定性可以被表示为很大程度上依赖于远端信号x(t)的能量。为了避免这一点，LMS算法通常以输入信号的能量来归一化C(t)C^T(t)。这给出广泛使用的归一化的最小均值平方算法。 $\left(15\right)--\hat{h}(t+1)=\hat{h}\left(t\right)+\frac{\mathrm{\μ}}{C\left(t\right)C^T\left(t\right)}{C}^T\left(t\right)\mathrm{\ϵ}{(t,h)}_{|h=\hat{h}\left(t\right)}$ $\hat{y}{(t+1,h)}_{|h=\hat{h}\left(t\right)}=C(t+1)\hat{h}\left(t\right)$ $\mathrm{\ϵ}{(t+1,h)}_{|h=\hat{h}\left(t\right)}=y(t+1)-\hat{y}{(t+1,h)}_{|h=\hat{h}\left(t\right)}.$

该算法是现在大多数回声消除器中应用的技术。

卡尔曼滤波器可以被理解为一种二阶牛顿算法的版本。它的导出式是有些复杂的[5]，并且不在次重新推导。为了这一目的，对于附录中假设A1中给出的模型，说明卡尔曼滤波器由(16)，(17)，(18)给出是充分的。 $\left(16\right),\hat{h}(t+1)=\hat{h}\left(t\right)+K\left(t\right)(y\left(t\right)-C\left(t\right)\hat{h}\left(t\right))$ $\left(17\right),K\left(t\right)=\frac{P\left(t\right)C^T\left(t\right)}{R_2\left(t\right)+C\left(t\right)P\left(t\right)C^T\left(t\right)}$ $\left(18\right)--P(t+1)=P\left(t\right)+R_1\left(t\right)-\frac{P\left(t\right)C^T\left(t\right)C\left(t\right)P\left(t\right)}{R_2\left(t\right)+C\left(t\right)P\left(t\right)C^T\left(t\right)}.$

其中，h^(t)和C(t)连同假设A1一起在附录中被解释。R2(t)是来自近端的附加噪声的方差(例如，背景噪声)。R1(t)是技术上所谓的系统噪声，它说明，在抽头是常数的假设中，人们可以有多大的自信度。见附录假设A1中第三个等式。P(0)是P等式的开始值(Riccati等式)。这个量等于抽头hi的先验协方差。所有这些量被假设为至少在时间t是可以知道。

给定p(0)，卡尔曼滤波器以(18)的一次迭代开始。这产生P(1)，然后可计算出K(1)。此后，假设可得到h^(1)(由于h^(0)是用户提供的，K(0)由P(0)得到)，h^(2)可被计算出。

注意到，如果FIR滤波器采用512抽头，卡尔曼滤波器P(t)就是一个512^*512阶的矩阵，这样与NLMS算法比较，该算法的复杂度是巨大的。但是它最佳快速收敛。这种快速收敛是由于增益矢量K(t)(对每一抽头可被看作是一个单个的μ)是时变的。这种时变性是针对小的t的，当对估计的抽头一无所知时，增益高是为了快速地获得一个粗略的估计。然后增益逐渐减小以获得高准确性并抑制扰动。

该发明的目标在于引入卡尔曼滤波器的简化，将它的计算复杂度降低到相似于NLMS算法复杂度的水平。

注意到计算复杂度的主要来源是Riccati等式(18)。该目标可以按下面实现。

简化必须被引入到卡尔曼滤波器，使得对所有t，P(t)简化为一个对角线矩阵。这样，对于一个512抽头长度的滤波器，复杂度降低了大约512倍。这些简化必须以这样一种方式来完成，就是在简化了的算法中，至少保留了一些时间最优性。

在附录中简要给出了这点的具体做法。其结果是具有近似最优时变性的块NLMS类型算法。

这使算法具有：

或者如同NLMS算法一样，对于同样的回声抑制具有降低的收敛时间，或者从初始起在相同的的收敛时间上，具有比NLMS算法更高的回声抑制。

调整可以根据有物理意义的量例如回声功率，回声波形，阶数和算法块来计算出的参数。

后一个特性给出下述该发明的第二个优点：

最初，用一小组调整参数使该算法收敛。在收敛之后，可得到回声能量和回声波形的更准确的估计。然后可得到该算法的新的，改进的调整，它继而产生改进的回声抑制。

导出式的启发性解释

推导是基于[3]和[4]的。在[3]中，介绍了一种块NLMS算法，其中，每一块的调整增益μ之间的关系遵循一种指数关系。该想法是：脉冲响应的早几抽头可能比后几抽头更大而且更重要。如果对重要的早几抽头施加更大的增益，而对不重要的抽头施加较小的增益，性能就可以改进。对不重要的小抽头施加较小增益的原因是，否则的话它们对由背景噪声造成的所谓的误调所起的作用与重要的早几抽头相等。但是，它们对回声作用不大，因此它们的等价SNR很低。因此，在模拟来自混合电路的回声时，它们在信号中的加权(坏信息)应该比较小。

为了继续附录中的描述，首先定义滤波器模型和卡尔曼滤波器的数量。如同注意到的，对(18)简化的一个主要障碍是C(t)预先不知道。因此，P(t)不能预先计算出。该思想是：为了获得一个可预先计算的量，在Riccati等式上简单的进行常规的统计平均。假设A2和A3表达了实现这一点需要什麽，并产生附录中的(20)。为了继续，需要关于流入远端信号白化的附加假设(A5)，以及技术假设A4。于是，如果用户选择R1(t)并根据A6和A7选择了P(0)，正如上面想要作的一样，Riccati矩阵等式(18)简化为n标量近似Riccati等式。然而，留下一个问题。精确的Riccati等式总是稳定的，并具有一些吸引人的特征[5]。由于附录中的(22)是近似的，引入了变换和另外的条件以得到精确的(23)。然后，对应于(26)，附录计算平均卡尔曼增益，并引入归一化，参照上面关于LMS和NLMS的讨论。

之后，附录给出用于图2中滤波器18中的详细的算法(AKFA)，并讨论怎样选取各个参数以获得好的性能。首先考虑算法，第一行初始化能量估计和滤波器估计。第二行引入标量Riccati等式的初始化。第三行是流入的远端信号功率的估计器，该功率是第五行中卡尔曼功率增益的归一化所必须的。卡尔曼增益计算也需要每一个抽头不同Riccati等式值。如果需要，这些可以被预先计算并存储在存储器中。对应于上面处理的NLMS算法，该算法最后一行描述了滤波器抽头的更改。

然后剩下选择α，pi(0)和γi的值。这在算法以下给予解释，结果是：

可以首先根据先验假设的最大回声功率(Power)来判定pi(0)，在算法中的块的数目(m)，先验假设的回声脉冲响应的包络波形(power，δ)及滤波器抽头数(n)。

然后根据p1(0)，滤波器抽头数及一个稳定因子(f1)判定α。

最后，可以根据另一个稳定因子α，滤波器抽头数和P1(0)判定γ1。然后根据判定Pi(0)时计算出的量得出γi。

一个非常重要的观察结果是所有的调整参数可以根据直接关联到的物理量如回声功率，回声波形，算法块，滤波器长度和两个容易判定的稳定因子的量值计算出。这使得上面描述的再优化过程可能成为该发明的一个特别的优点。

实验结果表明，与传统的NLMS算法相比，确定的抑制水平的建立时间可以减小5倍。图3表示了分别用于该发明的NLMS算法和卡尔曼算法的误差信号建立时间的比较。从图中可以看到，NLMS算法的误差在近似15秒后收敛到一个稳定不变的状态，而相应于该发明方法的误差在经过近似3秒后收敛到一个稳定不变的状态(该图根据15秒8000Hz采样的语音获得，但只表示出了每个第100个误差)。该例以针对给定回声抑制的收敛时间的形式来比较NLMS算法和AKFA算法。x(t)是一个零均值白信号，均匀分布且方差为1.00。高斯测量噪声被加在模拟回声信号上(由一个实际的回声冲击响应产生)。使得y(t)的SNR为10db。NLMS算法使用一个自适应增益0.01954且n＝512。通过利用P功率＝0.25(-60db)，n＝512，m＝16，δ＝0.001，f1＝1，fc＝0.07815及N＝512，AKFA被应用于同样的数据。

另一种方法，保持建立时间不变，回声抑制可以提高大约7db。

两种可选方法都可以以计算复杂度的少量百分比增加为代价来实现。

该领域的技术人员会明白，可以不叛离由附加权利要求定义的该发明的思想和范围而对该发明做各种的修正和改变。

                            附录

一个在回声消除方面广泛讨论的问题是回声消除器的最初收敛时间。这个问题涉及到一个事实，就是传统的NLMS算法的收敛时间受到回声产生信号的谱特性和估计FIR滤波器长度的限制。该附录的目的是获得具有减小的最初收敛时间的，低复杂度时变矩阵步长的NLMS类型算法。

平均的卡尔曼滤波

用h(t)来表示自适应FIR滤波器的未知参数矢量，其中t是离散时间。回声产生系统的输入x(t)被集中到矢量C(t)中。为了模拟滤波器初始化，假设x(t)＝0，t＜0，由y(t)来表示反射回声的测量值。这些测量值被掺杂了噪声e(t)。其中e(t)是具有零均值和协方差E[e(t)e(s)]＝δt，sR2R2＞0的高斯白噪声。这里E表示期望。该讨论是基于描述h(t)的展开的一种随机移动模型。系统噪声w(t)被假设为具有零均值及协方差矩阵E[w(t)w^T(s)]＝δt，sR1(t)，R1(t)＞＝0的高斯白噪声。卡尔曼滤波器的推导因此基于假设：

A1回声产生系统由下式描述h(t)＝(h₁(t)…h_n(t))^T

           C(t)＝(x(t-1)…x(t-n))，x(t)＝0，t＜0

                        h(t)＝h(t-1)+ω(t)

                        y(t)＝C(t)h(t)+e(t)

           E[e(t)e(s)]＝δ_t，sR₂，E[ω(t)ω^T(s)]＝δ_t，sR₁(t)

[5]现在给出下面的最佳一步预测器 $\left(19\right)--\hat{h}(t+1|t)=\hat{h}\left(t\right|t-1)+K\left(t\right)(y\left(t\right)-C\left(t\right)\hat{h}\left(t\right|t-1)).$

这里K(t)是根据Riccati等式[5]计算出的卡尔曼增益。C(t)预先不知道。因此有必要实时迭代全Riccati等式。然而，根据[4]，如果不同的信号实现之间P(t)的展开没有太大差别，应该可能以一个对应的集平均来代替P(t)。因此，假设：

A2：P(t)是一个具有现存集平均E[p(t)]＝ P(t)的矩阵值随机过程。

A3：ΔP(t)＝p(t)- P(t)很小。

利用A2计算Riccati等式左边和右边的平均值，并如同在[4]中一样准确地使用A3，给出下面的近似平均等式。 $\left(20\right)--\stackrel{-}{P}(t+1)=\stackrel{-}{P}\left(t\right)+R_1\left(t\right)-\stackrel{-}{P}\left(t\right)E\left[\frac{C^T\left(t\right)C\left(t\right)}{R_2+C\left(t\right)\stackrel{-}{P}\left(t\right)C^T\left(t\right)}\right]\stackrel{-}{P}\left(t\right).$

在[4]中使用了x(t)的常数系数特性以进一步简化(20)。在回声消除情况中，这一点是不可能的，需要一种新的方法。注意到R1(t)/R2→0意味着当t→∞时，P(t)→0。因此，如果R1(t)/R2很小，要求R2＞＞C(t) P(t)C^T(t)是合理的，如下式表示的：

A4：A4：R₂＞＞C(t) P(t)C^T(t)t.

这将(20)变换为一个期望表达式的分母被简化为标量R2的等式，该标量R2可以被从期望表达式中提出来。由于获得了低的计算复杂度，E[C^T(t)C(t)]的对角线型的可预计算近似式现在成为优选的。这一点由下式保证：

A5：x(t)是E[X(t)S(t)]＝S(t)δt，sσ²的零均值信号。

这里，S(t)是在t＝0时的单位步长函数，参照A1。当包括用于预白化和x(t)的偏差补偿的装置时，A5是合理的。这将(20)简化为： $\left(21\right)--\stackrel{-}{P}(t+1)=\stackrel{-}{P}\left(t\right)+R_1\left(t\right)-\frac{{\mathrm{\σ}}^2}{R_2}\stackrel{-}{P}\left(t\right)\mathrm{diag}{\left\{S(t-i)\right\}}_{i=1,...,n}\stackrel{-}{P}\left(t\right).$

为了降低复杂度，通过下式，(21)被简化为一系列的标量Riccati等式。

A6…R₁(t)＝diag{r_i(t)}_{i＝1，…，n}＝diag{S(t-i)r_i}_{i＝1，…，n}.

A7… P(0)＝diag{ p_i(0)}_{i＝1，…，n}＝diag(p_i(0)}_{i＝1，…，n}＝P(0).

A6和A7表明 P(t)＝diag{ Pi(t)}i＝1，…n，.t，其中 Pi(t)遵从标量等式： $\left(22\right)--{\stackrel{-}{p}}_i(t+1)={\stackrel{-}{p}}_i\left(t\right)+S{(t-i)}_{r_i}-S(t-i)\frac{{\mathrm{\σ}}^2}{R_2}{\stackrel{-}{p}}_i^2\left(t\right).$

由于(22)不是Riccati等式，可能不会保持 Pi(t)＞0，t。并且(22)会变的不稳定。然而，如果(22)中的R2可以由R2+S(t-i)σ² Pi(t)替换，(22)的结构将成为Riccati类型，它将保证稳定性及 Pi(t)＞0，t[4]。考虑A4。注意到由于 P(t)＝diag{ Pi(t)}i＝1…n，它保持C(t) P(t)C^T(t)＜sup_i，t| Pi(t)|nσ²。因此，如果R2＞＞sup_i，t｜pi(t)|nσ2，那麽A4为真并且R2＞＞S(t-i) Pi(t)nσ²。在进一步简化R₂＞＞sup_i，t| Pi(t)nσ²以得到A4’之前，产生的Riccati等式被表示为： $\left(23\right)--{\stackrel{-}{p}}_i(t+1)={\stackrel{-}{p}}_i\left(t\right)+S(t-i)r_i-\frac{S(t-i){\mathrm{\σ}}^2{\stackrel{-}{p}}_i^2\left(t\right)}{R_2+S(t-i){\mathrm{\σ}}^2{\stackrel{-}{p}}_i\left(t\right)}.$

通过注意到当t-＞∞时，(23)有稳定不变的解(24)，上面的假设现在可以被进一步的修改 ${\left(24\right)----\stackrel{-}{p}}_i^{+}(r_i,R_2,\mathrm{\σ})=\frac{r_i}{2}+\sqrt{\frac{r_i^2}{4}+\frac{r_i{R}_2}{{\mathrm{\σ}}^2}}>0.$

从(23)也可直接看出，对于t≥i $\left(25\right)--\mathrm{sign}({\stackrel{-}{p}}_i(t+1)-{\stackrel{-}{p}}_i\left(t\right))=\mathrm{sign}({\stackrel{-}{p}}_i^{+}(r_i,R_2,\mathrm{\σ})-{\stackrel{-}{p}}_i\left(t\right)).$

因此，Riccati等式(23)的稳定性和正值性，连同R2＞＞supi，t|pi(t)|nσ2，(24)和(25)一起促使了为什麽以下式替换A4是合理的。 $A{4}^{\′}:R_2>>{\max }_i{p}_i\left(0\right)n{\mathrm{\σ}}^2\mathrm{and}{\max }_i{p}_i\left(0\right)>{\max }_i{\stackrel{-}{p}}_i^{+}(r_i,R_2,\mathrm{\σ}).$

现在，平均的卡尔曼增益矢量可以根据A1而计算出(参照[2]，[4])。利用A2-A4’和(23)得到： $\left(26\right)--K\left(t\right)\≈\frac{\stackrel{-}{P}\left(t\right)C^T\left(t\right)}{R_2+C\left(t\right)\stackrel{-}{P}\left(t\right)C^T\left(t\right)}$ $\≈\frac{1}{R_2}{({\stackrel{-}{p}}_1\left(t\right)x(t-1)...{\stackrel{-}{p}}_n\left(t\right)x(t-n))}^T=\stackrel{-}{K}\left(t\right).$

由于涉及到几个近似，平均的卡尔曼滤波器要求稳定性检查。可以使用LMS分析中的技术来进行这种稳定性分析，在LMS分析中，t-＞∞并且得到了(23)的稳定不变的解。对于大的t，(24)和(26)给出： $\left(27\right)--\stackrel{-}{K}\left(t\right)\≈\frac{1}{R_2}(\left(\underset{t\→\∞}{\lim }{\stackrel{-}{p}}_1\left(t\right)\right)x(t-1)...\left(\underset{t\→\∞}{\lim }{\stackrel{-}{p}}_n\left(t\right)\right)x(t-n){)}^T.$

然而，通过利用(24)计算

及LMS稳定性理论，产生了条件A8 $A8:0<{\max }_i(\frac{r_i}{2R_2}+\sqrt{{\left(\frac{{r_i}^{\&prime;}}{2R_2}\right)}^2+\frac{r_i}{R_2{\mathrm{\&sigma;}}^2}})<\frac{2}{n{\mathrm{\&sigma;}}^2}.$

算法和调整

通过将R2看作为一个设计变量并根据(28)选取R2，现在可以得到一种NLMS类型的算法。

(28)       R₂＝ασ².

如果Nσ²被在线估计，根据A1，(19)，(23)，(26)可获得下面NLMS样的方案：

平均的卡尔曼滤波器算法(AKFA)该算法概括为以下各式。： ${\widehat{\mathrm{\&sigma;}}}^2(-1)=0;\hat{h}\left(0\right|-1)=E\left[h\left(0\right)\right|-1]$

         p_i(0)＝p_i(0)，i＝1，…，n

${\stackrel{-}{p}}_i(t+1)={\stackrel{-}{p}}_i\left(t\right)+S(t-i)r_i-\frac{S(t-i){\stackrel{-}{p}}_i^2\left(t\right)}{\mathrm{\&alpha;}+S(t-i){\stackrel{-}{p}}_i\left(t\right)}.$ $\stackrel{-}{K}\left(t\right)=\frac{N}{\mathrm{\&alpha;}{\widehat{\mathrm{\&sigma;}}}_N^2\left(t\right)}({\stackrel{-}{p}}_1\left(t\right)x(t-1)...{\stackrel{-}{p}}_n\left(t\right)x(t-n){)}^T$

        C(t)＝(x(t-1)…x(t-n))^T，x(t)＝0，t＜0 $\hat{h}(t+1|t)=\hat{h}\left(t\right|t-1)+\stackrel{-}{K}\left(t\right)(y\left(t\right)-C\left(t\right)\hat{h}\left(t\right|t-1)).$

如果需要， Pi(t)可以被预先计算并存储。对σ² _N(t)([ ]＞ε)的正值检测防止在卡尔曼增益计算中可能被零除的情况。

现在可以讨论pi(0)的选择。由于没有关于脉冲响应的详细的先验信息，传统的卡尔曼滤波器理论[5]，和A7表明E[h(0)-1]＝0并且pi(0)＝E[hi²(0)-1]，作为一个例子，作为i的函数，下面pi(0)的分段常指数分解式如在[3]中被引入。 $\left(29\right)--p_i\left(0\right)=\mathrm{\&beta;}{e}^{-r\left[\frac{i-1}{l}\right]}--l,\frac{n}{l}=m\&Element;Z_{+}.$

其中m等于分段常数间隔的数目。[(I-1)/l]表示(I-1)/l的整数部分。如果由Power来表示回声脉冲响应功率的一个先验上界，它遵循Power等于pi(0)，I＝1....n的和。将根据(29)产生的几何级数相加并解出β得到： $\left(30\right)\mathrm{\&beta;}=\frac{1}{n}\mathrm{Power},m=1;\mathrm{\&beta;}=\frac{1}{l}\mathrm{Power}\frac{1-e^{-r}}{1-e^{-\mathrm{mr}}},m>1.$

其中γ可以利用(29)通过确定与抽头1相比抽头n的残留功率来判定。如果该残留功率由因子确定，它遵循βe^-r(m-1)＝p_n(0)＝δp₁(0)＝δβ这给出： $\left(31\right)--\mathrm{\&delta;}=1,m=1;r=-\frac{\log \left(\mathrm{\&delta;}\right)}{m-1},m>1.$

通过选择δ，n，m和Power，相应的Pi(0)可以被计算出。这里有m个不同的Pi(0)值。如果根据与Pi(0)相同的原则来选取γi这意味着只有m个Riccati等式需要被迭代，见下面。通过注意到如果Pi(0)＝P_i+1(0)且如果γi＝γI+1则 Pi+1(t)＝ Pi(t-1)来处理乘以S(t-I)的操作。这可以将复杂度从与3n成比例降到与(2n+m)成比例(NLMS算法复杂度为2n)。由于通常n＞＞m，这表示接近50％的降低。考虑到 Pi(t)的存储，在规定的ROM中，这种降低更大。

现在(28)将第一条件A4，变换为 $\left(32\right)--\mathrm{\&alpha;}=f_{1}n\underset{i}{\max }{p}_i\left(0\right)=f_{1}n{p}_1\left(0\right).$

其中f1是稳定因子。为了判定γi，首先注意到 P⁺ _i(γI，R2，σ)是γi的单调递增函数。下面的必要条件用于判定允许的γi的范围，于是根据A4′，(28)和(29)产生： $\left(33\right)--\underset{r_i}{\max }(\frac{r_i}{2}+\sqrt{\frac{r_i^2}{4}+r_i\mathrm{\&alpha;}})<p_1\left(0\right).$

由于E[hi(t)-hi(0)]＝tri，通过随机移动模型A1，γi涉及第I个滤波器抽头的平均变化率。相应的时间因此是t＝pi(0)/γi。有理由假设对于所有抽头，pi(0)有同样的平均漂移时间。且(29)意味着 $\left(34\right)--{\mathrm{\&tau;}}_i={\mathrm{\&tau;}}_1{e}^{-\mathrm{\&gamma;}\left[\frac{i-1}{l}\right]},i=2,...,n.$

稳定条件A8也必须满足。利用(28)来推导A8很简单。其结果，连同(33)和(34)一起可以被集中在(35)。 $\left(35\right)--\frac{r_1}{2}+\sqrt{\frac{r_1^2}{4}+r_1\mathrm{\&alpha;}}=f_2\min (\frac{2\mathrm{\&alpha;}}{n},p_1\left(0\right))$ 其中f2是一个稳定因子。最终结果是： $\left(36\right)--r_1=f_2^2\frac{{\left(\min (\frac{2\mathrm{\&alpha;}}{n},p_1\left(0\right))\right)}^2}{\mathrm{\&alpha;}+f_2\min (\frac{2\mathrm{\&alpha;}}{n},p_1\left(0\right))}$

                        参考文献

[1]E.Hansler，“免提电话问题——带注释的参考书目”，信号处理，vol.27，no.3，pp.259-271，1992.

[2]H.Schutze和R.Zen，“快速递归最小平方正割自适应算法的数字特性——一种比较研究”信号处理，vol.27，no.3，pp.317-332，1992。

[3]S.Makino，Y.Kaneda和N.Koizumi，“基于一种房间脉冲响应统计的指数加权步长NLMS自适应滤波器”，IEEETrans.声学，语音，音频处理，vol.1，no.1，pp.101-108，1993.

[4]L.Lindbom，“衰减移动无线电通道的简化的卡尔曼估计：在LMS计算负载下的高性能”，inProc.ICASSP，Pt.III，pp.352-355，1993.

[5]S.Haykin，现代滤波器，Macmillan出版，1989.

Claims (2)

1.在电信系统的混合电路中消除由一个输入信号X(t)产生的回声信号的方法，包括步骤：在一个精细脉冲响应卡尔曼滤波器中对所述输入信号进行滤波，所述卡尔曼滤波器有一个时间变化卡尔曼增益矢量K(t)，所述矢量确定所述滤波器的时间变化滤波器抽头，所述增益矢量K(t)和以下矢量成比例： $\left[\begin{array}{c}{P}_1\left(t\right)x(t-1)\\ P_1\left(t\right)x(t-2)\\ \underset{.}{\overset{.}{.}}\\ P_n\left(t\right)x(t-n)\end{array}\right]$

其中n表示滤波器抽头数，t表示离散时间，并且P_i(t)，i＝1…n，表示矩阵P(t)的第i个对角线元素，所述矩阵P(t)是所有t的对角，独立于所述输入信号，并且每个所述对角元素满足以下形式的相应标量等式： ${\stackrel{-}{P}}_i(t+1)={\stackrel{-}{P}}_i\left(t\right)+S(t-i)r_i-\frac{S(t-i){\stackrel{-}{P}}_i^2\left(t\right)}{\mathrm{\&alpha;}+S(t-i){\stackrel{-}{P}}_i\left(t\right)},i=1...n$ 其中r_I和α是恒定参数，而s(t-i)是单位步长函数；并且从所述回声信号中减去所述滤波器信号。

2.根据权利要求1所述的方法，其中所述矩阵P(t)是由对角块矩阵形成的，在每个块矩阵中所述对角元素 P_i(t)相同。

Images