CN106788337A - 稳健的仿射投影符号自适应滤波算法 - Google Patents

Abstract

稳健的仿射投影符号自适应滤波算法，该方法引入步长函数，并将变步长的方法和比例矩阵的思想融合到一起，即变步长的改进比例仿射投影符号算法,本申请将变步长的方法和比例矩阵的思想融合到一起，引入步长函数，提出了一种稳健的仿射投影符号自适应滤波算法——变步长的改进比例仿射投影符号算法（VSS‑IPAPSA，Variable Step‑Size Improved Proportionate Affine Projection Sign Algorithm）。该算法不仅可以缓解收敛速度与稳态失调之间的矛盾，同时也可以增加其对系统的不同稀疏特性和噪声特性的适应性。理论分析和仿真结果验证了其稳健性和有效性。

Description

稳健的仿射投影符号自适应滤波算法

技术领域：本发明涉及一种稳健的仿射投影符号自适应滤波算法。

背景技术：非高斯的冲击噪声在现实世界广泛存在，严重影响了基于l₂范数优化准则的自适应滤波算法的性能。在各类自适应滤波算法中仿射投影符号算法(APSA，AffineProjection Sign Algorithm)结合了仿射投影算法(APA，Affine Projection Algorithm)良好的收敛特性和符号算法对非高斯冲击噪声干扰的抑制能力，因而其在非高斯冲击噪声条件下具有良好的性能。但是该算法的步长选择是固定的，且未考虑系统稀疏特性，因而在参数选择和收敛速度方面有一定的局限性。具体的说，自适应滤波算法在很多领域都有广泛的应用，如声学和网络回声消除、噪声抑制、信道估计等^[1-4]。文献[5]提出的仿射投影算法(APA，Affine Projection Algorithm)及其改进算法是一类重要的自适应滤波算法，该类算法在输入信号相关性较高的情况下仍具有良好的收敛性能。目前大多数算法是基于高斯噪声假设和l₂范数优化准则的。然而，现实世界中广泛存在着各类非高斯的冲击噪声，如图像中的椒盐噪声、电力开关中的脉动噪声等，这些噪声的存在破坏了基于l₂范数优化准则的自适应滤波算法的性能。为了提高算法对非高斯噪声的抑制性能，文献[6-8]相继提出了一些代替基于l₂范数优化准则的自适应滤波算法，如基于l₁范数优化准则的仿射投影符号算法(APSA，Affine Projection Sign Algorithm)，该算法在存在非高斯噪声干扰时具有良好的稳健性，但是该算法的步长选择是固定的，因而在参数选择和收敛速度方面有一定的局限性，因此，固定步长参数的选择势必要在收敛性能、跟踪性能、稳态误差等之间做出妥协。针对上述问题，文献[9]引入了变步长的概念，提出了基于最小均方偏差(MSD，MeanSquare Deviation)准则的变步长仿射投影符号算法(VSS-APSA，Variable Step-SizeAPSA)，该算法通过引入步长函数，采用随机逼近和移动平均的方法得到算法收敛的最优步长，该算法很大程度上减少了收敛速度与稳态失调之间的矛盾。

然而，上述算法均并未考虑系统的稀疏特性^[10](车载电话、电话会议等系统中声环境的温度和压力的变化或者电话持有者的位置变化会导致系统的稀疏性在一个较大的范围内变化)。

上述参考文献为：

[1]HAYKIN S.Adaptive filter theory[M].Upper Saddle River，NJ，USA：Prentice-Hall，2002.

[2]蔡宇，洪缨，原建平，等.语音系统中的子带自适应回声消除技术[J].仪器仪表学报，2013，34(7)：1448-1453.

[3]张红梅，韩万刚.一种新的变步长LMS自适应滤波算法研究及其应用[J].仪器仪表学报，2015，36(8)：1822-1830.

[4]倪锦根，马兰申.分布式子带自适应滤波算法[J].电子学报，2015，43(11)：2225-2231.

[5]OZEKI K，UMEDAT.An adaptive filtering algorithm using an orthogonalprojection to an affine subspace and its properties[J].ElectronicsCommunication，1984，67(A)：19-27.

[6]XIAO L S，WU M，YANG J.A new efficient filtered-x affine projectionsign algorithm for active control of impulsive noise[J].Signal Processing，2016，120(3)：456-461.

[7]NI J G，LI F.Efficient implementation of the affine projection signalgorithm[J]. IEEE Signal Processing Letters，2012，19(1)：24-26.

[8]SHAO TIANGE，ZHENG Y R，BENESTY J.An affine projection signalgorithm robust against impulsive interferences[J].IEEE Signal ProcessingLetters，2010，17(4)：327-330.

[9]SHIN J，YOO J，PARK P.Variable step-size affine projection signalgorithm[J].Electronics Letters，2012,48(9)：483-485.

[10]郭莹，侯明云.基于指数梯度和凸组合的稀疏自适应滤波算法[J].仪器仪表学报，2014，35(4)：894-900.

发明内容：

发明目的：发明提供一种稳健的仿射投影符号自适应滤波算法，其目的是解决以往所存在的问题。

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

