CN104076331B - 一种七元麦克风阵列的声源定位方法 - Google Patents

Abstract

本发明公开了一种七元麦克风阵列的声源定位方法，该方法建立了七元麦克风阵列模型与室内声场模型；利用该麦克风阵列接收室内声源语音信号，对接收的语音信号进行采集、滤波、加窗、谱减、倒谱处理；再利用相位变换加权广义互相关方法GCC-PHAT，得各麦克风之间的时间延迟；运用三维空间定位方法与坐标旋转数字式计算机方法，确定声源位置。本发明在混响、噪声等干扰信号存在的室内环境下，能准确确定声源位置，实时性强、准确性高，在语音信号处理领域有广泛的应用。

Description

一种七元麦克风阵列的声源定位方法

技术领域

本发明涉及声源定位领域，特别是一种七元麦克风阵列的声源定位方法。

背景技术

声源定位技术是阵列信号处理的重要技术之一。目前在声呐探测、视频电话会议、人工智能、语音追踪与识别、监控系统等多个领域有着广泛的应用。利用麦克风阵列对声源方位进行计算是声源定位的基本方法，它是将一组麦克风传感器按一定方式布置在空间不同位置上,形成麦克风阵列；利用麦克风阵列接收空间声源信号，再对阵列接收的信号进行处理，提取信号的有用特征，再通过一定计算方法得到声源的方位信息。

目前常用的声源定位方法主要有三类：基于可控波束形成器的声源定位法，主要是将各阵元采集到的信号加权求和，通过调控权值(权值取决于阵元信号的相位延迟，主要利用波达方向(DOA)(见文献：DespoinaPavlidi,AnthonyGriffin,MatthieuPuigt,AthanasiosMouchtaris.Real-timemultiplesoundsourcelocalizationandcountingusingacircularmicrophonearray[J].IEEETransactionsonAudio,Speech,andLanguageProcessing,2013,21(10):2193-2206.))等方法进行调控)，使阵列输出信号功率最大，从而进行声源定位；基于高分辨率谱估计的声源定位法，主要利用最小方差估计法、子空间法(如MUSIC法(见文献：JunliLiang,DingLiu.Passivelocalizationofmixednear-fieldandfar-fieldsourcesusingtwo-stageMUSICalgorithm[J].IEEETransactionsonSignalProcessing,2010,58(1):108-120.)、ESPRIT法(见文献：HaohaiSun,HeinzTeutsch,EdwinMabande,WalterKellermann.RoubstlocalizationofmultiplesourcesinreverberantenvironmentsusingEB-ESPRITwithsphericalmicrophonearrays[C].IEEEInternationalConferenceonAcousticsSpeechandSignalProcessing(ICASSP),Prague,2011:117-120.))等进行声源定位；基于到达时间差(TDOA)的声源定位法，主要结合自适应、广义互相关(GCC)等进行声源定位。利用多个麦克风构成一定形状的阵列，可在二维平面、三维立体空间中获得更为准确的声源目标方位。

现有定位技术中，主要是利用一维线性阵列、二维十字四麦克风阵列进行声源定位，与三维的阵列相比，这些阵列不能同时采集多个方向上的声源信号，难以精确定位。另外，现有定位技术中，仅考虑了环境噪声对定位精度的影响，而没有考虑混响的影响。因此，在噪声与混响同时干扰下，现有的声源定位方法是不能准确地进行声源定位的。

如何解决现有技术的不足已成为声源定位领域亟待解决的一大难题。

发明内容

本发明所要解决的技术问题是为了克服现有技术的不足，而提供一种七元麦克风阵列的声源定位方法。本发明方法在三维空间中实时、有效地解决了声源定位问题，在一定程度上降低了噪声、混响等信号对语音信号的干扰，提高了定位的精度与稳定性。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种七元麦克风阵列的声源定位方法，包括如下步骤：

步骤一、在三维坐标系中，建立麦克风阵列模型，确定各麦克风在三维空间的位置；

步骤二、对各麦克风接收到声源发出的语音信号s(t)进行采样，得到离散信号s(k)，t表示时间，k表示时间序列；

步骤三、对步骤二所述的离散信号s(k)进行滤波、加窗、谱减、倒谱处理，得到纯语音时域信号x(k)；具体如下：

(301)、将离散信号s(k)通过带通滤波器进行带通滤波后，得到带通信号；

