CN104158512A - 一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法 - Google Patents

Abstract

一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法，其步骤主要是：A、获取自适应FIR滤波器的输出信号和期望信号；B、将期望信号向量d(n)与输出信号向量相减得到自适应FIR滤波器的输出误差向量e(n)；C、计算自适应FIR滤波器的权值向量的比例因子列向量g(n)；D、计算下一时刻的自适应FIR滤波器的权值列向量w(n+1)；E、令n＝n+1，重复步骤A、B、C、D的操作，即连续得到稀疏系统在不同时刻的辨识值w(n)。该方法的收敛速度快，跟踪能力强，计算复杂度低。

Description

一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法

技术领域

本发明涉及一种在冲击噪声环境中的稀疏系统辨识方法，属于数字信号处理技术领域。

背景技术

近些年来，随着自适应滤波算法理论的飞速发展，自适应滤波器技术已广泛成功应用在系统辨识，信道均衡，有源噪声控制和回声消除(包括声学回声消除，网络回声消除和水下回声消除)等领域。本质上，这些应用都是一个自适应系统辨识问题，即利用一个自适应滤波器并根据未知系统的输入和输出信号的统计特性，去辨识未知系统的脉冲响应或传递函数。例如：在语音回声消除中，虽然是首先将扬声器的信号送给自适应滤波器产生回声的估计信号，然后与麦克风接收到的信号相减得到干净的信号并将其传送给远端，实现回声消除，但本质上是用自适应滤波器辨识扬声器与麦克风之间的回声路径的脉冲响应。

在实际中，常常假定未知系统具有有限脉冲响应(FIR)的结构，记为w_o＝[w_o1,w_o2,...,w_oM]^T，并构造一个具有相同长度M的自适应FIR滤波器，记为w(n)＝[w₁(n),w₂(n),...,w_i(n),...,w_M(n)]^T，其中w_i(n)为w(n)中的第i个权值，然后采用自适应滤波算法不断的更新w(n)使其逼近w_o，从而实现辨识未知系统的目的。众所周知，当发送给未知系统和自适应滤波器的输入信号是相关信号(或称有色信号)时，仿射投影算法(Affine projection algorithm，简称APA)比归一化最小均方(Normalized least mean square，简称NLMS)算法有着更快的收敛速度，是因为APA具有解相关的能力，并且这种能力随着仿射投影阶数的增加而增强。然而，在实际应用中，例如无线通信、声学回声消除、水下语音信道估计等，未知脉冲响应常常是稀疏的，即它的大部分系数的幅度都接近或等于零，并且这类系数被称为非活性系数而其它具有较大幅度的系数被称为活性系数。在这种情况下，APA和NLMS算法的收敛速度都会出现明显的下降。为了解决这个问题，一系列成比例的仿射投影算法被提出，例如成比例APA(简称PAPA)，改进的PAPA(简称IPAPA)，带记忆的IPAPA(简称MIP-APA)等。这些成比例算法的基本原理是：在每一次迭代过程中，为自适应FIR滤波器的每个系数分配一个独立的步长，并且较大的系数获得较大的步长，较小的系数获得较小的步长，从而加快全局收敛速度。

遗憾的是，当环境中存在冲击干扰(或噪声)时，上述算法的收敛速度被大大的降低，严重时将导致上述算法发散，这是因为上述算法是基于l2范数的最优化被提出的。为了适应冲击干扰环境，T.Shao等基于l1范数的最优化的优点以及结合仿射投影算法的解相关特性，于2010年提出了仿射投影符号算法(T.Shao,Y.R.Zheng,and J.Benesty,“An affine projection signalgorithm robust against impulsive interference,”IEEE Signal Process.Lett.,vol.17,no.4,pp.327–330,2010.)，简称APSA。随后，为了满足稀疏系统的要求，Z.Yang等将PAPA和IPAPA中的“成比例”思想扩展到APSA中(Z.Yang,Y.R.Zheng,and S.L.Grant,“Proportionate affine projection sign algorithms fornetwork echo cancellation,”IEEE Trans.Speech Audio Process.,vol.19,no.8,pp.2273–2284,2011.)，衍生出了成比例APSA(简称RP-APSA)和改进的RP-APSA(简称RIP-APSA)。近来，F.Albu和H.K.Kwan通过借鉴MIP-APA算法的思想，提出了带记忆的RIP-APSA(Memory improved proportionateaffine projection sign algorithm,”Electronics Letters,vol.48,no.20,pp.1279–1281,2012.)，简称MIP-APSA.

