基于鲁棒自适应盲系统辨识的多通道时延估计方法

摘要

本发明公开了一种基于鲁棒自适应盲系统辨识的多通道时延估计方法，初始化各传声器对应通道的自适应滤波器系数，然后各传声器分别对声信号进行采集并分帧保存，利用各传声器在同一时间段采集的信号帧进行多通道时延估计：首先由各传声器采集的时域信号与其前一帧时域信号、前一帧时域信号对应的自适应滤波器系数以及Cauchy估计器的一阶和二阶导数计算自适应滤波器系数的更新量，利用该更新量对前一次迭代获得的自适应滤波器系数进行更新，得到当前帧对应的自适应滤波器系数，然后根据计算得到的自适应滤波器系数进行时延估计。本发明使用Cauchy估计器定义新的频域代价函数，提升时延估计的性能。

说明书

技术领域

本发明属于声源定位技术领域，更为具体地讲，涉及一种基于鲁棒自适应盲系统辨识的多通道时延估计方法。

背景技术

时延估计，其目的是利用传感器阵列采集的信号估计声源与阵列间的到达时间差，在免提语音通信的声源定位和跟踪中扮演着重要角色。众多学者对这一主题进行了广泛的研究、开发了许多时延估计方法，包括著名的广义互相关方法、基于盲通道辨识的方法、多通道线性预测方法、基于信息论的方法等。

在时延估计方法中，文献“Y.Huang and J.Benesty,“A class of frequency-domain adaptive approaches to blind multichannel identification,”IEEE Trans.Signal Process.,vol.51,no.1,pp.11–24,Jan.2003.”和“J.Chen,Y.Huang,and J.Benesty,“An adaptive blind SIMO identification approach to joint multichannel time delay estimation,”in Proc.IEEE Int.Conf.Acoust.,Speech,Signal Process.(ICASSP),2004,pp.IV-53–IV-56.”中所公开的基于归一化多通道频域最小均方(NMCFLMS)自适应盲辨识算法受到极大关注，该方法的核心思想是先对声源到各传声器间的通道脉冲响应进行盲辨识，接着比较不同通道脉冲响应直达成分的到达时间差来确定时延。经验证，该算法在混响环境下获得了良好的鲁棒性，然而，它对噪声非常敏感，容易受到噪声影响而降低性能。

在文献“M.A.Haque and M.K.Hasan,“Noise robust multichannel frequency-domain LMS algorithms for blind channel identification,”IEEE Signal Process.Lett.,vol.15,pp.305–308,2008.”和“H.He,J.Chen,J.Benesty,and T.Yang,“Multichannel time delay estimation for acoustic source localization via robust adaptive blind system identification,”in Proc.Int.Workshop Acoust.Signal Enhancement(IWAENC),2016,to appear.”中公开了一种鲁棒归一化多通道频域最小均方(RNMCFLMS)自适应算法中，对通道脉冲响应引入了一种平坦谱约束，对NMCFLMS算法进行了改进，提高了噪声环境下的时延估计性能，但强噪声环境下的时延估计鲁棒性仍然是一个极大的挑战。而在文献“H.He,J.Lu,J.Chen,X.Qiu,and J.Benesty,“Robust blind identification of room acoustic channels in symmetric alphastable distributed noise environments,”J.Acoust.Soc.Amer.,vol.136,no.8,pp.693–704,Aug.2014.”和“H.He,J.Chen,J.Benesty,and T.Yang,“Multichannel time delay estimation for acoustic source localization via robust adaptive blind system identification,”in Proc.Int.Workshop Acoust.Signal Enhancement(IWAENC),2016,to appear.”中公开的另一种鲁棒归一化多通道频域最小平均M估计(RNMCFLMM)算法中，利用Huber估计器构造了一种鲁棒时域代价函数，由此导出的自适应滤波器应用于声通道盲辨识对时延进行估计，获得了对非高斯和高斯噪声的鲁棒性。然而，在低信噪比(SNR)环境下，语音的非平稳和有色特性仍然对时延估计精度产生不利影响。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于鲁棒自适应盲系统辨识的多通道时延估计方法，使用Cauchy估计器定义新的频域代价函数，利用该估计器对语音信号的脉冲谱的鲁棒性降低自适应滤波器的失调，从而提升时延估计的性能。

为实现上述发明目的，本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法包括以下步骤：

S1：初始化自适应滤波器系数其中为长度为L的系数列向量，k＝1,2,…,M，设置自适应步长μ_f、遗忘因子λ、Cauchy估计器的参数c；

S2：M个传声器分别对声信号进行持续采集，每个传声器将采集的时域信号以每L个样本为一帧进行保存，第i个传声器的第m帧时域信号记为x_i,L(m)＝[x_i(mL-L)>i(mL-L+1)…x_i(mL-1)]^T，其中i＝1,2,…,M，m＝1,2,…；当第1帧样本采集并保存完毕后，每当完成一帧样本采集，则进入步骤S3；

S3：根据每个传声器的第m帧、第m+1帧时域信号及其对应的自适应滤波器系数迭代计算得到第m+1帧时域信号对应的自适应滤波器系数其具体步骤包括：

S3.1：计算自适应滤波器系数向量的傅里叶变换向量

${\widehat{\underset{‾}{h}}}_k^{10}\left(m\right)=F_{2L\times 2L}{\left(\begin{array}{cc}{\hat{h}}_k^T\left(m\right)& 0_{1\times L}\end{array}\right)}^T$

其中，k＝1,2,…,M，F_2L×2L表示维度为2L×2L的傅里叶矩阵，0_1×L表示长度为L的所有元素为0的列向量；

S3.2：计算时域信号的傅里叶变换向量x_i,2L(m)：

x_i,2L(m)＝F_2L×2Lx_i,2L(m)

其中，i＝1,2,…,M，x_i,2L(m)＝[x_i(mL-L)>i(mL-L+1)…x_i(mL+L-1)]^T；

S3.3：按照以下步骤依次计算得到时域先验误差信号的傅里叶变换向量e_ik(m)：

${\bar{e}}_{ik}\left(m\right)=F_{2L\times 2L}^{-1}{\underset{‾}{\bar{e}}}_{ik}\left(m\right)$