(302)、对(301)中所得带通信号用窗函数进行加窗处理，得到加窗信号s_w(k)；

(303)、对(302)所得加窗信号s_w(k)用谱减法进行去噪处理，得到谱减信号X_w(f)，对X_w(f)作傅里叶反变换得到去噪时域信号x_w(k)；

(304)、对去噪时域信号x_w(k)由倒谱法进行消除混响处理，得到纯语音时域信号x(k)；

步骤四、对步骤三中的纯语音时域信号x(k)，利用相位变换加权广义互相关方法GCC-PHAT计算第i个麦克风接收的经步骤三处理所得纯语音时域信号x_i(k)与三维坐标原点处第0个麦克风接收的经步骤三处理所得纯语音时域信号x₀(k)间的时延值τ_i，1≤i≤6且i为自然数，则互相关函数为

$R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right)={\∫}_{-\∞}^{+\∞}{\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)X_0\left(f\right){X_i}^{*}\left(f\right)e^{{j2\mathrm{\πf\τ}}_i}\mathrm{df};$

式中，f表示频率，e表示自然指数，X₀(f)表示纯语音时域信号x₀(k)的傅里叶变换，X_i ^*(f)表示纯语音时域信号x_i(k)傅里叶变换的共轭，j为虚数单位，ψ_PHAT(f)为加权函数，且

${\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)=\frac{1}{\left|G_{x_0{x}_i}\left(f\right)\right|};$

式中，为相关函数的功率谱，且

$G_{x_0{x}_i}\left(f\right)={\∫}_{-\∞}^{+\∞}{R}_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)e^{{j2\mathrm{\πf\τ}}_i}{\mathrm{d\τ}}_i;$

式中， $R_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)=E\left[x_0\left(k\right)x_i(k+{\mathrm{\τ}}_i)\right],$ 表示x₀(k)与x_i(k)的互相关函数，则其中argmax表示使取得最大值时即为第0个麦克风与第i个麦克风所接收到信号的时延估计值；

步骤五、利用三维空间定位方法，结合步骤四中的得声源到三维坐标原点的距离r：

$r=\frac{{6L}^2-c^2\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\left({\widehat{\mathrm{\τ}}}_i\right)}^2}{2c\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\widehat{\mathrm{\τ}}}_i};$

式中，c表示声速，L表示第i个麦克风位置M_i与三维坐标原点处第0个麦克风位置M₀间的距离；

步骤六、利用坐标旋转数字式计算机方法，结合步骤四得到的得声源位置的俯仰角θ、方位角为

步骤七、根据直角坐标与球面坐标间的关系，得声源位置的直角坐标为

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(302)中窗函数采用汉明窗函数w(k)，即

式中，N表示时间序列的取值范围，且N为大于零的正整数。

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(303)中谱减法的公式为

$\left|X_w\left(f\right)\right|=\sqrt{{\left|S_w\left(f\right)\right|}^2-{\left|N_w\left(f\right)\right|}^2};$

式中，S_w(f)为加窗信号s_w(k)的傅里叶变换，N_w(f)为加窗信号s_w(k)中噪声信号n_w(k)的傅里叶变换。

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(304)中倒谱法处理过程为：

(a)、对去噪时域信号x_w(k)进行分帧处理，得到分帧信号；

(b)、对分帧信号用窗函数进行预加窗处理，得到预加窗信号x_ww(k)；

(c)、对预加窗信号x_ww(k)计算每帧的复倒谱信号x_c(k)，所述复倒谱计算公式为

x_c(k)＝IFFT{ln[FFT[x_ww(k)]]}；

式中，FFT表示傅里叶变换，IFFT表示傅里叶反变换，ln表示自然对数；

(d)、对复倒谱信号x_c(k)，利用复倒谱域低通滤波器进行滤波，得到滤波信号；

(e)、对滤波信号，在时域条件下将每帧信号叠接相加，得到纯语音时域信号x(k)。

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(301)中带通滤波器的上截止频率为f_H＝3500Hz，下截止频率为f_L＝250Hz。

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(a)中分帧信号的帧长为512点，帧移为128点。

