基于tdoa的麦克风阵列单声源定位(附详细代码)

Greener_Pat

已于 2024-11-14 11:01:03 修改

阅读量3.5k

点赞数 18

分类专栏：声纪元文章标签：硬件工程实时互动音视频

于 2024-08-23 12:37:51 首次发布

本文链接：https://blog.csdn.net/a2023816/article/details/141462701

版权

声纪元专栏收录该内容

2 篇文章

订阅专栏

基于tdoa的麦克风阵列单声源定位(附详细代码)

一、前言

人一生中接收到的所有信息中，听觉共占11%。仅仅通过双耳的听觉，人便可以分辨声源的方位，识别说话的对象，同时听取多个声音且互不干扰。而计算机自从1939年诞生起，经过摩尔定律的疯狂增长，如今的计算能力，无论是速度还是精度上都已经远超人类。那原则上讲，如果没有灵魂或者上帝作用到了人对世界的理解中，只要给了计算机以与人类所接收到的相同的信息，那么计算机也同样能完成上述这些奇妙的任务，甚至达到更高的的水平。能不能创造神迹不好说，但做一些有趣的工作不在话下~

但为什么不用视觉呢？视网膜高清所承载的信息难道不比听到充斥着杂音的听觉音频信息量含量高、信噪比低吗？某种意义上确实，但正是因为这些原因，人对听觉的信息利用率远远低于视觉，而这一个通道所承载的信息又是视觉所无法取代的。抛开这些，难道这种用无法“眼见为实”的方法，近似于魔法的方式来理解世界的方式不足以吸引你吗？话不多说

Let’s step into the Sound Era

二、信息论声学

从信息论的角度出发，声音所承载的信息可以简单地分为两种——信源和信道，或者简单点来说，声源和环境。有人可能要说了，环境不就是噪声吗？哪来的信息？这里需要说明一下，此处的环境是一种泛在的概念，指的是从发声源到接收器过程中所有参与到最终结果的合成中的信息，包括其他的非主要声源，也就是所谓的噪声。除此之外环境信息其实还包括如房间形状，物品陈设的空间信息 (自然可以想到人在小房间的回声和室外的空旷感是不同的，里面包含着所处地点的空间信息);

另外，这里还需要对刚才指一种所谓的“噪声”来一次平反，所谓噪声，都是没有被充分利用的信息，你怎么敢定义什么是噪声？假若我就是想听听窗外的喇叭声来判断交通情况而不是听电视里男女主深情告白呢？从信息论的角度来看。从信息论的角度来看，信源和信道是对立统一且本质相同的。这看起来或许有些不切实际，毕竟物理世界中有那么多参数，根本就是一个复杂系统，你又如何从进行了纷繁复杂的叠加后的复合信息中体取出原来的各个分量呢？

不得不说，物理世界诚然复杂，但不难发现它的复杂背后也有其秩序，不同变量对一个物理过程的影响程度(也就是权重)是不同的，我们只要抓住其中的主因子便能做一系列有趣的分析了，比如接下来的声源定向

三、声源定位简介

声源定位顾名思义指的就是定位声源的位置，最常见的便是我们人耳的声源定位，而在实际中，我们利用多个麦克风（麦克风阵列）在测量不同位置点对声源进行测量，而由于声信号到达不同麦克风的时间有不同程度的延迟（也被称为时延），利用算法对测量得到的声信号进行处理，由此获得声源点相对于麦克风的到达方向（包括方位角，俯仰角）和距离等

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

四、单声源定位技术路线

1. 使用设备

这里使用的都是ReSpeaker Mic Array v2.0，如图有四个麦克风，采样率为16000，采样深度为16位。录音会得到6声道数据，通道0为asr处理声道，通道5是合并播放，通道1-4则为4个麦克风所接收到的数据(分别对应图中右上、左上、左下、右下)

由于该设备只有平面上的四个摄像头，难以采集竖直方向的信息，所以后续的声源定位均是针对二维的

这是它的介绍页面，也有购买链接(只是给感兴趣的朋友，真的没有收钱啊😢)

2. 数学基础

根据到达时间差的声源定位根据声源距离麦克风距离可以分为远场和近场两种模型，这里会对两者的特点以及进行tdoa的分析做详细的说明(但最终使用的是远场模型)

远场模型

远场即麦克风与声源之间的距离远远大于麦克风之间的距离，此时只能通过不同麦克风接收到声源的相对时间差计算出声源的方向。但声源的具体位置无法确定，因为此时麦克风阵列接收到的近似于平面波，不足以区分同一方向不同距离的声源

远场声源数学模型大致可以简化为此，即确定一个麦克风为原点并建立直角坐标系，编号为 $R_0(0, 0)$ ; 其余的n-1个麦克风编号为 $R_i(x_i, y_i), i = 1, 2, ... n-1$ ，而声源则记为 $S (a, b)$

