“鸡尾酒会效应”^[1]告诉我们，人耳可以在嘈杂的环境下，从多个声源中选取并定位跟踪自己感兴趣的声源，人耳的这一能力很大程度上取决于“双耳效应”^[2]。但是对于单侧耳上佩戴电子耳蜗^[3]或助听器^[4]的听力障碍者，其感知声源方向的能力非常弱，因此需要通过安装在其上的近距离双麦克风对声源方向进行分析^[5-6]。

“双耳效应”的产生源于声音到达双耳时的耳间时间差(interaural time difference，ITD)以及耳间声级差(interaural intensity difference，IID)。在一般的研究中，信号的低频部分用时间延迟差(time difference，TD)进行分析，但对于信号的高频部分，用TD进行估计时会由于相位的多值性产生检测不准确，因此高频信号一般用能量衰减比(intensity difference，ID)进行分析^[7]。然而对于助听器或电子耳蜗体外机，安装在其上的双麦克风间的距离远小于双耳间的距离，因此在语音的频带范围内计算TD时不会产生相位估计的多值结果，故本文利用TD来估计声源方向而不需要考虑信号到达双麦克风的能量衰减比。

在通过TD值计算声源方向的研究中，文^[8]通过TD和IID的联合估计，在无混响和低混响情形下得到了对声源的准确定位，该算法的优点在于可以实时处理，即可以估计运动声源，但是该算法的双麦克风间距为25 cm，对于近距离双麦克风，该算法的估计会变得不准确。文^[9]提出了一种适用于近距离双麦克风的角度估计方法，该方法结合语音帧判断和功率谱估计技术，可对声源的方向参数进行实时估计，并且该方法对各频带进行了归一化波束调整，使其更适用于电子耳蜗植入者，但是该方法只能估计2个声源，并假设其中一个声源位于双麦克风正前方，而对于更多的声源，该方法失去有效性。文^[10]基于语音信号稀疏性的假设，实现了对多个声源的估计定位，但是该方法不具有实时性，即不能对运动声源的轨迹进行估计。文^[11]通过计算过零点时间差(zero-crossing time difference，ZCTD)并结合信噪比估计，在含有白噪声的低信噪比环境下得到了较好的TD估计值。文^[12]对文^[11]的算法作了进一步的改进，提出利用Teager算子提取语音包络的方法，从而选取有效的语音段估计TD，免去了文^[11]中估计信噪比的繁琐步骤，缩短了估计TD的时间。但是文^[11-12]的方法适用的情形均为声源在较长时间内(不小于1 s)静止不动，而在实际的听觉场景中，当多个声源不断运动时，上述算法估计的准确性会下降。

为了解决上述算法中存在的问题，本文提出一种基于近距离双麦克风的多个运动声源方向的估计方法，并分别对多声源静止、多声源运动的情形做了实验分析，结果均证明了算法的有效性。

1 算法介绍 1.1 算法基本流程

为了提高角度估计的准确性，本算法在以下4个方面作了改进： 首先，将传统的Gammatone滤波器^[13]进行改进，以增加麦克风1和2经过滤波器输出得到的ZCTD的数据量，减小ZCTD的计算标准差； 其次，在计算ZCTD时，分别利用上升ZCTD和下降ZCTD对结果的准确性进行判别，舍去会引起结果估计不准确的时频单元； 再次，利用得到的ZCTD的均值和标准差，生成一个更为准确地反映角度分布区域可能性的Gauss函数，并对所有这样的Gauss函数进行带权重的累加； 最后，经过平滑处理，去掉累加函数中无效的峰值，使得累加函数的峰值更可靠，并通过峰值检测得到准确的角度估计值。

算法的基本流程图如图 1所示。算法流程主要分为： 改进的Gammatone滤波器、过零点检测和ZCTD计算、 Gauss函数累加、平滑处理、峰值检测。

1.2 改进的Gammatone滤波器