作为本发明的一种七元麦克风阵列的声源定位方法进一步的优化方案，所述(d)中复倒谱域低通滤波器由通带、过渡带、阻带三部分组成，其中通带的最高截止点P＝4,过渡带的带宽h为64。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：本发明建立了七元麦克风阵列模型与室内声场模型，在室内声场中采集声源发出的一段语音信号；对采集的语音信号进行滤波、加窗、谱减、倒谱处理，利用相位变换加权广义互相关方法GCC-PHAT计算各麦克风接收信号间的时间延迟；再由三维空间定位方法与坐标旋转数字式计算机方法，确定声源方位角，最终确定声源位置。本发明在混响与噪声干扰都存在的环境下，准确地确定了声源位置，定位精度高且稳定性强，在语音信号处理领域，具有较强的实用性。

附图说明

图1是本发明的流程图。

图2是本发明设计的七元麦克风阵列模型。

图3是本发明中倒谱法的步骤。

图4是本发明中复倒谱域低通滤波器原理图。

图5是本发明中麦克风阵列模型：(a)七元麦克风阵列，(b)四元麦克风阵列。

图6是本发明麦克风间距与方位角、俯仰角的关系图：(a)麦克风间距为0.1m时与方位角的关系图，(b)麦克风间距为0.1m时与俯仰角的关系图。

图7是本发明麦克风间距与方位角、俯仰角的关系图：(a)麦克风间距为0.2m时与方位角的关系图，(b)麦克风间距为0.2m时与俯仰角的关系图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

一种七元麦克风阵列的声源定位方法，包括如下步骤：

步骤一、在三维坐标系中，建立麦克风阵列模型，确定各麦克风在三维空间的位置；

步骤二、对各麦克风接收到声源发出的语音信号s(t)进行采样，得到离散信号s(k)，t表示时间，k表示时间序列；

步骤三、对步骤二所述的离散信号s(k)进行滤波、加窗、谱减、倒谱处理，得到纯语音时域信号x(k)；具体如下：

(301)、将离散信号s(k)通过带通滤波器进行带通滤波后，得到带通信号；

(302)、对(301)中所得带通信号用窗函数进行加窗处理，得到加窗信号s_w(k)；

(303)、对(302)所得加窗信号s_w(k)用谱减法进行去噪处理，得到谱减信号X_w(f)，对X_w(f)作傅里叶反变换得到去噪时域信号x_w(k)；

(304)、对去噪时域信号x_w(k)由倒谱法进行消除混响处理，得到纯语音时域信号x(k)；

步骤四、对步骤三中的纯语音时域信号x(k)，利用相位变换加权广义互相关方法GCC-PHAT计算第i个麦克风接收的经步骤三处理所得纯语音时域信号x_i(k)与三维坐标原点处第0个麦克风接收的经步骤三处理所得纯语音时域信号x₀(k)间的时延值τ_i，1≤i≤6且i为自然数，则互相关函数为

$R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right)={\∫}_{-\∞}^{+\∞}{\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)X_0\left(f\right){X_i}^{*}\left(f\right)e^{{j2\mathrm{\πf\τ}}_i}\mathrm{df};$

式中，f表示频率，e表示自然指数，X₀(f)表示纯语音时域信号x₀(k)的傅里叶变换，X_i ^*(f)表示纯语音时域信号x_i(k)傅里叶变换的共轭，j为虚数单位，ψ_PHAT(f)为加权函数，且

${\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)=\frac{1}{\left|G_{x_0{x}_i}\left(f\right)\right|};$

式中，为相关函数的功率谱，且

$G_{x_0{x}_i}\left(f\right)={\∫}_{-\∞}^{+\∞}{R}_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)e^{{j2\mathrm{\πf\τ}}_i}{\mathrm{d\τ}}_i;$

式中， $R_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)=E\left[x_0\left(k\right)x_i(k+{\mathrm{\τ}}_i)\right],$ 表示x₀(k)与x_i(k)的互相关函数，则其中argmax表示使取得最大值时即为第0个麦克风与第i个麦克风所接收到信号的时延估计值；

步骤五、利用三维空间定位方法，结合步骤四中的得声源到三维坐标原点的距离r：

$r=\frac{{6L}^2-c^2\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\left({\widehat{\mathrm{\τ}}}_i\right)}^2}{2c\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\widehat{\mathrm{\τ}}}_i};$

式中，c表示声速，L表示第i个麦克风位置M_i与三维坐标原点处第0个麦克风位置M₀间的距离；