一种稳健的仿射投影符号自适应滤波算法，其特征在于：该方法引入步长函数，并将变步长的方法和比例矩阵的思想融合到一起，即变步长的改进比例仿射投影符号算法(VSS-IPAPSA,Variable Step-Size Improved Proportionate Affine Projection SignAlgorithm)。具体的，本申请主要进行了三方面的工作：1)引入步长函数的概念，采用随机逼近的方法推导得最优步长；2)在实际情况分析过程中，采用移动平均法处理期望项，得出准确的可变步长更新公式；3)将比例矩阵的思想引入传统的仿射投影算法中，提高算法在不同稀疏系统中的适用性，同时利用符号算法对非高斯噪声干扰的抑制能力。

该方法中将比例归一化最小均方算法(PNLMS，Proportionate Normalized LeastMean Square)中成比例的方法引入仿射投影符号算法(APSA，Affine Projection SignAlgorithm)算法中，得到比例仿射投影符号算法(PAPSA, Proportionate AffineProjection Sign Algorithm)算法的更新公式：

其中，sgn(·)表示符号运算，μ是步长因子，L是滤波器长度，n是时间系数，表示所估计的自适应滤波器权系数向量，用于辨识未知系统w₀，输入信号向量为x(n)＝[x(n)x(n-1)…x(n-L+1)]^T，未知系统的期望响应信号为d(n)＝x^T(n)w₀+v(n)(v(n)是噪声)，则先验输出误差向量和后验输出误差向量分别为和d(n)＝[d(n)d(n-1)…d(n-M+1)]^T是期望输出信号向量，M为投影阶数，X(n)是APA滤波器结构中包含M维输入向量的输入矩阵，即X(n)＝[x(n)x(n-1)…x(n-M+1)]，L×L对角矩阵G(n)表示成比例矩阵，g_i(n)是其对角矩阵G(n)的对角元素，表示滤波器第i个抽头参数在n时刻的成比例比重，G(n)和g_i(n)的计算方法如下:

G(n)＝diag{g₀(n),…g_L-1(n)} (2)

其中，max(·)表示求最大值，||·||_∞表示求无穷范数，||·||₂表示求2范数，参数δ_p是一个很小的正数，它在所有滤波器抽头权值为零的时候启动更新；ρ和f(n)能够防止远小于最大抽头的权值停止更新。

同理，将IPNLMS算法的思想应用于APSA算法中，得到本申请的IPAPSA算法，其与PAPSA算法的权系数更新公式类似，不同之处在于对角元素上：

式中，||·||₁表示求1范数，θ是与系统稀疏度有关的参数，稀疏度越大表示回声路径越稀疏，ε是一个比较小的正数，以防止分母为0。

(1)最优步长推导

本申请的算法推导是基于最小均方偏差(MSD，Mean Square Deviation)准则的。定义滤波器系数误差向量并将固定步长μ替代为可变步长μ(n)。为了得到最优步长，与传统方法不同，本申请不作噪声向量与误差向量相关性的假设，即不删除噪声项，采用随机逼近的方法推导出步长，具体推导过程如下：

对式两边取2范数平方的期望得：

其中，f(μ(n))是关于步长的一个函数，且v(n)是噪声向量，符号表示“定义为”。

为了使MSD的值达到最小化，本申请的改进算法从n到n+1的迭代过程中通过选择最优步长来使步长函数最小化，即步长函数写为：

当自适应滤波器收敛到最佳状态时，

此时，sgn(e^T(n))≈sgn(v(n))，且在稳定状态时：

由于准确地计算出sgn(e^T(n))v(n)的值比较困难，所以不能直接获得f(μ(n))；因此，本申请根据sgn(e^T(n))v(n)的上限，运用随机逼近方法得到下式：

但是，v(n)的绝对值不是精确的测量值，本申请将v(n)的绝对值近似为其期望值，v(n)有半正态分布的属性，且即将||v(n)||₁改为：

由式(12)得，

其中，M为投影阶数，σ_v为噪声信号v(n)的方差。

运用随机逼近的方法，步长函数的上限重写为：

从n到n+1的迭代过程中，最小化步长函数f(μ(n))大幅度降低了MSD的值；因此，对式(14)关于μ(n)求偏导，得，

令上式求导结果为0，得到的派生步长为：

(2)最优步长分析

通过式(16)推导出了步长最小化步长函数的上限，然而由于式(16)中期望项的存在，要得到准确的步长是非常困难的；需要考虑以下两种情况，首先定义：

情况1：即β(n)＞0时，步长因子按本申请中的迭代公式进行更新；

情况2:当即β(n)≤0时，此时μ(n)停止更新，即μ(n)＝μ(n-1)。

本申请采用移动平均法来处理上式中的期望项，得到如下：

其中，α(0≤α＜1)为平滑因子，本申请平滑因子的值设为k为一个恒定的值(本申请中k取4)；

因此，本申请所提出的VSS-IPAPSA算法的滤波器系数更新公式为：

G(n)＝diag{g₀(n),…g_L-1(n)}

本申请提出的算法总结如表1所示：

表1 VSS-IPAPSA算法总结

由于上述迭代过程中需要估算噪声信号v(n)的方差σ_v，考虑到冲击噪声干扰对噪声方差估算方法的影响，本申请通过使用中值滤波器来提高输出噪声方差对冲击噪声的稳健性，过程如下：

if mod(n,N_w)＝0

end

