CN103561185A - 一种稀疏路径的回声消除方法 - Google Patents

Abstract

本发明涉及一种针对稀疏系统的回声消除技术，属于信号处理技术领域。包括步骤：一、采样产生一个采样序列

Description

一种稀疏路径的回声消除方法

技术领域：本发明涉及一种针对稀疏系统的回声消除技术，基于凸组合和梯度下降的原理、采用自适应滤波方法实现对稀疏路径的辨识，属于信号处理技术领域。

背景技术：在车载电话、视频会议等免提通话系统中，远端话音会被受话器捕获并传回给对方，形成回声，干扰正常通话。在回路中采用自适应滤波器，对远端语音信号传输的路径进行建模，则其输出信号就是回声信号的一个复本，从近端信号中将这个复本减掉，即可实现回声消除。因此，回声消除技术的核心是采用何种自适应滤波方法有效地实现回声路径的建模，其实质是系统辨识问题。

目前回声消除技术面临的问题之一是很多系统的时间响应很长，相应地要求自适应滤波器也要具有足够的长度，从而导致其收敛速度、稳态失调等性能有所下降。但该类系统往往具有一种特殊的稀疏性，即只有一小部分的系统冲激响应分量具有较大幅值，其它的分量均很小甚至为零。如VoIP(Vioce over IP)网络中，传输时延是不确定的，且存在较大的抖动范围，从而使得网络回声路径的长度变得特别大，但其回声路径中大系数却很少；又如一般办公室或会议室环境里,由于地板墙壁等设施的吸音系数小，而导致混响时间达到64ms甚至128ms，却只有很少的非零系数。这些系统的这种特殊性有力地推动了系数比例自适应滤波方法的发展。该类方法通过对滤波器的不同权系数分配不同的步长，来保证每个权均独立地进行更新。同时要求该步长与所估计的滤波器权值成比例，即冲激响应中占主导地位的、较大的系数获得较大的步长，以加快整体收敛速度。第一个系数比例自适应滤波方法是由Duttweiler针对网络回声消除提出的PNLMS（Proportionate Normalized Least Mean Square），此后其改进方法被相继提出。例如，针对车载电话、电话会议等系统中声环境的温度和压力变化或者电话持有者的位置变化会导致系统的稀疏性在一个较大的范围内变化，人们提出了对不同系统稀疏度均适用的IPNLMS（Improved PNLMS）方法。将IPNLMS中比例系数的思想直接应用于APA（Affine Projection Algorithm）中，得到了IPAPA(Improved Proportionate Affine Projection Algorithm)，这种方法同时具有APA的快速收敛速度和IPNLMS对系统稀疏度的广泛适用性，且对于相关性很高的语音信号尤其有效。但该类方法一方面需在收敛速度和收敛精度这两个方面进行慎重地折中选择，另一方面针对不同稀疏度的系统也必须进行谨慎的参数选择，从而限制了方法的灵活性。

发明内容：

发明目的：本发明一种稀疏路径的回声消除方法，其目的是为了克服现有IPAPA方法的不足，提供一种基于凸组合和EG（EG:Exponentiated Gradient）法的自适应滤波方法，以简化参数选择和计算量，并同时提高收敛和跟踪性能。

技术方案：本发明是通过以下技术方案实施的：

一种稀疏路径的回声消除方法，其特征在于：该方法包括以下步骤：

步骤一、获得一个采样序列x(n)，其中n＝1,2,…,N，N为采样序列x(n)的采样点数量；所述的采样序列x(n)是一维信号，且其中包含N个采样点；

步骤二、获得滤波器输出信号：采用所得到的采样序列x(n)，根据误差信号e(n)，应用比例仿射投影算法，在每一离散时间点上进行迭代，对未知的回声路径h进行估计和更新调整，分别计算出两个自适应滤波器的输出

和

步骤三、获得组合滤波器的输出信号：利用组合系数λ(n)，将两个输出信号

和

进行组合，得到

且λ(n)∈[0,1]；

步骤四、组合系数的更新：采用EG法计算新的组合系数；这样利用组合参数和

得到组合自适应滤波器的输出；当自适应滤波器收敛后，该输出信号即为回声信号的一个复制，将其从期望信号d(n)中减去便可消除回声。

步骤二中应用比例仿射投影算法是在每一离散时间点上进行迭代时，其更新过程为：