将语音信号进行时间分块处理，每 0.05 s 为一个时间块。语音信号的采样率设置为 48 kHz，因此每个时间块有2 400个采样点。块与块之间的交叠为50%时间块长度，时间块信号加汉明窗。对于每个时间块内的数据，用Gammatone滤波器进行频带划分，因此每个频带内的信号近似为一个幅度调制的正弦。滤波器频带范围为80~5 000 Hz，通带个数为128。每个通带内滤波器的冲击响应如下：

${{g}_{{{f}_{c}}}}\left( t \right)={{t}^{N-1}}\text{exp}[-2\pi tb({{f}_{c}})]\text{cos}(2\pi {{f}_{c}}t+\phi )u\left( t \right).$

(1)

其中：N为滤波器阶数，本文设置为4阶；f_c为频带的中心频率；$\phi $为相位，由于人耳对相位信息不敏感，因此$\phi $设置为0；b为衰减速度， b的计算在文^[14]中通过等效矩形带宽(ERB)的计算得出，如下：

$b({{f}_{c}})=1.019\text{ERB}\left( f \right)=1.019\times (24.7+0.108{{f}_{c}}).$

(2)

在后续的处理过程中，本文希望对每个时间块内每个通带中的信号的所有过零点进行检测，并进行统计分析，信号内的周期数越多，这样的统计便越准确。但是由于每个时间块内信号时长只有 0.05 s，因此对于低频的通带，该通带每个时间块内的输出信号包含的周期数较少，这样会在后续计算ZCTD的均值和标准差时产生不精确的结果。因此对于Gammatone滤波器的低频通带，应尽量选取更多的数据点以获取更准确的过零点序列的统计值，即计算此时间块的Gammatone滤波器输出时将此时间块信号作扩展，增加时间块长度。但是过大的时间块长度同样会带来结果的不准确，这是由于时间块长度过大破坏了时间块长度选取的短时性原则，即意味着当前时间块内声源的位置可能发生的改变过大，这会造成最终估计声源角度时的不准确，因此对于低频通带，应同时注意上述分析，选取合适的时间块长度。

基于上述分析，将传统的Gammatone滤波器稍作改进，即将时间块长度进行调整，调整方式如下：