然而，RP-APSA的性能依赖于两个预定义的参数，换句话说，RP-APSA要获得比较好的性能，必须合适的选择这两个参数。此外，在RP-APSA中，由于分配给自适应FIR滤波器的每个权值的活性因子都是相同的，因此导致了所有非活性滤波器权值都被分配了相同步长，从而放慢了RP-APSA在整个自适应过程中的收敛速度。

发明内容

本发明的目的是提供一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法，该方法的收敛速度快和跟踪能力强，计算复杂度低。

本发明实现其发明目的所采用的技术方案是，一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法，包括以下步骤：

A、获取滤波器的期望信号和输出信号

发送输入信号x(n)给稀疏系统w_o，得到稀疏系统w_o的输出信号d(n)即自适应FIR滤波器的期望信号；同时发送输入信号x(n)给自适应FIR滤波器得到自适应FIR滤波器的输出信号 其中：

n表示当前时刻，上标T表示转置运算，w(n)＝[w₁(n),w₂(n),...,w_i(n),...,w_M(n)]^T为自适应FIR滤波器在当前时刻的权值列向量，也即稀疏系统w_o在当前时刻的列向量辨识值，其长度为M,；w_i(n)为w(n)中的第i个权值,i＝1,2,…M,各个权值w_i(n)的初始值为零；

X(n)＝[x(n),x(n-1),...,x(n-M+1)]^T是当前时刻的输入向量，该向量X(n)由当前时刻的输入信号x(n)与之前M-1个时刻的输入信号x(n-1),...,x(n-M+1)构成；

B、计算输出误差向量

B1、将最近P个时刻的输入向量X(n),X(n-1),...,X(n-P+1)构成输入信号矩阵U(n)，U(n)＝[X(n),X(n-1),...,X(n-P+1)]，将最近P个时刻的期望信号d(n),d(n-1),...,d(n-P+1)，构成期望信号向量D(n)，D(n)＝[d(n),d(n-1),...,d(n-P+1)]^T，其中P为仿射投影阶数，P＝2～20；

B2、将输入信号矩阵U(n)通过自适应FIR滤波器后得到其输出信号向量即 $\hat{Y}\left(n\right)=U^T\left(n\right)w\left(n\right);$

B3、将期望信号向量D(n)与输出信号向量相减得到自适应FIR滤波器的输出误差向量e(n)，即

C、比例因子列向量的计算；

C1、若当前时刻n不是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)与其前一时刻的值相等，即q_i(n)＝q_i(n-1)；活性因子q_i(n)的初始值为10^-3与10^-4之间的常数；

若当前时刻n是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)按下式得出：

$q_i\left(n\right)=\frac{1}{2}\left|w_i\left(n\right)\right|+\frac{1}{2}\max \{q_i(n-1),|w_i(n-1)\left|\right\}$

其中，|·|为求绝对值运算，max{·,·}为求最大值运算；

C2、根据自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)，计算出当前时刻自适应FIR滤波器的第i个权值w_i(n)的比例因子g_i(n)， $g_i\left(n\right)=\frac{\max \{q_i\left(n\right),|w_i\left(n\right)\left|\right\}}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\max \{q_i\left(n\right),|w_i\left(n\right)\left|\right\}};$