其中，A_e(n)＝[e(n)e(n-1)…e(n-N_w+1)]，B_e(n)＝[e²(n)e²(n-1)…e²(n-N_w+1)]。表示n时刻估计的误差方差的平方，表示n时刻估计的输入信号方差的平方，表示n时刻估计的噪声方差的平方，表示n时刻所估计的中间向量，模mod(n,N_w)表示整数n和N_w之间的余数，N_w为估计窗的长度，为有限样值校正因子，目的是减少输出噪声方差估计算法在每n次迭代时计算复杂性。

优点效果：

本申请将变步长的方法和比例矩阵的思想融合到一起，引入步长函数，提出了一种稳健的仿射投影符号自适应滤波算法——变步长的改进比例仿射投影符号算法(VSS-IPAPSA，Variable Step-Size Improved Proportionate Affine Projection SignAlgorithm)。该算法不仅可以缓解收敛速度与稳态失调之间的矛盾，同时也可以增加其对系统的不同稀疏特性和噪声特性的适应性。理论分析和仿真结果验证了其稳健性和有效性。

附图说明：

图1为不同特征指数条件下的概率密度函数。

图2(a)、(b)、(c)、(d)分别为α＝0.5,1.0,1.5,2.0时的α-稳定分布；

图3为自适应滤波器系统辨识原理图；

图4为仿真中用的信道，图4中(a)和(b)分别为仿真中所用的稀疏信道和非稀疏信道，(a)稀疏度为0.85349的稀疏信道，(b)稀疏度为0.32883的非稀疏信道；

图5为高斯噪声条件下各类算法在非稀疏系统下性能比较；

图6为非高斯噪声条件下各个算法在非稀疏系统下性能比较；

图7为非高斯噪声条件下各个算法在稀疏系统下性能比较；

图8为非高斯噪声条件下各个算法在不同系统下的跟踪性能。

具体实施方式：

如图1所示，发明提供一种稳健的仿射投影符号自适应滤波算法，

传统算法分析

α-稳定分布

α-稳定分布作为非高斯冲击噪声模型既可以满足随机噪声产生过程的合理假设又计算方便，从而可以当作信号处理中噪声的理想模型来使用。α-稳定分布有很多种定义方式^[11-12]，本申请主要介绍基于特征函数的定义^[13]。

α-稳定分布的概率密度函数并没有统一的闭合形式，但它的特征函数具有统一的表达式，这是表示α-稳定分布最简便的方法。若一个随机变量X服从α-稳定分布，那么其特征函数可以这样描述：

式中，-∞＜δ＜∞，γ＞0，0＜α≤2，-1≤β≤1，四个参数δ，γ，α，β分别代表不同的物理意义，δ表示位置参数，γ表示尺度参数，α表示特征值数，用来表征α-稳定分布概率密度函数拖尾的厚重程度(当α＝2时为高斯分布)，β为对称参数(当β＝0时，稳定分布是关于δ的对称分布，称为对称α-稳定分布(SαS))。

本申请所用的非高斯冲击噪声模型是用α-稳定分布来描述，图1所示为标准SαS分布在不同特征指数α下的概率密度函数曲线图，图2分别给出了不同α值的SαS稳定噪声。

由图1，图2可以看到，当α＝2时，SαS概率密度函数曲线实质上为均值为零且方差为2的高斯分布，由此可见，SαS概率密度函数曲线拥有许多的特征与高斯分布相同，例如曲线光滑、单峰分布、关于中值对称、钟形等特征。而与高斯分布不同之处在于SαS概率密度函数的拖尾要比高斯分布的更厚重，且特征指数的值越小，其拖尾越厚重，即大幅度样本发生的概率越大。因此，SαS分布非常适合来描述那些类似于高斯分布，却具有很强冲击性的非高斯分布。由于非高斯冲击噪声条件下的稳定分布不存在二阶距，因此服从SαS分布的随机变量的方差也不存在，即该分布破坏了基于l₂范数优化准则的自适应滤波算法的性能。

APSA

本申请以系统辨识^[1]为例来分析各类自适应滤波算法的性能，图3为自适应滤波器系统辨识原理图。

图3中，表示所估计的自适应滤波器权系数向量，用于辨识未知系统w₀，N是滤波器长度，n是时间系数，输入信号向量为x(n)＝[x(n)x(n-1)…x(n-N+1)]^T，则

未知系统的期望响应信号为：

d(n)＝x^T(n)w₀+v(n) (21)

设先验输出误差向量和后验输出误差向量分别为：

式(22)和(23)中，d(n)＝[d(n)d(n-1)…d(n-M+1)]^T是期望输出信号向量，M为投影阶数，X(n)是APA滤波器结构中包含M维输入向量的输入矩阵，即X(n)＝[x(n)x(n-1)…x(n-M+1)]。

传统APSA的滤波器系数更新公式可根据式(24)和(25)推导：

式中，||.||₁表示某矢量的l₁范数，τ²是保证滤波器系数在一次迭代中不会发生太大的变化的参数，一般τ²取较小值。利用拉格朗日乘数法，APSA算法的权系数更新公式为：

其中，μ为步长参数，式(26)中用先验误差向量e(n)代替后验误差向量e_p(n)。这是由于e_p(n)在实践中无法获得，而先验误差向量e(n)是对其较好的近似。