$e_{ik}\left(m\right)=\left(\begin{array}{cc}{0}_{L\times L}& I_{L\times L}\end{array}\right){\bar{e}}_{ik}\left(m\right)$

e_ik(m)＝F_L×Le_ik(m)＝[e_ik(mL)>e_ik(mL+1)…e_ik(mL+L-1)]^T

其中，⊙表示两个向量点乘，上标-1表示矩阵求逆，0_L×L是维度为L×L的零矩阵，I_L×L是维度为L×L的单位矩阵，F_L×L是维度为L×L的傅里叶矩阵；

S3.4：计算Cauchy估计器的一阶导数ρ′[|e_ik(n)|]：

${\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]=\frac{c^2\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}{c^2+{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2},n=mL,mL+1,...,mL+L-1$

S3.5：计算先验误差信号的Cauchy变换向量ψ[e_ik(m)]及其变换向量ψ⁰¹[e_ik(m)]：

$\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=\left(\begin{array}{c}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(mL\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}\left(mL\right)\left]\right\}\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+1)\left]\right\}\\ .\\ .\\ .\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+L-1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+L-1)\left]\right\}\end{array}\right)$

${\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=G_{2L\times L}^{01}\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]$

其中，exp{·}表示以常数e为底的指数函数，arg{·}表示求辐角主值；

S3.6：计算传声器采集的信号的谱矩阵

$D_{x_i}\left(m\right)=diag\left\{F_{2L\times 2L}{x}_{i,2L}\left(m\right)\right\}$

S3.7：计算Cauchy估计器的二阶导数ρ″[|e_ik(n)|]：