$e\left(n\right)=d\left(n\right)-X^T\left(n\right)\hat{h}(n-1)---\left(1\right)$

$\hat{h}\left(n\right)=\hat{h}(n-1)+\mathrm{\μG}(n-1)X\left(n\right){[{\mathrm{\δI}}_K+X\left(n\right)G(n-1)X\left(n\right)]}^{-1}e\left(n\right)---\left(2\right)$

其中，

是对未知的回声路径h进行估计和更新调整，0＜μ＜1是每个独立自适应滤波器的步长，δ是一个比较小的正数，e(n)是长度为K的误差信号向量，K表示投影阶数，d(n)＝[d(n),d(n-1),…,d(n-K+1)]^T是包含最近的K个样本值的目标信号向量，d(n)是期望信号，L是滤波器长度，X(n)＝[x(n),x(n-1),…,x(n-K+1)]^T是输入信号矩阵,x(n)＝[x(n),x(n-1),…,x(n-L+1)]^T，G(n-1)是一个L×L维对角矩阵，用于给滤波器每个权分配独一的步长，大权值获得大的步长以加快该系数的收敛速度，它由g_l(n-1)（0≤l≤L-1）确定：

$g_l(n-1)=\frac{1-\mathrm{\κ}}{2L}+(1+\mathrm{\κ})\frac{\left|{\hat{h}}_l(n-1)\right|}{2\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\left|{\hat{h}}_i(n-1)\right|+\mathrm{\ϵ}}---\left(3\right)$

其中，-1＜κ＜1是与系统稀疏度有关的参数，ε是一个比较小的正数，以防止分母为0，I_K是K×K的单位矩阵。

步骤四中进行凸组合自适应滤波方法的组合系数更新时，更新过程定义为：

$\begin{array}{c}{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\arg \min \left\{d\right[w,w_{\mathrm{\ρ}}\left(n\right)]+\frac{{\mathrm{\μ}}_{\mathrm{\ρ}}}{2}\×\\ {[e^2\left(n\right)+\frac{\∂[{d\left(n\right)-w^{T}y\left(n\right)]}^2}{\∂\mathrm{\λ}}|}_{\mathrm{\λ}={\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\\ \×(\mathrm{\λ}-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right))\left]\right\}\end{array}---\left(4\right)$

式中，w＝[λ,1-λ]^T，w_ρ(t)＝[λ_ρ(n),1-λ_ρ(n)]^T， $e\left(n\right)=d\left(n\right)-w_{\mathrm{\ρ}}^T\left(n\right)y\left(n\right),y\left(n\right)={[{\hat{d}}_1\left(n\right),{\hat{d}}_2\left(n\right)]}^T,$ $d[w,w_{\mathrm{\ρ}}\left(n\right)]=\mathrm{\λ}\ln \left(\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)+(1-\mathrm{\λ})\ln \left(\frac{1-\mathrm{\λ}}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)$ 是新的估计量λ和旧的估计量λ_ρ(n)之间的相对熵距离，μ_ρ＞0是组合参数进行更新时的步长因子，决定着组合参数迭代前后的变化大小。

组合系数更新过程等价于求满足下式的λ值：

$\ln (\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\·\frac{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{\mathrm{\λ}})+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)[{\hat{d}}_2\left(n\right)-{\hat{d}}_1\left(n\right)]=0---\left(5\right)$

所以，

${\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\mathrm{\λ}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(6\right)$

即

$1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\frac{[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(7\right)$

由此得到：

$\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\right]---\left(8\right)$

从而

$\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))---\left(9\right)$

这是组合系数进行更新的等价公式。

所述的组合系数更新过程，其特征在于通过定义即 ${\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)=\frac{1}{1+\exp [-\mathrm{\ρ}\left(n\right)]},$ 得到：

$\begin{array}{c}\mathrm{\ρ}(n+1)=\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\\ =\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}\{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)e_{a,1}\left(n\right)+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]e_{a,2}\left(n\right)+n\left(n\right)\}\\ \×[e_{a,1}\left(n\right)-e_{a,2}\left(n\right)]\end{array}---\left(10\right)$

这就是组合系数完整的更新迭代过程。

优点及效果：本发明与现有技术相比，具有如下优点：

（1）采用凸组合的方式，将两个不同参数的IPAPA滤波器进行结合，提高了本发明对系统不同稀疏特性的适应性，减少了参数选择的局限性，也使本发明在收敛速度和估计精度方面达到了更好的折中。