步骤六、利用坐标旋转数字式计算机方法，结合步骤四得到的得声源的俯仰角θ、方位角为

步骤七、根据直角坐标与球面坐标间的关系，得声源位置的直角坐标为

所述(302)中窗函数采用汉明窗函数w(k)，即

式中，N表示时间序列的取值范围，且N为大于零的正整数。

所述(303)中谱减法的公式为

$\left|X_w\left(f\right)\right|=\sqrt{{\left|S_w\left(f\right)\right|}^2-{\left|N_w\left(f\right)\right|}^2};$

式中，S_w(f)为加窗信号s_w(k)的傅里叶变换，N_w(f)为加窗信号s_w(k)中噪声信号n_w(k)的傅里叶变换。

所述(304)中倒谱法处理过程为：

(a)、对去噪时域信号x_w(k)进行分帧处理，得到分帧信号；

(b)、对分帧信号用窗函数进行预加窗处理，得到预加窗信号x_ww(k)；

(c)、对预加窗信号x_ww(k)计算每帧的复倒谱信号x_c(k)，所述复倒谱计算公式为

x_c(k)＝IFFT{ln[FFT[x_ww(k)]]}；

式中，FFT表示傅里叶变换，IFFT表示傅里叶反变换，ln表示自然对数；

(d)、对复倒谱信号x_c(k)，利用复倒谱域低通滤波器进行滤波，得到滤波信号；

(e)、对滤波信号，在时域条件下将每帧信号叠接相加，得到纯语音时域信号x(k)。

所述(301)中带通滤波器的上截止频率为f_H＝3500Hz，下截止频率为f_L＝250Hz；所述(a)中分帧信号的帧长为512点，帧移为128点；所述(d)中复倒谱域低通滤波器由通带、过渡带、阻带三部分组成，其中通带的最高截止点P＝4,过渡带的带宽h为64。

本方法的流程参见图1，本发明是的一种七元麦克风阵列的声源定位方法，利用三维七元麦克风阵列，结合空间定位方法进行声源定位，其具体实施步骤如下:

步骤一、在三维坐标系中，建立麦克风阵列模型，确定各麦克风在三维空间的位置；

建立麦克风阵列模型，如图2所示，坐标原点O(0,0,0)，七个麦克风的位置坐标分别为M₀(0,0,0)、M₁(L,0,0)，M₂(0,L,0)，M₃(-L,0,0)，M₄(0,-L,0)，M₅(0,0,-L)，M₆(0,0,L)；声源位置的直角坐标为S(x,y,z)，方位角为俯仰角为θ(0°≤θ≤90°)。声源到坐标原点(就是第0个麦克风的位置M₀)的距离为r，在远场条件下阵列接收的信号为平面波。d_i(1≤i≤6)为声源到第0个麦克风的位置M₀与声源到第i个麦克风的位置M_i的距离差值。

步骤二、对各麦克风接收到声源语音信号s(t)进行采样，得离散信号s(k)，t表示时间，k表示时间序列。为使仿真实验更加接近真实环境，所建立的房间模型中，房间的反射系数为0.5，声源为一段含噪声的语音信号，持续时间为8s，信号采样频率为16kHz，64kHz。

步骤三、对步骤二所述的离散信号s(k)进行：滤波、加窗、谱减、倒谱处理，处理过程为：

(301)、对步骤二中的离散信号s(k)通过带通滤波器进行带通滤波后，得到带通信号。滤波可以有效地抑制低频、高频信号对语音信号的污染。由于语音信号为有限带宽信号，一般情况下其频率主要集中在300Hz～3000Hz之间，因此我们选择带通滤波器，带通滤波器的上截止频率为f_H＝3500Hz，下截止频率为f_L＝250Hz。

(302)、对(301)所得带通信号用窗函数进行加窗处理，得到加窗信号s_w(k)。因为语音信号一般在10ms到30ms之间，是平稳的。为了更好的处理语音信号，对语音信号进行加窗处理，所得加窗信号s_w(k)有效地减少了频谱泄漏，保持语音信号的平稳性；所述窗函数w(k)采用汉明窗函数(在汉明窗函数中，其主瓣加宽并降低，使得信号的泄漏与波动都减小，汉明窗优于矩形窗；其加权系数不同，使得旁瓣达到更小，旁瓣衰减速度减慢，汉明窗优于汉宁窗；其在去除混响效果中，比矩形窗、汉宁窗、布莱克曼窗的效果更为稳定、均衡、优异。故选择汉明窗函数)，即

