CN105652243B - 多通道群稀疏线性预测时延估计方法 - Google Patents

多通道群稀疏线性预测时延估计方法 Download PDF

Info

Publication number
CN105652243B
CN105652243B CN201610142030.5A CN201610142030A CN105652243B CN 105652243 B CN105652243 B CN 105652243B CN 201610142030 A CN201610142030 A CN 201610142030A CN 105652243 B CN105652243 B CN 105652243B
Authority
CN
China
Prior art keywords
mrow
msub
matrix
time
msubsup
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610142030.5A
Other languages
English (en)
Other versions
CN105652243A (zh
Inventor
何宏森
陈景东
杨涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest University of Science and Technology
Original Assignee
Southwest University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest University of Science and Technology filed Critical Southwest University of Science and Technology
Priority to CN201610142030.5A priority Critical patent/CN105652243B/zh
Publication of CN105652243A publication Critical patent/CN105652243A/zh
Application granted granted Critical
Publication of CN105652243B publication Critical patent/CN105652243B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/22—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/26—Position of receiver fixed by co-ordinating a plurality of position lines defined by path-difference measurements

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

本发明公开了一种多通道群稀疏线性预测时延估计方法，采用各个可能的时延值计算多通道互相关系数，将多通道互相关系数的平方最大值所对应的时延作为时延估计值；多通道互相关系数的计算方法为：对每只传声器拾取的声信号分别截取一个长度为L的信号帧，对这些信号帧进行时移后叠放成信号向量以及信号矩阵，根据F/l1,2范数优化准则建立时延估计模型，求解群稀疏预测系数矩阵，然后计算预测误差相关矩阵，根据预测误差相关矩阵计算得到多通道互相关系数。本发明利用预测系数矩阵列向量的群稀疏特性构建一个F/l1,2范数优化准则，统一了多通道互相关系数方法和多通道空时预测方法，联合应用空间和时间线性预测的白化能力，以提高时延估计性能。

Description

多通道群稀疏线性预测时延估计方法
技术领域
本发明属于声源定位技术领域，更为具体地讲，涉及一种多通道群稀疏线性预测时延估计方法。
背景技术
时延估计是利用传声器阵列拾取的声信号估计声源到各阵元间的到达时间差，在雷达、声呐以及免提语音通信系统的声源定位和跟踪中扮演着极其重要的角色。常见的时延估计方法有广义互相关方法、盲通道辨识方法、信息论方法、基于语音信号的某些特征进行时延估计等。然而在室内声环境下，由于噪声和混响的不利影响，基于传声器阵列的时延估计是一项具有挑战性的工作。
为了提高时延估计对噪声和混响的鲁棒性，一种多通道互相关系数方法被提出，其具体方法可参考文献“J.Chen,J.Benesty,and Y.Huang,“Robust time delayestimation exploiting redundancy among multiple microphones,”IEEETrans.Speech Audio Process.,vol.11,no.6,pp.549-557,Nov.2003.”和“J.Benesty,J.Chen,and Y.Huang,“Time-delay estimation via linear interpolation and cross-correlation,”IEEE Trans.Speech Audio Process.,vol.12,no.5,pp.509-519,Sep.2004.”。相对于传统的双通道方法，多通道互相关系数方法利用多只传声器间的空间冗余信息抑制噪声和混响的影响，极大地提高了时延估计对噪声的鲁棒性。然而，这种多通道方法对混响仍然敏感。多通道空时预测方法对多通道互相关系数方法进行了推广，其具体方法可参考文献“H.He,L.Wu,J.Lu,X.Qiu,and J.Chen,“Time difference of arrivalestimation exploiting multichannel spatio-temporal prediction,”IEEETrans.Audio Speech Lang.Process.,vol.21,pp.463-475,Mar.2013.”。这种方法以目前最优的方式利用空间和时间信息预白化传声器信号，由此获得时延估计对混响的鲁棒性。然而，这种方法在低信噪比条件下性能降低。
线性预测广泛应用于语音处理。传统的线性预测器由长时预测器和短时预测器级联组成。对于语音信号，这种结构的预测系数向量具有高度的稀疏性。该特性目前已成功应用于语音编码领域。然而，当语音信号受到噪声污染时，这种稀疏性降低甚至消失，线性预测器的性能急剧下降。在最近的研究中“H.He,T.Yang,and J.Chen,“On time delayestimation from a sparse linear prediction perspective,”J.Acoust.Soc.Amer.,vol.137,no.2,pp.1044-1047,Feb.2015.”一种稀疏线性预测时延估计方法被提出。这种方法在最小二乘基础上引入一个稀疏正则化项，构成一种l2/l1范数优化方法，提高了时延估计对噪声的鲁棒性。但该方法是在双通道情形下对两个通道的数据分别进行预白化，也即仅仅对数据进行时间预白化，没有综合利用多通道空间和时间预测的白化能力，从而降低该方法的估计性能。
发明内容
本发明的目的在于克服现有技术的不足，提供一种多通道群稀疏线性预测时延估计方法，利用预测系数矩阵列向量的群稀疏特性构建一个F/l1,2范数优化准则，统一了多通道互相关系数方法和多通道空时预测方法，联合应用空间和时间线性预测的白化能力，以提高时延估计性能。
为实现上述发明目的，本发明多通道群稀疏线性预测时延估计方法包括以下步骤：
S1：M只传声器分别对声源信号进行持续采集，第m只传声器采集的时域信号记为Xm＝[xm(0),xm(1),…,xm(L-1)]，其中m＝1,2,…,M，xm(n)表示第m只传声器在时刻n的采集样本，n＝0,1,…,L-1，L表示每只传声器采集的样本数量；
S2：令时移序号d＝1，声源信号到达第1只和第2只传声器间的初始时移p1＝-pmax，pmax表示时延的最大可能值；
S3：分别将第m只传声器采集的信号Xm按照时延fm(pd)进行时移，fm(pd)表示声源信号到达第1只和第m只传声器间的相对时延，该时延是关于时移pd的函数；时移后的信号中每个样本记为xm(n,pd)，将M只传声器在时刻n的时移样本xm(n,pd)进行叠放，得到信号向量x(n,pd)；
S4：求解以下公式，得到预测系数矩阵A(pd)：
其中，||·||F表示矩阵的F范数，表示矩阵的l2范数，λ是正则化参数，取值范围为λ＞0；
X(0,pd)＝[x(0,pd) x(1,pd) … x(K+L-1,pd)]T
Y(-1,pd)＝[y(-1,pd) y(0,pd) … y(K+L-2,pd)]T
其中：
x(w,pd)＝[x1(w,pd),x2(w,pd),…xM(w,pd)]T
xm(w-1,pd)＝[xm(w-1,pd),xm(w-2,pd),…xm(w-K,pd)]T
w＝0,1,…,K+L-1；K表示预测器阶数，其取值范围为K＜L；对于xm(q,pd)，如果q＜0或q＞L-1，则令xm(q,pd)＝0；
S5：计算预测误差矩阵E(0,pd)：
E(0,pd)＝X(0,pd)-Y(-1,pd)A(pd)
S6：计算预测误差相关矩阵R(pd)：
S7：计算多通道互相关系数ρ(pd)：
其中det(·)表示方阵的行列式，rm,m(pd)是矩阵R(pd)的第m个对角元素；
S8：如果pd＜pmax，令pd＝pd+1，返回步骤S3，否则根据以下公式求得时延估计值
本发明多通道群稀疏线性预测时延估计方法，采用各个可能的时延值计算多通道互相关系数，将多通道互相关系数的平方最大值所对应的时延作为时延估计值；多通道互相关系数的计算方法为：对每只传声器拾取的声信号分别截取一个长度为L的信号帧，对这些信号帧进行时移后叠放成信号向量以及信号矩阵，根据F/l1,2范数优化准则建立时延估计模型，求解群稀疏预测系数矩阵，然后计算预测误差相关矩阵，根据预测误差相关矩阵计算得到多通道互相关系数。
本发明将基于稀疏线性预测的l2/l1范数优化方法推广到多通道，提出了一种多通道群稀疏线性预测时延估计方法。本发明利用预测系数矩阵列向量的群稀疏特性构建一个F/l1,2范数优化准则，从时延估计性能角度统一了多通道互相关系数方法和多通道空时预测方法。本发明还可以通过调节正则化参数来构建一组在传声器信号预白化和非预白化间可作折中处理的时延估计器，以适应不同需求。
附图说明
图1是基于最小二乘的多通道空时预测器对纯净语音预测时预测系数矩阵的一个列向量；
图2是基于最小二乘的多通道空时线性预测器对含噪语音预测时预测系数矩阵的一个列向量；
图3是采用本发明对含噪语音预测时预测系数矩阵的一个列向量；
图4是本发明多通道群稀疏线性预测时延估计方法的流程图；
图5是基于增广拉格朗日乘子交替方向法求解预测系数矩阵的流程图；
图6是噪声和混响环境下各方法时延估计性能随传声器个数变化的曲线图；
图7是噪声环境下各方法时延估计性能随混响时间变化的曲线图；
图8是轻度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图；
图9是中度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图；
图10是轻度混响环境下本发明时延估计性能随参数δ变化的曲线图；
图11是中度混响环境下本发明时延估计性能随参数δ变化的曲线图。
具体实施方式
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
实施例
为了更好地说明本发明的技术内容，首先对本发明所使用的信号模型及方法推导进行说明。
假设一宽带声源辐射平面波信号，M只传声器构成的阵列用于接收声信号。如果选择第1只传声器作为参考点，则第m只传声器在时刻n拾取的信号样本xm(n)可表示为：
xm(n)＝αms[n-t-fm(τ)]+wm(n) (1)
其中，m＝1,2,…,M，αm是声信号传输效应导致的衰减因子，s(n)是未知的零均值因果宽带声源信号，t是从声源到第1只传声器间的传输时间，wm(n)是第m只传声器上的加性噪声，假设该噪声与声源信号以及其它传声器的噪声不相关，τ是声源信号到达第1只和第2只传声器间的时延，fm(τ)是声源信号到达第1只和第m只传声器间的相对时延，它是时延τ的函数。本实施例为了易于推导与实验，假设使用等间距的线性阵列，则在远场条件下fm(τ)＝(m-1)τ。
在上述信号模型基础上，时延估计的目的是给定M只传声器采集的样本信号来估计时延τ。对于一个假定的时延p，利用时移信号xm(n+fm(p))对齐传声器信号。当p＝τ时，M只传声器信号对齐，此时M只传声器拾取的信号间相似性最大。为简化符号，将xm(n+fm(p))记为xm(n,p)。
将M只传声器在时刻n拾取的信号样本表示成一个向量：
x(n,p)＝[x1(n,p),x2(n,p),…xM(n,p)]T (2)
其中，xm(n,p)表示第m只传声器在时刻n的信号样本的时移信号，(·)T表示向量或矩阵的转置。
为了便于推导，对第m只传声器的信号在时刻n-1定义另一个向量：
xm(n-1,p)＝[xm(n-1,p),xm(n-2,p),…xm(n-K,p)]T (3)
其中K表示预测器阶数，其取值范围为K＜L。
现在利用M个通道过去的信号向量x1(n-1,p)、x2(n-1,p)、…、xM(n-1,p)来预测x(n,p)，即预测值表示为：
其中是多通道前向预测器的系数矩阵。
于是，可得预测误差向量e(n,p)的表达式为：
e(n,p)可以记为：
e(n,p)＝[e1(n,p),e2(n,p),…eM(n,p)]T (6)
A(p)是KM×M的多通道前向预测系数矩阵，可以表示为：
A(p)＝[A1(p) A2(p) … AM(p)]T (7)
y(n-1,p)是M只传声器在时刻n-1拾取的信号向量，可以表示为：
根据以上推导过程可以看出，公式(5)是采用M个通道时刻n-K到时刻n-1的样本预测M个通道时刻n的样本，以期获得对应的预测系数矩阵。在实际应用中，为了利用预测稀疏矩阵的群稀疏特性以提高时延估计的准确度，需要将(5)写成矩阵方程。因此利用M只传声器采集的长度分别为L的数据帧，根据公式(5)扩展得到：
E(n,p)＝X(n,p)-Y(n-1,p)A(p) (9)
其中：
E(n,p)＝[e(n,p) e(n+1,p) … e(n+K+L-1,p)]T (10)
X(n,p)＝[x(n,p) x(n+1,p) … x(n+K+L-1,p)]T (11)
Y(n-1,p)＝[y(n-1,p) y(n,p) … y(n+K+L-2,p)]T (12)
其中，X(n,p)是(K+L)×M的矩阵，Y(n-1,p)是(K+L)×KM的矩阵。本发明中，将传声器采集的长度为L的样本中第一个样本序号作为时刻n来构建X(n,p)和Y(n-1,p)，即n＝0。由于数据帧的长度为L，那么在构建X(n,p)和Y(n-1,p)时，对于时刻q的信号样本xm(q,p)，如果q＜0或q＞L-1，则令xm(q,p)＝0。
传统线性预测器由长时预测器和短时预测器级联组成，由此预测器系数向量具有高度的稀疏性。因此在多通道线性预测中，预测器系数矩阵也具有稀疏性。以预测器阶数为80、传声器采集的语音帧长度为1024个采样点、4只传声器用于时延估计为例，对多通道预测器的系数矩阵进行分析。图1是基于最小二乘的多通道空时预测器对纯净语音预测时预测系数矩阵的一个列向量。从图1可以看出，基于最小二乘的多通道空时线性预测系数矩阵的列向量是稀疏的，而且预测系数矩阵不同列向量的大能量样本基本上处于相同的位置，这说明所有列向量具有类似的群稀疏特性。图2是基于最小二乘的多通道空时线性预测器对含噪语音预测时预测系数矩阵的一个列向量。图2中语音信号的信噪比SNR＝5dB。从图2可以看出，当存在噪声时，预测器系数矩阵的稀疏特性降低。既然对于纯净语音信号，预测器系数矩阵是群稀疏的，则可利用这一特性来提高噪声环境下线性预测器的鲁棒性。为此，本发明在最小二乘基础上引入一个预测系数矩阵的群稀疏正则化约束，强制预测系数矩阵的列向量具有相同的稀疏结构。于是，提出如下F/l1,2范数优化准则来预处理传声器信号：
其中，||·||F表示矩阵的F范数，λ是正则化参数，取值范围为λ＞0。可表示为：
其中，表示矩阵或向量的l2范数，Ai(p)表示预测系数矩阵A(p)的第i行。
本发明中，令n＝0，求解公式(13)即可得到预测系数矩阵A(p)，再根据预测系数矩阵A(p)即可计算出时延。
图3是采用本发明对含噪语音预测时预测系数矩阵的一个列向量。图3中语音信号的信噪比SNR＝5dB，从图3可以看出，在噪声环境下通过本发明所提出的如下F/l1,2范数优化准则可以获得稀疏的预测系数矩阵。
基于以上推导，可以得到本发明多通道群稀疏线性预测时延估计方法。图4是本发明多通道群稀疏线性预测时延估计方法的流程图。如图4所示，本发明多通道群稀疏线性预测时延估计方法的具体步骤包括：
S401：信号采集：
M只传声器分别对声信号进行持续采集，第m只传声器采集的时域信号记为Xm＝[xm(0),xm(1),…,xm(L-1)]，其中m＝1,2,…,M，xm(n)表示第m只传声器在时刻n的采集样本，n＝0,1,…,L-1，L表示每只传声器采集的样本数量。
S402：初始化时移序号：
令时移序号d＝1，声源信号到达第1只和第2只传声器间的初始时移p1＝-pmax，pmax表示时延的最大可能值。
S403：信号时移：
分别将第m只传声器采集的信号Xm按照时延fm(pd)进行时移，fm(pd)表示声源信号到达第1只和第m只传声器间的相对时延，该时延是关于时移pd的函数。时移后的信号样本即为xm(t+fm(pd))，为简化表达，记第m只传声器在时刻n的时移信号为xm(n,pd)。将M只传声器在时刻n的时移信号xm(n,pd)进行叠放，得到信号向量x(n,pd)。
S404：求解预测系数矩阵：
求解以下公式，得到预测系数矩阵A(pd)：
其中，A(pd)是KM×M的矩阵，||·||F表示矩阵的F范数，表示矩阵或向量的l2范数，λ是正则化参数，取值范围为λ＞0，
X(0,pd)＝[x(0,pd) x(1,pd) … x(K+L-1,pd)]T (16)
Y(-1,pd)＝[y(-1,pd) y(0,pd) … y(K+L-2,pd)]T (17)
其中：
x(w,pd)＝[x1(w,pd),x2(w,pd),…xM(w,pd)]T (18)
xm(w-1,pd)＝[xm(w-1,pd),xm(w-2,pd),…xm(w-K,pd)]T (20)
w＝0,1,…,K+L-1；K表示预测器阶数，其取值范围为K＜L。对于构建矩阵X(0,pd)和Y(-1,pd)的时刻q的信号样本xm(q,pd)，如果q＜0或q＞L-1，则令xm(q,pd)＝0。显然根据公式(16)至公式(20)可知，q的取值范围为[-K,K+L-1]。
可表示为：
其中，表示矩阵或向量的l2范数，Ai(pd)表示预测系数矩阵A(pd)的第i行。
S405：计算预测误差矩阵：
将步骤S404得到的预测系数矩阵，代入(9)式，即可得到预测误差矩阵E(0,pd)，即：
E(0,pd)＝X(0,pd)-Y(-1,pd)A(pd) (22)
S406：计算预测误差相关矩阵：
根据预测误差矩阵E(0,pd)计算预测误差相关矩阵R(pd)：
S407：计算多通道互相关系数：
利用预测误差相关矩阵R(pd)计算多通道互相关系数ρ(pd)：
其中det(·)表示方阵的行列式，rm,m(pd)是矩阵R(pd)的第m个对角元素。
S408：判断是否pd＜pmax，如果是，进入步骤S409，否则进入步骤S410。
S409：令pd＝pd+1，返回步骤S403。
S410：计算得到时延估计值：
根据以下公式求得时延估计值
公式(25)的含义是，在[-pmax,pmax]范围内的各个时移pd中，当其对应的ρ2(pd)最大，即将该pd作为时延估计值。
显然，在步骤S404中，(15)式是一个凸优化问题，可用多种方法求解，例如线性规划、内点法、原始-对偶内点法等。增广拉格朗日乘子交替方向法能有效利用多个变量可分离(可解耦)的特性，因此本实施例中采用该方法求解这一问题。引入一个辅助矩阵Z(pd)，(13)式可等价地表示为：
其增广拉格朗日子问题可表示为
式中是线性约束拉格朗日乘子矩阵；β是惩罚参数，其取值范围为β＞0；<θ(pd),A(pd)-Z(pd)>是矩阵内积，其表达式为：
<θ(pd),A(pd)-Z(pd)>＝tr{θT(pd)[A(pd)-Z(pd)]} (28)
其中，tr(·)表示矩阵的迹；
(27)式中大括号内的第四项、即增广项，用于确保目标函数是严格凸的。(27)式的迭代求解过程中，给定Zk(pd)，θk(pd)，可利用(27)式交替保持一个未知矩阵固定对另一个矩阵最小化求得Ak+1(pd)，Zk+1(pd)，θk+1(pd)。
首先，当Z(pd)＝Zk(pd)，θ(pd)＝θk(pd)固定，(27)式对A(pd)最小化等价于：
其解为：
式中I表示KM×KM的单位矩阵。
当A(pd)＝Ak+1(pd)，θ(pd)＝θk(pd)固定，(27)式对Z(pd)最小化等价于：
Z(pd)的第i行Zi(pd)的解为：
其中表示辅助矩阵Zk+1(pd)的第i行，i＝1,2,…,MK；
soft函数定义为：
最后，线性约束拉格朗日乘子矩阵θ(pd)更新如下：
θk+1(pd)＝θk(pd)+β(Ak+1(pd)-Zk+1(pd)) (34)
因此通过迭代计算(30)、(32)和(34)式获得(13)式的解。图5是基于增广拉格朗日乘子交替方向法求解预测系数矩阵的流程图。如图5所示，本实施例中求解预测系数矩阵A(pd)的具体过程为：
S501：初始化参数：
令迭代次数k＝1，初始化KM×M的辅助矩阵Z1(pd)和拉格朗日乘子矩阵θ1(pd)。
S502：计算预测系数矩阵：
S503：更新辅助矩阵：
S504：更新拉格朗日乘子矩阵：
θk+1(pd)＝θk(pd)+β(Ak+1(pd)-Zk+1(pd)) (37)
S505：判断是否k＜Q，Q表示最大迭代次数，如果是，进入步骤S506，否则进入步骤S507。
S506：令k＝k+1，返回步骤S502。
S507：得到预测系数矩阵：
将当前的预测系数矩阵Ak+1(pd)作为最终的预测系数矩阵A(pd)，即A(pd)＝Ak+1(pd)。
此外，从(15)式可看出，参数λ在控制预测器系数矩阵的稀疏程度方面扮演着重要角色。该参数主要受传声器信号影响，也即受矩阵X(0,p)和Y(-1,p)的影响，因此可以通过如下方式确定参数λ：
其中表示求取矩阵的无穷范数，δ是一个正常数，根据需要进行设置。
在实际应用中，通常需要持续实时地对时延进行估计，因此采用本发明对M只传声器分别采集的长度为L的数据帧进行时延估计后，M只传声器需要分别采集长度为L的新数据帧，再利用本发明进行下一次时延估计。每次进行时延估计时，都以第1个样本(即时刻0的样本)作为基准来构建矩阵，以求解预测系数矩阵。
为了更好地说明本发明的技术效果，采用本发明(记为MCSTGSP)对一个具体的实施例进行实验验证，并以多通道互相关系数(MCCC)方法、具有预白化的多通道互相关系数(预白化MCCC)方法、多通道空时预测(MCSTP)方法作为对照算法对实验结果进行对照。
实验场地为一个7m×6m×3m的房间。使用6只等间距全指向传声器构成的直线阵列拾取传声器信号，阵元间距为0.1m。假设房间地面西南角为坐标原点，房间内任意位置的坐标表示为(x,y,z)。阵列的第1只和第6只传声器分别放置于(3.25,3.00,1.40)和(3.75,3.00,1.40)。声源位于(2.49,1.27,1.40)。
声源信号是一段预先录制的采样率为16kHz的纯净语音信号，其长度约为1分钟。声源到6只传声器间的脉冲响应通过Image模型产生，Image模型可参考文献“J.B.Allenand D.A.Berkley,“Image method for efficiently simulating small-roomacoustics,”J.Acoust.Soc.Amer.,vol.65,pp.943-950,Apr.1979.”脉冲响应长度为2048个采样点。声源信号与对应通道脉冲响应进行卷积获得传声器信号，在其中加入零均值高斯白噪声以控制信噪比SNR。
在实验中，将传声器信号分割成长度为64ms互不重叠的信号帧。对每一帧信号用Hamming窗函数加窗后用于时延估计。采用畸变估计概率和非畸变估计的根均方误差这两种准则评价时延估计方法的性能(有关这两种准则的定义和分类方法参见文献“H.He,L.Wu,J.Lu,X.Qiu,and J.Chen,“Time difference of arrival estimation exploitingmultichannel spatio-temporal prediction,”IEEE Trans.Audio SpeechLang.Process.,vol.21,pp.463-475,Mar.2013.”，“J.P.Ianniello,“Time delayestimation via cross-correlation in the presence of large estimation errors,”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 estimationin the presence of room reverberation,”IEEE Trans.Speech Audio Process.,vol.4,pp.148-152,Mar.1996.”)。用于性能统计的语音帧总数为936帧(帧长为1024个采样点)。声源到前两只传声器间的真实时延为2.0个采样间隔。
图6是噪声和混响环境下各方法时延估计性能随传声器个数变化的曲线图。图6展示了噪声(SNR＝10dB)和混响(混响时长T60＝300ms)环境下，本发明方法与对比方法的时延估计性能随传声器个数变化的关系。其中对于本发明方法，δ＝0.001。可以看出，随着传声器个数增加，四种时延估计方法的畸变估计概率和非畸变估计的根均方误差基本上都降低，说明适当增加传声器个数能有效改善时延估计的鲁棒性。当使用两只传声器进行时延估计时，所有具有预白化能力的时延估计方法都未获得明显的优势，而原始的MCCC方法具有一定的鲁棒性。当使用多只传声器时，尽管MCCC方法获得较小的畸变估计概率，但对应的根均方误差较大。尽管预白化MCCC和多通道空时预测MCSTP方法获得适当的根均方误差，然而其畸变估计概率较大。相比之下，本发明提出的MCSTGSP方法在畸变估计概率和非畸变估计的根均方误差间获得了良好的折中，这表明在噪声和混响环境下本发明用于时延估计是有效的。
图7是在噪声环境下各方法时延估计性能随混响时间变化的曲线图。图7中传声器信号的信噪比SNR＝10dB。其中对于本发明方法，δ＝0.001。如图7所示，当混响时长T60＝0ms时，MCCC获得了最好的性能，表明该方法对噪声具有鲁棒性。随着T60增加，尽管MCCC方法的畸变估计概率较小，但根均方误差较大，表明MCCC对混响鲁棒性较差。尽管MCSTP方法对混响具有良好的鲁棒性，但对噪声不具有鲁棒性。对于预白化MCCC、MCSTP和MCSTGSP这三种具有预白化能力的时延估计方法而言，MCSTGSP方法获得了最好的性能，表明提出的多通道群稀疏预测时延估计方法在不同混响环境是有效的。
图8是轻度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图。图9是中度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图。图8中混响时长T60＝120ms，图9中混响时长T60＝300ms。其中对于本发明方法，δ＝0.001。根据图8和图9可以看出原始的MCCC算法对噪声具有最好的鲁棒性，尤其在低信噪比条件下，然而在高信噪比条件下MCCC算法对混响最敏感。相对于MCCC，预白化MCCC对混响获得了更好的鲁棒性。MCSTP和预白化MCCC在低信噪比条件下获得相当的性能，然而，由于其最优的预白化能力，MCSTP对混响更具有鲁棒性。尽管这两种具有预白化能力的时延估计方法的性能在混响条件下获得了较大的提高，但在噪声影响下其性能降低。本发明MCSTGSP在MCCC和MCSTP间获得了良好的折中，尤其是MCSTGSP增强了对噪声的免疫力，这证实了利用预测系数矩阵的群稀疏特性可提高多通道空时预测对噪声的鲁棒性。
图10是轻度混响环境下本发明时延估计性能随参数δ变化的曲线图。图11是中度混响环境下本发明时延估计性能随参数δ变化的曲线图。图10中混响时长T60＝120ms，图11中混响时长T60＝300ms。从图10和图11可以看出，一方面正则化参数δ越小，MCSTGSP越接近MCSTP方法，极限情况是δ＝0，即MCSTGSP退化成MCSTP。另一方面，随着δ增加，预测系数矩阵越来越稀疏，因此MCSTGSP对噪声越鲁棒，对混响越敏感；极限情况是传声器信号未进行预白化处理直接用于计算MCCC，因此对应的MCSTGSP退化成原始的MCCC。因此调节正则化参数δ的不同取值，MCSTGSP方法可构成一组具有不同程度预白化能力的时延估计器。
以上实验结果表明，提出的MCSTGSP方法在噪声和混响环境下获得了良好的鲁棒性，在MCCC和MCSTP两种方法间获得了有效的折中。而且，正则化参数的不同取值使MCSTGSP方法构成一组具有不同预白化能力的时延估计器，可以根据用户的实际需要来进行调节。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种多通道群稀疏线性预测时延估计方法，其特征在于，包括以下步骤：
S1：M只传声器分别对声源信号进行持续采集，第m只传声器采集的时域信号记为Xm＝[xm(0),xm(1),…,xm(L-1)]，其中m＝1,2,…,M，xm(n)表示第m只传声器在时刻n的采集样本，n＝0,1,…,L-1，L表示每只传声器采集的样本数量；
S2：令时移序号d＝1，声源信号到达第1只和第2只传声器间的初始时移p1＝-pmax，pmax表示时延的最大可能值；
S3：分别将第m只传声器采集的信号Xm按照时延fm(pd)进行时移，fm(pd)表示声源信号到达第1只和第m只传声器间的相对时延，该时延是关于时移pd的函数；时移后的信号中每个样本记为xm(n,pd)，将M只传声器在时刻n的时移样本xm(n,pd)进行叠放，得到信号向量x(n,pd)；
S4：求解以下公式，得到预测系数矩阵A(pd)：
<mrow> <munder> <mi>min</mi> <mrow> <mi>A</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow> </munder> <mo>{</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>X</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>Y</mi> <mrow> <mo>(</mo> <mo>-</mo> <mn>1</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mi>A</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>A</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mo>|</mo> <msub> <mi>l</mi> <mrow> <mn>1</mn> <mo>,</mo> <mn>2</mn> </mrow> </msub> </msub> <mo>}</mo> </mrow>
其中，||·||F表示矩阵的F范数，λ是正则化参数，取值范围为λ＞0；
X(0,pd)＝[x(0,pd) x(1,pd) … x(K+L-1,pd)]T
Y(-1,pd)＝[y(-1,pd) y(0,pd) … y(K+L-2,pd)]T
其中：
x(w,pd)＝[x1(w,pd),x2(w,pd),…xM(w,pd)]T
<mrow> <mi>y</mi> <mrow> <mo>(</mo> <mi>w</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>x</mi> <mn>1</mn> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>w</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <msubsup> <mi>x</mi> <mn>2</mn> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>w</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <msubsup> <mi>x</mi> <mi>M</mi> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>w</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> </mrow>
xm(w-1,pd)＝[xm(w-1,pd),xm(w-2,pd),…xm(w-K,pd)]T
w＝0,1,…,K+L-1；K表示预测器阶数，其取值范围为K＜L；对于xm(q,pd)，如果q＜0或q＞L-1，则令xm(q,pd)＝0；
表示矩阵的l2范数，Ai(pd)表示预测系数矩阵A(pd)的第i行；
S5：计算预测误差矩阵E(0,pd)：
E(0,pd)＝X(0,pd)-Y(-1,pd)A(pd)
S6：计算预测误差相关矩阵R(pd)：
<mrow> <mi>R</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>K</mi> <mo>+</mo> <mi>L</mi> </mrow> </mfrac> <msup> <mi>E</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mi>E</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow>
S7：计算多通道互相关系数ρ(pd)：
<mrow> <msup> <mi>&amp;rho;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <mi>det</mi> <mrow> <mo>(</mo> <mi>R</mi> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msubsup> <mo>&amp;Pi;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msub> <mi>r</mi> <mrow> <mi>m</mi> <mo>,</mo> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
其中det(·)表示方阵的行列式，rm,m(pd)是矩阵R(pd)的第m个对角元素；
S8：如果pd＜pmax，令pd＝pd+1，返回步骤S3，否则根据以下公式求得时延估计值
<mrow> <mover> <mi>&amp;tau;</mi> <mo>^</mo> </mover> <mo>=</mo> <mi>arg</mi> <munder> <mi>max</mi> <msub> <mi>p</mi> <mi>d</mi> </msub> </munder> <msup> <mi>&amp;rho;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>.</mo> </mrow> 1
2.根据权利要求1所述的时延估计方法，其特征在于，所述步骤S4中预测系数矩阵的求解方法为：
S4.1：令迭代次数k＝1，初始化大小为KM×M的辅助矩阵Z1(pd)和拉格朗日乘子矩阵θ1(pd)；
S4.2：计算预测系数矩阵：
Ak+1(pd)＝[YT(-1,pd)Y(-1,pd)+βI]-1
×[YT(-1,pd)X(0,pd)+βZk(pd)-θk(pd)]
S4.3：更新辅助矩阵：
<mrow> <msubsup> <mi>Z</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>s</mi> <mi>o</mi> <mi>f</mi> <mi>t</mi> <mrow> <mo>(</mo> <msubsup> <mi>A</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> <mo>+</mo> <msubsup> <mi>&amp;theta;</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>(</mo> <msub> <mi>p</mi> <mi>d</mi> </msub> <mo>)</mo> <mo>/</mo> <mi>&amp;beta;</mi> <mo>,</mo> <mi>&amp;lambda;</mi> <mo>/</mo> <mi>&amp;beta;</mi> <mo>)</mo> </mrow> </mrow>
其中表示辅助矩阵Zk+1(pd)的第i行，i＝1,2,…,MK；soft函数定义为：
S4.4：更新拉格朗日乘子矩阵：
θk+1(pd)＝θk(pd)+β(Ak+1(pd)-Zk+1(pd))
S4.5：如果k＜Q，Q表示最大迭代次数，令k＝k+1，返回步骤S4.12，否则令预测系数矩阵A(pd)＝Ak+1(pd)。
3.根据权利要求1所述的时延估计方法，其特征在于，所述正则化参数λ根据以下公式计算：
<mrow> <mi>&amp;lambda;</mi> <mo>=</mo> <mi>&amp;delta;</mi> <mo>|</mo> <mo>|</mo> <msup> <mi>X</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mi>p</mi> <mo>)</mo> </mrow> <mi>Y</mi> <mrow> <mo>(</mo> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mo>|</mo> <msub> <mi>l</mi> <mi>&amp;infin;</mi> </msub> </msub> </mrow>
其中，表示求取矩阵的无穷范数，δ是一个正常数。
CN201610142030.5A 2016-03-14 2016-03-14 多通道群稀疏线性预测时延估计方法 Expired - Fee Related CN105652243B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610142030.5A CN105652243B (zh) 2016-03-14 2016-03-14 多通道群稀疏线性预测时延估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610142030.5A CN105652243B (zh) 2016-03-14 2016-03-14 多通道群稀疏线性预测时延估计方法

Publications (2)

Publication Number Publication Date
CN105652243A CN105652243A (zh) 2016-06-08
CN105652243B true CN105652243B (zh) 2017-12-05

Family

ID=56493323

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610142030.5A Expired - Fee Related CN105652243B (zh) 2016-03-14 2016-03-14 多通道群稀疏线性预测时延估计方法

Country Status (1)

Country Link
CN (1) CN105652243B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106249204B (zh) * 2016-09-09 2018-05-29 西南科技大学 基于鲁棒自适应盲系统辨识的多通道时延估计方法
CN108418718B (zh) * 2018-03-06 2020-07-10 曲阜师范大学 一种基于边缘计算的数据处理延迟优化方法及系统
CN110187306B (zh) * 2019-04-16 2021-01-08 浙江大学 一种应用于复杂室内空间的tdoa-pdr-map融合定位方法
CN113655440B (zh) * 2021-08-09 2023-05-30 西南科技大学 一种自适应折中预白化的声源定位方法
CN113655441B (zh) * 2021-08-11 2023-05-30 西南科技大学 一种低复杂度折中预白化的鲁棒声源定位方法
CN115825870B (zh) * 2023-02-17 2023-05-12 北京理工大学 基于群稀疏的离网格压缩匹配场处理声源定位方法
CN117825898B (zh) * 2024-03-04 2024-06-11 国网浙江省电力有限公司电力科学研究院 一种gis分布式振声联合监测方法、装置及介质

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101867853B (zh) * 2010-06-08 2014-11-05 中兴通讯股份有限公司 基于传声器阵列的语音信号处理方法及装置
CN102131288B (zh) * 2011-03-30 2013-06-19 北京交通大学 基于超分辨率uwb信号传播时延估计的室内定位方法及系统

Also Published As

Publication number Publication date
CN105652243A (zh) 2016-06-08

Similar Documents

Publication Publication Date Title
CN105652243B (zh) 多通道群稀疏线性预测时延估计方法
CN106251877B (zh) 语音声源方向估计方法及装置
CN107039045B (zh) 用于语音增强的全局优化最小二乘后滤波
US8098842B2 (en) Enhanced beamforming for arrays of directional microphones
Benesty et al. Time-delay estimation via linear interpolation and cross correlation
US8565446B1 (en) Estimating direction of arrival from plural microphones
US7626889B2 (en) Sensor array post-filter for tracking spatial distributions of signals and noise
CN110085247B (zh) 一种针对复杂噪声环境的双麦克风降噪方法
CN111044973A (zh) 一种用于麦克风方阵的mvdr目标声源定向拾音方法
CN109884591B (zh) 一种基于麦克风阵列的多旋翼无人机声信号增强方法
EP1899954A1 (en) System and method for extracting acoustic signals from signals emitted by a plurality of sources
Li et al. Reverberant sound localization with a robot head based on direct-path relative transfer function
CN111239686B (zh) 一种基于深度学习的双通道声源定位方法
CN112363112B (zh) 一种基于线性麦克风阵列的声源定位方法及装置
CN113903353B (zh) 一种基于空间区分性检测的定向噪声消除方法及装置
Dietzen et al. Joint multi-microphone speech dereverberation and noise reduction using integrated sidelobe cancellation and linear prediction
Huang et al. Time delay estimation and source localization
CN110111802A (zh) 基于卡尔曼滤波的自适应去混响方法
KR100940629B1 (ko) 잡음 제거 장치 및 방법
Kumatani et al. Microphone array post-filter based on spatially-correlated noise measurements for distant speech recognition
Šarić et al. Supervised speech separation combined with adaptive beamforming
Triki et al. Delay and predict equalization for blind speech dereverberation
CN111157949A (zh) 一种语音识别及声源定位方法
CN106249204B (zh) 基于鲁棒自适应盲系统辨识的多通道时延估计方法
CN113655440A (zh) 一种自适应折中预白化的声源定位方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171205