大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

一. 关于麦克风阵列

麦克风阵列： 麦克风阵列是由一定数目的声学传感器（麦克风）按照一定规则排列的多麦克风系统，而基于麦克风阵列的声源定位是指用麦克风拾取声音信号，通过对麦克风阵列的各路输出信号进行分析和处理，得到一个或者多个声源的位置信息。

麦克风阵列系统的声源定位技术研究意义在于： 输入的信息只有两个方向难以确定声源的位置，人类的听觉系统主要取决于头和外耳气压差声波实现声源定位。假使没有这个压力差，只能定位在平面上声源的位置，但就无法知道声音是从前面，或从后面传来的。因此，由人的听觉系统，科技研发人员得到了灵感，使用多个麦克风系统可以实现在三维空间中的声源位置的定位，麦克风的数量越多，所接收到的信息量也越多。声源的声源定位和声源增强是实现智能处理的两个关键问题，而声源定位是实现语音增强的前提和基础。一个麦克风的信息量较少，使得声源定位所需的信息缺乏，而麦克风阵列克服了上述缺点，充分利用每个麦克风信号之间的数据相关性，并加以融合，可以实现声源定位。

麦克风阵列声源定位技术的应用： 广泛应用于国防、智能机器人、视频会议及语音增强等众多领域,尤其在当下以智能办公和智能家居为主要室内场景的远场语音交互系统中。

二. 关于声源定位

目前基于麦克风阵列的声源定位方法主要有三种：基于最大输出功率的可控波束成形的定位方法、基于高分辨谱估计的定位方法、基于到达时延差估计的定位方法(Time Difference of Arrival，TDOA)。

1. 基于最大输出功率的可控波束成形定位法。 波束形成法的原理是将麦克风接收到的信号进行滤波加权求和来形成波束，按照一定的规律对声源位置进行搜索，当麦克风达到最大输出功率时，为时搜索到的声源位置即为真实的声源方位。波束形成可分为常规的波束形成CBF（Conventional Beam Forming）、CBF+Adaptive Filter和自适应波束形成ABF（Adaptive Beam Forming）。
2. 基于高分辨谱估计的定位法。 基于高分辨率谱估计的定位方法通过分解协方差矩阵估计声源方位。适合多个声源的情况，且声源的分辨率与阵列尺寸无关，突破了物理限制。该方法的优点是不受采样频率限制，且在一定程度下可以实现任意程度的定位，但是该方法计算复杂度较高，抗噪和抗混响性能较差，因此该方法适合在一些特定的环境下使用。这类方法可以拓展到宽带处理，但是对误差十分敏感（如麦克风单体误差，通道误差），适合远场模型，且矩阵运算量巨大。
3. 基于到达时延差估计的定位法。 TDOA（time difference of arrival）是先后估计声源到达不同麦克风的时延差，通过时延来计算距离差，再利用距离差和麦克风阵列的空间几何位置来确定声源的位置。可分为TDOA估计（估计信号到达各麦克风的时间差）和TDOA定位（运用几何关系确定声源位置）两步。
   

三. 基于广义互相关(GCC)计算时延

时延估计有很多种，比较经典就是广义互相关函数 (Generalized Cross Correlation, GCC) 估计时延，这里简单介绍基于广义互相关函数估计时延的方法。

在噪声存在情况下，一个由远处声源发出的，并且被两个不同空间中的麦克风监听的信号可以数学建模为：

$x_1=s_1(t)+n_1(t)$                          (1)

$x_2=\alpha s_1(t-D)+n_2(t)$               (2)

其中，$s(t)$ 是声音信号，$n_1(t)、n_2(t)$是两个声音传感器检测噪声。 三者是稳定的随机过程，且互不相关。

计算 $x_1$ 与 $x_2$ 的互相关函数：

$R_{x_1{x}_2}(\tau )=E[x_1(t)\ast x_2(t-\tau )]$               (3)

${\hat{R}}_{x_1{x}_2}(\tau )=\frac{1}{T-\tau }{\int }_{\tau }^T{x}_1(t)x_2(t-\tau )\mathrm{d}t$       (4)