（2）在步骤四中，采用EG准则进行组合参数的更新，该更新过程与现有技术相比，不需要再乘以[λ_ρ(n)(1-λ_ρ(n))]。这样不仅减少了自适应回声消除器的计算量，也避免了当λ_ρ(n)接近0或1时，因子[λ_ρ(n)(1-λ_ρ(n))]会导致组合系数的更新减慢甚至停止的问题。

附图说明：

图1是本发明的自适应回声消除器的原理框图；

图2是本发明的应用流程图；

图3是本发明实施例中用到的回声路径，(a)稀疏度为0.69的实测室内信道(b)稀疏度为0.49的仿真信道；

图4是两个不同步长的IPAPA方法及其组合方法的性能比较图，(a)输入信号为高斯白噪声(b)输入信号为USASI信号；

图5是在不同稀疏度条件下，IPAPA方法及其组合方法的性能比较图。

具体实施方式：

本发明的技术解决方案是：利用凸组合思想，将两个IPAPA滤波器通过组合系数形成一个整体，更新过程中每个滤波器均性能不同且独立地对未知系统进行建模，而权系数、误差函数通过组合系数进行统一调整。同时，采用在稀疏自适应滤波理论中非常有效的EG法代替传统凸组合自适应滤波方法中的梯度下降法，得到新的组合参数迭代公式。具体地，本发明包括如下步骤：

步骤一、获得一个采样序列x(n)，其中n＝1,2,…,N，N为采样序列x(n)的采样点数量；所述的采样序列x(n)是一维信号，且其中包含N个采样点；

步骤二、应用改进的比例仿射投影算法（IPAPA，Improved ProportionateAffine Projection Algorithm）方法在每一离散时间点上进行迭代，根据误差信号e(n)且按照公式 $\hat{h}\left(n\right)=\hat{h}(n-1)+\mathrm{\μG}(n-1)X\left(n\right){[{\mathrm{\δI}}_K+X\left(n\right)G(n-1)X\left(n\right)]}^{-1}e\left(n\right),$ 对未知的回声路径h进行估计和更新调整，其中δ是一个比较小的正数，0＜μ＜1是每个独立自适应滤波器的步长，这样分别计算出两个自适应滤波器的输出

和

e(n)是长度为K的误差信号向量（K表示投影阶数），

d(n)＝[d(n),d(n-1),…,d(n-K+1)]^T是包含最近的K个样本值的目标信号向量，d(n)是期望信号，X(n)＝[x(n),x(n-1),…,x(n-K+1)]^T是输入信号矩阵，其中x(n)＝[x(n),x(n-1),…,x(n-L+1)]^T，L是滤波器长度；G(n-1)是一个L×L维对角矩阵，用于给滤波器每个权分配独一的步长，大权值获得大的步长以加快该系数的收敛速度，它由g_l(n-1)（0≤l≤L-1）确定：

-1＜κ＜1是与系统稀疏度有关的参数，ε是一个比较小的正数，以防止分母为0，I_K是K×K的单位矩阵；

步骤三、利用组合系数λ(n)，（λ(n)∈[0,1]）将两个输出信号

和按照公式进行组合得到 $\hat{d}\left(n\right)=\mathrm{\λ}\left(n\right){\hat{d}}_1\left(n\right)+[1-\mathrm{\λ}\left(n\right)]{\hat{d}}_2\left(n\right);$

步骤四、按照公式

进行组合系数的更新，其中

$\begin{array}{c}\mathrm{\ρ}(n+1)=\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\\ =\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}\{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)e_{a,1}\left(n\right)+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]e_{a,2}\left(n\right)+n\left(n\right)\}\\ \×[e_{a,1}\left(n\right)-e_{a,2}\left(n\right)]\end{array}.$

下面通过附图对本发明加以具体说明：

参照图1本发明实施例中，自适应滤波器

用于建模未知的网络回声路径h，其中L是滤波器长度，n是时间系数，e(n)是误差信号，

是对h的估计。当远端语音信号x(n)通过h时，会产生回声信号y(n)，其叠加到与近端信号s(n)（一般受到加性噪声v(n)的污染）一起作为自适应滤波器的期望信号d(n)。当自适应滤波器的收敛后，该输出信号