${\rho }^{\prime \prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]=\frac{c^2[c^2-{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2]}{{[c^2+{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2]}^2},n=mL,mL+1,...,mL+L-1$

S3.8：计算加权系数φ_ik(m)：

φ_ik(m)＝max{η_ik(mL),η_ik(mL+1),…η_ik(mL+L-1)}

其中：

${\underset{‾}{\eta }}_{ik}\left(n\right)={\rho }^{\prime \prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]+\frac{{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]}{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|},n=mL,mL+1,...,mL+L-1$

S3.9：计算传声器采集的信号的加权功率谱矩阵P_k(m)：

$P_k\left(m\right)=\frac{1}{2}\underset{i=1,i\ne k}{\overset{M}{\Sigma }}{\phi }_{ik}\left(m\right)D_{x_i}^{*}\left(m\right)D_{x_i}\left(m\right)$

S3.10：计算递推功率谱矩阵Q_k(m)：

Q_k(m)＝λQ_k(m-1)+(1-λ)P_k(m)

其中，λ表示遗忘因子，0＜λ＜1；

S3.11：计算第m+1帧时域信号对应的自适应滤波器系数增长向量

${\hat{h}}_k^{10}(m+1)={\hat{h}}_k^{10}\left(m\right)-F_{2L\times 2L}^{-1}[{\mu }_f▿J_{NF,k}^{01}\left(m\right)-{\mu }_f\beta \left(m\right)▿J_{SC,k}^{10}\left(m\right)]$

其中：

$▿J_{NF,k}^{01}\left(m\right)=Q_k^{-1}\left(m\right)\underset{i=1}{\overset{M}{\Sigma }}{D}_{x_i}^{*}\left(m\right){\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]$

$\beta \left(m\right)=\left|\frac{{[▿J_{SC}^{10}\left(m\right)]}^H▿J_{NF}^{01}\left(m\right)}{\left|\right|▿J_{SC}^{10}\left(m\right)|{|}_2^2}\right|$

$▿J_{SC}^{10}\left(m\right)={\left(\begin{array}{cccc}{[▿J_{SC,1}^{10}\left(m\right)]}^T& {[▿J_{SC,2}^{10}\left(m\right)]}^T& \mathrm{...}& {[▿J_{SC,M}^{10}\left(m\right)]}^T\end{array}\right)}^T$

$▿J_{NF}^{01}\left(m\right)={\left(\begin{array}{cccc}{[▿J_{NF,1}^{01}\left(m\right)]}^T& {[▿J_{NF,2}^{01}\left(m\right)]}^T& \mathrm{...}& {[▿J_{NF,M}^{01}\left(m\right)]}^T\end{array}\right)}^T$

其中，表示两个长度相同的向量按元素相除，1_2L×1表示长度为2L的所有元素为1的向量，|·|²表示按元素进行模的平方运算，||·||₂表示对向量求二范数；

S3.12：计算第m+1帧时域信号对应的自适应滤波器系数向量

${\hat{h}}_k(m+1)=\left(\begin{array}{cc}{I}_{L\times L}& 0_{L\times L}\end{array}\right){\hat{h}}_k^{10}(m+1)$

S4：根据自适应滤波器系数进行当前时刻的时延估计：

${\hat{\tau }}_{k,k^{\prime }}(m+1)=\arg \underset{l}{max}\left|{\hat{h}}_{k,l}(m+1)\right|-\arg \underset{l}{max}\left|{\hat{h}}_{k^{\prime },l}(m+1)\right|$

其中，和分别表示和中的第l个元素，l＝1,2,…,L。

本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法，初始化各传声器对应通道的自适应滤波器系数，然后各传声器分别对声信号进行采集并分帧保存，利用各传声器在同一时间段采集的信号帧进行多通道时延估计：首先由各传声器采集的时域信号与其前一帧时域信号、前一帧时域信号对应的自适应滤波器系数以及Cauchy估计器的一阶和二阶导数计算自适应滤波器系数的更新量，利用该更新量对前一次迭代获得的自适应滤波器系数进行更新，得到当前帧对应的自适应滤波器系数，然后根据计算得到的自适应滤波器系数进行时延估计。

本发明使用Cauchy估计器定义一个频域代价函数，利用该估计器对具有脉冲谱特征的声源信号的鲁棒性以降低自适应滤波器的失调。经实验证明，对于具有脉冲谱特性的语音信号，本发明可以实现较为准确的多通道时延估计。

附图说明

图1是频域Cauchy估计器与二次函数的比较曲线图；

图2是本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法的流程图；

图3是本发明和四种对比算法对通道1和通道2间的时延估计结果对比图；

图4是本发明和四种对比算法对通道1和通道3间的时延估计结果对比图；

图5是本发明和四种对比算法对通道2和通道3间的时延估计结果对比图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

为了更好说明本发明的技术方案，对鲁棒自适应盲系统辨识的原理进行推导说明。

假设一个单输入多输出声系统由一个声源和M个传声器组成，则该系统有如下输入输出关系：

x_k(n)＝s(n)*h_k+v_k(n),k＝1,2,…,M>

式中，x_k(n)表示第k个传声器的输出信号，s(n)表示声源信号，*表示线性卷积，h_k是一个有限冲击响应滤波器，代表声源到第k个传声器间的通道脉冲响应，v_k(n)是对应第k个传声器的加性噪声。如果忽略(1)式中的加性噪声项，则可获得两个传声器采集的信号和通道脉冲响应之间具有如下关系：

$\begin{array}{c}{x}_i\left(n\right)*h_j=s\left(n\right)*h_i*h_j\\ =x_j\left(n\right)*h_i,i,j=1,2,\mathrm{...},M,i\ne j\end{array}---\left(2\right)$

上式可写成向量形式：

$x_i^T\left(n\right)h_j-x_j^T\left(n\right)h_i=0---\left(3\right)$

式中x_i(n)表示第i个通道长度为L的观测信号向量：

x_i(n)＝[x_i(n)>i(n-1)…x_i(n-L+1)]^T>

h_i是第i个通道长度为L的脉冲响应向量：

h_i＝[h_i,0>i,1…h_i,L-1]^T>

其中(·)^T表示向量或矩阵的转置。

当噪声存在或者脉冲响应的估计值偏离真实值时，(3)式不再等于零，于是在时刻n时通道i和通道j之间便产生一先验误差信号e_ij(n)：

$e_{ij}\left(n\right)=x_i^T\left(n\right){\hat{h}}_j\left(n\right)-x_j^T\left(n\right){\hat{h}}_i\left(n\right)---\left(6\right)$

式中是h_i在时刻n的一个估计值，是h_j在时刻n的一个估计值。该先验误差信号可用于定义一个代价函数，通过对代价函数最小化可获得脉冲响应的最优估计。

本发明将先验误差信号e_ij(n)变换到频域，开发一种鲁棒自适应频域盲系统辨识算法。对于典型的声通道，其脉冲响应的幅度谱具有平坦性。然而，对于语音激励信号，其幅度谱具有明显的脉冲状特性。在低信噪比(SNR)环境，通常和不能精确估计，意味着不等于0。由此，强调了语音声源信号s(n)的脉冲谱在先验误差信号e_ij(n)的幅度谱中的存在。结果，语音信号中的脉冲谱在基于二次函数的频域代价函数中占据主导地位，这严重影响通道辨识的精确性和鲁棒性。为提高语音激励下的声通道辨识性能，本发明在频域定义一个新的递推最小二乘代价函数如下：

$J\left(m\right)=(1-\lambda )\underset{r=0}{\overset{m}{\Sigma }}{\lambda }^{m-r}{J}_{\rho }\left(r\right)---\left(7\right)$

式中m表示块序号，λ表示遗忘因子，0＜λ＜1，J_ρ(r)的表达式为：

$J_{\rho }\left(r\right)=\underset{i=1}{\overset{M-1}{\Sigma }}\underset{j=i+1}{\overset{M}{\Sigma }}\underset{n=rL}{\overset{rL+L-1}{\Sigma }}\rho \left[\right|{\underset{‾}{e}}_{ij}\left(n\right)\left|\right]---\left(8\right)$

其中，当r＝m时，n＝mL,mL+1,…,mL+L-1，e_ij(n)表示时域先验误差块信号的傅里叶变换频域分量，其傅里叶变换向量e_ij(m)定义为：

e_ij(m)＝[e_ij(mL)>e_ij(mL+1)…e_ij(mL+L-1)]^T>

根据(6)式，(9)式可进一步表示成：

${\underset{‾}{e}}_{ij}\left(m\right)=G_{L\times 2L}^{01}[D_{x_i}\left(m\right)G_{2L\times L}^{10}{\underset{‾}{\hat{h}}}_j\left(m\right)-D_{x_j}\left(m\right)G_{2L\times L}^{10}{\underset{‾}{\hat{h}}}_i\left(m\right)]---\left(10\right)$

其中：

$D_{x_i}\left(m\right)=\mathrm{diag}\left\{F_{2L\times 2L}{x}_{i,2L}\left(m\right)\right\}---\left(11\right)$

$D_{x_j}\left(m\right)=\mathrm{diag}\left\{F_{2L\times 2L}{x}_{j,2L}\left(m\right)\right\}---\left(12\right)$

x_i,2L(m)＝[x_i(mL-L)>i(mL-L+1)…x_i(mL+L-1)]^T>

x_j,2L(m)＝[x_j(mL-L)>j(mL-L+1)…x_j(mL+L-1)]^T>

$G_{L\times 2L}^{01}=F_{L\times L}\left(\begin{array}{cc}{0}_{L\times L}& I_{L\times L}\end{array}\right)F_{2L\times 2L}^{-1}---\left(15\right)$

$G_{2L\times L}^{10}=F_{2L\times 2L}{\left(\begin{array}{cc}{I}_{L\times L}& 0_{L\times L}\end{array}\right)}^T{F}_{L\times L}^{-1}---\left(16\right)$

${\widehat{\underset{‾}{h}}}_i\left(m\right)=F_{L\times L}{\hat{h}}_i\left(m\right)---\left(17\right)$

${\underset{‾}{\hat{h}}}_j\left(m\right)=F_{L\times L}{\hat{h}}_j\left(m\right)---\left(18\right)$

F_L×L和分别是维度为L×L的傅里叶矩阵和逆傅里叶矩阵，0_L×L是维度为L×L的零矩阵，I_L×L是维度为L×L的单位矩阵，diag{·}表示由一个向量构成的对角矩阵，ρ[·]是一个鲁棒频域M估计器。上标01表示矩阵的构建过程中有乘矩阵[0_L×L>L×L]的操作，上标10表示矩阵的构建过程中有乘矩阵[I_L×L>L×L]^T的操作。

本发明使用Cauchy估计器作为M估计器，于是通道i和通道j间的Cauchy估计器可写作：

$\rho \left[\right|{\underset{‾}{e}}_{ij}\left(n\right)\left|\right]=\frac{c^2}{2}log[1+{\left(\frac{\left|{\underset{‾}{e}}_{ij}\left(n\right)\right|}{c}\right)}^2]---\left(19\right)$

式中，参数c是一个正常数。

图1是频域Cauchy估计器与二次函数的比较曲线图。如图1所示，对于平坦的语音谱导致的较小的先验误差，这两种代价函数具有类似的变化率，使得两种代价函数对这种谱产生相似的自适应能力。对于脉冲语音谱导致的较大的先验误差，图1中二次代价函数对这种较大的先验误差进行平方，产生更大的误差，然而图1中的Cauchy估计器并不强调这个大的误差。本发明正是利用频域Cauchy估计器的这种鲁棒性来提高频域自适应滤波器的性能。

本发明使用牛顿迭代法推导自适应滤波算法，那么基于(7)式的代价函数，可获得该最优化问题的递推解如下：

${\underset{‾}{\hat{h}}}_k(m+1)={\underset{‾}{\hat{h}}}_k\left(m\right)-2\mu (1-\lambda )S_k^{-1}\left(m\right)▿J_{\rho }\left(m\right),k=1,2,\mathrm{...},M---\left(20\right)$

$S_k\left(m\right)={\mathrm{\lambda S}}_k(m-1)+(1-\lambda )▿{▿}^H{J}_{\rho }\left(m\right)---\left(21\right)$

式中μ表示自适应步长，λ表示遗忘因子，0＜λ＜1，是J_ρ(m)关于的梯度，表示J_ρ(m)对应的Hessian矩阵，(·)^H表示向量或矩阵的共轭转置，(·)^-1表示求逆矩阵。

为了方便进一步推导，可将e_ij(n)的共轭表示成：

$\begin{array}{c}{\underset{‾}{e}}_{ij}^{*}\left(n\right){\underset{‾}{e}}_{ij}^H\left(m\right)u_{n-mL+1}\\ =[{\underset{‾}{\hat{h}}}_j^H\left(m\right){\left(G_{2L\times L}^{10}\right)}^H{D}_{x_i}^H\left(m\right)-{\underset{‾}{\hat{h}}}_i^H\left(m\right){\left(G_{2L\times L}^{10}\right)}^H{D}_{x_j}^H\left(m\right)]{\left(G_{L\times 2L}^{10}\right)}^H{u}_{n-mL+1}\end{array}---\left(22\right)$

式中u_n-mL+1是单位矩阵I_L×L的第n-mL+1列，由于n＝mL,mL+1,…,mL+L-1，那么n-mL+1＝1,2,…,L，(·)^*表示求共轭矩阵。于是，对的偏导数可推导为：

$\begin{array}{c}\frac{\partial }{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\underset{i=1}{\overset{M-1}{\Sigma }}\underset{j=i+1}{\overset{M}{\Sigma }}\rho \left[\right|{\underset{‾}{e}}_{ij}\left(n\right)\left|\right]\\ =\frac{\partial }{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\left\{\underset{i=1}{\overset{k-1}{\Sigma }}\rho \right[\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|]+\underset{j=i+1}{\overset{M}{\Sigma }}\rho [\left|{\underset{‾}{e}}_{kj}\left(n\right)\right|\left]\right\}\\ =\underset{i=1}{\overset{k-1}{\Sigma }}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]·\frac{\partial \left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}·\frac{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}+\underset{j=i+1}{\overset{M}{\Sigma }}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{kj}\left(n\right)\left|\right]·\frac{\partial \left|{\underset{‾}{e}}_{kj}\left(n\right)\right|}{\partial {\underset{‾}{e}}_{kj}^{*}\left(n\right)}·\frac{\partial {\underset{‾}{e}}_{kj}^{*}\left(n\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\\ =\frac{1}{2}\underset{i=1}{\overset{k-1}{\Sigma }}{G}_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{10}{u}_{n-mL+1}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}\left(n\right)\left]\right\}\\ -\frac{1}{2}\underset{j=i+1}{\overset{M}{\Sigma }}{G}_{L\times 2L}^{10}{D}_{x_j}^{*}\left(m\right)G_{2L\times L}^{10}{u}_{n-mL+1}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{kj}\left(n\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{kj}\left(n\right)\left]\right\}\\ =\frac{1}{2}\underset{i=1}{\overset{M}{\Sigma }}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|\left]\right\}{G}_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{10}{u}_{n-mL+1}\end{array}---\left(23\right)$