由远场模型中声波可视为平面波的特点，声源S到 $R_0, R_i$ 的距离差 $R_0 S| - |R_i S|可视为|R_0 P|$ ，所以有公式

$|R_0 R_i|\cos{\theta} = |R_0 P| = v_{sound}\Delta t_i$

为更好地使用坐标信息，这里采用向量化的计算方式

$|R_0 R_i| \cdot \frac{\vec{R_0 R_i} \cdot \vec{R_0 R_i}}{|R_0 R_i|\cdot|R_0 R_i|} = v_{sound}\Delta t_i$

具体将坐标带入即(其中只有a, b是未知量):

$\frac{a}{\sqrt{a^2 + b^2}} x_i + \frac{b}{\sqrt{a^2 + b^2}} y_i = v_{sound}\Delta t_i$

由于存在多个麦克风，所以做矩阵方程:

其实求解这个二元一次的变量组，只需要两个线性无关的方程即可，也就是三个麦克风。但如果只是用三个声道的话，首先对资源的利用，也就是这里的四个麦克风(甚至在别的地方可能更多)是不到位的; 另一面，考虑到实际对声源到不同麦克风的时间差的测量可能是存在误差的，综合多个麦克风得到的信息，才能减少误差的影响，也就是这里用到的最小二乘法。

最小二乘法在矩阵上的计算可以参考这篇博客😄

最终得到的结果:

$\vec{x} = (M^TM)^{-1}M^T \cdot v_{sound}\vec{b}$

由此可以计算出声源的方向:

$\angle SR_0x = \arctan2(x[0], x[1])$

近场模型

近场即麦克风与声源的距离接近麦克风之间的距离较小，此时声源发出的波可视为球面波。借助麦克风阵列，通过三角计算等方式，可以确定声源的具体位置(包括方向和距离)

同理做数学推导:

$|R_i S|^2 = |R_0 R_i|^2 + |R_0 S|^2 - 2 \vec{R_0 S}\cdot\vec{R_0 R_i}$

化简得到:

$\vec{R_0 S}\cdot\vec{R_0 R_i} - 2m\Delta d_i + (\Delta d_i)^2 - r_i^2$

整理为矩阵形式:

请添加图片描述

同样用最小二乘解出 $\vec{x}$ ，再取出前两项即是声源的坐标

远近场的区分主要是比较声源与麦克风之间距离L与 $\frac{2d^2}{\lambda}$ 之间的关系，其中d为阵列间距， $\lambda$ 为声波波长。若L更大，则为原厂，反之为近场。但是其实也是一个模糊的选择方法，这里考虑到计算简单，就选择了远场模型(实际效果也还不错~)

3. 代码实现(远场模型)

这里先对其中求交叉相关的部分做一些简单的介绍，这里交叉相关的目的即是要求声源到不同麦克风之间的距离差

已知声速的情况下，要求距离差也就是求声音到达的时间差。所谓交叉相关就是通过将一个信号或序列与另一个信号或序列进行滑动比较来计算的，它衡量了一个信号在不同时间延迟下与另一个信号的相似程度

但是要应用到实时检测的场景中，还要解决一些问题，一个是实时录取的音频中由于不是静态完成的音频段，所以必须对该音频进行分段后再处理。分段的长度的选择即是要在保证能正确比较出两个信号的先后距离的条件下，尽可能的短(会使得交叉相关的复杂度小)

但是，如果只是对每一小段进行声源分析，很有可能因为一些噪声或计算误差导致结果的不稳定。但如果是对多段求平均的话，在实际测试后发现，又有可能因为异常的大偏离值而出现问题。所以这里采取一种“投票”的方式，只是因为发现了由于麦克风阵列之间的距离较小，导致交叉相关实际计算出的时间差由于计算精度的影响呈现出离散的现象，相对正确的值(根据每一段计算出)往往连续有多个，所以便取一定数量的段(维护一个队列)中计数量最大的时间差作为结果代入TDOA的计算

整体代码如下

# 可以先用这段代码确定麦克风阵列对应的index(即是下文的RESPEAKER_INDEX)
import pyaudio
p = pyaudio.PyAudio()
info = p.get_host_api_info_by_index(0)
numdevices = info.get('deviceCount')

for i in range(0, numdevices):
        if (p.get_device_info_by_host_api_device_index(0, i).get('maxInputChannels')) > 0:
            print ("Input Device id ", i, " - ", p.get_device_info_by_host_api_device_index(0, i).get('name'))

# 具体的实时声源检测程序如下
import pyaudio
import numpy as np
from collections import Counter
import matplotlib.pyplot as plt

# 参数设置(根据麦克风阵列的具体参数做调整)
RESPEAKER_RATE = 16000
RESPEAKER_CHANNELS = 6  # change base on firmwares, default_firmware.bin as 1 or i6_firmware.bin as 6
USED_CHANNELS = 4       # 有多少个声道是具体使用的
RESPEAKER_WIDTH = 2
RESPEAKER_INDEX = 1  # refer to input device id
CHUNK = 1024

