CN101568058B - 一种基于加权次梯度投影的数字助听器回声路径估计方法 - Google Patents

Abstract

本发明提出了一种用于数字助听器回声路径估计的加权次梯度投影自适应系统估计方法。本发明通过推导得出加权次梯度投影公式，并将其用于数字助听器自适应回声消除系统之中。其具体方法是利用数字助听器连续输出信号序列和输入信号序列，通过加权次梯度投影算法使自适应滤波器系数h逼近实际回声路径h’。当算法的权值矩阵设置合理的情况下，算法的收敛速度和精度相比较以前的算法有很大提高。并且，无论在白噪声环境还是有色噪声环境下，算法都取得了卓越的失准性能。

Description

一种基于加权次梯度投影的数字助听器回声路径估计方法

一.技术领域

本发明涉及一种数字助听器回声路径自适应估计方法。 

二.背景技术

回声是困扰数字助听器使用的一个普遍问题。通常回声消除的主要思想是先估计回声路径，然后根据回声路径估计回声信号，并将其从输入信号中减去。自适应滤波器法是最常见的回声消除方法。这种方案的核心是设计自适应滤波器估计回声路径，进而估计出回声信号，并将其从输入信号中减去。因为回声路径随环境而变化，故自适应滤波器应能快速跟踪这种变化。 

数字助听器回声路径的自适应估计算法有很多种。在本发明之前，最常用的是是LMS(Least Mean Square，最小均方误差)算法、NLMS(Normal Least Square，归一化最小均方误差)算法和RLS(Recursive Least Square，迭代最小二乘)算法。2000年J.M.Kates在U.S.patent US6072884(Feedback cancellation apparatus and methods)中介绍了LMS算法和NLMS算法在数字助听器回声消除系统的应用，这类算法通过最小化均方误差目标函数的方法估计回声路径，主要优点是计算量较小，有利于实时实现，但是算法的收敛速度和收敛精度都不高，并且在低信噪比和高有色噪声的环境下性能明显下降。S.Haykin在Adaptive Filter Theory.(Upper Saddle River，NJ：Prentice-Hall，1996)一书中介绍了RLS自适应估计算法，该算法虽然能够取得较高的收敛性能，但计算复杂性随滤波器阶数快速上升，难以在数字助听器中实时实现。近年来AP(Affine Projection，仿射投影)算法和POCS(Projection of Onto Convex Sets，凸集投影)算法等投影算法也被应用于通信系统回声路径估计，在这类应用中，AP算法改善了经典NLMS算法在有色信号输入时的收敛性能，而POCS算法相比AP算法的优点在于POCS利用投影算子向特定的闭合凸集多次投影，在低信噪比情况下取得了比AP算法更快的收敛速度。但是，这类投影算法的缺点在于它们所使用的投影算子常常难以定义，而且也没有在数字助听器回声消除系统中应用的研究。 

三.发明内容

1.所要解决的技术问题 

本发明的目的就在于克服上述缺陷，针对数字助听器回声消除模型，设计实现了一种新的加权次梯度投影回声路径估计方法。 

2.技术方案 

一种数字助听器回声路径自适应估计方法，其主要技术步骤如下： 

(1)初始化自适应滤波器h₀，该滤波器系数在算法中将迭代更新，收敛于实际反馈路径h^*； 

(2)构造数字助听器输出信号矩阵U_k＝[u_k，u_k-1，..，u_k-r+1]^T∈R^N×r，这里k表示离散信号时间序号，r为算法的步长，u_k＝[u(k)，u(k-1)，...，u(k-N+1)]^T为助听器输出语音信号向量，N为自适应滤波器长度； 

(3)获取数字助听器输入信号向量 $d_k:=U_k^T{h}^{*}+n_k={[d\left(k\right),d(k-1),...,d(k-r+1)]}^T,$ 其中h^*表示实际反馈路径，n_k表示助听器中原始输入信号s_k和叠加噪声的总和； 

(4)计算凸函数g(h_k)：＝‖U^Th_k-d‖²-ρ，其中参数ρ使用公式ρ：＝m_ξ＝rσ²计算，这里σ为n_k的方差； 

(5)计算次梯度 $t=\▿g\left(h_k\right)={2U}_k(U_k^T{h}_k-d_k);$

(6)根据已有的数字助听器系统先验信息设定指数衰减权重矩阵 