式中ρ′[·]是ρ[·]的一阶导数，exp{·}表示以常数e为底的指数函数，arg{·}表示求辐角主值，

$G_{L\times 2L}^{10}=F_{L\times L}\left(\begin{array}{cc}{I}_{L\times L}& 0_{L\times L}\end{array}\right)F_{2L\times 2L}^{-1}---\left(24\right)$

$G_{2L\times L}^{01}=F_{2L\times 2L}{\left(\begin{array}{cc}{0}_{L\times L}& I_{L\times L}\end{array}\right)}^T{F}_{L\times L}^{-1}---\left(25\right)$

F_2L×2L和分别是维度为2L×2L的傅里叶矩阵和逆傅里叶矩阵。

因此，J_ρ(m)关于的梯度为：

$\begin{array}{c}▿J_{\rho }\left(m\right)=2\frac{\partial J_{\rho }\left(m\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\\ =2\underset{i=1}{\overset{M-1}{\Sigma }}\underset{j=i+1}{\overset{M}{\Sigma }}\underset{n=mL}{\overset{mL+L-1}{\Sigma }}\frac{\partial \rho \left[\right|{\underset{‾}{e}}_{ij}\left(n\right)\left|\right]}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\\ =\underset{i=1}{\overset{M}{\Sigma }}\underset{n=mL}{\overset{mL+L-1}{\Sigma }}{G}_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}{u}_{n-mL+1}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}\left(n\right)\left]\right\}\\ =\underset{i=1}{\overset{M}{\Sigma }}{G}_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(n\right)\right]\end{array}---\left(26\right)$

式中

$\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=\left(\begin{array}{c}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(mL\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}\left(mL\right)\left]\right\}\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+1)\left]\right\}\\ .\\ .\\ .\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+L-1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+L-1)\left]\right\}\end{array}\right)---\left(27\right)$

