我是后背
要想实现一个声源的定位，我们需要：
- 声源的到达方向 (Direction of Arrival，DOA) ，主要计算其的方位角和俯仰角
- 声源的距离估计
目前我们的重点先放在DOA上，对于DOA估计，主要有三个方法：
- 基于到达时间差（TDOA）：计算不同阵元之间接收到的信号的时延差，根据该时延差即可推断出声源位置
- 基于波束形成的方法：进行不同角度遍历扫描，分别进行波束形成（即假设时延差），功率最大处即为声源位置
- 基于高分辨谱估计：然后采用特征分解操作，从协方差矩阵中提取出信号子空间。之后，利用空间谱 估计技术对该子空间进行分析，最终得出声源方向的估计结果。
根据定义我们就可以看出，基于相对时延估计的技术点主要在于如何计算出时延差；而波束形成方法则更简单暴力，技术点在于使用什么波束形成方法更好
本章我们将对基于相对时延估计的方法进行展开，其可以分为两部分：
其中，时间差计算为核心
声场模型
根据麦克风阵列和声源距离的远近，我们可以将声场模型分为：近场模型和远场模型，在近场模型下我们将声波看作为球面波，其主要考虑麦克风阵列各阵元接收到信号的幅度差; 而远场模型将声波看作为平面波，忽略阵元接收信号间的幅度差，近似认为各阵元接收信号之间为简单的时延关系
评判模型的标准为声源到阵列参考阵元的判据：
当绝对距离大于该值时即为远场模型
由公式可知，当选择的频率越大，远场判据越大，即该频率越容易表现出近场特性，因此若要使用更简单的远场模型，我们选择的频率应该偏小一点，工程中常选取主要频带上的最大频率值（关于距离这一块后面的章节会细说）
例子：当我们阵列直径为0.2m时，瑞利距离则约为0.47m
如果对近场模型使用了基于远场模型的算法，则会导致导向矢量的失配，最终导致结果的不准确
传统互相关法
假设参考阵元信号为x1，则我们可以直接通过互相关函数计算阵元x2与x1的时延差：
当互相关值最大时，对应的tau即为时延差。然而，在实际室内环境中，由于声波反射（混响）和环境噪声的影响，传统互相关函数的峰值往往会变得非常平缓，甚至出现多个伪峰，导致估计彻底失效。
广义互相关法GCC
由维纳-辛钦定理可知，随机信号的相关函数和功率谱密度函数服从一对傅里叶变换的关系（即互相关时域函数在频域的表达就是互功率谱密度），功率谱密度函数可以表示为：
这里我们就建立了频域表达与互相关函数（时域）的等式。为了增强互相关函数峰值，我们引入一个加权函数ψ来增强信号中有用的频率部分，在频域增强完后再转到时域：
对于这个加权函数的设计，则就是GCC各种分支的由来了
互功率谱密度是指在f处的共有的能量密集度和在f下两信号的相位差，而我们现在只需要通过相位差来得到时延值，因此，可以直接令加权函数为功率谱密度的模（幅值）的倒数：
这样的话，互相关函数就是一个仅与功率谱密度的相位相关的函数了：
也就是说，PHAT本身就是指通过相位差来计算时延值，当t取到t0时，互相关函数就是最大值
py实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
| def gcc_get_tdoa(self, freq_range=None):
if freq_range is None:
freq_range = [1000, 2000]
mask = (self.freqbin >= freq_range[0]) & (self.freqbin <= freq_range[1])
X_ref = self.data[0, 0:1, :]
G = np.conjugate(X_ref) * self.data[0,:,:]
eps = 1e-12
G_phat = G / (np.abs(G) + eps)
G_filtered = G_phat * mask[None, None, :]
G_avg = np.mean(G_filtered, axis=0)
gcc_time = np.fft.irfft(G_avg, n=self.num_freqs, axis=-1)
gcc_time_shifted = np.fft.fftshift(gcc_time, axes=-1)
lags = np.arange(-self.num_freqs // 2, self.num_freqs // 2) / self.sampling_frequency
tdoa_indices = np.argmax(gcc_time_shifted, axis=-1)
tdoa = lags[tdoa_indices]
return tdoa, gcc_time_shifted, lags
|
这里简单实现了基于麦克风第一个通道的GCC-PHAT计算，这里需要注意几点：
np.conjugate(X_ref)：对参考通道进行共轭转置，及对应上面的公式中的-t0，这样晚于参考通道的信号的t就是正值- 带通滤波：若要进行窄带运算，必须要使用带通滤波，而不是简单的仅提取对应通道，因为后面要使用ifft
gcc_time_shifted：当ifft对象为普通信号时，其会正常还原时域波形；但如果对象为互功率谱，其对应的时域为互相关函数，包含正负时延，因此ifft的输出会首先输出正时延值，再输出负时延值（与常规顺序相反），gcc_time_shifted会根据0值位置，将正负互换lags = np.arange...刚刚提到，互相关函数的完整物理定义域是有正负的，具体其实就是[-0.5T,0.5T]（因为任何傅里叶变换得到的都是一个周期函数/序列，ifft得到的就是各频段周期函数的叠加，而一个长度为N的周期序列，任何索引n都可以映射到唯一区间[-0.5T,0.5T]），因此gcc_time_shifted对应的时间索引是一定对称的，而lags就是把这个索引值算出来，因为gcc_time_shifted没有包含每个索引对应的时延值，它代表的是每个时延值上对应的互相关值
优缺点
优点：由于仅通过相位差来得到时延值，因此该方法拥有极强的抗混响能力（混响对相位影响小）
缺点：当信号中存在噪声时，当某些频段上噪声为主要成分时，也就意为着此时得到的相位也是噪声的相位，而PHAT舍弃了幅值，因此也将放大噪声，最终导致互相关函数上真峰降低，并出现大量伪峰（因为互相关函数代表的是不同t上信号的相关程度，伪峰即为噪声与信号产生的相关）
后记，由于是对每个频段分别进行相位差计算，因此在噪声源存在的情况下几乎没办法准确定位声源（即使噪声比较微弱）
准确来说这个方法已经不算是直接计算DOA，但是跟GCC-PHAT强相关，所以也放这里了
SRP-PHAT先在空间中划分网格（候选声源点）。对于空间中任意一个假设的声源点，计算它到达所有麦克风对的理论时延。然后将所有麦克风对在这些理论时延处的 GCC-PHAT 值进行空间累加。
这和基于波束形成的方法非常像，只是波束形成的输出是声压，而SRP-PHAT的输出是互相关函数值：
显而易见，该方法要计算任意两个阵元的互相关值，计算量极大，即使预先计算出所有的时延，仍要计算大量次互相关，因此实际使用时一般配合一些简化计算量的优化：
由粗到精的搜索（Coarse-to-Fine Search）：先用极稀疏的网格扫描空间，锁定高能量的局部区域；再在这些候选区域内细化网格进行微观搜索。
随机区域收缩法（Stochastic Region Contraction, SRC）：基于动态优化思想，每次随机抽取空间样本点，根据响应功率不断收缩搜索空间收敛到最优解，能将计算量降低几个数量级。
GPU 并行加速：由于每个空间网格点的功率计算 彼此完全独立，属于典型的数据并行（SIMD）任务。在现代工程中，利用 CUDA 将网格映射分发给 GPU 执行，可以极其轻松地实现高分辨率的三维实时 SRP-PHAT 定位。
代码实现：
建立麦克风对tdoa
由于pair tdoa是srp-phat的主要计算量处，而其结果仅与阵列坐标和扫描角度有关，因此通常是预计算，在后续计算中直接读取保存的toda结果
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
| def build_pair_tdoa_table(
mic_struct,
phi_scan,
theta_scan,
sound_speed=343.0):
M = mic_struct.shape[0]
pairs = list(combinations(range(M), 2))
Npair = len(pairs)
pair_tdoa = np.zeros(
(
len(theta_scan),
len(phi_scan),
Npair
),
dtype=np.float32
)
for it, theta in enumerate(tqdm(theta_scan,desc="计算各麦克风pairs时延中")):
sin_t = np.sin(theta)
cos_t = np.cos(theta)
for ip, phi in enumerate(phi_scan):
u = np.array([
sin_t*np.cos(phi),
sin_t*np.sin(phi),
cos_t
])
tau = -mic_struct @ u / sound_speed
for k, (i, j) in enumerate(pairs):
pair_tdoa[it, ip, k] = (
tau[i] - tau[j]
)
pair_lag = np.round(pair_tdoa * 51200).astype(np.int32)
np.savez_compressed(
"pair_tdoa.npz",
pair_tdoa=pair_tdoa,
pairs=pairs,
pair_lag = pair_lag,
phi_scan=phi_scan,
theta_scan=theta_scan,
)
return 1
|
我们后续使用的主要是pair_lag（每对时延值对应的互相关函数采样索引，后面再说详细情况）和pairs（麦克风对索引）
计算麦克风对的互相关函数
对每个麦克风对的信号进行互相关计算
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
| def compute_gcc_phat(spec, pairs, N_fft):
eps = 1e-12
Npair = len(pairs)
gcc = np.zeros(
(Npair, N_fft),
dtype=np.float32
)
for k, (i, j) in enumerate(tqdm(pairs, desc="gcc计算中")):
cross = spec[:,i,:] * np.conj(spec[:,j,:])
cross /= np.abs(cross) + eps
avg_cross = np.mean(cross,axis=0)
gcc_time = np.fft.irfft(avg_cross)
gcc_time = np.real(
np.fft.fftshift(gcc_time)
)
gcc[k] = gcc_time
return gcc
|
得到是每个麦克风对对应的互相关函数
这个函数需要注意，输入的spec如果是单边谱，则使用irfft还原，如果是双边谱，就用ifft。如果这里没有对应，最后输出的功率谱上声源方向会出现缩放情况
SRP-PHAT扫描
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
| def srp_phat_fast(
gcc,
pair_lag):
Npair, Nlag = gcc.shape
center = Nlag // 2
idx = center + pair_lag
idx = np.clip(
idx,
0,
Nlag - 1
)
gcc_expand = gcc[
np.arange(Npair)[None, None, :],
idx
]
power_map = np.sum(
gcc_expand,
axis=2
)
return power_map
|
之前得到的pair_lag是根据0点计算的索引，也就是一个带正负的索引，所以在用之前要先进行idx = center + pair_lag转化为真正的[0,N_fft]的索引，用于索引gcc
SRP-PHAT就是直接用假想方向的tdoa来索引gcc，使用客观索引避免了噪声干扰，再通过每个麦克风对求和来增加稳定性；gcc-phat就是直接找gcc的最大值，因此更容易被噪声干扰
个人实际使用下来，SPR-PHAT通过查表方式和DAS的计算时长都还不错，接近的计算时间下SPR-PHAT精度甚至更高
DAS成像：
SRP-PHAT成像：
从TDOA得到DOA
本节默认在远场模型下计算（远场模型只能得到DOA，无法得到距离，想象声源在参考阵元射出的一条射线方向，而波就是以该射线为法线的平面）
先从二维说起，对于两个距离为d的阵元，对于平面波的入射角（与阵元连线的夹角），有：
c为声速，t为TODA，即可通过反三角求解θ
现在来到三维，方向就变为了方位角Φ和俯仰角θ(与xy平面夹角），与球坐标系定义相同，此时声源方向单位向量就可以表示为：
于是就有：
理论上只需要3阵元（也就是3对TDOA）就可以求解这个方程，当阵元大于3时，一般使用最小二乘法求最优解，对于平面阵列，有：
线性方程组可以紧凑地表示为：
最小二乘法就是让损失函数（Am-b)^2最小，即让该函数对m的偏导为0，化简后则有：
（上式为最小二乘法通用解），得到方向向量后根据三角函数关系即可得到具体角度
总结
基于TDOA的声源定位方法实现原理较为简单，速度也还算可以（SRP不算该方法），核心就是互相关函数和如何设计加权函数，普遍存在的缺点就是极容易被噪声干扰。