C3、将当前时刻自适应FIR滤波器的所有权值w_i(n)的比例因子g_i(n)组成当前时刻自适应FIR滤波器的比例因子列向量G(n)即，G(n)＝[g₁(n),g₂(n),...,g_M(n)]^T

D、下一时刻滤波器的权值列向量的计算

D1、根据当前时刻自适应FIR滤波器的比例因子列向量G(n)和当前时刻自适应FIR滤波器的输入向量X(n)，得到M行P列的中间矩阵变量H(n)，其中，表示点乘运算，H_-1(n)由前一时刻的中间矩阵变量H(n-1)的前面P-1个列向量构成；初始时刻n＝0的中间矩阵变量H(0)则由初始时刻n＝0的自适应FIR滤波器的比例因子列向量G(0)生成对角矩阵diag{G(0)}再和初始时刻n＝0的输入信号矩阵U(0)相乘得到，即H(0)＝diag{G(0)}U(0)；

D2、根据B1步的输入信号矩阵U(n)、B3步的输出误差向量e(n)和D1步的中间矩阵变量H(n)，得到下一时刻n+1的自适应FIR滤波器的权值列向量w(n+1)，也即稀疏系统w_o在下一时刻n+1的列向量辨识值

$w(n+1)=w\left(n\right)+\frac{\mathrm{\μH}\left(n\right)\mathrm{sgn}\left[e\left(n\right)\right]}{\sqrt{\mathrm{\δ}+{\left[H\left(n\right)\mathrm{sgn}\right[e\left(n\right)\left]\right]}^{T}H\left(n\right)\mathrm{sgn}\left[e\left(n\right)\right]}}$

其中，sgn[e(n)]表示求输出误差向量e(n)中每个元素的符号运算，δ为正则化参数，其取值为0.01，μ为步长，其取值范围为0＜μ＜＜1；

E、令n＝n+1，重复步骤A、B、C、D的操作，即连续得到稀疏系统w_o在不同时刻的辨识值w(n)。

与现有技术相比，本发明的有益效果是：

一、收敛速度快，跟踪能力强

在本发明中，自适应FIR滤波器每个权值都对应了一个独立的活性因子，即第i个滤波器权值w_i(n)对应的活性因子是q_i(n)，并且q_i(n)的值与w_i(n)的绝对值成正比，使得每个滤波器权值都有一个与其幅度成正比的独立步长。因此，本发明具有更快的收敛速度，尤其当稀疏系统w_o发生突变后，本发明有着更强的跟踪能力。

二、计算复杂度低

中间矩阵变量是通过递归的方式进行计算，从而考虑了比例因子的记忆性，降低了计算复杂度。在自适应FIR滤波器的长度是M，仿射投影阶数是P的情况下，本发明完成一次更新需要(2P+1)M+1次加法运算和(P+3)M+2次乘法运算，而RP-APSA需要(2P+2)M－1次加法运算和(2P+2)M+1次乘法运算；很明显，本发明比RP-APSA减少了M－2次加法运算和(P－1)M－1次乘法运算，并且这个优点会随着仿射投影阶数P的增加而变得更加明显。

下面结合附图和具体实施方式对本发明进行详细说明

附图说明

图1a是本发明与MIP-APA和APSA的权值均方误差归一化曲线。

图1b是本发明与RP-APSA和MIP-APSA的权值均方误差归一化曲线。

图2a是本发明与MIP-APA和APSA的权值均方误差归一化曲线。

图2b是本发明与RP-APSA和MIP-APSA的权值均方误差归一化曲线。

具体实施方式

实施例

一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法，包括以下步骤：

A、获取滤波器的期望信号和输出信号

发送输入信号x(n)给稀疏系统w_o，得到稀疏系统w_o的输出信号d(n)即自适应FIR滤波器的期望信号；同时发送输入信号x(n)给自适应FIR滤波器得到自适应FIR滤波器的输出信号 其中：