接着，对于Hessian矩阵，可计算如下：

$\begin{array}{c}▿{▿}^H{J}_{\rho }\left(m\right)=2\frac{\partial }{\partial {\widehat{\underset{‾}{h}}}_k^{*}\left(m\right)}{[▿J_{\rho }\left(m\right)]}^H\\ =2\frac{\partial }{\partial {\widehat{\underset{‾}{h}}}_k^{*}\left(m\right)}\left\{\underset{i=1,i\ne k}{\overset{M}{\Sigma }}{\underset{‾}{\psi }}^H\right[{\underset{‾}{e}}_{ik}\left(m\right)\left]G_{L\times 2L}^{01}{D}_{x_i}{G}_{2L\times L}^{10}\right\}\\ =2\underset{i=1,i\ne k}{\overset{M}{\Sigma }}\frac{\partial {\underset{‾}{\psi }}^H\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]}{\partial {\widehat{\underset{‾}{h}}}_k^{*}\left(m\right)}{G}_{L\times 2L}^{01}{D}_{x_i}{G}_{2L\times L}^{10}\end{array}---\left(28\right)$

可以推导得到：

$\begin{array}{c}\frac{\partial }{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}{\left\{{\rho }^{\prime }\right[\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|\left]\exp \right\{\sqrt{-1}\arg \left[{\underset{‾}{e}}_{ik}\left(n\right)\right]\left\}\right\}}^{*}\\ =\frac{\partial \left\{{\rho }^{\prime }\right[\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|\left]\exp \right\{-\sqrt{-1}\arg \left[{\underset{‾}{e}}_{ik}\left(n\right)\right]\left\}\right\}}{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}·\frac{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\\ =\left\{{\rho }^{\prime \prime }\right[\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|\left]\frac{\partial \left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}\exp \right\{-\sqrt{-1}\arg \left[{\underset{‾}{e}}_{ik}\left(n\right)\right]\}\\ +{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]\frac{\partial \left\{\exp \right\{\sqrt{-1}\arg \left[{\underset{‾}{e}}_{ik}\left(n\right)\right]\left\}\right\}}{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}\}·\frac{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\\ =\frac{1}{2}{\underset{‾}{\eta }}_{ik}\left(n\right)\frac{\partial {\underset{‾}{e}}_{ik}^{*}\left(n\right)}{\partial {\underset{‾}{\hat{h}}}_k^{*}\left(m\right)}\end{array}---\left(29\right)$

式中，ρ″[·]是ρ[·]的二阶导数，

${\underset{‾}{\eta }}_{ik}\left(n\right)={\rho }^{\prime \prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]+\frac{{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]}{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}---\left(30\right)$

并由此推导得到

$\begin{array}{c}\frac{\partial {\underset{‾}{\psi }}^H\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]}{\partial {\hat{h}}_k^{*}\left(m\right)}=\frac{1}{2}\left(\begin{array}{cc}{\underset{‾}{\eta }}_{ik}\left(mL\right)\frac{\partial {\underset{‾}{e}}_{ik}^{*}\left(mL\right)}{\partial {\hat{h}}_k^{*}\left(m\right)}& {\underset{‾}{\eta }}_{ik}(mL+1)\frac{\partial {\underset{‾}{e}}_{ik}^{*}(mL+1)}{\partial {\hat{h}}_k^{*}\left(m\right)}\end{array}\right)\\ \left(\begin{array}{cc}\mathrm{...}& {\underset{‾}{\eta }}_{ik}(mL+L-1)\frac{\partial {\underset{‾}{e}}_{ik}^{*}(mL+L-1)}{\partial {\hat{h}}_k^{*}\left(m\right)}\end{array}\right)\\ =\frac{1}{2}\left(\begin{array}{cc}{\underset{‾}{\eta }}_{ik}\left(mL\right)G_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{10}{u}_1& {\underset{‾}{\eta }}_{ik}(mL+1)G_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}{u}_2\end{array}\right)\\ \left(\begin{array}{cc}\mathrm{...}& {\underset{‾}{\eta }}_{ik}(mL+L-1)G_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}{u}_L\end{array}\right)\\ =\frac{1}{2}{G}_{L\times 2L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}{T}_{ik}\left(m\right)\end{array}---\left(31\right)$

式中

T_ik(m)＝diag{η_ik(mL)>η_ik(mL+1)…η_ik(mL+L-1)}>

将(31)代入(28)可得：

$▿{▿}^H{J}_{\rho }\left(m\right)=G_{L\times 2L}^{10}{P}_k\left(m\right)G_{2L\times L}^{10}---\left(33\right)$

式中P_k(m)表示传声器采集的信号的加权功率谱矩阵：

$P_k\left(m\right)=\underset{i=1,i\ne k}{\overset{M}{\Sigma }}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}{T}_{ik}\left(m\right)G_{L\times 2L}^{01}{D}_{x_i}\left(m\right)---\left(34\right)$

根据(21)式可以计算得到：

$\begin{array}{c}{S}_k\left(m\right)={\mathrm{\lambda S}}_k(m-1)+(1-\lambda )▿{▿}^H{J}_{\rho }\left(m\right)\\ =\lambda (1-\lambda )\underset{i=0}{\overset{m-1}{\Sigma }}{\lambda }^{m-1-i}▿{▿}^H{J}_{\rho }\left(i\right)+(1-\lambda )▿{▿}^H{J}_{\rho }\left(m\right)\\ =(1-\lambda )\underset{i=0}{\overset{m-1}{\Sigma }}{\lambda }^{m-i}▿{▿}^H{J}_{\rho }\left(i\right)\end{array}---\left(35\right)$

将(33)代入(35)可得：

$\begin{array}{c}{S}_k\left(m\right)=(1-\lambda )\underset{i=0}{\overset{m}{\Sigma }}{\lambda }^{m-i}{G}_{L\times 2L}^{10}{P}_k\left(i\right)G_{2L\times L}^{10}\\ =G_{L\times 2L}^{10}{Q}_k\left(m\right)G_{2L\times L}^{10}\end{array}---\left(36\right)$

式中Q_k(m)表示递推功率谱矩阵：