即为回声信号的一个复制，将其从期望信号d(n)中减去便可消除回声。

参照图2，本实施例的具体实现步骤如下：

步骤一、获得一个采样序列x(n)，其中n＝1,2,…,N，N为采样序列x(n)的采样点数量；所述的采样序列x(n)是一维信号，且其中包含N个采样点。

步骤二、应用IPAPA方法在每一离散时间点上进行迭代，其更新过程为：

$e\left(n\right)=d\left(n\right)-X^T\left(n\right)\hat{h}(n-1)---\left(1\right)$

$\hat{h}\left(n\right)=\hat{h}(n-1)+\mathrm{\μG}(n-1)X\left(n\right){[{\mathrm{\δI}}_K+X\left(n\right)G(n-1)X\left(n\right)]}^{-1}e\left(n\right)---\left(2\right)$

其中,

是对未知的回声路径h进行估计和更新调整，0＜μ＜1是每个独立自适应滤波器的步长，因子δ是一个比较小的正数，e(n)是长度为K的误差信号向量（K表示投影阶数）；d(n)＝[d(n),d(n-1),…,d(n-K+1)]^T是包含最近的K个样本值的目标信号向量，d(n)是期望信号；L是滤波器长度；X(n)＝[x(n),x(n-1),…,x(n-K+1)]^T是输入信号矩阵,x(n)＝[x(n),x(n-1),…,x(n-L+1)]^T，G(n-1)是一个L×L维对角矩阵，用于给滤波器每个权分配独一的步长，大权值获得大的步长以加快该系数的收敛速度。它由g_l(n-1)（0≤l≤L-1）确定：

$g_l(n-1)=\frac{1-\mathrm{\κ}}{2L}+(1+\mathrm{\κ})\frac{\left|{\hat{h}}_l(n-1)\right|}{2\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\left|{\hat{h}}_i(n-1)\right|+\mathrm{\ϵ}}---\left(3\right)$

式中，-1＜κ＜1是与系统稀疏度有关的参数，ε是一个比较小的正数，μ_ρ＞0以防止分母为0，I_K是K×K的单位矩阵。

步骤三、利用组合系数λ(n)，（λ(n)∈[0,1]）将两个输出信号

和

进行按照公式进行组合得到 $\hat{d}\left(n\right)=\mathrm{\λ}\left(n\right){\hat{d}}_1\left(n\right)+[1-\mathrm{\λ}\left(n\right)]{\hat{d}}_2\left(n\right).$

步骤四、进行组合系数的更新：按照统计信号处理理论，几乎所有已知自适应滤波方法的代价函数都能够写为式(4)的形式：

G(w)＝d(w,w_ρ)+μ_ρL(d,w·x)     (4)

其中，w_ρ表示旧的向量估计，w是当前向量的估计，x表示输入向量，μ_ρ是组合参数进行更新时的步长因子，决定着组合参数迭代前后的变化大小，d是目标输出。等式(4)右侧的第一项函数d(·)是新、旧权向量估计的距离函数，用于表征代价函数的稳定性，第二项L(·)是误差的损失函数，表征代价函数的矫正性。对G(w)进行最小化，即求G(w)对w的导数，并令导数为零：

$\frac{\∂d(w,w_{\mathrm{\ρ}})}{\∂w}+\mathrm{\μ}\frac{\∂L(d,w\·x)}{\∂w}x=0---\left(5\right)$

可以用泰勒公式对损失函数进行表达：

L(d,w·x)＝L(d,w_ρ·x)+L'(d,w_ρ·x)·x·(w-w_ρ)     (6)

在基于梯度下降的算法中（如NLMS类算法），等式(4)中的距离函数d(·)是欧氏距离，损失函数L(·)是平方损失函数。因此，对应的凸组合自适应滤波方法中组合系数的更新也是基于欧氏距离和平方损失函数的。基于指数梯度的自适应滤波方法是一种十分适合回声消除系统的方法，在该类方法中距离函数d(·)不是欧氏距离而是相对熵距离。因此，可以将凸组合自适应滤波方法中组合系数按照等式(4)、(5)、(6)将更新过程定义为：

