基于到达时延差估计的定位法。 TDOA(time difference of arrival)是先后估计声源到达不同麦克风的时延差,通过时延来计算距离差,再利用距离差和麦克风阵列的空间几何位置来确定声源的位置。可分为TDOA估计(估计信号到达各麦克风的时间差)和TDOA定位(运用几何关系确定声源位置)两步。
x
1
=
s
1
(
t
)
+
n
1
(
t
)
x_1=s_1(t)+n_1(t)
x1=s1(t)+n1(t) (1)
x
2
=
α
s
1
(
t
−
D
)
+
n
2
(
t
)
x_2=αs_1(t-D)+n_2(t)
x2=αs1(t−D)+n2(t) (2)
其中,
s
(
t
)
s(t)
s(t) 是声音信号,
n
1
(
t
)
、
n
2
(
t
)
n_1(t)、n_2(t)
n1(t)、n2(t)是两个声音传感器检测噪声。 三者是稳定的随机过程,且互不相关。
计算
x
1
x_1
x1 与
x
2
x_2
x2 的互相关函数:
R
x
1
x
2
(
τ
)
=
E
[
x
1
(
t
)
∗
x
2
(
t
−
τ
)
]
R_{x_1x_2}(τ)=E [ x_1(t)*x_2(t-τ) ]
Rx1x2(τ)=E[x1(t)∗x2(t−τ)] (3)
R
^
x
1
x
2
(
τ
)
=
1
T
−
τ
∫
τ
T
x
1
(
t
)
x
2
(
t
−
τ
)
d
t
\hat R_{x_1x_2}(τ)=\frac{1}{T-τ}\int_τ^T {x_1(t)x_2(t-τ)} \,{\rm d}t
R^x1x2(τ)=T−τ1∫τTx1(t)x2(t−τ)dt (4)
其中估计的时延
D
D
D 为互相关函数值达到最大值时取得的
τ
τ
τ 值,即:
D
^
=
a
r
g
m
a
x
τ
R
x
1
x
2
(
τ
)
\hat D=\underset {\ τ}{argmax} R_{x_1x_2}(τ)
D^=τargmaxRx1x2(τ) (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);
一个三维麦克风阵列,麦克风分别为
m
i
(
i
=
0
、
1
⋅
⋅
⋅
n
)
m_i(i=0、1···n)
mi(i=0、1⋅⋅⋅n),声源
S
S
S 符合近场模型。现以麦克风
m
0
m_0
m0 为原点,如图建立直角坐标系。推导公式之前,需要先确定以下概念:
序号
概念
字符
1
麦克风的坐标
m
i
m_i
mi
i
∈
[
0
,
n
]
i \in[0,n]
i∈[0,n]
2
声源的估计坐标
S
S
S
3
声源
s
s
s 到
m
i
m_i
mi 与
m
1
m_1
m1 的估计距离差
d
^
i
\hat d_i
d^i
4
麦克风
m
i
m_i
mi 到原点的的距离
R
i
R_i
Ri
5
声源
s
s
s 到原点的的距离
R
s
R_s
Rs
如上图,根据三角形
m
0
m
i
s
m_0 m_i s
m0mis,由余弦定理有:
(
R
s
+
d
i
)
2
=
R
i
2
−
2
m
i
T
s
+
R
s
2
(R_s+d_i)^2=R_i^2-2m_i^Ts+R_s^2
(Rs+di)2=Ri2−2miTs+Rs2 (12)
公式 (12) 展开后得到:
R
i
2
−
2
m
i
T
s
−
2
R
s
d
i
−
d
i
2
=
0
R_i^2-2m_i^Ts-2R_sd_i-d_i^2=0
Ri2−2miTs−2Rsdi−di2=0 (13)
由于
d
i
d_i
di 是通过估计的时延得到的,与实际值之间会有一个偏差,因此公式 (10) 不为零,可写为:
ε
=
(
R
i
2
−
d
i
2
)
−
2
m
i
T
s
−
2
R
s
d
i
ε=(R_i^2-d_i^2)-2m_i^Ts-2R_sd_i
ε=(Ri2−di2)−2miTs−2Rsdi (14)
此时已经得到目标值的误差函数,使用最小二乘法求解,问题可以转化为:估计声源坐标
s
(
x
,
y
,
z
)
s(x,y,z)
s(x,y,z),使得最终的误差平方和最小,即:
a
r
g
m
i
n
s
(
x
,
y
,
z
)
∑
i
=
1
n
[
(
R
i
2
−
d
i
2
)
−
2
m
i
T
s
−
2
R
s
d
i
]
2
\underset {\ s(x,y,z)}{argmin}\displaystyle\sum_{i=1}^{n} [(R_i^2-d_i^2)-2m_i^Ts-2R_sd_i]^2
s(x,y,z)argmini=1∑n[(Ri2−di2)−2miTs−2Rsdi]2 (15)
将公式 (14) 写成矩阵形式:
ε
=
2
R
s
D
^
+
2
M
s
−
δ
ε=2R_s\hat D+2Ms-\delta
ε=2RsD^+2Ms−δ (16)
其中,
M
=
[
m
1
m
2
m
3
.
.
.
.
m
n
]
=
[
x
1
y
1
z
1
x
2
y
2
z
2
x
3
y
3
z
3
.
.
.
.
.
.
.
.
.
x
n
y
n
z
n
]
M=\begin{bmatrix}m_1\\m_2\\m_3\\....\\m_n\end{bmatrix}=\begin{bmatrix}x_1&y_1&z_1\\x_2&y_2&z_2\\x_3&y_3&z_3\\...&...&...\\x_n&y_n&z_n\end{bmatrix}
M=⎣⎡m1m2m3....mn⎦⎤=⎣⎡x1x2x3...xny1y2y3...ynz1z2z3...zn⎦⎤
D
^
=
[
d
^
1
d
^
2
d
^
3
.
.
.
d
^
n
]
\hat D=\begin{bmatrix}\hat d_1\\\hat d_2\\\hat d_3\\\ ... \\\hat d_n\end{bmatrix}
D^=⎣⎡d^1d^2d^3...d^n⎦⎤
δ
=
[
R
1
2
−
d
^
1
2
R
2
2
−
d
^
2
2
R
3
2
−
d
^
3
2
.
.
.
.
.
.
.
.
.
.
R
n
2
−
d
^
n
2
]
\delta=\begin{bmatrix}R_1^2-\hat d_1^2\\R_2^2-\hat d_2^2\\R_3^2-\hat d_3^2\\..........\\R_n^2-\hat d_n^2\end{bmatrix}
δ=⎣⎡R12−d^12R22−d^22R32−d^32..........Rn2−d^n2⎦⎤
公式 (16) 可以简化为:
ε
=
A
μ
−
b
ε=A\mu-b
ε=Aμ−b (17)
其中,
A
=
[
M
D
^
]
A=\begin{bmatrix}M&\hat D\end{bmatrix}
A=[MD^]
μ
=
[
s
R
s
]
\mu=\begin{bmatrix}s\\R_s\end{bmatrix}
μ=[sRs]
b
=
1
2
δ
b=\frac{1}{2}\delta
b=21δ
公式 (17) 的最小二乘解可以表示为:
μ
^
=
(
A
T
A
)
−
1
A
T
b
\hat\mu=(A^TA)^{-1}A^Tb
μ^=(ATA)−1ATb (18)
关于矩阵形式下最小二乘解推导过程可以参考下面的链接(点击跳转),此处不详细讲解。
公式 (18) 即为计算结果,下面对结果进行进一步化简:
定义沿
A
A
A 到
D
^
\hat D
D^ 的投影矩阵为:
P
A
=
D
^
D
^
T
D
^
T
D
^
P_A=\frac{\hat D \hat D^T}{\hat D^T \hat D }
PA=D^TD^D^D^T (19)
沿
D
^
\hat D
D^ 到
A
A
A 的投影矩阵即为:
P
D
=
I
−
P
A
=
I
−
D
^
D
^
T
D
^
T
D
^
P_D=I-P_A=I-\frac{\hat D \hat D^T}{\hat D^T \hat D }
PD=I−PA=I−D^TD^D^D^T(
I
I
I 为单位矩阵) (20)
根据投影矩阵的性质,可以得到:
A
=
P
D
[
M
0
]
A=P_D\begin{bmatrix}M&0\end{bmatrix}
A=PD[M0] (21)
将公式 (21) 带入 (17) 中得到:
ε
=
P
D
M
s
−
b
ε=P_D M s-b
ε=PDMs−b (22)
公式(22)的最小二乘解即为最终简化结果。
最终简化结果:
s
=
(
M
T
P
D
M
)
−
1
M
T
P
D
b
s=( M^T P_D M )^{-1}M^TP_Db
s=(MTPDM)−1MTPDb
M
M
M只与麦克风阵列的个数与排列有关,
P
D
P_D
PD 与
b
b
b 在计算出时延估计的距离差 (
d
i
d_i
di) 后也可以求得。