由式(26)可见，APSA利用输出误差向量e(n)的符号进行滤波器系数更新，因此，它可有效地抑制非高斯冲激噪声干扰。此外，APSA通过重复利用过去的输入向量提高了算法在输入信号严重相关时的收敛性能。但是，式(26)中的固定步长使得APSA算法不能同时满足高速收敛与低稳态误差的性能需求，VSS-APSA算法^[9]可有效改善这一缺陷。

VSS-APSA

该算法的推导是基于最小MSD准则的。定义滤波器系数误差向量并将固定步长μ替代为可变步长μ(n)。据此，式(26)可改写为

对式(27)两边取2范数平方的期望，并对所得等式关于μ(n)求导并令其为零，得变步长表达式:

其中，式(28)可进一步变换为

该方法首先假定噪声向量与误差向量不相关即其次用μ(n)的时间平均替代统计平均，从而进一步得到最优步长：

其中，α(0≤α＜1)为平滑因子，最小函数min(.)是为保证算法在冲激噪声干扰下的稳定性，即保证步长在迭代过程中一直减小。

由此得出变步长的仿射投影符号算法权系数的更新表达式为：

新算法的提出及分析

上述VSS-APSA算法虽然提高了APSA的收敛性能与稳态误差性能，但其对噪声项进行了删除，导致其步长的计算并非最优的。同时，式(30)中最小函数min(.)的应用使得算法在非稳定环境下(如待辨识系统冲激响应在某时刻从w₀变为-w₀)性能严重衰减。

本申请引入步长函数，并将变步长的方法和比例矩阵的思想融合到一起，提出了一种新的算法，即变步长的改进比例仿射投影符号算法(VSS-IPAPSA)。该算法不仅降低了IPAPSA算法在非稀疏系统下稳态误差，而且在存在冲击噪声环境下稳健性和收敛速度都要优于传统的仿射投影算法，之后的理论分析和仿真结果验证了其有效性和稳健性。

本申请算法中成比例思想的引入

将PNLMS算法^[14]中的成比例的方法引入APSA算法中，可得到PAPSA算法^[15]的更新公式

其中，L×L对角矩阵G(n)表示成比例矩阵，g_i(n)表示滤波器第i个抽头参数在n时刻的成比例比重。其对角矩阵G(n)的对角元素g_i(n)计算方法如下:

G(n)＝diag{g₀(n),…g_L-1(n)}

其中,参数δ_p是一个很小的正数,它在所有滤波器抽头权值为零的时候启动更新；ρ和f(n)能够防止远小于最大抽头的权值停止更新。

同理，将IPNLMS算法^[16]的思想应用于APSA算法中，得到本申请的IPAPSA算法，其与PAPSA算法的权系数更新公式相同，不同之处在于对角元素上：

式中:θ是与系统稀疏度(稀疏度越大表示回声路径越稀疏)有关的参数，ε是一个比较小的正数，以防止分母为0。

3.2本申请算法的最优步长推导及分析

(1)最优步长推导

由于式(29)中噪声项的删除，使得传统的滤波算法中所得步长不是最优的，因此，为了得到最优步长，本申请不作噪声向量与误差向量相关性的假设，即不删除噪声项，采用随机逼近的方法推导出步长，定义滤波器系数误差向量并将固定步长μ替代为可变步长μ(n)。具体推导过程如下：

对式(26)两边取2范数平方的期望得:

其中，f(μ(n))是关于步长的一个函数，且v(n)是噪声向量，符合表示定义为。

为了使MSD的值达到最小化，本申请的改进算法从n到n+1的迭代过程中通过选择最优步长来使步长函数最小化，即步长函数写为：

当自适应滤波器收敛到最佳状态时，

此时，sgn(e^T(n))≈sgn(v(n))，且在稳定状态时：

由于准确地计算出sgn(e^T(n))v(n)的值比较困难，所以不能直接获得f(μ(n))；因此，本申请根据sgn(e^T(n))v(n)的上限，运用随机逼近方法得到：

但是，v(n)的绝对值不是精确的测量值，本申请将v(n)的绝对值近似为其期望值，v(n)有半正态分布的属性，且即将||v(n)||₁改为：

由式(39)得，

其中，M为投影阶数，σ_v为噪声信号v(n)的方差。

运用随机逼近的方法，步长函数的上限重写为：

从n到n+1的迭代过程中，最小化步长函数f(μ(n))大幅度降低了MSD的值；因此，对式(41)关于μ(n)求偏导，得，

令上式求导结果为0，得到的派生步长为：

(2)最优步长分析

由式(43)推导出的是步长最小化步长函数的上限。然而，由于期望项的存在，要得到准确的步长是非常困难的。需要考虑以下两种情况：

情况1：即β(n)＞0时，步长因子按本申请中的迭代公式进行更新；

情况2:当即β(n)≤0时，此时μ(n)停止更新，即μ(n)＝μ(n-1)。本申请采用移动平均法来处理上式中的期望项，得到如下：

其中，α(0≤α＜1)为平滑因子，本申请平滑因子的值设为k为一个恒定的值(本申请中k取4)；

因此，本申请所提出的VSS-IPAPSA算法的滤波器系数更新公式为：

G(n)＝diag{g₀(n),…g_L-1(n)}

3.3算法总结

综上所述，本申请提出的算法总结如表2所示。

表2 VSS-IPAPSA算法总结