$\begin{array}{c}{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\arg \min \left\{d\right[w,w_{\mathrm{\ρ}}\left(n\right)]+\frac{{\mathrm{\μ}}_{\mathrm{\ρ}}}{2}\×\\ {[e^2\left(n\right)+\frac{\∂[{d\left(n\right)-w^{T}y\left(n\right)]}^2}{\∂\mathrm{\λ}}|}_{\mathrm{\λ}={\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\\ \×(\mathrm{\λ}-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right))\left]\right\}\end{array}---\left(7\right)$

式中，w＝[λ,1-λ]^T，w_ρ(t)＝[λ_ρ(n),1-λ_ρ(n)]^T， $e\left(n\right)=d\left(n\right)-w_{\mathrm{\ρ}}^T\left(n\right)y\left(n\right),y\left(n\right)={[{\hat{d}}_1\left(n\right),{\hat{d}}_2\left(n\right)]}^T,$ $d[w,w_{\mathrm{\ρ}}\left(n\right)]=\mathrm{\λ}\ln \left(\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)+(1-\mathrm{\λ})\ln \left(\frac{1-\mathrm{\λ}}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)$ 是新的估计量λ和旧的估计量λ_ρ(n)之间的相对熵距离。

式(7)右侧是λ的函数，因此分别对每一项求关于λ的导数。第一项： $\frac{\∂\left[d(w,w_{\mathrm{\ρ}}\left(n\right))\right]}{\∂\mathrm{\λ}}=\ln (\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\·\frac{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{\mathrm{\λ}});$ 第二项中不含有λ，所以， $\frac{\∂\left[\frac{{\mathrm{\μ}}_{\mathrm{\ρ}}}{2}{e}^2\left(n\right)\right]}{\∂\mathrm{\λ}}=0;$ 将第三项展开，得到

$\frac{\∂\left\{{\mathrm{\μ}}_{\mathrm{\ρ}}\right\{d\left(n\right)-[{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right){\hat{d}}_1\left(n\right)+(1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)){\hat{d}}_2\left(n\right)]\}\·[-{\hat{d}}_1\left(n\right)+{\hat{d}}_2\left(n\right)]\·(\mathrm{\λ}-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right))\}}{\∂\mathrm{\λ}}={\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)[{\hat{d}}_2\left(n\right)-{\hat{d}}_1\left(n\right)].$ 这里用到了e(n)的定义式。因此，由式(5)求λ的最小值，就等价于求满足下式的λ值：

$\ln (\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\·\frac{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{\mathrm{\λ}})+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)[{\hat{d}}_2\left(n\right)-{\hat{d}}_1\left(n\right)]=0---\left(8\right)$

所以，

${\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\mathrm{\λ}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(9\right)$

为将本方法与基于梯度下降得到的组合系数更新方法进行比较，这里需要对上述公式做一定的变形和定义。

首先，由式(9)得到：

$1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\frac{[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(10\right)$

式(9)和式(10)相除，得到：

$\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\right]---\left(11\right)$

对式(11)的左右两侧取对数，得到：

$\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))---\left(12\right)$

定义 $\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}=\mathrm{\ρ}\left(n\right),$ 即

${\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)=\frac{1}{1+\exp [-\mathrm{\ρ}\left(n\right)]}---\left(13\right)$

进一步，式(12)写为：

$\begin{array}{c}\mathrm{\ρ}(n+1)=\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\\ =\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}\{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)e_{a,1}\left(n\right)+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]e_{a,2}\left(n\right)+n\left(n\right)\}\\ \×[e_{a,1}\left(n\right)-e_{a,2}\left(n\right)]\end{array}---\left(14\right)$

这样分别计算出两个自适应滤波器的输出

和

利用式(13)和式(14)所示的组合参数得到组合自适应滤波器的输出。当自适应滤波器收敛后，该输出信号即为回声信号的一个复制，将其从期望信号d(n)中减去便可消除回声。

本实施例中，步骤一中的信号远端输入信号分别选用的是零均值、方差为1的WGN（White Gausssian Noise）信号和一段与人类语音信号的频谱类似的USASI(USA Standards Institute)信号。噪声为外加的信噪比为25dB的WGN（White Gausssian Noise），近端信号s(t)＝0。因此，目标信号分别为两种不同的远端输入信号叠加WGN。

本实施例中，回声路径采用两种形式。一种是实测的大小为4m×3m，墙壁为硬墙的室内环境的回声路径，即选择激励信号为幅频特性为平坦谱的白噪声信号，记录由扬声器输出的激励信号在麦克风接收端的响应，然后通过解卷积的方法求得室内冲激响应。在此过程中，白噪声的采样频率为8KHz，精度是16位；麦克风放在1m左右高的桌子上，且与扬声器的距离为1m左右。由于测量结果会随时间而不同，因此为获得较为准确的冲激响应，这里我们记录了15次，