其中估计的时延 $D$ 为互相关函数值达到最大值时取得的 $\tau$ 值，即：

$\hat{D}=\underset{\tau }{\mathrm{argmax}}{R}_{x_1{x}_2}(\tau )$                                 (5)

MATLAB 例程：

% 导入两个麦克风的音频数据
[y_0,Fs] = audioread('音轨-0.wav'); 
[y_1] = audioread('音轨-1.wav');

fprintf('采样频率：%d\n ······\n', Fs);

% 取出两段音频前2048个采样点，并作互相关处理，绘制曲线。
A = y_0(1:2048);
B = y_1(1:2048);

[value,delay] = xcorr(A,B);

subplot(1,2,1);
plot(delay, value);

D = zeros(1,926); 

% 以2048个点为一帧，计算互相关得到的时延，绘制出两段音频的时延变化。
for a = 1:2048:size(y_0,1)-2048
    A = y_0(a:a+2048);
    B = y_1(a:a+2048);
    [value,delay]=xcorr(A,B);
    value_max_idx = find(value==max(value));
    D1 = delay(value_max_idx);
    D((a-1)/2048+1) = D1;
end

subplot(1,2,2);
plot(D);

使用ReSpeaker的树莓派六麦克风套件，拿手机放歌作为声源围绕麦克风阵列移动，录制音频，并保存为6个wav音频文件，采用率为44100Hz。对其中两个麦克风音频文件做处理，如果取初始的2048个采样点做互相关处理就能得到下面图一，图一的峰值的横坐标即为估计时延。如果以2048个采样点作为一帧，对两段音频的每一帧进行互相关，并找到每帧的估计时延，绘制出来就能得到时延的变化曲线，如下面图二。

注意： 这里的时延并不是指具体的时间，而是采样点数。想要得到具体时间还需要除以采样频率。以2048采样点为一帧，就是以0.464s为一帧（$\frac{2048}{44100}\approx 0.464$s）。

四. 基于时延差的声源定位法

在对麦克风阵列进行建模之前，我们需要分清楚什么近场与远场。顾名思义，离麦克风近则符合近场模型 ，离得远则符合远场模型 。

假设$L$ 为阵列间距，$\lambda$ 为声波波长，$M$为声源与麦克风的距离，我们定义 $\frac{2L^2}{\lambda }$ 为远近场临界值。

当 $M<\frac{2L^2}{\lambda }$ 时，符合近场模型，此时声源到达麦克风阵列的波形视为球面波。

当 $M<\frac{2L^2}{\lambda }$ 时，符合远场模型，此时声源到达麦克风阵列的波形视为平面波。

可听声的机械波频带为20Hz ~20000Hz，机械波波长大约在1.7cm ~ 17m（声速取340m/s）。然而在现实生活中，过高频率或过低频率的声波都是非常少量的，以人的声音为例，人语音频带的范围大概为300Hz至3400Hz，波长范围为0.1m~ 1.12m，当阵列间隔取4cm时，远近场临界值范围为 0m~3.2m。若取中间值，则可以认为1.6m内符合近场模型，1.6m外符合远场模型。

1. 近场模型

当 $M<\frac{2L^2}{\lambda }$ 时，符合 近场模型，此时声源到达麦克风阵列的波形视为 球面波。

近场模型至少需要3个麦克风，以最简单的3麦克风模型为例（如图：$y_1、y_2、y_3$）。假设 ${\tau }_{12}、{\tau }_{13}$分别表示第一个麦克风与第二和第三个麦克风之间的时延，那么有：

${\tau }_{12}=\frac{r_2-r_1}{c}$         (6)

${\tau }_{13}=\frac{r_3-r_1}{c}$         (7)

其中，$c$ 为声速。标准大气压、$15°$ 条件下，声速为 $340m/s$。

根据麦克风阵列的的几何关系，由余弦定理，可以得到：

$r_2^2=r_1^2+d^2+2r_{1}dcos{\theta }_1$     (8)

$r_3^2=r_1^2+4d^2+4r_{1}dcos{\theta }_1$    (9)