$\text{帧长=}\left\{ \begin{array}{*{35}{l}} \frac{20}{{{f}_{\text{c}}}\left( i \right)\times 0.05}\times 2400, & 1\le i\le 15; \\ \frac{40}{{{f}_{\text{c}}}\left( i \right)\times 0.05}\times 2400, & 16\le i\le 55; \\ 2400, & 56\le i\le 128, \\ \end{array} \right.$

(3)

其中： f_c(i)为通带的中心频率; i为通带序数。而对于任意的时间块，低频部分和高频部分相比，由于其时间块长度变长，因此其块与块之间的交叠也会相应变长，故低频部分与高频部分的时间差不会随着时间推移而变大。

1.3 过零点检测和ZCTD计算

对于麦克风1和2采集到语音信号每个时间块内的每个通带滤波器输出，分别检测其相应的过零点并计算ZCTD。将上述改进后的Gammatone滤波器的输出记为s_i,j(k)，其中i表示时间块序数，j表示通带序数，k表示通带输出的采样点序数。计算所有使得s_i,j正负性发生改变的k值，并将k值按照信号从正到负和从负到正2种情况分别定义为上升过零点和下降过零点。因为采样点很有可能无法正好采样到零点，此时无法准确地得到k值，所以应该对零点附近的点进行插值，而由于正弦函数在零点附近的二阶导数值也为0，即正弦函数在零点附近斜率的变化很小，因此用线性插值便可快速地计算出准确的k值。具体的上升过零点和下降过零点的计算方式如下：

$\begin{align} & \text{上升过零点=}k-\frac{{{s}_{i,j}}\left( k \right)}{{{s}_{i,j}}\left( k+1 \right)-{{s}_{i,j}}\left( k \right)}|{{s}_{i,j}}\left( k \right)＜0\text{且}{{s}_{i,j}}\left( k+1 \right)>0; \\ & \text{下降过零点=}k+\frac{{{s}_{i,j}}\left( k \right)}{{{s}_{i,j}}\left( k \right)-{{s}_{i,j}}\left( k+1 \right)}|{{s}_{i,j}}\left( k \right)>0\text{且}{{s}_{i,j}}\left( k+1 \right)＜0. \\ \end{align}$

(4)

由于上述的操作对麦克风1和2是同时进行的，因此根据式(4)可以分别得到2组上升过零点和下降过零点序列，并将这些序列记为P_up1、 P_up2、 P_down1、 P_down2，其中脚标1和2分别对应麦克风1和2。 对P_up1中的每个数，在P_up2中分别寻找和P_up1最接近的数，并将结果记为P_up2，通过式(5)计算2个向量的差向量ΔP_up，并按照同样的方式求出ΔP_down，ΔP_up和ΔP_down代表通过上升ZCTD序列和下降ZCTD序列。

$\Delta {{P}_{up}}={{P}_{up1}}-{{{\bar{P}}}_{up2}}$

(5)

然后分别计算ΔP_up和ΔP_down内所有数的均值和标准差，并分别记作E_up、 σ_up、 E_down、 σ_down。 接着对求出的ZCTD的均值进行准确与否的判别，如果E_up落在区间[E_down－σ_down,E_down+σ_down]内，同时E_down落在区间[E_up－σ_up,E_up+σ_up]内，则表示ZCTD序列检测准确; 反之，检测错误，停止此时频单元的计算，并开始计算下一时频单元。

若判断出ZCTD序列正确，那么将上升ZCTD序列ΔP_up和下降ZCTD序列ΔP_down进行合并，得到ΔP_up-down，并计算其均值E_up-down和标准差σ_up-down。 通过上述过程，便求出了麦克风1和2采集到的语音信号在每个时间块内的每个通带对应的时间延迟差。

1.4 Gauss函数累加

为了更准确地对当前时频单元内时延点数的值进行判断，应该用计算出的ZCTD序列的均值和标准差共同确定其范围。真实的时延点数靠近均值的可能性最大，而距离均值越远，可能性越小。此外，较小的标准差和较大的标准差相比时，真实的时延点数在均值处的可能性更大，并且所有可能的时延点数分布的概率密度函数的积分应该为一常数。由于Gauss函数正好满足上述特点，因此本文用上面得到的每个时间块内的每个通带对应的时间延迟差ZCTD的均值E_up-down和标准差σ_up-down生成一个Gauss函数g_i,j(x)，如式(6)所示，在式(6)中，分别用E和σ代表E_up-down和σ_up-down，

${{g}_{i,j}}\left( x \right)=\left\{ \begin{array}{*{35}{l}} \frac{1}{\sigma }\text{exp}\left[ \frac{-{{\left( x-E \right)}^{2}}}{2{{\sigma }^{2}}} \right] & E-3\sigma \le x\le E+3\sigma ; \\ 0 & x＜E-3\sigma \text{或}x>E+3\sigma . \\ \end{array} \right.$

(6)

其中： i表示时间块序数； j表示通带序数； x代表时延点数(可以为小数)。由于语音信号有语音段也有静音段，因此对于当前选中的时频区域，若此区域为语音区，则进行式(6)的Gauss函数累加，并乘以相应的权重系数，语音越强，权重系数越大； 若此区域为背景噪声区或静音区，则Gauss函数累加时的权重很小或不进行累加。这样更能反映声源在语音段的方位信息，而不会将噪声段或静音段计算在内。因此对于每个固定的时间块内(如第i₀块)的信号，计算该时间块内所有通带的Gauss函数，并将这些Gauss函数进行带权重的求和累加，累加权重为2个麦克风在当前时频区域能量的乘积，累加后的序列记为h_i0(x)，如式(7)所示，h_i0(x)反映了当前时间块内时间延迟差的分布情况，

${{h}_{{{i}_{0}}}}\left( x \right)=\sum\limits_{j}{|MI{{C}_{1}}({{i}_{0}},j)\cdot MI{{C}_{2}}({{i}_{0}},j)|\cdot {{g}_{{{i}_{0}},j}}\left( x \right)}.$

(7)

其中，MIC₁和MIC₂表示2个麦克风的短时傅里叶变换。为了更直观地体现声源的角度，通过式(8)，将式(7)中的h_i0(x)的自变量按照时延点数和角度之间的对应关系映射为h_i0(θ)，

$\theta =\arccos \left( \frac{x\cdot c}{d\cdot {{f}_{\text{s}}}} \right).$

(8)

其中： θ代表声源的角度； c为声速； d为双麦克风之间的距离； f_s为采样频率。 h_i0(θ)函数的峰值即可反映当前时间块(第i₀块)内声源的方向信息，其峰值个数代表了当前时间块内正在发声的声源的个数，每个峰所对应的位置代表了该声源所在的角度值。

1.5 平滑处理

对上述得到的当前时间块内的角度累积函数h_i0(θ)进行平滑，以便进行后续的峰值检测。先求出h_i0(θ)所有的极大值点，然后将所有极大值点进行3次样条插值。将插值得到的平滑后函数记为g_i0(θ)，g_i0(θ)即为当前时间块对h_i0(θ)平滑后的结果。

1.6 峰值检测

为了最终输出角度估计结果，对上面得到的平滑后的角度累计函数g_i0(θ)进行峰值检测。求出g_i0(θ)中的最大值g_i0(θ)_max和极大值序列 g_i0(θ)_max(l)以及这些极大值附近的点，其中l为极大值序数。极大值附近的点的判别标准为： 若某点到某极大值点之间的序列为单调序列，并且该点的值不小于极大值的95%，则该点到极大值之间的所有点均视为极大值附近的点。保留极大值及其附近的点中不小于最大值1/4的部分，并将保留的极大值记为peak_i0(θ)，peak_i0(θ)的计算式如下：

$\text{pea}{{\text{k}}_{{{i}_{0}}}}\left( \theta \right)=\left\{ \begin{array}{*{35}{l}} 1 & \overline{{{g}_{{{i}_{0}}}}{{\left( \theta \right)}_{max}}}\left( l \right)\ge \frac{1}{4}{{g}_{{{i}_{0}}}}{{\left( \theta \right)}_{max}}; \\ 0 & \overline{{{g}_{{{i}_{0}}}}{{\left( \theta \right)}_{max}}}\left( l \right)＜\frac{1}{4}{{g}_{{{i}_{0}}}}{{\left( \theta \right)}_{max}}. \\ \end{array} \right.$

(9)

通过式(9)求出peak_i0(θ)曲线的峰值个数和每个峰所在位置，便可确定当前时间块内发声声源个数以及这些声源各自所处的方向。方向一旦得到，则当前时间块内的所有计算结束，转向下一时间块的分析计算，计算方式和当前时间块的计算方式一致，重复此过程，直到求出所有时间块内的峰值函数peak_i(θ)，并将这些峰值函数peak_i(θ)按照时间块的顺序合并为二维函数peak(i,θ)，peak(i,θ)的二维灰度图即可反映声源角度随时间变化的结果。

2 实验处理结果 2.1 实验基本设置

本文设计实验来验证算法对于角度估计的有效性。根据助听器和电子耳蜗体外机的实际尺寸，实验中，双麦克风距离设置为15 mm，实验环境为无混响环境。实验中选取的语音来源于Neospeech语音库(http://www.neospeech.com)，分别为低频男声Liang、 女声Lily和高频男声Wang朗读的3段语音，3段语音的语义见表 1。

2.2 对多个静止声源角度的估计结果

为了测试算法对于多个静止声源角度估计的效果，在第1次测试中，让3个不同的朗读者发出的声音固定在3个不同的方位，分别为45°、 80°和110°，同时进行循环播放，声音的持续时间为30 s。 为了进一步验证算法是否能区分较为接近的角度，在第2次测试中，让3个不同的朗读者发出的声音固定的方位分别为89°、 90°和91°，声音的持续时间仍为30 s。

对多个静止声源角度的估计结果如图 2所示，具体的角度估计值如表 2所示。

从图 2可以看出，算法估计出的多声源角度和真实的角度几乎一致，而当多声源角度非常接近时，这一估计结果仍然和真实角度保持一致。从表 2可以看出，在这2次测试中，估计值与真实值的差距均不超过0.4°，误差值较小。通过实验证明： 算法对于多个静止声源角度的估计结果较为准确。

2.3 对多个运动声源角度的估计结果

为了测试算法对于多个运动声源角度轨迹估计的效果，实验中让3个不同的朗读者发出的声音同时循环播放，并同时进行移动。在30 s的时间内，播放朗读者Liang的声源的位置从60°均匀地移动到120°，Lily的位置从70°均匀地移动到130°，Wang的位置从160°均匀的移动到40°。 对多个运动声源角度的估计结果如图 3所示。

从图 3可以看出，实验用到的3个声源中既有较快速运动的声源，又有较慢速运动的声源； 既有运动轨迹完全不同的声源，又有运动轨迹比较接近的声源。对比图 3b—3d可以看出： 经过平滑处理和峰值提取，角度的轨迹线去掉了一些干扰点，变得更加准确。对比图 3a和3d可以看出： 算法对于混合在一起的不同类型声源的运动方向轨迹的检测能力都是较好的。

3 分析与讨论

本文之所以可以通过近距离双麦克风对多运动声源的方向进行准确的估计，是由于在以下4个方面的创新研究： 1) 算法根据通带大小对时间块长度的调节，使得低频通带ZCTD的数据量大大增加，从而得到了更准确的ZCTD的统计值。 2) 算法联合上升ZCTD和下降ZCTD，对ZCTD统计值的准确性做判别，使得算法可以识别出错误的估计结果。 3) 算法利用Gauss函数的特性，更准确地体现了角度分布在各个区域的可能性。 4) 算法通过平滑处理和峰值检测，进一步得到了更为准确的角度估计值。

