当前位置：首页 > news >正文

从麦克风阵列到声源坐标：手把手实现Python版SRP-PHAT定位（含代码）

news 2026/4/22 2:49:07

从麦克风阵列到声源坐标：手把手实现Python版SRP-PHAT定位（含代码）

在智能音箱、会议系统甚至机器人听觉领域，声源定位技术正悄然改变人机交互的方式。想象一下，当你说出"打开客厅灯"时，设备不仅能理解指令，还能准确判断声音来自沙发左侧——这种空间感知能力的核心，正是我们今天要实现的SRP-PHAT算法。不同于传统理论综述，本文将带你在Python环境中构建完整的声源定位流水线，从麦克风信号采集到三维空间坐标输出，每个步骤都配有可运行的代码片段。

1. 环境搭建与数据准备

1.1 硬件选型与配置

市面上的USB麦克风阵列如ReSpeaker 4-Mic Array（约$59）或MiniDSP UMA-8（约$199）都是理想的实验设备。以ReSpeaker为例，其四麦克风呈环形排列，直径10cm，适合桌面级应用。连接电脑后，通过pyaudio库可快速验证设备是否正常工作：

import pyaudio p = pyaudio.PyAudio() for i in range(p.get_device_count()): dev = p.get_device_info_by_index(i) print(f"{i}: {dev['name']} (输入通道:{dev['maxInputChannels']})"

1.2 软件依赖安装

建议使用conda创建独立环境，避免库版本冲突：

conda create -n srp_env python=3.8 conda activate srp_env pip install numpy scipy librosa matplotlib pyaudio

对于需要GPU加速的场景，可额外安装cupy库替代numpy进行矩阵运算。

1.3 模拟数据生成

当没有物理麦克风阵列时，可以通过仿真数据验证算法。以下代码生成一个位于(1.5, 0.8, 0.6)米处的声源信号，被四个虚拟麦克风接收：

import numpy as np from scipy.signal import chirp fs = 44100 # 采样率 duration = 2 # 秒 t = np.linspace(0, duration, int(fs*duration), endpoint=False) source_signal = chirp(t, f0=100, f1=8000, t1=duration, method='logarithmic') # 麦克风位置（单位：米） mic_positions = np.array([[0, 0, 0], [0.1, 0, 0], [0, 0.1, 0], [0.1, 0.1, 0]]) source_pos = np.array([1.5, 0.8, 0.6])

2. SRP-PHAT核心算法实现

2.1 广义互相关计算（GCC-PHAT）

时延估计是定位的基础，PHAT加权能有效抑制混响影响：

def gcc_phat(sig1, sig2, fs, max_tau=None): n = len(sig1) + len(sig2) - 1 fft_size = 2**np.ceil(np.log2(n)).astype(int) # 计算互功率谱 S1 = np.fft.fft(sig1, fft_size) S2 = np.fft.fft(sig2, fft_size) R = S1 * np.conj(S2) # PHAT加权 R_phat = R / (np.abs(R) + 1e-10) # 避免除以零 # 逆变换得到时域相关 cc = np.fft.ifft(R_phat) cc = np.fft.fftshift(cc) # 时延轴 if max_tau: max_shift = int(max_tau * fs) cc = cc[fft_size//2 - max_shift : fft_size//2 + max_shift + 1] delays = np.arange(-max_shift, max_shift + 1) / fs else: delays = np.arange(-(fft_size//2), fft_size//2 + 1) / fs return cc, delays

2.2 空间网格构建策略

搜索空间的划分直接影响计算效率和定位精度。对于3米×3米×2米的房间，建议初始采用0.1米分辨率：

def build_search_grid(x_range=(-1.5, 1.5), y_range=(-1.5, 1.5), z_range=(0, 2), step=0.1): x = np.arange(x_range[0], x_range[1] + step, step) y = np.arange(y_range[0], y_range[1] + step, step) z = np.arange(z_range[0], z_range[1] + step, step) xx, yy, zz = np.meshgrid(x, y, z) grid = np.vstack((xx.flatten(), yy.flatten(), zz.flatten())).T return grid

实际应用中可采用多级网格策略：先用0.3米粗网格定位大致区域，再在目标区域使用0.05米精细网格

2.3 功率映射与峰值检测

计算每个网格点的SRP值并可视化结果：

def srp_phat(mic_signals, mic_positions, grid, fs, speed_of_sound=343): n_mics = mic_positions.shape[0] srp = np.zeros(grid.shape[0]) for i in range(n_mics): for j in range(i+1, n_mics): # 计算理论时延 tau_ij = (np.linalg.norm(grid - mic_positions[i], axis=1) - np.linalg.norm(grid - mic_positions[j], axis=1)) / speed_of_sound # 计算GCC-PHAT cc, delays = gcc_phat(mic_signals[i], mic_signals[j], fs) # 插值获取对应时延的相关系数 cc_interp = np.interp(tau_ij, delays, np.real(cc)) srp += cc_interp return srp # 可视化结果 def plot_srp_map(grid, srp_values, elev=30, azim=45): import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') sc = ax.scatter(grid[:,0], grid[:,1], grid[:,2], c=srp_values, cmap='viridis', s=20) ax.view_init(elev=elev, azim=azim) plt.colorbar(sc, label='SRP Value') plt.show()

3. 性能优化技巧

3.1 计算加速方案

原始SRP-PHAT的计算复杂度为O(N²M)，其中N是麦克风对数，M是网格点数。通过以下方法可显著提升速度：

矩阵化运算：替换嵌套循环为广播操作

# 优化后的时延矩阵计算 diff_vectors = grid[:, np.newaxis, :] - mic_positions # (M,4,3) distances = np.linalg.norm(diff_vectors, axis=2) # (M,4) tau_matrix = (distances[:, :, np.newaxis] - distances[:, np.newaxis, :]) / speed_of_sound # (M,4,4)

Numba加速：对关键函数添加JIT编译

from numba import jit @jit(nopython=True) def compute_srp_numba(tau_matrix, cc_values, delays): # 实现numba优化版本 ...

3.2 混响环境下的改进

强混响会导致定位偏移，可通过以下策略缓解：

频带选择：只保留500Hz-4kHz语音主要能量频段

def bandpass_filter(signal, fs, lowcut=500, highcut=4000): from scipy.signal import butter, filtfilt nyq = 0.5 * fs b, a = butter(4, [lowcut/nyq, highcut/nyq], btype='band') return filtfilt(b, a, signal)

峰值增强：对GCC-PHAT结果进行非线性放大

cc_enhanced = np.sign(cc) * np.abs(cc)**1.5

3.3 实时处理框架

构建可处理音频流的类结构：

class RealTimeSRP: def __init__(self, mic_positions, fs=44100, chunk_size=2048): self.mic_positions = mic_positions self.fs = fs self.chunk_size = chunk_size self.buffer = np.zeros((len(mic_positions), chunk_size * 3)) def update(self, new_frames): # 更新环形缓冲区 self.buffer = np.roll(self.buffer, -self.chunk_size, axis=1) self.buffer[:, -self.chunk_size:] = new_frames def locate(self, grid): # 使用缓冲区最新数据计算位置 current_data = self.buffer[:, -self.chunk_size*2:] srp = srp_phat(current_data, self.mic_positions, grid, self.fs) return grid[np.argmax(srp)]

4. 实战案例与问题排查

4.1 典型问题解决方案

问题现象	可能原因	解决方案
定位点随机跳动	背景噪声过大	增加VAD语音活动检测
定位偏向墙面	强混响干扰	应用峰值增强和频带过滤
计算速度过慢	网格分辨率过高	采用两级网格搜索策略
垂直方向误差大	麦克风共面布置	增加高度方向麦克风