n表示当前时刻，上标T表示转置运算，w(n)＝[w₁(n),w₂(n),...,w_i(n),...,w_M(n)]^T为自适应FIR滤波器在当前时刻的权值列向量，也即稀疏系统w_o在当前时刻的列向量辨识值，其长度为M；M的取值根据稀疏系统的情况而定，越复杂的稀疏系统其取值越大，通常用于电话回声消除的稀疏系统辨识时，M取值为128、256、512或1024。w_i(n)为w(n)中的第i个权值,i＝1,2,…M,各个权值w_i(n)的初始值为零。

X(n)＝[x(n),x(n-1),...,x(n-M+1)]^T是当前时刻的输入向量，该向量X(n)由当前时刻的输入信号x(n)与之前M-1个时刻的输入信号x(n-1),...,x(n-M+1)构成。

B、计算输出误差向量

B1、将最近P个时刻的输入向量X(n),X(n-1),...,X(n-P+1)构成输入信号矩阵U(n)，U(n)＝[X(n),X(n-1),...,X(n-P+1)]，将最近P个时刻的期望信号d(n),d(n-1),...,d(n-P+1)，构成期望信号向量D(n)，D(n)＝[d(n),d(n-1),...,d(n-P+1)]^T，其中P为仿射投影阶数，P＝2～20；

B2、将输入信号矩阵U(n)通过自适应FIR滤波器后得到其输出信号向量即 $\hat{Y}\left(n\right)=U^T\left(n\right)w\left(n\right);$

B3、将期望信号向量D(n)与输出信号向量相减得到自适应FIR滤波器的输出误差向量e(n)，即

C、比例因子列向量的计算；

C1、若当前时刻n不是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)与其前一时刻的值相等，即q_i(n)＝q_i(n-1)；活性因子q_i(n)的初始值为10^-3与10^-4之间的常数；

若当前时刻n是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)按下式得出：

$q_i\left(n\right)=\frac{1}{2}\left|w_i\left(n\right)\right|+\frac{1}{2}\max \{q_i(n-1),|w_i(n-1)\left|\right\}$

其中，|·|为求绝对值运算，max{·,·}为求最大值运算；

C2、根据自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)，计算出当前时刻自适应FIR滤波器的第i个权值w_i(n)的比例因子g_i(n)， $g_i\left(n\right)=\frac{\max \{q_i\left(n\right),|w_i\left(n\right)\left|\right\}}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\max \{q_i\left(n\right),|w_i\left(n\right)\left|\right\}};$

C3、将当前时刻自适应FIR滤波器的所有权值w_i(n)的比例因子g_i(n)组成当前时刻自适应FIR滤波器的比例因子列向量G(n)即，G(n)＝[g₁(n),g₂(n),...,g_M(n)]^T

D、下一时刻滤波器的权值列向量的计算

D1、根据当前时刻自适应FIR滤波器的比例因子列向量G(n)和当前时刻自适应FIR滤波器的输入向量X(n)，得到M行P列的中间矩阵变量H(n)，其中，表示点乘运算，H_-1(n)由前一时刻的中间矩阵变量H(n-1)的前面P-1个列向量构成；初始时刻n＝0的中间矩阵变量H(0)则由初始时刻n＝0的自适应FIR滤波器的比例因子列向量G(0)生成对角矩阵diag{G(0)}再和初始时刻n＝0的输入信号矩阵U(0)相乘得到，即H(0)＝diag{G(0)}U(0)；

D2、根据B1步的输入信号矩阵U(n)、B3步的输出误差向量e(n)和D1步的中间矩阵变量H(n)，得到下一时刻n+1的自适应FIR滤波器的权值列向量w(n+1)，也即稀疏系统w_o在下一时刻n+1的列向量辨识值

$w(n+1)=w\left(n\right)+\frac{\mathrm{\μH}\left(n\right)\mathrm{sgn}\left[e\left(n\right)\right]}{\sqrt{\mathrm{\δ}+{\left[H\left(n\right)\mathrm{sgn}\right[e\left(n\right)\left]\right]}^{T}H\left(n\right)\mathrm{sgn}\left[e\left(n\right)\right]}}$

