本发明涉及一种凸约束自适应回声消除方法。
背景技术：
近年来，通信系统(比如免提电话、电视电话会议系统等)的广泛应用使得通信用户对通话质量的要求越来越高。语音通话质量主要受到回声的影响，是用户在通信过程中反复听到自己声音的现象，具有延时短、多路径及时变特性。目前，回声消除仍然是一个热门的研究课题。自适应回声消除方法造价低，是被公认的最具有前景的回声消除技术，也是现在回声消除的主流技术。
所谓自适应回声消除技术，是利用自适应滤波器的一种回声消除方法。对于一些未知的环境来说，自适应滤波器能够在工作过程中逐步学习出所需的统计特性，并以此为依据自动调整滤波器系数，以达到最佳滤波的效果；一旦输入信号的统计特性发生变化，它又可以跟踪这种变化，自动调整滤波器参数，使滤波性能重新达到最佳。通信中的回声信道通常属于高阶稀疏信道，传统自适应滤波器在这种信道下的回声消除收敛速度慢，稳态误差大。针对系统稀疏的特性，文献1(eksioglu,e.m.&tanc,a.k.(2011).rlsalgorithmwithconvexregularization.ieeesignalprocessingletters,18(8),470-473)，提出了一种凸约束递归最小二乘算法，简称为crrls算法。该算法在rls算法的代价函数后添加滤波器权值参数的凸约束项，在权值参数更新方程中产生了一个零吸引子，进而使权值参数在更新过程中趋于0，一定程度上减小了稳态失调。然而，当系统的背景噪声被冲击噪声所干扰时，crrls算法会面临不收敛的问题。
技术实现要素：
本发明的目的是提供一种凸约束自适应回声消除方法，该方法具有抗冲击能力，更低的稳态失调，回声消除效果更好。
本发明实现其发明目的所采用的技术方案是，一种凸约束自适应回声消除方法，其步骤如下：
a、远端信号采样
将当前时刻n和前l-1个时刻的远端采样信号x(n),x(n-1),…,x(n-l 1)，组成当前时刻n的输入信号向量x(n)，x(n)＝[x(n),x(n-1),…,x(n-l 1)]t；其中，上标t表示转置，l为自适应滤波器的抽头长度，取值为32、64或128；
b、回声信号估计
将当前时刻n的输入信号向量x(n)通过滤波器得到当前时刻n的回声信号的估计值y(n)，即y(n)＝wt(n)x(n)，其中，w(n)＝[w1(n),w2(n),...,wi(n),...wl(n)]t为当前时刻n的滤波器抽头权系数向量，其初始值为零，即w(0)＝0，wi(n)为当前时刻n的滤波器的第i个抽头权系数，i＝1,2,3,...,l；
c、回声信号消除
将当前时刻n的近端采样信号d(n)，减去步骤b获得的当前时刻n的回声信号的估计值y(n)，得到当前时刻n的误差信号e(n)，e(n)＝d(n)-y(n)；
d、误差信号阈值的确定
d1、方差纠正量的计算
由当前时刻n的误差信号e(n)；得到当前时刻n的误差信号平方e2(n)，将当前时刻n和前nw-1个时刻的误差信号平方e2(n),e2(n-1),...,e2(n-nw 1)，组成误差信号平方时间窗向量e(n)，e(n)＝[e2(n),e2(n-1),...,e2(n-nw 1)]，其中，nw是时间窗长度，取值为9；
计算出当前时刻n的方差纠正量α(n)，α(n)＝c(1-ζ)med(e(n))，其中，c为方差纠正量比例系数，c＝1.483{1 5/(nw-1)}；ζ为方差纠正量延迟参数，其取值为0.9～0.99；med(·)表示取中位数运算；
d2、估计方差的计算
根据当前时刻n的方差纠正量α(n)和前一时刻n-1的估计方差σ2(n-1)，计算出当前时刻n的估计方差σ2(n)，σ2(n)＝σ2(n-1) α(n)；估计方差σ2(n)的初始值为0，即σ2(0)＝0；
d3、误差信号阈值的确定
计算出当前时刻n的误差信号阈值ξ(n)，
e、误差信号标识值的计算
根据当前时刻n的误差信号e(n)及当前时刻n的误差信号阈值ξ(n)，得出当前时刻n的m估计器误差信号标识值q(n)：
f、权向量更新
f1、增益向量的计算
由下式得到当前时刻n的增益向量k(n)，
式中，λ表示遗忘因子，其取值为0.9～0.99999；p(n-1)为l×l的方阵，是前一时刻n-1的递归矩阵；其初始值p(0)为对角元素均为0.0001的对角矩阵；
f2、递归矩阵的计算
由下式得到当前时刻n的递归矩阵p(n)，
p(n)＝λ-1p(n-1)-λ-1k(n)xt(n)p(n-1)
f3、权系数惩罚项向量的计算
由当前时刻n的递归矩阵p(n)和当前时刻n的权系数向量w(n)，得到权系数惩罚项向量γ(n)，
γ(n)＝γ(1-λ)p(n)sgn(w(n))
式中，γ为惩罚力度因子，其取值为1～2；sgn(·)表示取符号运算；
f4、滤波器权向量更新
计算出下一时刻n 1的权向量w(n 1)，w(n 1)＝w(n) k(n)e(n)-γ(n)；
g、迭代
令n＝n 1，重复步骤a、b、c、d、e、f的操作，直至通话结束。
与现有技术相比，本发明的有益效果是：
一、抗冲击性能强
本发明通过对误差信号进行递归m估计器运算，实现误差信号的阈值实时更新，误差信号标识值q(n)使本发明具有抗冲击特性，即在误差信号幅值大于阈值时，判定系统被冲击噪声干扰，误差信号标识值q(n)取值为0，进而使增益向量k(n)取值为0，滤波器权值更新不受误差调整量k(n)e(n)影响，即滤波器不受冲击噪声的干扰，系统鲁棒性好，抗冲击能力强；当误差信号幅值小于阈值时，判定系统中未出现冲击噪声，滤波器权值受到误差调整量k(n)e(n)的影响，实现对误差的实时追踪，加快了收敛速度，减小了稳态失调。
二、稳态失调小
权系数向量更新公式中的最后一项权系数惩罚项向量γ(n)，是一个使权系数趋近于零的量，抽头权系数向量在更新时减去该惩罚项，有利于权系数迅速向零靠近，加快了收敛速度，减小了稳态失调；提高了回声消除的质量。
下面结合附图和具体实施方式对本发明做进一步的详细说明。
附图说明
图1是文献1的cr-rls算法和本发明方法的仿真实验的归一化稳态失调曲线。
具体实施方式
实施例
本发明的一种具体实施方式是，一种凸约束自适应回声消除方法，其步骤如下：
a、远端信号采样
将当前时刻n和前l-1个时刻的远端采样信号x(n),x(n-1),…,x(n-l 1)，组成当前时刻n的输入信号向量x(n)，x(n)＝[x(n),x(n-1),…,x(n-l 1)]t；其中，上标t表示转置，l为自适应滤波器的抽头长度，取值为32、64或128；
b、回声信号估计
将当前时刻n的输入信号向量x(n)通过滤波器得到当前时刻n的回声信号的估计值y(n)，即y(n)＝wt(n)x(n)，其中，w(n)＝[w1(n),w2(n),...,wi(n),...wl(n)]t为当前时刻n的滤波器抽头权系数向量，其初始值为零，即w(0)＝0，wi(n)为当前时刻n的滤波器的第i个抽头权系数，i＝1,2,3,...,l；
c、回声信号消除
将当前时刻n的近端采样信号d(n)，减去步骤b获得的当前时刻n的回声信号的估计值y(n)，得到当前时刻n的误差信号e(n)，e(n)＝d(n)-y(n)；
d、误差信号阈值的确定
d1、方差纠正量的计算
由当前时刻n的误差信号e(n)；得到当前时刻n的误差信号平方e2(n)，将当前时刻n和前nw-1个时刻的误差信号平方e2(n),e2(n-1),...,e2(n-nw 1)，组成误差信号平方时间窗向量e(n)，e(n)＝[e2(n),e2(n-1),...,e2(n-nw 1)]，其中，nw是时间窗长度，取值为9；
计算出当前时刻n的方差纠正量α(n)，α(n)＝c(1-ζ)med(e(n))，其中，c为方差纠正量比例系数，c＝1.483{1 5/(nw-1)}；ζ为方差纠正量延迟参数，其取值为0.9～0.99；med(·)表示取中位数运算；
d2、估计方差的计算
根据当前时刻n的方差纠正量α(n)和前一时刻n-1的估计方差σ2(n-1)，计算出当前时刻n的估计方差σ2(n)，σ2(n)＝σ2(n-1) α(n)；估计方差σ2(n)的初始值为0，即σ2(0)＝0；
d3、误差信号阈值的确定
计算出当前时刻n的误差信号阈值ξ(n)，
e、误差信号标识值的计算
根据当前时刻n的误差信号e(n)及当前时刻n的误差信号阈值ξ(n)，得出当前时刻n的m估计器误差信号标识值q(n)：
f、权向量更新
f1、增益向量的计算
由下式得到当前时刻n的增益向量k(n)，
式中，λ表示遗忘因子，其取值为0.9～0.99999；p(n-1)为l×l的方阵，是前一时刻n-1的递归矩阵；其初始值p(0)为对角元素均为0.0001的对角矩阵；
f2、递归矩阵的计算
由下式得到当前时刻n的递归矩阵p(n)，
p(n)＝λ-1p(n-1)-λ-1k(n)xt(n)p(n-1)
f3、权系数惩罚项向量的计算
由当前时刻n的递归矩阵p(n)和当前时刻n的权系数向量w(n)，得到权系数惩罚项向量γ(n)，
γ(n)＝γ(1-λ)p(n)sgn(w(n))
式中，γ为惩罚力度因子，其取值为1～2；sgn(·)表示取符号运算；
f4、滤波器权向量更新
计算出下一时刻n 1的权向量w(n 1)，w(n 1)＝w(n) k(n)e(n)-γ(n)；
g、迭代
令n＝n 1，重复步骤a、b、c、d、e、f的操作，直至通话结束。
仿真实验
为了验证本发明方法的有效性，我们进行了仿真实验，并与cr-rls算法做了性能对比。
仿真实验中自适应滤波器抽头长度l为64，远端的输入信号采用一阶自回归(ar(1))信号，在房间为长6.25m，宽3.75m，高2.5m，温度20℃，湿度50％的安静密闭房间内，将接收到的远端信号经扬声器播放后，在房间中用麦克风按采样频率为8000hz，共拾取出1000个时刻点的近端信号d(n)。
实验中各算法的参数具体取值如下表:
图1是文献1和本发明方法的实验的归一化稳态失调曲线。
从图1中可以看出在含有冲击噪声的稀疏系统中，文献1基于误差二阶矩的cr-rls算法在时刻50发生冲击噪声后，系统稳态误差始终高于-5db，系统不收敛，即cr-rls算法在冲击噪声的环境中工作时，其收敛效果很差。本发明在在时刻50发生冲击噪声后，系统在时刻100左右起，系统收敛、稳态误差降低到-40db以下，具有较快的收敛速度和较低的稳态失调，表明其能够在含冲击噪声的环境中正常工作。
1.一种凸约束自适应回声消除方法，其步骤如下：
a、远端信号采样
将当前时刻n和前l-1个时刻的远端采样信号x(n),x(n-1),…,x(n-l 1)，组成当前时刻n的输入信号向量x(n)，x(n)＝[x(n),x(n-1),…,x(n-l 1)]t；其中，上标t表示转置；l为自适应滤波器的抽头长度，取值为32、64或128；
b、回声信号估计
将当前时刻n的输入信号向量x(n)通过滤波器得到当前时刻n的回声信号的估计值y(n)，即y(n)＝wt(n)x(n)，其中，w(n)＝[w1(n),w2(n),...,wi(n),...wl(n)]t为当前时刻n的滤波器抽头权系数向量，其初始值为零，即w(0)＝0，wi(n)为当前时刻n的滤波器的第i个抽头权系数，i＝1,2,3,...,l；
c、回声信号消除
将当前时刻n的近端采样信号d(n)，减去步骤b获得的当前时刻n的回声信号的估计值y(n)，得到当前时刻n的误差信号e(n)，e(n)＝d(n)-y(n)；
d、误差信号阈值的确定
d1、方差纠正量的计算
由当前时刻n的误差信号e(n)；得到当前时刻n的误差信号平方e2(n)，将当前时刻n和前nw-1个时刻的误差信号平方e2(n),e2(n-1),...,e2(n-nw 1)，组成误差信号平方时间窗向量e(n)，e(n)＝[e2(n),e2(n-1),...,e2(n-nw 1)]，其中，nw是时间窗长度，取值为9；
计算出当前时刻n的方差纠正量α(n)，α(n)＝c(1-ζ)med(e(n))，其中，c为方差纠正量比例系数，c＝1.483{1 5/(nw-1)}；ζ为方差纠正量延迟参数，其取值为0.9～0.99；med(·)表示取中位数运算；
d2、估计方差的计算
根据当前时刻n的方差纠正量α(n)和前一时刻n-1的估计方差σ2(n-1)，计算出当前时刻n的估计方差σ2(n)，σ2(n)＝σ2(n-1) α(n)；估计方差σ2(n)的初始值为0，即σ2(0)＝0；
d3、误差信号阈值的确定
计算出当前时刻n的误差信号阈值ξ(n)，
e、误差信号标识值的计算
根据当前时刻n的误差信号e(n)及当前时刻n的误差信号阈值ξ(n)，得出当前时刻n的m估计器误差信号标识值q(n)：
f、权向量更新
f1、增益向量的计算
由下式得到当前时刻n的增益向量k(n)，
式中，λ表示遗忘因子，其取值为0.9～0.99999；p(n-1)为l×l的方阵，是前一时刻n-1的递归矩阵；其初始值p(0)为对角元素均为0.0001的对角矩阵；
f2、递归矩阵的计算
由下式得到当前时刻n的递归矩阵p(n)，
p(n)＝λ-1p(n-1)-λ-1k(n)xt(n)p(n-1)
f3、权系数惩罚项向量的计算
由当前时刻n的递归矩阵p(n)和当前时刻n的权系数向量w(n)，得到权系数惩罚项向量γ(n)，
γ(n)＝γ(1-λ)p(n)sgn(w(n))
式中，γ为惩罚力度因子，其取值为1～2；sgn(·)表示取符号运算；
f4、滤波器权向量更新
计算出下一时刻n 1的权向量w(n 1)，w(n 1)＝w(n) k(n)e(n)-γ(n)；
g、迭代
令n＝n 1，重复步骤a、b、c、d、e、f的操作，直至通话结束。
技术总结