式中，N表示时间序列的取值范围，且N为大于零的正整数。

(303)、对(302)所得加窗信号s_w(k)用谱减法进行去噪处理，得到的谱减信号X_w(f)最大程度地消除现场噪音对语音信号的干扰，使得分离出来的语音信号纯净、稳定，减少了后续数据处理量，提高了定位精度。

所述谱减法原理为：在噪声信号与语音信号相互独立的前提下，通过快速傅里叶变换(FFT),从含噪语音信号中去除噪声功率谱，得到去除噪声的语音频谱。

所述谱减法处理过程为：

设x_w(k)为去噪时域信号，n_w(k)为噪声信号，k表示时间序列，则有

s_w(k)＝x_w(k)+n_w(k)(2)

对式(2)进行傅里叶变换，得

S_w(f)＝X_w(f)+N_w(f)(3)

将式(3)取模的平方，由于x_w(k)与n_w(k)相互独立，有

|S_w(f)|²＝|X_w(f)|²+|N_w(f)|²(4)

则谱减信号X_w(f)为

$\left|X_w\left(f\right)\right|=\sqrt{{\left|S_w\left(f\right)\right|}^2-{\left|N_w\left(f\right)\right|}^2}---\left(5\right)$

对X_w(f)作傅里叶反变换IFFT，得去噪时域信号x_w(k)。

(304)、对(303)所得去噪时域信号x_w(k)由倒谱法进行消除混响处理，得到纯语音时域信号x(k)。由于语音信号的复倒谱通常集中在复倒谱域的原点附近，而混响则由复倒谱域扩展至远离原点的位置，故在复倒谱域中用低通滤波法去除尖峰的办法来削减混响，然后再通过逆向操作，获得混响受到削减的声源信号。图3为倒谱法处理流程图，其中x_w(k)为(303)所得去噪时域信号，x(k)为纯语音时域信号；图4为复倒谱域低通滤波器原理图，其中P为通带的最高截止点，h为过渡带的带宽。

所述倒谱法处理过程为：

(a)、对去噪时域信号x_w(k)进行分帧处理，得到分帧信号，这里选定帧长为512点，帧移为128点；

(b)、对分帧信号进行预加窗处理，得到预加窗信号x_ww(k)，所述窗函数仍采用(302)中的汉明窗函数；

(c)、对预加窗信号x_ww(k)计算每帧的复倒谱x_c(k)，所述复倒谱计算公式为

x_c(k)＝IFFT{ln[FFT[x_ww(k)]]}(6)

式中，FFT表示傅里叶变换，IFFT表示傅里叶反变换，ln表示自然对数；

(d)、对复倒谱信号x_c(k)，利用复倒谱域低通滤波器进行滤波，得到滤波信号；所述复倒谱域低通滤波器由通带、过渡带、阻带三部分组成，其中通带的最高截止点P＝4，过渡带的带宽h为64；

(e)、对滤波信号，在时域条件下将每帧信号叠接相加，得到纯语音时域信号x(k)。

步骤四、对步骤三中的纯语音时域信号x(k)，利用相位变换加权广义互相关方法GCC-PHAT(GeneralizedCrossCorrelation-ThePhaseTransform,GCC-PHAT)计算时延值。

所述的GCC-PHAT计算时延的过程为：

τ_i表示第i个麦克风接收的经步骤三处理所得纯语音时域信号x_i(k)与坐标原点处第0个麦克风接收的经步骤三处理所得纯语音时域信号x₀(k)间的时延值，1≤i≤6且i为自然数。信号x_i(k)与x₀(k)的互相关函数为

$R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right)={\∫}_{-\∞}^{+\∞}{\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)X_0\left(f\right){X_i}^{*}\left(f\right)e^{{j2\mathrm{\πf\τ}}_i}\mathrm{df}---\left(7\right)$

式中，f表示频率，e表示自然指数，X₀(f)表示时域信号x₀(k)的傅里叶变换，X_i ^*(f)表示时域信号x_i(k)傅里叶变换的共轭，为虚数单位，ψ_PHAT(f)为加权函数，且

${\mathrm{\ψ}}_{\mathrm{PHAT}}\left(f\right)=\frac{1}{\left|G_{x_0{x}_i}\left(f\right)\right|}---\left(8\right)$