这里A和B满秩，N为自适应滤波器阶数，α_i＝α₀γ^i-1，i＝1，......，N，γ为指数衰减因子，满足0＜γ＜1，γ的值通过先验知识获得或经过测量得到； 

(7)计算加权次梯度投影 

$P_{H^{-}\left(h_k\right)}\left(h_k\right)=\left\{\begin{array}{cc}{h}_k,& g\left(h_k\right)\≤0\\ h_k-\frac{\mathrm{Ag}\left(h_k\right)}{{(\▿g\left(h_k\right))}^{T}A\▿g\left(h_k\right)}\▿g\left(h_k\right),& g\left(h_k\right)>0\end{array}\right.;$ (式1) 

(8)更新滤波器系数 

$h_{k+1}=h_k+{\mathrm{\λ}}_k(P_{H^{-}\left(h_k\right)}\left(h_k\right)-h_k),{\mathrm{\λ}}_k\∈\left[\mathrm{0,2}\right];$ (式2) 

k＝k+1，返回步骤(2)，直至信号序列处理完毕。 

3.有益效果 

本发明的优点和有益效果在于： 

(1)相比较LMS类以均方误差为目标函数的估计算法，本方法收敛速度与收敛精度更高，且在有色噪声环境下性能鲁棒性更强； 

(2)相比较RLS类递归最小二乘算法，本方法计算量更小； 

(3)相比较AP(Affine Projection，仿射投影)算法和POCS(Projection Onto Convex Sets，凸集投影)算法等投影类估计算法，本方法通过次梯度投影技术将投影的可行域扩展为一个半空间并得到投影的一个简单的闭合表达式，不涉及矩阵求逆，且不需要精确的对凸集的正交投影公式，因此实现简单、计算量小； 

(4)相比较一般的自适应次梯度投影(APSM，Adaptive Projection of Subgradient Method)估计算法，本方法由于引入了权值矩阵以利用系统的先验知识，所以算法的收敛速度更快、收敛精度更高。 

本发明的其他优点和效果将在下面继续描述。 

四.附图说明

图1——数字助听器回声消除系统模型框图。 

图2——加权次梯度投影的几何意义示意图。 

图3——h1回声路径脉冲响应图。 

图4——g1回声路径脉冲响应图。 

图5——系统失准系数性能比较图。 

五.具体实施方式

下面结合附图和实施例，对本发明所述的技术方案作进一步的阐述。 

1.数字助听器回声消除系统模型 

本发明使用的数字助听器回声消除系统模型如图1所示。图中G(z)为前向路径信号处理系统，例如响度补偿系统，主要实现声音的放大功能。H^*(z)是待估计的从扬声器到麦克风的外界真实回声路径，其参数为h^*，H(z)是估计回声路径，它的参数h由自适应估计算法产生。Δ是前向路径引入的延时单元，用来对输入信号和反馈信号解相关。s(n)为助听器输入信号，u(n)为助听器输出信号，d(n)是麦克风拾取的全部信号，f(n)是真实反馈信号， 

是估计出的反馈信号，ε(n)是减去估计反馈信号后的残差信号，是前向路径的真正输入。 

2.加权次梯度自适应投影算法 

本方法考虑在图1所示的回声消除系统中，利用连续语音信号序列估计外界回声路径的方法，目的是用自适应滤波器系数h逼近实际回声路径h^*。本方法采用加权次梯度投影算法实现对外界回声路径的估计，其具体实施方式是： 

(1)首先，设定自适应滤波器阶数N，并初始化自适应滤波器系数h₀，该系数可以初始化为常数向量。 

(2)然后，在数字助听器语音信号接收的每个离散时刻k，构造语音信号矩阵U_k＝[u_k，u_k-1，...，u_k-r+1]^T∈R^N×r，其中u_k＝[u(k)，u(k-1)，...，u(k-N+1)]^T为助听器输出信号向量，r为算法的步长。 

(3)获取数字助听器输入端信号向量d_k：＝[d(k)，d(k-1)，..，d(k-r+1)]^T，根据图1所示的模型，该向量 $d_k:=U_k^T{h}^{*}+n_k,$ 其中h^*表示实际反馈路径，n_k表示助听器中原始输入信号s_k 和叠加噪声的总和。 

(4)为了在第k个时刻通过U_k、d_k，使自适应滤波器系数h向实际反馈路径h^*逼近，采取加权次梯度投影的逼近方法，其公式推导如下。 

定义凸集 

$C_k:=\arg {\min }_{h\∈H}\left|\right|U_k^{T}h-d_k\left|\right|,$ (式3) 

{C_k，k＝1，...，m}满足如下约束 

(式4) 

这样对实际反馈路径h^*的估计问题就转变成了向凸集C的投影问题。 

定义凸函数 

(式5) 

以及梯度算子 

(式6) 

那么可以定义基于该凸函数的半空间 

(式7) 

当 $\∀x\∈C_k\left(\mathrm{\ρ}\right)$ 有g(x)≤0，并且可以得到 

             (x-y)^Ts+g(y)≤g(x)≤0，                 (式8) 

(式8)意味着x∈H^-(y)，因此 $C_k\left(\mathrm{\ρ}\right)\⋐H^{-}\left(h_k\right),$ 即对C_k(ρ)的投影可以扩展为对半空间H^-(h_k)的投影，此为次梯度投影，其投影公式为 

$P_{H^{-}\left(h_k\right)}\left(h_k\right):=\left\{\begin{array}{cc}{h}_k& h_k{\∈H}^{-}\left(h_k\right)\\ h_k+\frac{-g\left(h_k\right)}{{\left|\right|\▿g_k\left(h_k\right)\left|\right|}^2}\▿g_k\left(h_k\right)& h_k\∉H^{-}\left(h_k\right)\end{array},.\right.$  

参数ρ的取值直接影响了h^*∈C(ρ)的概率。较大的ρ的取值意味着约束集合包含真实反馈路径概率较大，但收敛速度较慢；相反，较小的ρ的取值意味着收敛速度较快，但约束集合包含真实反馈路径的概率较小，可能导致自适应系统收敛到错误的反馈路径。因此ρ的取值应与相应的噪声环境相吻合。假定随机变量ξ＝‖n_k‖²满足χ² p.d.f.，可以使用简单的公式ρ：＝m_ξ＝rσ²计算ρ，这里σ表示n_k的方差。 

由于在助听器回声路径估计时泄漏在耳道内的脉冲响应是遵循指数衰减规律的，因此充分考虑这样的衰减特性可以提高算法的收敛速度和收敛精度。下面介绍利用回声路径的先验信息设定权值矩阵，使自适应滤波器的系数按照权值矩阵的规律收敛到理想路径的加权次梯度投影方法。这里设定指数衰减权重矩阵 

A和B满秩，N为自适应滤波器阶数，α_i＝α₀γ^i-1，i＝1，......，N，γ为指数衰减因子，满足0＜γ＜1，γ的值通过先验知识获得或经过测量得到。 

当考虑投影的权重矩阵A后，定文 $\tilde{h}=B^{-1}h$ 和 $\tilde{U}=\mathrm{BU},$ 则构造凸函数 

定义为 

$g\left(\tilde{h}\right):={\left|\right|{\tilde{U}}^T\tilde{h}-d\left|\right|}^2-\mathrm{\ρ},$ (式10) 

并且注意到 

$g\left(\tilde{h}\right):={\left|\right|{\left(\mathrm{BU}\right)}^T{B}^{-1}h-d\left|\right|}^2-\mathrm{\ρ}$

$={\left|\right|U^{T}h-d\left|\right|}^2-\mathrm{\ρ}=g\left(h\right),$ (式11) 

因此凸集 $C_k\left(\mathrm{\ρ}\right):=\{\tilde{h}\∈H:g\left(\tilde{h}\right)\≤0\}$ 实际上和(式4)定义的凸集完全相同。 

重新定义半空间 

使其能够分离 ${\tilde{h}}_k=B^{-1}{h}_k$ 和凸集C_k(ρ)，根据次梯度投影原理可以得到 

(式12) 

这里 

表示在 点凸函数 

对 

的次梯度，并且 

(式13) 

这里注意到 

$\widetilde{\▿}g\left(\tilde{h}\right)=2\mathrm{BU}(U^T{B}^T{B}^{-1}h-d)=B\▿g\left(h\right)$ (式14) 

和 

$=H^{-}\left(h_k\right)$ (式15) 

根据(式9)从 

向半空间 的投影可以定义为 

$P_{{\tilde{H}}^{-}\left({\tilde{h}}_k\right)}\left({\tilde{h}}_k\right)=\left\{\begin{array}{cc}{\tilde{h}}_k& {\tilde{h}}_k{\∈\tilde{H}}^{-}\left({\tilde{h}}_k\right)\\ {\tilde{h}}_k+\frac{-g\left({\tilde{h}}_k\right)}{{\left|\right|\widetilde{\▿}{g}_k\left({\tilde{h}}_k\right)\left|\right|}^2}\widetilde{\▿}{g}_k\left({\tilde{h}}_k\right)& {\tilde{h}}_k\∉{\tilde{H}}^{-}\left({\tilde{h}}_k\right)\end{array}\right..$ (式16) 

根据定义 $\tilde{h}=B^{-1}h$ 和 $\tilde{U}=\mathrm{BU},$ 我们可以将(式16)写为 

$P_{{\tilde{H}}^{-}\left({\tilde{h}}_k\right)}\left({B^{-1}h}_k\right)=\left\{\begin{array}{cc}{B}^{-1}{h}_k& g\left({B^{-1}h}_k\right)\≤0\\ B^{-1}{h}_k+\frac{-g\left(B^{-1}{h}_k\right)}{{\left|\right|\widetilde{\▿}{g}_k\left({\tilde{h}}_k\right)\left|\right||}^2}\widetilde{\▿}{g}_k\left({\tilde{h}}_k\right)& g\left({B^{-1}h}_k\right)>0\end{array}\right.,$ (式17) 

利用(式11)、(式14)和(式16)，(式17)等价于 

$P_{H^{-}\left(h_k\right)}\left(h_k\right)=\left\{\begin{array}{cc}{h}_k& g\left(h_k\right)\≤0\\ h_k+\frac{-\mathrm{BBg}\left(h_k\right)}{{(\▿g_k\left(h_k\right))}^T\mathrm{BB}\▿g_k\left(h_k\right)}\▿g_k\left(h_k\right)& g\left(h_k\right)>0\end{array}\right.,$ (式18) 

即 

$P_{H^{-}\left(h_k\right)}\left(h_k\right)=\left\{\begin{array}{cc}{h}_k& g\left(h_k\right)\≤0\\ h_k+\frac{-\mathrm{Ag}\left(h_k\right)}{{(\▿g_k\left(h_k\right))}^{T}A\▿g_k\left(h_k\right)}\▿g_k\left(h_k\right)& g\left(h_k\right)>0\end{array}\right..$ (式19) 

(式19)定义了加权次梯度投影公式，表示为WAPSM(Weighted Adaptive ProjectionSubgradient Method，加权次梯度投影算法)。它意味着h_k的更新方向由设定的权矩阵A修正，该投影的几何意义解释(简单起见，设M＝2)由图2所示。 

本方法根据(式19)计算滤波器系数h_k的加权次梯度投影。首先计算凸函数g(h_k)：＝‖U^Th_k-d‖²-ρ，其中参数ρ使用公式ρ：＝m_ξ＝rσ²计算，这里σ为n_k的方差；然后计算次梯度 $t=\▿g\left(h_k\right)={2U}_k(U_k^T{h}_k-d_k);$ 最后根据(式19)计算计算加权次梯投影P_H ^- _(hk)(h_k)。 

(5)更新滤波器系数 

$h_{k+1}=h_k+{\mathrm{\λ}}_k(P_{H^{-}\left(h_k\right)}\left(h_k\right)-h_k),{\mathrm{\λ}}_k\∈\left[\mathrm{0,2}\right];$ (式20) 

(6)k＝k+1，返回步骤(2)，直至信号序列处理完毕。 

3.性能评价分析 

以下的实验比较了本发明提出的WAPSM(Weighted Adaptive Projection SubgradientMethod，加权次梯度投影)算法与NLMS算法和APSM算法的性能差异。实验使用了两种回声路径，第一种回声路径h1(图3)为空间尺寸5m×4m×3m的真实房间作为回声路径，第二种回声路径g1(图4)由下式给出 

$g1\left(n\right)=\left\{\begin{array}{cc}0& 0\≤n\≤10\\ q\left(n\right)e^{-0.05n}& 11\≤n\≤200\end{array}\right.,$ (式21) 

这里q(n)为均值0方差0.5的高斯白噪声，之所以选用这种路径是因为它非常类似于助听器扬声器经耳道和耳罩泄露到麦克风的真实测量回声路径。 

实验考虑两种类型的助听器输入信号，第一种使用高斯白噪声作为输入信号，第二种使用TIMIT数据库中的“an101-mtms-senn”语音段作为输入信号，该语音片段采样率为16KHz。 

实验系统的前向路径采用固定增益G₀，并具有单一延迟单元，其传递函数为G(z)＝G₀z^-1。实验针对G₀＝3，5，10多种取值进行，不考虑自动增益控制(AGC)单元。在实验中NLMS算法取归一化收敛因子μ＝0.5，在h1路径估计时γ＝0.99，在g1路径估计时γ＝0.95123。 

回声路径估计方法的性能由自适应滤波器失准系数(Misalignment)，其定义为 

(式22) 

这里ω和 

为真实反馈路径的参数及估计路径的参数，‖·‖为范数算子。 

比较WAPSM算法和NLMS算法以及APSM算法在固定增益G₀情况下的失准性能。实际结果如图5所示。图中实验结果显示，不论是白噪声作为输入信号还是语音片段作为输入信号，WAPSM算法都取得了比NLMS算法和APSM算法更稳定的收敛性能和更小的失准系数。对h1反馈路径，由于路径的指数衰减并不明显，因此WAPSM算法和APSM算法性能相似，二者在白噪声输入的情况下，都取得了比NLMS算法低15dB的最终失准系数；对语音片段输入，二者比NLMS算法失准系数降低了6dB。对g1反馈路径，由于路径有明显的指数衰减特征，WAPSM算法比APSM和NLMS算法的性能更优越。在白噪声输入的情况下，WAPSM算法的失准系数比APSM算法下降3dB，比NLMS算法下降4dB；在语音片段作为输入的情况下，WAPSM算法的失准系数比APSM下降了1dB，比NLMS算法下降了5dB。另外，WAPSM算法的失准系数下降的收敛速度和收敛稳定性得到显著提升。特别从图6(c)和图6(d)可以看出，在实际语音片段作为助听器输入信号的情况下，实验结果显示WAPSM算法取得了卓越的失准性能，而绝大多数常用算法，如NLMS算法，在实际语音信号作为输入时性能都有显著下降。 

Claims (1)

1.一种数字助听器回声路径自适应估计方法，其特征是使用加权次梯度投影算法估计回声路径，其主要步骤为：

(1)初始化自适应滤波器h_o，该自适应滤波器系数在算法中将迭代更新，收敛于实际反馈路径h^*；

(2)构造数字助听器输出信号矩阵U_k＝[u_k，u_k-1，...，u_k-r+1]^T∈R^N×r，这里k表示离散信号时间序号，r为算法的步长，u_k＝[u(k)，u(k-1)，...，u(k-N+1)]^T为助听器输出语音信号向量，N为自适应滤波器长度；

(3)获取数字助听器输入信号向量 $d_k:=U_k^T{h}^{*}+n_k={[d\left(k\right),d(k-1),...,d(k-r+1)]}^T,$ 其中

h^*表示实际反馈路径，n_k表示数字助听器中原始输入信号s_k和叠加噪声的总和；

(4)计算凸函数g(h_k)：＝||U^Th_k-d||²-ρ，其中参数ρ使用公式ρ：＝m_ζ＝rσ²计算，这里σ为n_k的方差；

(5)计算次梯度 $t=\▿g\left(h_k\right)={2U}_k(U_k^T{h}_k-d_k);$

(6)根据已有的数字助听器系统先验信息设定指数衰减权重矩阵

这里A和B满秩，N为自适应滤波器阶数，

α_i＝α₀γ^i-1，i＝1，......，N，γ为指数衰减因子，满足0＜γ＜1；

(7)计算加权次梯度投影

$P_{H^{-}\left(b_k\right)}\left(h_k\right)=\left\{\begin{array}{cc}{h}_k,& g\left(h_k\right)\≤0\\ h_k-\frac{\mathrm{Ag}\left(h_k\right)}{{(\▿g\left(h_k\right))}^{T}A\▿g\left(h_k\right)}\▿g\left(h_k\right),& g\left(h_k\right)>0\end{array}\right.;$ (式1)

(8)更新滤波器系数

$h_{k+1}=h_k+{\mathrm{\λ}}_k(p_{H^{-}\left(h_k\right)}\left(h_k\right)-h_k),{\mathrm{\λ}}_k\∈\left[\mathrm{0,2}\right];$ (式2)

(9)k＝k+1，返回步骤(2)，直至信号序列处理完毕。

Images