[ 猩猩也会用 ] 如何使用Python实现4x4麦克风阵列的二维波束形成及实时Heatmap呈现

主要重点

• 麦克风阵列配置与信号採集 - 确保音讯数据的正确收集与前处理
• 二维波束形成算法实现 - 使用延迟与总和方法增强目标方向信号
• 实时Heatmap生成与可视化 - 利用Matplotlib动态呈现声源分布

1. 麦克风阵列配置与信号採集

1.1 硬体配置

您的麦克风阵列为4x4布局，共16个麦克风，每个麦克风在水平与垂直方向上的间隔均为42mm。这样的配置允许在多个方向上精确地捕捉音讯信号，为后续的波束形成奠定基础。麦克风的位置对于波束形成算法的準确性至关重要，因此必须正确设定每个麦克风的坐标。

1.2 音讯数据採集

使用Python进行音讯数据的实时採集，可以利用

pyaudio

或

sounddevice

等库。以下范例展示如何初始化音频流并读取数据：

import pyaudio
import numpy as np

# 音讯参数设定
CHUNK = 1024 # 每次读取的样本数
RATE = 44100 # 採样率 (Hz)
CHANNELS = 16 # 麦克风通道数

# 初始化pyaudio
p = pyaudio.PyAudio()
stream = p.open(format=pyaudio.paInt16,
channels=CHANNELS,
rate=RATE,
input=True,
frames_per_buffer=CHUNK)

1.3 麦克风阵列几何结构

定义麦克风的位置对于计算延迟至关重要。以下是麦克风阵列的几何结构设定：

mic_positions = np.array([
[0, 0], [0, 0.042], [0, 0.084], [0, 0.126],
[0.042, 0], [0.042, 0.042], [0.042, 0.084], [0.042, 0.126],
[0.084, 0], [0.084, 0.042], [0.084, 0.084], [0.084, 0.126],
[0.126, 0], [0.126, 0.042], [0.126, 0.084], [0.126, 0.126]
]) # 单位：米

1.4 麦克风性能与校準

为了确保波束形成的準确性，所有麦克风必须经过校準。校準过程包括检查每个麦克风的灵敏度、频率响应和相位一致性。任何麦克风的不一致性都可能导致波束形成结果的偏差。

1.5 资料前处理

在进行波束形成之前，需要对音讯数据进行前处理，包括去除直流分量、应用高通或低通滤波器以消除不需要的频率成分，以及对多通道信号进行标準化处理，以确保各通道信号的均衡。

2. 二维波束形成算法实现

2.1 波束形成基本原理

波束形成（Beamforming）是一种信号处理技术，用于通过麦克风阵列中的多个麦克风信号，通过调整各个信号的相位与幅度，来增强特定方向的声音信号，同时抑制其他方向的噪声。这样可以精确地定位声源方向并提升信噪比。

2.2 延迟与总和（Delay-and-Sum）波束形成

延迟与总和是一种基本的波束形成方法。其步骤如下：

2.3 二维波束形成实现

实现二维波束形成需要在水平方向（Azimuth）和垂直方向（Elevation）上对每个目标角度进行计算。以下程式码展示了如何计算延迟并进行波束形成：

def beamforming(audio_data, mic_positions, angle_h, angle_v, speed_of_sound=343):
"""
执行二维波束形成
:param audio_data: (通道数, 样本数) 的音讯数据
:param mic_positions: (通道数, 2) 的麦克风位置数组
:param angle_h: 水平角度 (度)
:param angle_v: 垂直角度 (度)
:param speed_of_sound: 声速 (默认343 m/s)
:return: 波束形成后的信号
"""
delays = np.zeros(len(mic_positions))
for i, pos in enumerate(mic_positions):
# 计算每个麦克风的延迟时间
delay = (pos[0] * np.sin(np.radians(angle_h)) + pos[1] * np.sin(np.radians(angle_v))) / speed_of_sound
delays[i] = delay

# 对音讯数据进行延迟补偿
delayed_data = np.zeros_like(audio_data)
for i in range(len(mic_positions)):
sample_shift = int(delays[i] * RATE)
delayed_data[i] = np.roll(audio_data[i], sample_shift)

# 合併延迟后的音讯数据
beamformed_signal = np.sum(delayed_data, axis=0)
return beamformed_signal

2.4 信号前处理与噪声抑制

为了提升波束形成的效果，可以在前处理阶段进行噪声抑制，例如使用滤波器来去除背景噪声，或者应用自适应波束形成算法来动态调整波束以应对不同的噪声环境。下方程式码展示了如何应用高通滤波器去除低频噪声：

from scipy.signal import butter, lfilter