其中， ${\tau }_{12}、{\tau }_{13}$可以通过互相关GCC求得（下面会细说），且 $c、d$ 为已知。结合(6) ~ (9) 即可求出未知量 $r_1、r_2、r_3、\theta$ ，结合坐标系可以求出$s(k)$ 的坐标。

2. 远场模型

当 $M<\frac{2L^2}{\lambda }$ 时，符合 远场模型，此时声源到达麦克风阵列的波形视为 平面波。

对于两个麦克风的情况，只能计算出方位角，无法计算出方位距离。假设 $\tau$ 为声波到达两个麦克风之间的时延，则有：

$\tau =\frac{dcos\theta }{c}$                      (10)

$\theta =arccos\frac{c\tau }{d}$              (11)

五. 三维空间阵列的声源定位系统实现

计算方法：几何求解、最小二乘法

上面所举的近场、远场的定位模型例子，分析都是最简模型下，平面声源的定位。下面以三维麦克风阵列为例，通过几何的方式求解三维空间下声源的定位。

1. 推导过程

一个三维麦克风阵列，麦克风分别为 $m_i(i=0、1\cdot \cdot \cdot n)$，声源 $S$ 符合近场模型。现以麦克风 $m_0$ 为原点，如图建立直角坐标系。推导公式之前，需要先确定以下概念：

序号	概念	字符
1	麦克风的坐标	$m_i$   $i\in [0,n]$
2	声源的估计坐标	$S$
3	声源 $s$ 到 $m_i$ 与 $m_1$ 的估计距离差	${\hat{d}}_i$
4	麦克风 $m_i$ 到原点的的距离	$R_i$
5	声源 $s$ 到原点的的距离	$R_s$

如上图，根据三角形 $m_0{m}_{i}s$，由余弦定理有：

$(R_s+d_i{)}^2=R_i^2-2m_i^{T}s+R_s^2$                      (12)

公式 (12) 展开后得到：

$R_i^2-2m_i^{T}s-2R_s{d}_i-d_i^2=0$                         (13)

由于 $d_i$ 是通过估计的时延得到的，与实际值之间会有一个偏差，因此公式 (10) 不为零，可写为：

$\epsilon =(R_i^2-d_i^2)-2m_i^{T}s-2R_s{d}_i$                      (14)

此时已经得到目标值的误差函数，使用最小二乘法求解，问题可以转化为：估计声源坐标 $s(x,y,z)$，使得最终的误差平方和最小，即：

$\underset{s(x,y,z)}{\mathrm{argmin}}\sum _{i=1}^n[(R_i^2-d_i^2)-2m_i^{T}s-2R_s{d}_i{]}^2$      (15)

将公式 (14) 写成矩阵形式：

$\epsilon =2R_s\hat{D}+2Ms-\delta$                                          (16)

其中，$M=\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\\ \dots .\\ m_n\end{array}\right]=\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ x_3& y_3& z_3\\ \dots & \dots & \dots \\ x_n& y_n& z_n\end{array}\right]$      $\hat{D}=\left[\begin{array}{c}{\hat{d}}_1\\ {\hat{d}}_2\\ {\hat{d}}_3\\ \dots \\ {\hat{d}}_n\end{array}\right]$         $\delta =\left[\begin{array}{c}{R}_1^2-{\hat{d}}_1^2\\ R_2^2-{\hat{d}}_2^2\\ R_3^2-{\hat{d}}_3^2\\ \dots \dots \dots .\\ R_n^2-{\hat{d}}_n^2\end{array}\right]$

公式 (16) 可以简化为：

$\epsilon =A\mu -b$                                   (17)

其中，$A=\left[\begin{array}{cc}M& \hat{D}\end{array}\right]$       $\mu =\left[\begin{array}{c}s\\ R_s\end{array}\right]$        $b=\frac{1}{2}\delta$

公式 (17) 的最小二乘解可以表示为：

$\hat{\mu }=(A^{T}A{)}^{-1}{A}^{T}b$                       (18)

关于矩阵形式下最小二乘解推导过程可以参考下面的链接（点击跳转），此处不详细讲解。