然后取平均，每次记录1秒（8000个权值）的数据。为简化自适应滤波方法的测试，这里取160个权，即20ms的系统冲激响应，详见图3（a）。另一种是根据G.168标准产生的非稀疏的信道，详见图3（b）。图3中，横轴表示单位是ms的时间变量，纵轴是幅度变量。

实际进行取值时，ε和δ是一个比较小的正数；μ_ρ＞0是组合参数进行更新时的步长因子，决定着组合参数迭代前后的变化大小，它一般与信噪比、两个独立自适应滤波器的步长因子等因素有关，可根据具体需要进行相应调整；参数κ是与系统稀疏度有关的参数，是对APA（κ＝-1）和IPAPA（κ＝1）的折中，所以范围是-1＜κ＜1，可根据具体需要进行选取；0＜μ＜1是每个独立自适应滤波器的步长因子，与所有基于梯度下降的自适应滤波方法一样，IPAPA需要权衡收敛速度与估计精度，大的步长μ具有较快的收敛速度，但却会带来较大的失调量。在本实施例，步骤二中的自适应滤波器长度均设为160，ε＝10^-6，δ＝0.001，组合系数的步长μ_a＝100，每个实验均是50次的平均结果。考虑如下两种情况：（1）当选择κ₁＝κ₂＝-0.5，μ₁＞μ₂即μ₁＝0.9,μ₂＝0.1时，组合方法可以实现由μ₁带来的较快的收敛速度和由μ₂带来的较好的稳态性能，从而实现收敛速度与稳态失调的折中，具体结果见图4。（2）当选择κ₁＜0即κ₁＝-0.5，κ₂≈1即κ₂＝0.9，μ₁＝μ₂＝0.1时，组合方法可以同时实现对稀疏度较小(由κ₁保证)和稀疏度较大(由κ₂保证)的信道的良好收敛能力，即提高了本发明对信道不同稀疏度的适应能力，具体结果见图5。图4,5中，横轴表示迭代次数，纵轴表示单位为dB的失调量。

Claims (5)

1.一种稀疏路径的回声消除方法，其特征在于：该方法包括以下步骤：

步骤一、获得一个采样序列x(n)，其中n＝1,2,…,N，N为采样序列x(n)的采样点数量；所述的采样序列x(n)是一维信号，且其中包含N个采样点；

步骤二、获得滤波器输出信号：采用所得到的的采样序列x(n)，根据误差信号e(n)，应用比例仿射投影算法，在每一离散时间点上进行迭代，对未知的回声路径h进行估计和更新调整，分别计算出两个自适应滤波器的输出

和

步骤三、获得组合滤波器的输出信号：利用组合系数λ(n)，将两个输出信号

和

进行组合，得到 $\hat{d}\left(n\right)=\mathrm{\λ}\left(n\right){\hat{d}}_1\left(n\right)+[1-\mathrm{\λ}\left(n\right)]{\hat{d}}_2\left(n\right),$ 且λ(n)∈[0,1]；

步骤四、组合系数的更新：采用自适应滤波法计算新的组合系数；这样利用组合参数和得到组合自适应滤波器的输出；当自适应滤波器的收敛后，该输出信号即为回声信号的一个复制，将其从期望信号d(n)中减去便可消除回声。

2.根据权利要求1所述的稀疏路径的回声消除方法，其特征在于：步骤二中应用比例仿射投影算法是在每一离散时间点上进行迭代时，其更新过程为：

$e\left(n\right)=d\left(n\right)-X^T\left(n\right)\hat{h}(n-1)---\left(1\right)$

$\hat{h}\left(n\right)=\hat{h}(n-1)+\mathrm{\μG}(n-1)X\left(n\right){[{\mathrm{\δI}}_K+X\left(n\right)G(n-1)X\left(n\right)]}^{-1}e\left(n\right)---\left(2\right)$

其中