由于上述迭代过程中需要估算噪声信号v(n)的方差σ_v，考虑到冲击噪声干扰对噪声方差估算方法的影响，本申请通过使用中值滤波器来提高输出噪声方差对冲击噪声的稳健性，过程如下：

if mod(n,N_w)＝0

end

其中，A_e(n)＝[e(n)e(n-1)…e(n-N_w+1)]，B_e(n)＝[e²(n)e²(n-1)…e²(n-N_w+1)]。表示n时刻估计的误差方差的平方，表示n时刻估计的输入信号方差的平方，表示n时刻估计的噪声方差的平方，表示n时刻所估计的中间向量，模mod(n,N_w)表示整数n和N_w之间的余数，N_w为估计窗的长度，为有限样值校正因子，目的是减少输出噪声方差估计算法在每n次迭代时计算复杂性。

仿真实验

仿真条件

本申请在MATLAB环境下进行仿真实验。输入的有色信号均由零均值高斯白噪声通过一阶AR系统产生，这里假设自适应滤波器的长度和未知系统的长度相等，均为120。实验中变步长算法的初始步长统一设置为μ(0)＝0.1，各个算法的投影阶数均为4，且算法的比较都在公平原则下进行。每个仿真均是20次实验的平均结果。信道h由式(47)产生:

其中，k表示自适应滤波器的权系数序号，本申请中k取1至120。

性能指标及参数设置

实验是使算法在等效步长的条件下，以权误差向量(WEVN，Weight Error Vector)和归一化均方偏差(NMSD，Normalized MSD)的收敛曲线来评价算法的收敛性能。

WEVN(权误差向量)的值愈小意味着自适应滤波器愈逼近未知系统。

NMSD(归一化的均方偏差)用于度量自适应滤波器与目标系统的逼近程度。

仿真中的参数是在等效步长^[17]的条件下设定的，如表3所示。

表3实验中的参数设置

仿真结果与分析

(1)高斯噪声条件下各类算法在非稀疏系统中的性能比较

图5为APA、APSA、VSS-APSA和本申请改进的新算法(VSS-IPAPSA)在高斯噪声条件下的NMSD收敛曲线，各算法参数设置如表3所示。

由图5可以看到，步长参数对APA算法和APSA算法收敛速度有明显的影响，VSS-APSA算法利用变步长的方法使得其在保证收敛速度加快的同时，稳态误差降低。本申请提出的VSS-IPAPSA算法结合了变步长方法和比例矩阵的优势，既具有快的收敛速度又具有小的稳态误差。

(2)非高斯噪声条件下各类算法在非稀疏系统中的性能比较

图6为APA、APSA、VSS-APSA、IPAPSA和本申请改进的新算法(VSS-IPAPSA)在非高斯噪声条件下的NMSD收敛曲线，各算法参数设置如表3所示。

由图6可以看到，当干扰噪声为非高斯噪声时，APA算法失效，APSA算法的收敛速度变慢，VSS-APSA算法的收敛速度要优于APSA算法，IPAPSA算法的稳健性和收敛速度要优于APSA算法和其VSS-APSA算法，但是IPAPSA算法的稳态性能比APSA和VSS-APSA差，而本申请提出的新的算法，即变步长IPAPSA算法(VSS-IPAPSA)，在保证其他性能优于传统算法的同时稳态性能得到了大大的提高。

(3)非高斯噪声条件下各类算法在稀疏系统中的性能比较

图7为APSA、VSS-APSA、IPAPSA和本申请改进的新算法(VSS-IPAPSA)在稀疏系统且存在非高斯噪声条件下的NMSD收敛曲线。

由图7可以看到，在稀疏系统下IPAPSA算法的性能要优于APSA算法和其VSS-APSA算法，本申请提出的VSS-IPAPSA算法结合了变步长算法和比例矩阵对收敛速度的优势，使其收敛速度、稳健性及稳态误差都要优于传统的仿射投影算法。

(4)非高斯噪声条件下各类算法的跟踪性能比较

图8中(a)和(b)为非高斯噪声条件下的APSA、VSS-APSA、IPAPSA和新算法(VSS-IPAPSA)分别在稀疏系统与非稀疏系统下的跟踪性能比较图。当迭代次数n≤3000时，所使用的信道是稀疏度为0.85349的稀疏信道(见图4(a))，当迭代次数n≥3000时，所使用的信道是稀疏度为0.32883的非稀疏信道(见图4(b))。

由图8可以看到，本申请改进的新算法(VSS-IPAPSA)在系统突变时的跟踪性能稳定，且优于其他算法。

(5)不同信噪比下，各类算法的权误差向量

在以上仿真实验中，测量噪声v(n)加到输出信号(X^T(n)w₀)上，其信噪比为(SNR)25dB。SNR的定义如下：

然而，不同的信噪比对自适应算法的影响也是不一样的，表4为各类算法在不同信噪比下的权误差向量对比表，其中，其他实验条件同实验3一致。

表4各类算法在不同信噪比下的权误差向量(WEVN)

由表4可得，通过对比，发现本申请改进的新算法(VSS-IPAPSA)在不同信噪比条件下的稳态性能要优于同一条件下的其他算法。

综上：