$\begin{array}{c}{Q}_k\left(m\right)=(1-\lambda )\underset{i=0}{\overset{m}{\Sigma }}{\lambda }^{m-i}{P}_k\left(i\right)\\ ={\mathrm{\lambda Q}}_k(m-1)+(1-\lambda )P_k\left(m\right)\end{array}---\left(37\right)$

将(26)和(36)代入(20)，并在等式两端左乘可获得多通道频域自适应滤波器的更新方程如下：

$\begin{array}{c}{G}_{2L\times L}^{10}{\underset{‾}{\hat{h}}}_k(m+1)=G_{2L\times L}^{10}{\underset{‾}{\hat{h}}}_k\left(m\right)-2\mu (1-\lambda )G_{2L\times L}^{10}{\left[G_{L\times 2L}^{10}{Q}_k\left(m\right)G_{2L\times L}^{10}\right]}^{-1}\\ \times \underset{i=1}{\overset{M}{\Sigma }}{G}_{2L\times L}^{10}{D}_{x_i}^{*}\left(m\right)G_{2L\times L}^{01}\underset{‾}{\psi }\left[e_{ik}\left(m\right)\right],k=1,2,\mathrm{...},M\end{array}---\left(38\right)$

根据如下等式关系：

$D_{2L\times L}^{10}{\left[D_{L\times 2L}^{10}{Q}_k\left(m\right)D_{2L\times L}^{10}\right]}^{-1}{D}_{L\times 2L}^{10}=D_{2L\times 2L}^{10}{Q}_k^{-1}\left(m\right)---\left(39\right)$

式中：

$\begin{array}{c}{G}_{2L\times 2L}^{10}=G_{2L\times L}^{10}{G}_{L\times 2L}^{10}\\ =F_{2L\times 2L}\left(\begin{array}{cc}{I}_{L\times L}& 0_{L\times L}\\ 0_{L\times L}& 0_{L\times L}\end{array}\right)F_{2L\times 2L}^{-1}\end{array}---\left(40\right)$

可获得简化的更新方程：

${\widehat{\underset{‾}{h}}}_k^{10}(m+1)={\widehat{\underset{‾}{h}}}_k^{10}\left(m\right)-{\mu }_f{Q}_k^{-1}\left(m\right)\underset{i=1}{\overset{M}{\Sigma }}{D}_{x_i}^{*}\left(m\right){\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right],k=1,2,\mathrm{...},M---\left(41\right)$

式中μ_f＝μ/2是新的自适应步长，

${\widehat{\underset{‾}{h}}}_k^{10}\left(m\right)=G_{2L\times L}^{10}{\widehat{\underset{‾}{h}}}_k\left(m\right)---\left(42\right)$

${\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=G_{2L\times L}^{01}\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]---\left(43\right)$

当L较大时，可将(32)式近似表示成：

T_ik(m)≈φ_ik(m)I_L×L>

式中φ_ik(m)表示加权系数：

φ_ik(m)＝max{η_ik(mL),η_ik(mL+1),…η_ik(mL+L-1)}>

注意到

$\begin{array}{c}{G}_{2L\times L}^{01}{T}_{ik}\left(m\right)G_{L\times 2L}^{01}={\phi }_{ik}\left(m\right)G_{2L\times 2L}^{01}\\ \approx \frac{1}{2}{\phi }_{ik}\left(m\right)I_{2L\times 2L}\end{array}---\left(46\right)$

式中

$G_{2L\times 2L}^{01}=F_{2L\times 2L}\left(\begin{array}{cc}{0}_{L\times L}& 0_{L\times L}\\ 0_{L\times L}& I_{L\times L}\end{array}\right)F_{2L\times 2L}^{-1}---\left(47\right)$

于是，(34)式的矩阵P_k(m)可简化成：

$P_k\left(m\right)=\frac{1}{2}\underset{i=1,i\ne k}{\overset{M}{\Sigma }}{\phi }_{ik}\left(m\right)D_{x_i}^{*}\left(m\right)D_{x_i}\left(m\right)---\left(48\right)$

据上式可知，矩阵P_k(m)化简成一个对角阵，这大大降低了对P_k(m)求逆的计算复杂度。

类似于鲁棒归一化多通道频域最小均方(RNMCFLMS)自适应算法，为了利用声通道脉冲响应幅度谱的平坦性，本发明引入一个谱约束到上述算法。因此本发明最终使用的鲁棒自适应滤波算法如下：

${\widehat{\underset{‾}{h}}}_k^{10}(m+1)={\widehat{\underset{‾}{h}}}_k^{10}\left(m\right)-{\mu }_f▿J_{NF,k}^{01}\left(m\right)+{\mu }_f\beta \left(m\right)▿J_{SC,k}^{10}\left(m\right),k=1,2,\mathrm{...},M---\left(49\right)$

式中

$▿J_{NF,k}^{01}\left(m\right)=Q_k^{-1}\left(m\right)\underset{i=1}{\overset{M}{\Sigma }}{D}_{x_i}^{*}\left(m\right){\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]---\left(50\right)$

$\beta \left(m\right)=\left|\frac{{[▿J_{SC}^{10}\left(m\right)]}^H▿J_{NF}^{01}\left(m\right)}{\left|\right|▿J_{SC}^{10}\left(m\right)|{|}_2^2}\right|---\left(52\right)$

$▿J_{SC}^{10}\left(m\right)={\left(\begin{array}{cccc}{[▿J_{SC,1}^{10}\left(m\right)]}^T& {[▿J_{SC,2}^{10}\left(m\right)]}^T& \mathrm{...}& {[▿J_{SC,M}^{10}\left(m\right)]}^T\end{array}\right)}^T---\left(53\right)$

$▿J_{NF}^{01}\left(m\right)={\left(\begin{array}{cccc}{[▿J_{NF,1}^{01}\left(m\right)]}^T& {[▿J_{NF,2}^{01}\left(m\right)]}^T& \mathrm{...}& {[▿J_{NF,M}^{01}\left(m\right)]}^T\end{array}\right)}^T---\left(54\right)$

表示两个长度相同的向量按元素相除，1_2L×1表示长度为2L的所有元素为1的向量，|·|²表示按元素进行模的平方运算，||·||₂表示对向量求二范数。