是对未知的回声路径h进行估计和更新调整，0＜μ＜1是每个独立自适应滤波器的步长，δ是一个比较小的正数，e(n)是长度为K的误差信号向量，K表示投影阶数，d(n)＝[d(n),d(n-1),…,d(n-K+1)]^T是包含最近的K个样本值的目标信号向量，d(n)是期望信号，L是滤波器长度，X(n)＝[x(n),x(n-1),…,x(n-K+1)]^T是输入信号矩阵,x(n)＝[x(n),x(n-1),…,x(n-L+1)]^T，G(n-1)是一个L×L维对角矩阵，用于给滤波器每个权分配独一的步长，大权值获得大的步长以加快该系数的收敛速度，它由g_l(n-1)（0≤l≤L-1）确定：

$g_l(n-1)=\frac{1-\mathrm{\κ}}{2L}+(1+\mathrm{\κ})\frac{\left|{\hat{h}}_l(n-1)\right|}{2\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\left|{\hat{h}}_i(n-1)\right|+\mathrm{\ϵ}}---\left(3\right)$

其中，-1＜κ＜1是与系统稀疏度有关的参数，ε是一个比较小的正数，以防止分母为0，I_K是K×K的单位矩阵。

3.根据权利要求1所述的稀疏路径的回声消除方法，其特征在于：步骤四中进行凸组合自适应滤波方法的组合系数更新时，更新过程定义为：

$\begin{array}{c}{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\arg \min \left\{d\right[w,w_{\mathrm{\ρ}}\left(n\right)]+\frac{{\mathrm{\μ}}_{\mathrm{\ρ}}}{2}\×\\ {[e^2\left(n\right)+\frac{\∂[{d\left(n\right)-w^{T}y\left(n\right)]}^2}{\∂\mathrm{\λ}}|}_{\mathrm{\λ}={\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\\ \×(\mathrm{\λ}-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right))\left]\right\}\end{array}---\left(4\right)$

式中，w＝[λ,1-λ]^T，w_ρ(t)＝[λ_ρ(n),1-λ_ρ(n)]^T， $e\left(n\right)=d\left(n\right)-w_{\mathrm{\ρ}}^T\left(n\right)y\left(n\right),y\left(n\right)={[{\hat{d}}_1\left(n\right),{\hat{d}}_2\left(n\right)]}^T,$ $d[w,w_{\mathrm{\ρ}}\left(n\right)]=\mathrm{\λ}\ln \left(\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)+(1-\mathrm{\λ})\ln \left(\frac{1-\mathrm{\λ}}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\right)$ 是新的估计量λ和旧的估计量λ_ρ(n)之间的相对熵距离，μ_ρ＞0是组合参数进行更新时的步长因子，决定着组合参数迭代前后的变化大小。

4.根据权利要求1或3所述的稀疏路径的回声消除方法，其特征在于：组合系数更新过程等价于求满足下式的λ值：

$\ln (\frac{\mathrm{\λ}}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\·\frac{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{\mathrm{\λ}})+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)[{\hat{d}}_2\left(n\right)-{\hat{d}}_1\left(n\right)]=0---\left(5\right)$

所以，

${\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\mathrm{\λ}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(6\right)$

即

$1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)=\frac{[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_1\left(n\right)\right]+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right){\hat{d}}_2\left(n\right)\right]}---\left(7\right)$

由此得到：

$\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}\exp \left[{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\right]---\left(8\right)$

从而

$\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}(n+1)}=\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))---\left(9\right)$

这是组合系数进行更新的等价公式。

5.根据权利要求1或3所述的稀疏路径的回声消除方法，其特征在于：所述的组合系数更新过程，其特征在于通过定义 $\ln \frac{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}{1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)}=\mathrm{\ρ}\left(n\right),$ 即 ${\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)=\frac{1}{1+\exp [-\mathrm{\ρ}\left(n\right)]},$ 得到：

$\begin{array}{c}\mathrm{\ρ}(n+1)=\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}e\left(n\right)({\hat{d}}_1\left(n\right)-{\hat{d}}_2\left(n\right))\\ =\mathrm{\ρ}\left(n\right)+{\mathrm{\μ}}_{\mathrm{\ρ}}\{{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)e_{a,1}\left(n\right)+[1-{\mathrm{\λ}}_{\mathrm{\ρ}}\left(n\right)]e_{a,2}\left(n\right)+n\left(n\right)\}\\ \×[e_{a,1}\left(n\right)-e_{a,2}\left(n\right)]\end{array}---\left(10\right)$

这就是组合系数完整的更新迭代过程。

Images