本申请提出的基于l₁范数的改进的变步长比例仿射投影符号自适应滤波算法(VSS-IPAPSA)是针对非高斯噪声干扰对自适应滤波算法性能的影响和系统的不同稀疏性所提出的。具体地，本申请主要进行了三方面的工作：1) 引入步长函数的概念，采用随机逼近的方法推导得最优步长；2)在实际情况分析过程中，采用移动平均法处理期望项，得出准确的步长；3)将比例矩阵的思想引入传统的仿射投影算法中，提高算法在不同稀疏系统中的适用性，同时利用符号算法对非高斯噪声干扰的抑制能力。实验结果和理论分析表明，本申请的算法对不同稀疏特性的系统具有更好的适应性，且对非高斯噪声的抑制能力、收敛速度、稳态性能、系统发生突变时的跟踪性能均优于其他的自适应滤波算法。

本申请将IPNLMS中比例系数的思想应用于APSA，得到了改进的比例仿射投影符号算法(IPAPSA，Improved Proportionate APSA)，同时结合变步长的方法，提出了一种变步长的改进比例仿射投影符号算法(VSS-IPAPSA，Variable Step-Size IPAPSA)，本算法一方面引入步长函数，采用随机逼近和移动平均的方法得到算法收敛的最优步长，另一方面结合比例矩阵的优势，使得算法对不同系统稀疏性具有广泛的适用性，仿真结果和理论分析表明本申请提出的算法在非高斯噪声和高斯噪声干扰下均具有良好的收敛性能、稳态性能和跟踪性能。

上述的描述中：

[11]单志明.α-稳定分布参数估计及自适应算法研究[D].哈尔滨：哈尔滨工程大学，2012.

[12]张凯，单甘霖.基于高斯和SCKF的非线性非高斯滤波算法[J].仪器仪表学报，2015,35(11)：2524-2530.

[13]SAMORODNITSKY G，TAQQU M S.Stable non-Gaussian random processes：stochastic models with infinite variance[D].New York：Chapman and Hall，1994.

[14]DUTTWEILER D L.Proportionate normalized least-mean-squaresadaptation in echo cancellers[J].IEEE Transactions on Speech and AudioProcessing，2000，8(5)：508-518.

[15]ALBU FELIX，KWAN H.K.New proportionate affine projection signalgorithms[J].IEEE International Symposium on Circuits and Systems，2013：521-524.

[16]BENESTY J，GAY S L.An improved PNLMS algorithm[C].IEEEInternational Conference On Acoustics，Speech，and Signal Processing，Orlando，FL，USA，2002：1881-1884.

[17]郭莹，高媛.几种稀疏自适应滤波算法的对比研究[J].微型机与应用，2014，33(8)：1-3。

Claims (4)

1.一种稳健的仿射投影符号自适应滤波算法，其特征在于：该方法引入步长函数，并将变步长的方法和比例矩阵的思想融合到一起，即变步长的改进比例仿射投影符号算法(VSS-IPAPSA,Variable Step-Size Improved Proportionate Affine Projection SignAlgorithm)。具体的，本申请主要进行了三方面的工作：1)引入步长函数的概念，采用随机逼近的方法推导得最优步长；2)在实际情况分析过程中，采用移动平均法处理期望项，得出准确的可变步长更新公式；3)将比例矩阵的思想引入传统的仿射投影算法中，提高算法在不同稀疏系统中的适用性，同时利用符号算法对非高斯噪声干扰的抑制能力。

2.根据权利要求1所述的稳健的仿射投影符号自适应滤波算法，其特征在于：该方法中将比例归一化最小均方算法(PNLMS，Proportionate Normalized Least Mean Square)中成比例的方法引入仿射投影符号算法(APSA，Affine Projection Sign Algorithm)算法中，得到比例仿射投影符号算法(PAPSA,Proportionate Affine Projection SignAlgorithm)算法的更新公式：

$\hat{w}(n+1)=\hat{w}\left(n\right)+\frac{\mathrm{\μ}G\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}{\sqrt{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)G^T\left(n\right)G\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}}---\left(1\right)$

其中，sgn(·)表示符号运算，μ是步长因子，L是滤波器长度，n是时间系数，表示所估计的自适应滤波器权系数向量，用于辨识未知系统w₀，输入信号向量为x(n)＝[x(n)x(n-1)…x(n-L+1)]^T，未知系统的期望响应信号为d(n)＝x^T(n)w₀+v(n)(v(n)是噪声)，则先验输出误差向量和后验输出误差向量分别为和d(n)＝[d(n)d(n-1)...d(n-M+1)]^T是期望输出信号向量，M为投影阶数，X(n)是APA滤波器结构中包含M维输入向量的输入矩阵，即X(n)＝[x(n)x(n-1)…x(n-M+1)]，L×L对角矩阵G(n)表示成比例矩阵，g_i(n)是其对角矩阵G(n)的对角元素，表示滤波器第i个抽头参数在n时刻的成比例比重，G(n)和g_i(n)的计算方法如下:

G(n)＝diag{g₀(n),…g_L-1(n)} (2)

$f\left(n\right)=\mathrm{\ρ}max({\mathrm{\δ}}_p,|\left|\hat{w}\right(n)\left|{|}_{\mathrm{\∞}}\right)---\left(5\right)$

其中,max(·)表示求最大值，||·||_∞表示求无穷范数，||·||₂表示求2范数，参数δ_p是一个很小的正数，它在所有滤波器抽头权值为零的时候启动更新；ρ和f(n)能够防止远小于最大抽头的权值停止更新；