其中，sgn[e(n)]表示求输出误差向量e(n)中每个元素的符号运算，δ为正则化参数，其取值为0.01，μ为步长，其取值范围为0＜μ＜＜1；

E、令n＝n+1，重复步骤A、B、C、D的操作，即连续得到稀疏系统w_o在不同时刻的辨识值w(n)。

仿真实验：

为了验证本发明的有效性，对已知的稀疏系统w_o进行了辨识仿真实验，并与算法MIP-APA、APSA、RP-APSA和MIP-APSA做了对比。

一、仿真条件

已知稀疏系统w_o有M＝100个系数，其中第1、30、35、85个系数的值分别是0.1、1.0、-0.5和0.1，其余系数的值为0。同时，为了比较各算法的跟踪能力，该稀疏系统w_o在时刻n＝10000处突然变为-w_o。通过稀疏度量化公式可知，该稀疏系统w_o的稀疏程度为φ＝0.9435，是高

度稀疏的，其中||·||₁和||·||₂分别表示求向量的l₁和l₂范数。期望信号d(n)可通过公式d(n)＝x^T(n)w_o+v(n)+η(n)计算获得，其中背景噪声v(n)是信噪比(Signal-to-noise ratio，SNR)为30dB的白色高斯随机序列，η(n)表示冲击干扰。通常，冲击干扰被建模为η(n)＝z(n)A(n)，其中z(n)是伯努利随机过程序列，其概率密度函数为p{z(n)＝1}＝P_r和p{z(n)＝0}＝1-P_r，A(n)是零均值白色高斯序列，其方差为在本实验中，我们设置冲击干扰发生的概率P_r＝0.01和其中y(n)＝x^T(n)w_o。通过使用辨识出的权值w(n)与其已知稀疏系统的权值w_o的均方误差归一化曲线来比较各算法的性能；均方误差归一化值定义为单位为分贝，并且所有的仿真曲线都是独立运行50次的平均结果。这些算法的参数取值如表1所示，是为了保证它们具有相同的稳态均方权值误差，然后公平地比较它们的收敛速度与跟踪能力。

二、仿真结果

2.1一阶自回归相关输入

输入信号x(n)是由零均值且方差为1的白色高斯序列通过传递函数为H₁(z)＝1/(1-0.9z^-1)的一阶自回归系统获得。图1a是本发明与MIP-APA和APSA方法的权值均方误差归一化曲线，图1b是本发明与RP-APSA和RIP-APSA方法的权值均方误差归一化曲线。

表1各算法的参数取值

从图1a和图1b中可以看到：1)在冲击干扰的环境中，仅仅只有MIP-APA发散，是因为它源自于l₂范数的最优化，而APSA，RP-APSA、MIP-APSA和本发明是收敛的，是因为它们源自于l₁范数的最优化；2)对于仿真的高度稀疏系统(稀疏程度为φ＝0.9435)，与RP-APSA和MIP-APSA相比，本发明有最快的收敛速度。此外，当稀疏系统在时刻n＝10000发生突变后，本发明也有更强的跟踪能力。

2.2二阶自回归相关输入

输入信号x(n)是由零均值且方差为1的白色高斯序列通过传递函数为H₂(z)＝1/(1-0.4z^-1+0.4z^-2)的二阶自回归系统获得。图2a是本发明与MIP-APA和APSA的权值均方误差归一化曲线，图2b是本发明与RP-APSA和RIP-APSA的权值均方误差归一化曲线。从图2a和图2b中同样可以得到与图1a和图1b中相同的结论，即对于高度稀疏的稀疏系统，本发明提供了更快的收敛速度以及更强的跟踪能力。

Claims (1)

1.一种基于独立活性因子的抗冲击干扰的自适应稀疏系统辨识方法，包括以下步骤： 