根据以上推导过程，可以得到本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法的具体实现过程。图2是本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法的流程图。如图2所示，本发明基于鲁棒自适应盲系统辨识的多通道时延估计方法的具体步骤包括：

S201：初始化参数：

初始化自适应滤波器系数其中是对应通道k长度为L的滤波器系数列向量的初始值，k＝1,2,…,M，设置自适应步长μ_f、遗忘因子λ、Cauchy估计器的参数c。

本实施例中设置μ_f＝0.8，Cauchy估计器的参数c＝5.0。

S202：信号采集：

M个传声器分别对声源信号进行持续采集，每个传声器将采集的时域信号以每L个样本为一帧进行保存，第i个传声器的第m帧时域信号记为x_i,_L(m)＝[x_i(mL-L)x_i(mL-L+1)…x_i(mL-1)]^T，其中i＝1,2,…,M，m＝1,2,…；当第1帧样本采集并保存完毕后，每当完成一帧样本采集，则进入步骤S203。也就是说，当m＝2,3,…时，才能进行时延估计。

S203：迭代更新自适应滤波器系数：

根据每个传声器的第m帧、第m+1帧时域信号及其对应的自适应滤波器系数迭代计算得到第m+1帧时域信号对应的自适应滤波器系数根据之前的理论推导过程，其具体步骤包括：

S3.1：计算自适应滤波器系数向量的傅里叶变换向量

${\widehat{\underset{‾}{h}}}_k^{10}\left(m\right)=F_{2L\times 2L}{\left(\begin{array}{cc}{\hat{h}}_k^T\left(m\right)& 0_{1\times L}\end{array}\right)}^T---\left(55\right)$

其中，k＝1,2,…,M，F_2L×2L表示维度为2L×2L的傅里叶矩阵，0_1×L表示长度为L的所有元素为0的列向量。

S3.2：计算传声器采集的时域信号的傅里叶变换向量x_i,2L(m)：

x_i,2L(m)＝F_2L×2Lx_i,2L(m)>

其中，i＝1,2,…,M，x_i,2L(m)＝[x_i(mL-L)>i(mL-L+1)…x_i(mL+L-1)]^T，显然x_i,2L(m)是由第m帧、第m+1帧时域信号拼接而成。

S3.3：计算时域先验误差信号的傅里叶变换向量e_ik(m)：

${\bar{e}}_{ik}\left(m\right)=F_{2L\times 2L}^{-1}{\underset{‾}{\bar{e}}}_{ik}\left(m\right)---\left(58\right)$

$e_{ik}\left(m\right)=\left(\begin{array}{cc}{0}_{L\times L}& I_{L\times L}\end{array}\right){\bar{e}}_{ik}\left(m\right)---\left(59\right)$

e_ik(m)＝F_L×Le_ik(m)＝[e_ik(mL)>e_ik(mL+1)…e_ik(mL+L-1)]^T>

其中，⊙表示两个向量点乘，上标-1表示矩阵求逆，0_L×L是维度为L×L的零矩阵，I_L×L是维度为L×L的单位矩阵，F_L×L是维度为L×L的傅里叶矩阵。

S3.4：计算Cauchy估计器的一阶导数ρ′[|e_ik(n)|]：

${\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]=\frac{c^2\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}{c^2+{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2},n=mL,mL+1,...,mL+L-1---\left(61\right)$

S3.5：计算先验误差的Cauchy变换向量ψ[e_ik(m)]及其变换向量ψ⁰¹[e_ik(m)]：

$\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=\left(\begin{array}{c}{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(mL\right)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}\left(mL\right)\left]\right\}\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+1)\left]\right\}\\ .\\ .\\ .\\ {\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}(mL+L-1)\left|\right]\exp \left\{\sqrt{-1}\arg \right[{\underset{‾}{e}}_{ik}(mL+L-1)\left]\right\}\end{array}\right)---\left(62\right)$

${\underset{‾}{\psi }}^{01}\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]=G_{2L\times L}^{01}\underset{‾}{\psi }\left[{\underset{‾}{e}}_{ik}\left(m\right)\right]---\left(63\right)$

其中，exp{·}表示以常数e为底的指数函数，arg{·}表示求辐角主值。

S3.6：计算传声器采集的信号的谱矩阵

$D_{x_i}\left(m\right)=diag\left\{F_{2L\times 2L}{x}_{i,2L}\left(m\right)\right\}---\left(64\right)$

S3.7：计算Cauchy估计器的二阶导数ρ″[|e_ik(n)|]：