同理，将IPNLMS算法的思想应用于APSA算法中，得到本申请的IPAPSA算法，其与PAPSA算法的权系数更新公式相同，不同之处在于对角元素上,

$g_i\left(n\right)=\frac{1-\mathrm{\&theta;}}{2L}+\frac{(1+\mathrm{\&theta;})\left|w_i\right(n\left)\right|}{2\left|\right|w\left(n\right)|{|}_1+\mathrm{\&epsiv;}},-1<\mathrm{\&theta;}<1---\left(6\right)$

式中:||·||₁表示求1范数，θ是与系统稀疏度有关的参数，稀疏度越大表示回声路径越稀疏，ε是一个比较小的正数，以防止分母为0。

3.根据权利要求2所述的稳健的仿射投影符号自适应滤波算法，其特征在于：

(1)最优步长推导

本申请的算法推导是基于最小均方偏差(MSD，Mean Square Deviation)准则的。定义滤波器系数误差向量并将固定步长μ替代为可变步长μ(n)。为了得到最优步长，与传统方法不同，本申请不作噪声向量与误差向量相关性的假设，即不删除噪声项，采用随机逼近的方法推导出步长，具体推导过程如下：

对式两边取2范数平方的期望得：

$\begin{array}{c}E\left|\right|\tilde{w}(n+1)|{|}_2^2=E|\left|\tilde{w}\left(n\right)\right|{|}_2^2-2\mathrm{\&mu;}\left(n\right)E\left(\frac{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)\tilde{w}\left(n\right)}{\sqrt{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}}\right)\\ +{\mathrm{\&mu;}}^2\left(n\right)\stackrel{\mathrm{\&Delta;}}{=}E\left|\right|\tilde{w}\left(n\right)|{|}_2^2+f\left(\mathrm{\&mu;}\right(n\left)\right)\end{array}---\left(7\right)$

其中，f(μ(n))是关于步长的一个函数，且v(n)是噪声向量，符合表示定义为。

为了使MSD的值达到最小化，本申请的改进算法从n到n+1的迭代过程中通过选择最优步长来使步长函数最小化，即步长函数写为：

$\begin{array}{c}f\left(\mathrm{\&mu;}\right(n\left)\right)=-2\mathrm{\&mu;}\left(n\right)E\left(\frac{sign\left(e^T\right(n\left)\right)\left(e\right(n)-v(n\left)\right)}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)+{\mathrm{\&mu;}}^2\left(n\right)=\\ -2\mathrm{\&mu;}\left(n\right)E\left(\frac{\left|\right|e\left(n\right)|{|}_1-sign\left(e^T\right(n\left)\right)v\left(n\right)}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)+{\mathrm{\&mu;}}^2\left(n\right)\end{array}---\left(8\right);$

当自适应滤波器收敛到最佳状态时，

$E\left(\frac{\mathrm{sgn}\left(e^T\right(n\left)\right)v\left(n\right)}{\sqrt{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}}\right)\&ap;E\left(\frac{\left|\right|v\left(n\right)|{|}_1}{\sqrt{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}}\right)---\left(9\right)$

此时，sgn(e^T(n))≈sgn(v(n))，且在稳定状态时：

$E\left(\right|\left|X^T\right(n)\tilde{w}\left(n\right)\left|{|}_1\right)<E\left(\right|\left|v\right(n)\left|{|}_1\right)---\left(10\right)$

由于准确地计算出sgn(e^T(n))v(n)的值比较困难，所以不能直接获得f(μ(n))；因此，本申请根据sgn(e^T(n))v(n)的上限，运用随机逼近方法得到下式：

$\mathrm{sgn}\left(e^T\right(n\left)\right)v\left(n\right)=\mathrm{sgn}\left({\tilde{w}}^T\right(n\left)X\right(n)+v^T(n\left)\right)v\left(n\right)\&le;\left|\right|v\left(n\right)|{|}_1---\left(11\right)$

但是，v(n)的绝对值不是精确的测量值，本申请将v(n)的绝对值近似为其期望值，v(n)有半正态分布的属性，且即将||v(n)||₁改为：

$\begin{array}{c}\left|\right|v\left(n\right)|{|}_1\&ap;E\left(\right|\left|v\right(n)|{|}_1)=E\left({\mathrm{\&Sigma;}}_{n=0}^{M-1}\right|v\left(n\right)|)=M\&CenterDot;E\left(\right|v\left(n\right)|)\\ =M\&CenterDot;{\&Integral;}_0^{M}v\&CenterDot;\frac{1}{{\mathrm{\&sigma;}}_v\sqrt{2\mathrm{\&pi;}}}{e}^{-\frac{v^2}{2{\mathrm{\&sigma;}}_v^2}}dv\end{array}---\left(12\right)$

由式(12)得：

$\left|\right|v\left(n\right)|{|}_1\&ap;E\left(\right|\left|v\right(n)|{|}_1)=M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v---\left(13\right)$

其中，M为投影阶数，σ_v为噪声信号v(n)的方差。

运用随机逼近的方法，步长函数的上限重写为：