式中，为相关函数的功率谱，且

$G_{x_0{x}_i}\left(f\right)={\∫}_{-\∞}^{+\∞}{R}_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)e^{{j2\mathrm{\πf\τ}}_i}{\mathrm{d\τ}}_i---\left(9\right)$

式中， $R_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)=E\left[x_0\left(k\right)x_i(k+{\mathrm{\τ}}_i)\right],$ 表示x₀(k)与x_i(k)的互相关函数，则

${\widehat{\mathrm{\τ}}}_i=\arg ,\max ,R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right)---\left(10\right)$

式中，argmax表示使取得最大值时即为第0个麦克风与第i个麦克风所接收到信号的时延估计值。

步骤五、利用三维空间定位方法，求得声源到坐标原点(就是第0个麦克风的位置M₀)的距离r。

空气中声速为

$c=20.05\sqrt{273.15+T}---\left(11\right)$

式中，T代表室温，单位为℃。根据室温为14.5℃，取声速c为340.053m/s。

结合图2，得麦克风阵列与声源的几何位置的关系式为

$\left\{\begin{array}{c}\mathrm{OS}:x^2+y^2+z^2=r^2\\ {\mathrm{SM}}_1:{(x-L)}^2+y^2+z^2={(r+d_1)}^2\\ {\mathrm{SM}}_2:x^2+{(y-L)}^2+z^2={(r+d_2)}^2\\ {\mathrm{SM}}_3:{(x+L)}^2+y^2+z^2={(r+d_3)}^2\\ {\mathrm{SM}}_4:x^2+{(y+L)}^2+z^2={(r+d_4)}^2\\ {\mathrm{SM}}_5:x^2+y^2+{(z+L)}^2={(r+d_5)}^2\\ {\mathrm{SM}}_6:x^2+y^2+{(z-L)}^2={(r+d_6)}^2\end{array}\right.---\left(12\right)$

式中，d_i(1≤i≤6)为声源到第0个麦克风的位置M₀与声源到第i个麦克风的位置M_i的距离差值。

结合步骤四中求出的时延计算值，得声源与坐标原点的距离r为：

$r=\frac{{6L}^2-c^2\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\left({\widehat{\mathrm{\τ}}}_i\right)}^2}{2c\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\widehat{\mathrm{\τ}}}_i}---\left(13\right)$

步骤六、利用坐标旋转数字式计算机方法，结合步骤四中求出的时延估计值，得声源的俯仰角θ、方位角

由于是远场条件下，式(14)可化简为

步骤七、根据直角坐标与球面坐标间的关系，得声源位置的直角坐标为

结合步骤五中的r，步骤六中的θ、即可确定声源位置的直角坐标S(x,y,z)。

图5是本发明中麦克风阵列模型，图5中的(a)为七元麦克风阵列，图5中的(b)为四元麦克风阵列。图5中的(a)所示阵列进行声源定位效果与图5中的(b)进行比较，以说明本发明定位效果；房间尺寸为6m×5m×3.5m。

七元麦克风阵列与四元麦克风阵列的声源定位比较：

麦克风间距为0.1m时，七元麦克风的位置依次为M₀(0.05,0.05,0.05),M₁(0.1,0.05,0.05),M₂(0.05,0.1,0.05),M₃(0,0.05,0.05),M₄(0.05,0,0.05),M₅(0.05,0.05,0),M₆(0.05,0.05,0.1)；四元麦克风的位置依次为：M₁(0.1,0.05,0.05),M₂(0.05,0.1,0.05),M₃(0,0.05,0.05),M₄(0.05,0,0.05)。

麦克风间距为0.2m时，七元麦克风的位置依次为M₀(0.1,0.1,0.1),M₁(0.2,0.1,0.1),M₂(0.1,0.2,0.1),M₃(0,0.1,0.1),M₄(0.1,0,0.1),M₅(0.1,0.1,0),M₆(0.1,0.1,0.2)，四元麦克风的位置依次为M₁(0.2,0.1,0.1),M₂(0.1,0.2,0.1),M₃(0,0.1,0.1),M₄(0.1,0,0.1)。