A、获取滤波器的期望信号和输出信号 

发送输入信号x(n)给稀疏系统w_o，得到稀疏系统w_o的输出信号d(n)即自适应FIR滤波器的期望信号；同时发送输入信号x(n)给自适应FIR滤波器得到自适应FIR滤波器的输出信号其中： 

n表示当前时刻，上标T表示转置运算，w(n)＝[w₁(n),w₂(n),...,w_i(n),...,w_M(n)]^T为自适应FIR滤波器在当前时刻的权值列向量，也即稀疏系统w_o在当前时刻的列向量辨识值，其长度为M,；w_i(n)为w(n)中的第i个权值,i＝1,2,…M,各个权值w_i(n)的初始值为零； 

X(n)＝[x(n),x(n-1),...,x(n-M+1)]^T是当前时刻的输入向量，该向量X(n)由当前时刻的输入信号x(n)与之前M-1个时刻的输入信号x(n-1),...,x(n-M+1)构成； 

B、计算输出误差向量 

B1、将最近P个时刻的输入向量X(n),X(n-1),...,X(n-P+1)构成输入信号矩阵U(n)，U(n)＝[X(n),X(n-1),...,X(n-P+1)]，将最近P个时刻的期望信号d(n),d(n-1),...,d(n-P+1)，构成期望信号向量D(n)，D(n)＝[d(n),d(n-1),...,d(n-P+1)]^T，其中P为仿射投影阶数，P＝2～20； 

B2、将输入信号矩阵U(n)通过自适应FIR滤波器后得到其输出信号向量 即

B3、将期望信号向量D(n)与输出信号向量相减得到自适应FIR滤波器的输出误差向量e(n)，即

C、比例因子列向量的计算； 

C1、若当前时刻n不是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)与其前一时刻的值相等，即q_i(n)＝q_i(n-1)；活性因子q_i(n)的初始值为10^-3与10^-4之间的常数； 

若当前时刻n是滤波器长度M的整数倍，则自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)按下式得出： 

其中，|·|为求绝对值运算，max{·,·}为求最大值运算； 

C2、根据自适应FIR滤波器在当前时刻的权值列向量w(n)的第i个权值w_i(n)的活性因子q_i(n)，计算出当前时刻自适应FIR滤波器的第i个权值w_i(n)的比例因子g_i(n)，

C3、将当前时刻自适应FIR滤波器的所有权值w_i(n)的比例因子g_i(n)组成当前时刻自适应FIR滤波器的比例因子列向量G(n)即，G(n)＝[g₁(n),g₂(n),...,g_M(n)]^T

D、下一时刻滤波器的权值列向量的计算 

D1、根据当前时刻自适应FIR滤波器的比例因子列向量G(n)和当前时刻自适应FIR滤波器的输入向量X(n)，得到M行P列的中间矩阵变量H(n)， 其中，表示点乘运算，H_-1(n)由前一时刻的中间矩阵变量H(n-1)的前面P-1个列向量构成；初始时刻n＝0的中间矩阵变量H(0)则由初始时刻n＝0的自适应FIR滤波器的比例因子列向量G(0)生成对角矩阵diag{G(0)}再和初始时刻n＝0的输入信号矩阵U(0)相乘得到，即H(0)＝diag{G(0)}U(0)； 

D2、根据B1步的输入信号矩阵U(n)、B3步的输出误差向量e(n)和D1步 的中间矩阵变量H(n)，得到下一时刻n+1的自适应FIR滤波器的权值列向量w(n+1)，也即稀疏系统w_o在下一时刻n+1的列向量辨识值 

其中，sgn[e(n)]表示求输出误差向量e(n)中每个元素的符号运算，δ为正则化参数，其取值为0.01，μ为步长，其取值范围为0＜μ＜＜1； 

E、令n＝n+1，重复步骤A、B、C、D的操作，即连续得到稀疏系统w_o在不同时刻的辨识值w(n) 。