$\begin{array}{c}f\left(\mathrm{\&mu;}\right(n\left)\right)=-2\mathrm{\&mu;}\left(n\right)E\left(\frac{\left|\right|e\left(n\right)|{|}_1-|\left|v\left(n\right)\right|{|}_1}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)+{\mathrm{\&mu;}}^2\left(n\right)\&ap;\\ -2\mathrm{\&mu;}\left(n\right)E\left(\frac{\left|\right|e\left(n\right)|{|}_1-M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)+{\mathrm{\&mu;}}^2\left(n\right)\stackrel{\mathrm{\&Delta;}}{=}{f}_u\left(\mathrm{\&mu;}\right(n)\end{array}---\left(14\right)$

从n到n+1的迭代过程中，最小化步长函数f(μ(n))大幅度降低了MSD的值；因此，对式(14)关于μ(n)求偏导，得：

$\frac{\&part;f_u\left(\mathrm{\&mu;}\right(n\left)\right)}{\&part;\mathrm{\&mu;}\left(n\right)}=-2E\left(\frac{\left|\right|e\left(n\right)|{|}_1-M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)+2\mathrm{\&mu;}\left(n\right)---\left(15\right)$

令上式求导结果为0，得到的派生步长为：

$\mathrm{\&mu;}\left(n\right)=E\left(\frac{\left|\right|e\left(n\right)|{|}_1-M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}\right)---\left(16\right)$

(2)最优步长分析

通过式(16)推导出了步长最小化步长函数的上限，然而由于式(16)中期望项的存在，要得到准确的步长是非常困难的；需要考虑以下两种情况，首先定义：

$\mathrm{\&beta;}\left(n\right)=\frac{\left|\right|e\left(n\right)|{|}_1-M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}---\left(17\right)$

情况1：即β(n)＞0时，步长因子按本申请中的迭代公式进行更新；

情况2:当即β(n)≤0时，此时μ(n)停止更新，即μ(n)＝μ(n-1)。

本申请采用移动平均法来处理上式中的期望项，得到如下：

其中，α(0≤α＜1)为平滑因子，本申请平滑因子的值设为k为一个恒定的值(本申请中k取4)；

因此，本申请所提出的VSS-IPAPSA算法的滤波器系数更新公式为：

$\hat{w}(n+1)=\hat{w}\left(n\right)+\frac{\mathrm{\&mu;}\left(n\right)G\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}{\sqrt{\mathrm{sgn}\left(e^T\right(n\left)\right)X^T\left(n\right)G^T\left(n\right)G\left(n\right)X\left(n\right)\mathrm{sgn}\left(e\right(n\left)\right)}}---\left(19\right)$

G(n)＝diag{g₀(n),…g_L-1(n)}

$g_i\left(n\right)=\frac{1-\mathrm{\&theta;}}{2L}+\frac{(1+\mathrm{\&theta;})\left|w_i\right(n\left)\right|}{2\left|\right|w\left(n\right)|{|}_1+\mathrm{\&epsiv;}},-1<\mathrm{\&theta;}<1$

$\mathrm{\&beta;}\left(n\right)=\frac{\left|\right|e\left(n\right)|{|}_1-M\&CenterDot;\sqrt{\frac{2}{\mathrm{\&pi;}}}{\mathrm{\&sigma;}}_v}{\sqrt{sign\left(e^T\right(n\left)\right)X^T\left(n\right)X\left(n\right)sign\left(e\right(n\left)\right)}}.$

4.根据权利要求3所述的稳健的仿射投影符号自适应滤波算法，其特征在于：

本申请提出的算法总结如表1所示：

表1 VSS-IPAPSA算法总结

由于上述迭代过程中需要估算噪声信号v(n)的方差σ_v，考虑到冲击噪声干扰对噪声方差估算方法的影响，本申请通过使用中值滤波器来提高输出噪声方差对冲击噪声的稳健性，过程如下：

if mod(n,N_w)＝0

$\hat{r}\left(n\right)=\mathrm{\&alpha;}\hat{r}(n-1)+C(1-\mathrm{\&alpha;})x\left(n\right)median\&lsqb;A_e\left(n\right)\&rsqb;$

${\widehat{\mathrm{\&sigma;}}}_x^2\left(n\right)=\mathrm{\&alpha;}{\widehat{\mathrm{\&sigma;}}}_x^2(n-1)+(1-\mathrm{\&alpha;}){\mathrm{\&mu;}}^2\left(n\right)$

${\widehat{\mathrm{\&sigma;}}}_e^2\left(n\right)=\mathrm{\&alpha;}{\widehat{\mathrm{\&sigma;}}}_e^2(n-1)+C(1-\mathrm{\&alpha;})median\&lsqb;B_e\left(n\right)\&rsqb;$

${\widehat{\mathrm{\&sigma;}}}_v^2\left(n\right)={\widehat{\mathrm{\&sigma;}}}_e^2\left(n\right)-\frac{{\hat{r}}^T\left(n\right)\hat{r}\left(n\right)}{{\widehat{\mathrm{\&sigma;}}}_x^2\left(n\right)}$

end

其中，A_e(n)＝[e(n)e(n-1)…e(n-N_w+1)]，B_e(n)＝[e²(n)e²(n-1)…e²(n-N_w+1)]。表示n时刻估计的误差方差的平方，表示n时刻估计的输入信号方差的平方，表示n时刻估计的噪声方差的平方，表示n时刻所估计的中间向量，模mod(n,N_w)表示整数n和N_w之间的余数，N_w为估计窗的长度，为有限样值校正因子，目的是减少输出噪声方差估计算法在每n次迭代时计算复杂性。