REMAIN_LENGTH = 10      # 保留的投票段数(考虑应对变化的速度，设小会更灵敏但可能造成不稳定)
SEG_NUM = 3             # 累计多少段之后做一次交叉相关(避免音频写入快于处理)
SEG_LEN = 500           # 用于交叉相关的采样点数量
OVERLAP = 50            # 每段的重叠长度
SOUND_V = 343           # 声速
ANGLE_UNIT = 0.00001       # 角度最小单位(用于离散化)

# 图形化界面相关参数
RADIUS = 3
WIDTH = 10

# 每个声道相对于声道0的交叉相关结果
cor_res_lists = [[] for _ in range(USED_CHANNELS - 1)]
# 存储每个声道音频数据的列表
audio_data = [[] for _ in range(USED_CHANNELS)]
# 用于交叉相关的切割数据
cor_data = [[] for _ in range(USED_CHANNELS)]
# 与声源的距离之差
relative_dis = np.zeros(3)
# 各声道坐标(坐标系实际是按照1号麦克风为原点，1-2号为x轴正方向，1-4号为y轴正方向建立的，但画图时还是按照上面设备图的方向作的上下左右)
M = np.array([
    [0.045, 0],
    [0.045, 0.045],
    [0, 0.045]
])
keep_theta = [0]

def tdoa():
    # 交叉相关处理
    for i in range(USED_CHANNELS-1):
        # 对其中n-1个麦克风求相对于0号麦克风早接受的时间
        for j in range(SEG_NUM):
            tmp_cor_res = np.correlate(cor_data[0][j*SEG_LEN:(j+1)*SEG_LEN], 
                                        cor_data[i][j*SEG_LEN:(j+1)*SEG_LEN], mode='full')
            
            discrete_res = (np.argmax(tmp_cor_res) - (SEG_LEN - 1)) / RESPEAKER_RATE
            # discrete_res = round((np.argmax(tmp_cor_res) - (SEG_LEN - 1)) / RESPEAKER_RATE / ANGLE_UNIT) * ANGLE_UNIT
            cor_res_lists[i].append(discrete_res)
        cor_res_lists[i] = cor_res_lists[i][-REMAIN_LENGTH:]

        # 计数器，用于投票
        counter = Counter(cor_res_lists[i])
        relative_dis[i] = counter.most_common(1)[0][0] * SOUND_V

    # TDOA的数学计算
    S = np.linalg.inv(M.T @ M) @ M.T @ relative_dis
    keep_theta[0] = np.arctan2(S[1], S[0])
    return keep_theta[0]

###############################################################

audio = pyaudio.PyAudio()

# 打开音频流
stream = audio.open(
            rate=RESPEAKER_RATE,
            format=audio.get_format_from_width(RESPEAKER_WIDTH),
            channels=RESPEAKER_CHANNELS,
            input=True,
            input_device_index=RESPEAKER_INDEX,)

# 初始化图形和点
plt.ion()  # 开启交互模式
fig, ax = plt.subplots(figsize=(5, 5))
point, = ax.plot([], [], 'o', color='gold')     # 圆点
ax.set_xlim(0, WIDTH)
ax.set_ylim(0, WIDTH)
circle = plt.Circle((WIDTH / 2, WIDTH / 2), RADIUS, color='lightgreen', fill=False)
ax.add_artist(circle)

try:
    print("Recording...")
    while True:
        # 读取音频数据
        data = stream.read(CHUNK, exception_on_overflow=False)
        # 将数据转换为 NumPy 数组
        audio_chunk = np.frombuffer(data, dtype=np.int16)
        
        # 将数据分配到各个声道
        for i in range(USED_CHANNELS):
            audio_data[i].extend(audio_chunk[i+1::RESPEAKER_CHANNELS].tolist())

        # 如果已经累积了足够多段了，则可以进行计算处理   
        if len(audio_data[0]) > SEG_LEN * SEG_NUM:
            for i in range(USED_CHANNELS):
                cor_data[i] = audio_data[i][0:SEG_LEN * SEG_NUM]
                audio_data[i] = audio_data[i][SEG_LEN * SEG_NUM - OVERLAP:]
            theta = tdoa()

            # 更新绘图
            x = WIDTH / 2 - RADIUS * np.cos(theta - np.pi / 2)
            y = WIDTH / 2 - RADIUS * np.sin(theta - np.pi / 2)
            point.set_data(x, y)
            plt.draw()
            plt.pause(0.1)  # 暂停以便于观察

            if not plt.fignum_exists(fig.number):
                print("Recording stopped.")
                break

except KeyboardInterrupt:
    print("Recording stopped.")

# 停止和关闭流
stream.stop_stream()
stream.close()
audio.terminate()

plt.ioff()  # 关闭交互模式
plt.show()  # 显示最终图形