${\rho }^{\prime \prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]=\frac{c^2[c^2-{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2]}{{[c^2+{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|}^2]}^2},n=mL,mL+1,\mathrm{...},mL+L-1---\left(65\right)$

S3.8：计算加权系数φ_ik(m)：

φ_ik(m)＝max{η_ik(mL),η_ik(mL+1),…η_ik(mL+L-1)}>

其中：

${\underset{‾}{\eta }}_{ik}\left(n\right)={\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]+\frac{{\rho }^{\prime }\left[\right|{\underset{‾}{e}}_{ik}\left(n\right)\left|\right]}{\left|{\underset{‾}{e}}_{ik}\left(n\right)\right|},n=mL,mL+1,...,mL+L-1---\left(67\right)$

S3.9：计算传声器采集的信号的加权功率谱矩阵P_k(m)：

$P_k\left(m\right)=\frac{1}{2}\underset{i=1,i\ne k}{\overset{M}{\Sigma }}{\phi }_{ik}\left(m\right)D_{x_i}^{*}\left(m\right)D_{x_i}\left(m\right)---\left(68\right)$

S3.10：计算递推功率谱矩阵Q_k(m)：

Q_k(m)＝λQ_k(m-1)+(1-λ)P_k(m)>

其中，λ表示遗忘因子，0＜λ＜1。

S3.11：计算第m+1帧时域信号对应的自适应滤波器系数增长向量

将式(49)的等号两侧分别乘以即可得到如下计算公式：

${\hat{h}}_k^{10}(m+1)={\hat{h}}_k^{10}\left(m\right)-F_{2L\times 2L}^{-1}[{\mu }_f▿J_{NF,k}^{01}\left(m\right)-{\mu }_f\beta \left(m\right)▿J_{SC,k}^{10}\left(m\right)]---\left(70\right)$

S3.12：计算第m+1帧时域信号对应的自适应滤波器系数向量

${\hat{h}}_k(m+1)=\left(\begin{array}{cc}{I}_{L\times L}& 0_{L\times L}\end{array}\right){\hat{h}}_k^{10}(m+1)---\left(71\right)$

S204：时延估计：

通过步骤S203实现对声源到各传声器间的通道脉冲响应进行自适应盲辨识后，则可通过比较不同通道脉冲响应直达成分间的到达时间差确定相对时延。因此，得到第m+1帧时域信号所对应时刻任意两个不同通道k,k′间的时延估计如下：

${\hat{\tau }}_{k,k^{\prime }}(m+1)=\arg \underset{l}{max}\left|{\hat{h}}_{k,l}(m+1)\right|-\arg \underset{l}{max}\left|{\hat{h}}_{k^{\prime },l}(m+1)\right|---\left(72\right)$

其中，和分别表示和中的第l个元素，l＝1,2,…,L。

实施例

为了说明本发明的技术效果，采用一个具体的实例对本发明进行实验验证，并将其实验效果与四种对比算法进行对比，四种对比算法分别为：经典的广义互相关-相位变换法(GCC-PHAT)(参见文献“Y.Huang,J.Benesty,and J.Chen,Acoustic MIMO Signal Processing.Berlin,Germany:Springer,2006.”)，归一化多通道频域LMS算法(NMCFLMS)(参见文献“J.Chen,Y.Huang,and J.Benesty,“An adaptive blind SIMO identification approach to joint multichannel time delay estimation,”in Proc.IEEE Int.Conf.Acoust.,Speech,Signal Process.(ICASSP),2004,pp.IV-53–IV-56.”)，鲁棒归一化多通道频域LMS算法(RNMCFLMS)(参见文献“M.A.Haque and M.K.Hasan,“Noise robust multichannel frequency-domain LMS algorithms for blind channel identification,”IEEE Signal Process.Lett.,vol.15,pp.305–308,2008.”)和鲁棒归一化多通道频域最小平均M估计算法(RNMCFLMM)(参见文献“H.He,J.Chen,J.Benesty,and T.Yang,“Multichannel time delay estimation for acoustic source localization via robust adaptive blind system identification,”in Proc.Int.Workshop Acoust.Signal Enhancement(IWAENC),2016,to appear.”)。

本次实验验证采用美国贝尔实验室可控混响室的测试数据对算法进行验证。该混响室长6.7m宽6.1m高2.9m。为方便表示，选择该室地面西北角作为坐标原点。使用三个全指向传声器构成的等距线性阵列采集声信号。这三个传声器的位置坐标分别为(2.437m,0.500m,1.400m)，(3.137m,0.500m,1.400m)，(3.837m,0.500m,1.400m)。声源位于(0.337m,3.938m,1.600m)。在打开可控混响室墙壁89％的嵌板情况下(对应的混响时间为280ms)，以48kHz的采样率测量声源与传声器间的声通道传递函数。对获得的脉冲响应以16kHz采样率对其进行降采样，并截取前1024个点作为多通道盲辨识实验中声通道的真实脉冲响应。即本实施例中每帧样本长度L＝1024。

声源信号是采样率为16kHz长度为1分钟的语音信号，前半段是男声语音，后半段是女声语音。将声源信号和测量的通道脉冲响应相卷积，在其基础上加入给定信噪比(SNR)条件下的白高斯噪声，以获得多通道系统的输出。在实验中，相关的参数设置和文献“H.He,J.Lu,J.Chen,X.Qiu,and J.Benesty,“Robust blind identification of room acoustic channels in symmetric alphastable distributed noise environments,”J.Acoust.Soc.Amer.,vol.136,no.8,pp.693–704,Aug.2014.”中的相同。

实验中，对传声器每拾取长度为64ms的一帧信号进行一次时延估计。使用两种性能测试准则，即畸变估计的概率和非畸变估计的根均方误差评价所提算法的性能(有关这两种准则的定义和方法参见文献“H.He,L.Wu,J.Lu,X.Qiu,and J.Chen,“Time difference of arrival estimation exploiting multichannel spatio-temporal prediction,”IEEE Trans.Audio Speech Lang.Process.,vol.21,pp.463-475,Mar.2013.”，“J.P.Ianniello,“Time delay estimation via cross-correlation in the presence of large estimationerrors,”IEEE Trans.Acoust.,Speech,Signal Process.,vol.ASSP-30,pp.998-1003,Dec.1982.”和“B.Champagne,S.Bedard,and A.Stephenne,“Performance of time-delay estimation in the presence of room reverberation,”IEEE Trans.Speech Audio Process.,vol.4,pp.148-152,Mar.1996.”)。声源到三个不同传声器对的真实时延分别是τ_1,2＝19个样本间隔，τ_1,3＝42个样本间隔，τ_2,3＝23个样本间隔。

图3是本发明和四种对比算法对通道1和通道2间的时延估计结果对比图；图4是本发明和四种对比算法对通道1和通道3间的时延估计结果对比图；图5是本发明和四种对比算法对通道2和通道3间的时延估计结果对比图。从这三幅图可以看出，在噪声和混响环境，NMCFLMS算法优于经典的GCC-PHAT算法，这是由于NMCFLMS算法能盲辨识声通道的直达成分。尽管RNMCFLMS算法引入一个声通道脉冲响应的谱能量约束改善了NMCFLMS在噪声环境下对系统盲辨识的性能，但该算法用于时延估计时未展现明显的鲁棒性。RNMCFLMM算法的时延估计性能优于上述三种算法，这是因为该算法在时域使用Huber估计器定义了一个代价函数，由此产生的自适应滤波器在均方误差准则和平均绝对值误差准则间交替运行，对脉冲样本具有识别能力。然而，在这五种时延估计算法中，本发明获得了最好的性能，尤其是在强噪声环境。这种优越性来自于如下事实：对于一帧语音信号的脉冲谱，所提鲁棒频域自适应滤波算法能有效地利用频域Cauchy估计器对这种脉冲谱的免疫性；同时，对于同一帧语音信号的平坦谱，本发明具有与基于二次代价函数的频域自适应滤波器相似的自适应能力。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。