在声源的实际位置为S(3.25，2.15,1.55)时，表1给出了麦克风间距为0.1m、信号采样率为16kHz时的声源定位结果比较；表2给出了麦克风间距为0.2m、信号采样率为16kHz时的声源定位结果比较；表3给出了麦克风间距为0.1m、信号采样率为64kHz时的声源定位结果比较。表中ΔX,ΔY,ΔZ分别表示声源位置直角坐标X,Y,Z绝对误差。

表1

表2

表3

表1表明，在麦克风间距为0.1m、信号采样率为16kHz条件下，用七元麦克风阵列所得声源位置X、Y、Z的绝对误差分别比四元麦克风阵列减少0.2751m、0.0520m、0.1625m。

表2表明，在麦克风间距为0.2m、信号采样率为16kHz条件下，用七元麦克风阵列所得声源位置X、Y、Z的绝对误差分别比四元麦克风阵列减少0.0186m、0.1058m、0.2244m。

表3表明，在麦克风间距为0.1m、信号采样率为64kHz条件下，用七元麦克风阵列所得声源位置X、Y、Z的绝对误差分别比四元麦克风阵列减少0.2004m、0.1325m、0.3044m。

因此，在信号采样率为16kHz不变、麦克风间距增加的条件下，或在麦克风间距为0.1m不变、信号采样率增加的条件下，七元麦克风阵列声源定位的绝对误差比四元麦克风阵列声源定位的绝对误差都小，定位精度都高。

声源与坐标原点间的距离r变化时，七元麦克风阵列的定位结果比较：

为了用七元麦克风阵列对处于不同位置声源进行定位，图6、图7给出了声源与坐标原点间的距离r变化时的定位结果。

图6是本发明麦克风间距与方位角、俯仰角的关系图，麦克风间距为0.1m，采样率为16kHz条件下；图6中的(a)为麦克风间距为0.1m时与方位角的关系图，图6中的(b)为麦克风间距为0.1m时与俯仰角的关系图。

图7是本发明麦克风间距与方位角、俯仰角的关系图，麦克风间距为0.2m，采样率为16kHz条件下；图7中的(a)为麦克风间距为0.2m时与方位角的关系图，图7中的(b)为麦克风间距为0.2m时与俯仰角的关系图。真实值为空心原点，测量值为实心原点。

由图6与图7所示的声源位置球面坐标数据值，利用球面坐标与直角坐标的关系，得到声源位置的直角坐标测量值与绝对误差，如表4与表5所示。表4给出了麦克风间距为0.1m、信号采样率为16kHz条件下的声源定位结果，表5给出了麦克风间距为0.2m、信号采样率为16kHz下的声源定位结果。

表4

表5

表4与表5表明：在距离r由0.5m增到8.3m、信号采样率为16kHz不变时，随着麦克风间距增加，七元麦克风阵列的声源定位平均绝对误差逐渐减小、定位精度逐步提高。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (7)

1.一种七元麦克风阵列的声源定位方法，其特征在于，包括如下步骤：

步骤一、在三维坐标系中，建立麦克风阵列模型，确定各麦克风在三维空间的位置；

步骤二、对各麦克风接收到声源发出的语音信号s(t)进行采样，得到离散信号s(k)，t表示时间，k表示时间序列；

步骤三、对步骤二所述的离散信号s(k)进行滤波、加窗、谱减、倒谱处理，得到纯语音时域信号x(k)；具体如下：

(301)、将离散信号s(k)通过带通滤波器进行带通滤波后，得到带通信号；

(302)、对(301)中所得带通信号用窗函数进行加窗处理，得到加窗信号s_w(k)；

(303)、对(302)所得加窗信号s_w(k)用谱减法进行去噪处理，得到谱减信号X_w(f)，对X_w(f)作反傅里叶变换得到去噪时域信号x_w(k)；

(304)、对去噪时域信号x_w(k)由倒谱法进行消除混响处理，得到纯语音时域信号x(k)；

步骤四、对步骤三中的纯语音时域信号x(k)，利用相位变换加权广义互相关方法GCC-PHAT计算第i个麦克风接收的经步骤三处理所得纯语音时域信号x_i(k)与三维坐标原点处第0个麦克风接收的经步骤三处理所得纯语音时域信号x₀(k)间的时延值τ_i，1≤i≤6且i为自然数，则互相关函数为

$R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{\mathrm{\ψ}}_{PHAT}\left(f\right)X_0\left(f\right)X_i^{*}\left(f\right)e^{j2{\mathrm{\πf\τ}}_i}df;$