公式 (18) 即为计算结果，下面对结果进行进一步化简：

1. 定义沿 $A$ 到 $\hat{D}$ 的投影矩阵为： $P_A=\frac{\hat{D}{\hat{D}}^T}{{\hat{D}}^T\hat{D}}$                                                     (19)

2. 沿 $\hat{D}$ 到 $A$ 的投影矩阵即为：$P_D=I-P_A=I-\frac{\hat{D}{\hat{D}}^T}{{\hat{D}}^T\hat{D}}$（ $I$ 为单位矩阵）     (20)

3. 根据投影矩阵的性质，可以得到：$A=P_D\left[\begin{array}{cc}M& 0\end{array}\right]$                       (21)

4. 将公式 (21) 带入 (17) 中得到：$\epsilon =P_{D}Ms-b$                              (22)

5. 公式（22）的最小二乘解即为最终简化结果。

最终简化结果：

$s=(M^T{P}_{D}M{)}^{-1}{M}^T{P}_{D}b$

$M$只与麦克风阵列的个数与排列有关，$P_D$ 与 $b$ 在计算出时延估计的距离差 ($d_i$) 后也可以求得。

六. 六元圆形麦克风阵列声源定位

硬件平台： Seed Studio 的六元麦克风阵列套件（官网)）

软件平台： Matlab 2019a

实际应用中，二维圆形阵列主要采用远场模型较多，因为可以在保证计算精度的同时，大大降低复杂性，减少运算量，且不需要较大的阵列间距，此处使用近场模型纯属因为要做毕设?。公式都是根据三维阵列来推导的，之所以使用二维平面阵列是因为：没钱搭三维阵列，但又需要点实验数据放进毕业论文中。

二维阵列缺少$z$轴的基元，使用第六章中的公式求解出的声源坐标值没有实际意义，但是可以通过 $y$ 与 $x$ 的比值得到声源的偏航角。

${\theta }_{yaw}=arctan\frac{y}{x}$        (23)

MATLAB 代码：

#代码写得太烂，自行脑补

现以树莓派Seed Studio麦克风阵列套件为硬件平台进行环境声音的获取，将各个麦克风所获取到的音频储存为wav格式文件，采样率为44100Hz，将音频文件在PC机中使用Matlab进行处理运算。麦克风阵列为六元均匀圆阵，属于二维平面阵列，阵元之间的间距为4cm。

如上图所示，为六个麦克风 收集到的声音曲线，其中横坐标表示时间，纵坐标表示声音幅值。麦克风录制到的信号是一段长时间的音频片段，因此需要进行分帧处理，此处以2048个采样点为一帧，由于时间=采样点数/采样频率，因此以2048个采样点为一帧即以0.0464秒为一帧。我们假设一帧时间内，声音信号是联合平稳的，声源相对于参考坐标系的相位运动足够小可以忽略不计。

如上图所示，左图为该阵列对静态声源的估计曲线，右图则为对动态声源的估计曲线。

左图： 在空旷且安静的室内环境，使用手机播放歌曲作为声源，在距离麦克风阵列0.5m处保持静止，同时记录下六个麦克风录制的音频信息，采样率为44100Hz，以2048个采样点为一帧，使用声达时延法计算每一帧声源的方位角。

右图： 在空旷且安静的室内环境，使用手机播放歌曲作为声源，在距离麦克风阵列中心0.5m处绕麦克风阵列中心匀速旋转，同时记录下六个麦克风录制的音频信息，采样率为44100Hz，以2048个采样点为一帧，使用声达时延法计算每一帧声源的方位角。

七. 相关链接

相关链接：

1. 矩阵形式下最小二乘解推导：https://zhuanlan.zhihu.com/p/87582571

2. 关于投影矩阵与最小二乘：https://www.cnblogs.com/bigmonkey/p/9897047.html

3. 关于投影矩阵概念及其性质：https://blog.csdn.net/niu_123ming/article/details/86371308

4. 参考文献：基于六元空间阵列的声源定位系统实现

5. 参考文献：基于麦克风阵列的室内三维声源定位优化算法