def highpass_filter(data, cutoff=1000, fs=RATE, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
# 设计高通滤波器
b, a = butter(order, normal_cutoff, btype=\'high\', analog=False)
y = lfilter(b, a, data)
return y

# 应用高通滤波器
filtered_data = highpass_filter(audio_data)

2.5 综合波束形成流程

以下流程图展示了二维波束形成的整体流程：

3. 实时Heatmap生成与可视化

3.1 Heatmap的基本概念

Heatmap（热图）是一种可视化工具，用于表示多维数据的强度或频率。在波束形成中，Heatmap可用于展示声源在不同方向上的强度分布，帮助直观理解声源定位结果。热图中，每个像素点代表一个特定的方向，颜色深浅则表示该方向上的信号强度。

3.2 使用Matplotlib生成实时Heatmap

Matplotlib是一个强大的Python图形库，适合用来生成和更新Heatmap。以下程式码展示如何初始化并实时更新Heatmap：

import matplotlib.pyplot as plt

# 初始化Heatmap
fig, ax = plt.subplots()
heatmap = ax.imshow(np.zeros((180, 180)), cmap=\'hot\', origin=\'lower\', extent=[-90, 90, -90, 90])
plt.colorbar(heatmap)
plt.xlabel(\'水平角度 (度)\')
plt.ylabel(\'垂直角度 (度)\')
plt.title(\'波束形成热图\')

# 实时更新Heatmap
while True:
# 读取音讯数据
audio_data = np.frombuffer(stream.read(CHUNK), dtype=np.int16).reshape(-1, CHANNELS).T

# 初始化Heatmap数据
heatmap_data = np.zeros((180, 180))

# 计算每个角度的讯号强度
for angle_h in range(-90, 90):
for angle_v in range(-90, 90):
beamformed_signal = beamforming(audio_data, mic_positions, angle_h, angle_v)
heatmap_data[angle_v + 90, angle_h + 90] = np.sum(np.abs(beamformed_signal))

# 更新Heatmap
heatmap.set_data(heatmap_data)
plt.pause(0.1) # 更新间隔为0.1秒
ax.cla()
heatmap = ax.imshow(heatmap_data, cmap=\'hot\', origin=\'lower\', extent=[-90, 90, -90, 90])
plt.colorbar(heatmap)
plt.xlabel(\'水平角度 (度)\')
plt.ylabel(\'垂直角度 (度)\')
plt.title(\'波束形成热图\')

3.3 提升Heatmap生成效率

由于波束形成计算量庞大，实时生成Heatmap时可能会出现延迟。以下是一些优化建议：

• 并行运算：利用多执行绪或多处理器对不同角度的计算进行并行化。可以使用Python的或库来实现。

  multiprocessing

  threading

• 算法优化：使用NumPy等高效数值运算库进行矩阵计算，或考虑使用Cython进行C语言级别的优化，提升运算速度。

• 减少角度范围或增加间隔：降低Heatmap的解析度（例如，每隔2度计算一次），以减少计算量，同时仍能保留足够的定位精度。

3.4 Heatmap美化与增强

为了让Heatmap更加清晰易读，可以进行以下美化处理：

• 使用对数刻度：将信号强度转换为对数刻度，以更好地展示强度差异。
• 添加标记与注释：在Heatmap上标记主要声源方向，或添加数值标注。
• 调整色彩映射：根据需求选择合适的色彩映射方案，例如从冷到热的渐变，以突出重要区域。

4. 完整的实现示例

4.1 完整程式码范例

以下是一个结合所有部分的完整Python程式码示例，展示如何从麦克风阵列的信号採集到实时生成Heatmap：

import numpy as np
import pyaudio
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter
import threading

# 音讯参数设定
CHUNK = 1024
RATE = 44100
CHANNELS = 16

# 麦克风阵列位置设定（单位：米）
mic_positions = np.array([
[0, 0], [0, 0.042], [0, 0.084], [0, 0.126],
[0.042, 0], [0.042, 0.042], [0.042, 0.084], [0.042, 0.126],
[0.084, 0], [0.084, 0.042], [0.084, 0.084], [0.084, 0.126],
[0.126, 0], [0.126, 0.042], [0.126, 0.084], [0.126, 0.126]
])

# 初始化pyaudio
p = pyaudio.PyAudio()
stream = p.open(format=pyaudio.paInt16,
channels=CHANNELS,
rate=RATE,
input=True,
frames_per_buffer=CHUNK)

# 高通滤波器
def highpass_filter(data, cutoff=1000, fs=RATE, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype=\'high\', analog=False)
y = lfilter(b, a, data)
return y

# 波束形成函数
def beamforming(audio_data, mic_positions, angle_h, angle_v, speed_of_sound=343):
delays = np.zeros(len(mic_positions))
for i, pos in enumerate(mic_positions):
delays[i] = (pos[0] * np.sin(np.radians(angle_h)) + pos[1] * np.sin(np.radians(angle_v))) / speed_of_sound
delayed_data = np.zeros_like(audio_data)
for i in range(len(mic_positions)):
sample_shift = int(delays[i] * RATE)
delayed_data[i] = np.roll(audio_data[i], sample_shift)
beamformed_signal = np.sum(delayed_data, axis=0)
return beamformed_signal

# Heatmap更新函数
def update_heatmap():
while not stop_event.is_set():
# 读取音讯数据
audio_data = np.frombuffer(stream.read(CHUNK), dtype=np.int16).reshape(-1, CHANNELS).T
# 应用高通滤波器
audio_data = highpass_filter(audio_data)

# 初始化Heatmap数据
heatmap_data = np.zeros((180, 180))

# 计算每个角度的讯号强度
for angle_h in range(-90, 90):
for angle_v in range(-90, 90):
beamformed_signal = beamforming(audio_data, mic_positions, angle_h, angle_v)
heatmap_data[angle_v + 90, angle_h + 90] = np.sum(np.abs(beamformed_signal))

# 更新Heatmap
ax.clear()
im = ax.imshow(heatmap_data, cmap=\'hot\', origin=\'lower\', extent=[-90, 90, -90, 90])
plt.colorbar(im, ax=ax)
plt.xlabel(\'水平角度 (度)\')
plt.ylabel(\'垂直角度 (度)\')
plt.title(\'波束形成热图\')
plt.pause(0.001) # 短暂暂停以更新图形

# 初始化Heatmap
fig, ax = plt.subplots()
heatmap = ax.imshow(np.zeros((180, 180)), cmap=\'hot\', origin=\'lower\', extent=[-90, 90, -90, 90])
plt.colorbar(heatmap, ax=ax)
plt.xlabel(\'水平角度 (度)\')
plt.ylabel(\'垂直角度 (度)\')
plt.title(\'波束形成热图\')

# 启动Heatmap更新执行绪
stop_event = threading.Event()
heatmap_thread = threading.Thread(target=update_heatmap)
heatmap_thread.start()

try:
plt.show(block=True)
except KeyboardInterrupt:
pass
finally:
stop_event.set()
heatmap_thread.join()
stream.stop_stream()
stream.close()
p.terminate()

4.2 执行步骤说明

4.3 执行与测试

在执行此程式前，请确保所有麦克风均已正确连接并能够稳定运行。运行程式后，应能实时看到波束形成的Heatmap呈现声源方向及强度分布。在测试过程中，您可以移动声源的位置，观察Heatmap的变化，以验证系统的準确性。

4.4 程式码优化建议

为了提升实时运行的效率，建议对程式码进行以下优化：

• **并行运算：**利用Python的或库，将波束形成计算和Heatmap更新分配到不同的执行绪或进程。

  multiprocessing

  threading

• **算法优化：**使用NumPy进行向量化运算，减少for迴圈的使用，以提高计算速度。

• **减少计算范围：**可以考虑将角度范围调整为[-90, 90)度，并增加角度间隔（例如，每2度计算一次），以减少计算量。

5. 优化与扩展建议

5.1 性能优化

为了提升实时运行的效率，可以考虑以下方法：

• **多执行绪处理：**将音讯数据採集与Heatmap生成分离到不同的执行绪，减少延迟并提高处理效率。
• **使用快速数值计算库：**利用NumPy的向量化运算功能，替代传统的for迴圈，提升数据处理速度。
• **减少计算范围或增加角度间隔：**通过降低Heatmap的解析度，或增大角度间隔（如每2度计算一次），以减少计算量，从而提高实时性能。

5.2 增强算法

可以进一步提升波束形成的準确性和对噪声的抑制能力：

• **自适应波束形成：**根据环境动态调整波束形状，以更有效地抑制噪声。这需要实现更复杂的算法，如最小方差无偏估计（MVDR）。
• **多频率处理：**结合不同频率的信号来提升定位準确度，特别是在多声源环境中。
• **信号预处理：**应用高通或低通滤波器，进行降噪处理，提升信号质量，从而提高波束形成效果。

5.3 扩展功能

除了Heatmap外，还可以考虑以下扩展功能：

• **音源跟踪：**在Heatmap基础上增加音源的移动轨迹追踪，实现动态声源定位。
• **多声源辨识：**提升算法以同时定位多个声源，适应复杂的声环境。
• **3D热图：**将Heatmap扩展到三维，以包含更多维度的数据，提供更加直观的声源分布视觉化。

6. 结论

本文详细介绍了如何使用Python实现4x4麦克风阵列的二维波束形成，并将结果以实时Heatmap呈现。通过步骤化的指导与范例程式码，您可以有效地设置麦克风阵列、进行信号处理及可视化，进而实现精确的声源定位。进一步的优化与扩展将使系统在更复杂的环境中表现更佳，满足不同应用需求。