式中，f表示频率，e表示自然指数，X₀(f)表示纯语音时域信号x₀(k)的傅里叶变换，X_i ^*(f)表示纯语音时域信号x_i(k)傅里叶变换的共轭，j为虚数单位，ψ_PHAT(f)为加权函数，且

${\mathrm{\ψ}}_{PHAT}\left(f\right)=\frac{1}{\left|G_{x_0{x}_i}\left(f\right)\right|};$

式中，为相关函数的功率谱，且

$G_{x_0{x}_i}\left(f\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{R}_{x_0{x}_i}\left({\mathrm{\τ}}_i\right)e^{j2{\mathrm{\πf\τ}}_i}{\mathrm{d\τ}}_i;$

式中， 表示x₀(k)与x_i(k)的互相关函数，则 ${\widehat{\mathrm{\τ}}}_i=\arg \max R_{x_0{x}_i}^{\left(g\right)}\left({\mathrm{\τ}}_i\right),$ 其中argmax表示使取得最大值时 ${\mathrm{\τ}}_i={\widehat{\mathrm{\τ}}}_i,$ 即为第0个麦克风与第i个麦克风所接收到信号的时延估计值；

步骤五、利用三维空间定位方法，结合步骤四中的得声源到三维坐标原点的距离r：

$r=\frac{6L^2-c^2\underset{i=1}{\overset{6}{\Σ}}{\left({\widehat{\mathrm{\τ}}}_i\right)}^2}{2c\underset{i=1}{\overset{6}{\Σ}}{\widehat{\mathrm{\τ}}}_i};$

式中，c表示声速，L表示第i个麦克风位置M_i与三维坐标原点处第0个麦克风位置M₀间的距离；

步骤六、利用坐标旋转数字式计算机方法，结合步骤四得到的得声源的俯仰角θ、方位角为

其中，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₁(L,0,0)的第1个麦克风所接收到信号的时延估计值，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₂(0,L,0)的第2个麦克风所接收到信号的时延估计值，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₃(-L,0,0)的第3个麦克风所接收到信号的时延估计值，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₄(0,-L,0)的第4个麦克风所接收到信号的时延估计值，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₅(0,0,-L)的第5个麦克风所接收到信号的时延估计值，表示位置坐标为M₀(0,0,0)的第0个麦克风与位置坐标为M₆(0,0,L)的第6个麦克风所接收到信号的时延估计值；

步骤七、根据直角坐标与球面坐标间的关系，声源位置的直角坐标为

2.根据权利要求1所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(302)中窗函数采用汉明窗函数w(k)，即

式中，N表示时间序列的取值范围，且N为大于零的正整数。

3.根据权利要求1所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(303)中谱减法的公式为

$\left|X_w\left(f\right)\right|=\sqrt{|S_w\left(f\right){|}^2-|N_w\left(f\right){|}^2};$

式中，S_w(f)为加窗信号s_w(k)的傅里叶变换，N_w(f)为加窗信号s_w(k)中噪声信号n_w(k)的傅里叶变换。

4.根据权利要求1所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(304)中倒谱法处理过程为：

(a)、对去噪时域信号x_w(k)进行分帧处理，得到分帧信号；

(b)、对分帧信号用窗函数进行预加窗处理，得到预加窗信号x_ww(k)；

(c)、对预加窗信号x_ww(k)计算每帧的复倒谱信号x_c(k)，所述复倒谱计算公式为

x_c(k)＝IFFT{ln[FFT[x_ww(k)]]}；

式中，FFT表示傅里叶变换，IFFT表示傅里叶反变换，ln表示自然对数；

(d)、对复倒谱信号x_c(k)，利用复倒谱域低通滤波器进行滤波，得到滤波信号；

(e)、对滤波信号，在时域条件下将每帧信号叠接相加，得到纯语音时域信号x(k)。

5.根据权利要求1所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(301)中带通滤波器的上截止频率为f_H＝3500Hz，下截止频率为f_L＝250Hz。

6.根据权利要求4所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(a)中分帧信号的帧长为512点，帧移为128点。

7.根据权利要求4所述的一种七元麦克风阵列的声源定位方法，其特征在于，所述(d)中复倒谱域低通滤波器由通带、过渡带、阻带三部分组成，其中通带的最高截止点P＝4,过渡带的带宽h为64。