本文提到在探究声源定位时，低频声音主要用TD进行分析，而高频声音在用TD分析的时候会由于相位的多值性而产生计算的不准确。对于不会使TD的估计发生多值性的频率范围，文^[11]中给出的频率上限为625 Hz。 而在本文中，为了适合助听器和电子耳蜗体外机，将麦克风间距设置为 15 mm，因此按照上面的分析方法进行计算，不会发生多值结果的频率上限为11.3 kHz，远大于语音的频带上限，因此本文的算法不会因为ZCTD计算的多值性而产生计算误差。

另外，在对传统的Gammatone滤波器进行改进的时候，节1.2提及了一个注意点，即改进后的时间块长度不能过大。在式(3)中，对于每一种情况，时间块长度均不会随着频率增加而变大，因此只需要知道式(3)的3种情形中最低通带的时间块长度即可，其他通带的时间块长度一定小于最低通带的时间块长度。在3种情形的最低通带中，通带1在扩展后的时间块长度为0.25 s，通带16扩展后的时间块长度为0.198 s，而通带56不进行扩展，时间块长度仍为0.05 s，因此不会出现因扩张后时间块长度过长而导致角度估计不准确的情况。

4 结 论

本文针对电子耳蜗或助听器上双麦克风距离小、声源不断运动所导致的声源方向估计不准的问题，提出了一种改进的Gammatone滤波器分析方法和一种准确的过零点检测判别方法，并利用Gauss函数进行角度的累加估计。通过实验证明了在麦克风距离很近的情况下，对于多个静止声源，在其方向远离或接近时均能得到很好的角度估计值； 对于多个运动声源，在其运动轨迹相似或完全不同时也均能得到很好的角度估计轨迹。但是，本文提出的算法的计算量和传统算法相比较为复杂，每个时间块的处理时间约为时间块长度的1.5倍，因此在减小算法复杂度、 提高程序运行效率方面，本文还需做进一步努力。