在昨天的推文中，讨论了利用声音测距的基本方法。在具体应用中，具体会碰到哪些问题，下面给出今天的一些测试结果。
实验装置和数据采集
数据时通过基于STM32F103RE的AD、DA1采样板来控制音响发送Chirp声音信号和采集MIC接收到的声音信号的。
实验中测试面包板
1. 直线滑轨定位
演示中所使用的滑轨长度大约为1米，滑块有效移动距离为0.9米。通过ZIGBEE无线控制命令控制滑块做直线均匀运动。
设置在滑轨上的声音传感器
实验中的Chirp声音信号是由单片机产生并通过DA输出给蓝牙音箱的播放。信号从数字缓冲区通过DA转换成模拟信号，声音信号经过AD转换为数字信号的速率为都设置为10kHz。输出和输入声音的缓存区的长度为2500，声波时间长度为0.25秒。
实际声音的长度为2048,。因此在输出和采样前后各有250个数据点（大约25ms）左右时静音的声音。下图给出了一个典型的发送声音信号和接收到的声音采集信号波形。
每个采集数据包里的数据波形
采集数据处理
在每个位置点都采集到相应的发送和接受数据，通过相关计算获得时间延迟。
1. 发送和接收数据之间的互相关系数计算
计算每个数据包里的发送信号与接收信号之间的相关系数，从而确定最大的峰值位置。下图给出了一个典型的发送和接受信号的相关结果曲线：
第一个数据包中发送和接收数据之间的相关系数
2. 相关系数的最大、最小位置与数据采样之间的关系
下图将每个数据包中相关系数的最大值点和最小值点的位置与数据采样点之间的绘制出来。随着数据序号的增多，对应接收声音的MIC位置原理音箱，相关系数的极值位置也在发生变化。基本上都是随着采用序列的增多而线性下降，反映了MIC的位置是线性原理声源的位置。
相关系数最大值最小值的位置与数据采用之间的关系
数据中最大值的位置基本上都是单调下降，只是在后期出现了略微的抖动。但是最小值的位置出现过较大的波动。这说明最大值的位置用于确定延迟时间比较稳定。
下图显示了相关系数最大值和最小值随着采样点位置之间的关系。随着采样点增多，音箱和麦克风之间的距离增加。相关系数的绝对值因为噪声的原因都会下降。整体上极大值点的相关系数呈现单调下降的趋势，而极小值点的数值变化呈现非单调的形态。这也从另外一个角度反映了为什么前面相关系数极小值点的位置出现抖动的原因。
相关系数最大值和最小值与采样位置之间的关系
3. 分析相关峰值位置的变化与采样距离之间的关系
对于相关系数最大值峰值点的变化从开始的 p 1 = p_1 = p1= 2497，一直变化到最后的 p 2 = p_2 = p2= 2473。变化的数值： p 1 − p 2 = p_1 - p_2 = p1−p2= 24。
由于实验中声音信号的采集频率
f
s
=
10
,
000
f_s = 10,000
fs=10,000Hz，所以上述峰值位置的变化所对应的时间延迟为：
Δ
t
=
p
1
−
p
2
f
s
=
24
10
×
1
0
3
=
2.4
(
m
s
)
\Delta t = {{p_1 - p_2 } \over {f_s }} = {{24} \over {10 \times 10^3 }} = 2.4\,\,\left( {ms} \right)
Δt=fsp1−p2=10×10324=2.4(ms)
如果考虑到当时的室温为25℃左右，对应的空气中声音传播的速度：
v
a
i
r
=
331
+
0.6
×
T
=
331
+
0.6
×
25
=
346
(
m
/
s
)
v_{air} = 331 + 0.6 \times T = 331 + 0.6 \times 25 = 346 \:\:(m/s)
vair=331+0.6×T=331+0.6×25=346(m/s)
那么前面 Δ t \Delta t Δt所对应的距离为： D = Δ t × v a i r = 2.4 × 1 0 − 3 × 346 = 0.83 ( m ) D = \Delta t \times v_{air} = 2.4 \times 10^{ - 3} \times 346 = 0.83\:\:(m) D=Δt×vair=2.4×10−3×346=0.83(m)
测量接收声音的MIC实际移动的距离为0.671米。这中间出现了 χ = 0.83 − 0.671 0.671 × 100 % = 23.7 % \chi = {{0.83 - 0.671} \over {0.671}} \times 100\% = 23.7\% χ=0.6710.83−0.671×100%=23.7%的误差。
具体为什么？现在还不可得而知。
测量实际接收声音的MIC移动的距离
声音延迟的分辨率与空间分辨率
通过相关运算可以获得声音传播的延迟，近而获得声源与接收麦克之间的距离。
由于声音信号是通过离散时间采样，因此对于时间 延迟的分辨率就会受到采样时间
T
s
T_s
Ts的影响。在前面给出了测量结果中可以看出时间延迟曲线呈现明显的台阶，这是由于采样时间所引起的时间分辨率引起的。
根据实验延迟计算出空间距离，同样也会具有一个分辨率下限。
Δ
d
=
T
s
⋅
v
a
i
r
\Delta d = T_s \cdot v_{air}
Δd=Ts⋅vair。其中的
v
a
i
r
v_{air}
vair是空气中的升速，
T
s
T_s
Ts是声音信号采样时间。
下面讨论如何提高声音延迟的分辨率、计算效率、以及麦克的不同空间指向对于测量结果的影响。
快速相关运算
1. 利用FFT加速计算相关运算
（1）相关运算的复杂度
对于两个时间信号
x
(
t
)
,
y
(
t
)
x\left( t \right),y\left( t \right)
x(t),y(t)，它们的相关运算结果
R
x
y
(
t
)
R_{xy} \left( t \right)
Rxy(t)定义如下：
R
x
y
(
t
)
=
∫
−
∞
∞
x
(
τ
)
⋅
y
(
τ
−
t
)
∗
⋅
d
τ
R_{xy} \left( t \right) = \int_{ - \infty }^\infty {x\left( \tau \right) \cdot y\left( {\tau - t} \right)^* \cdot d\tau }
Rxy(t)=∫−∞∞x(τ)⋅y(τ−t)∗⋅dτ如果这两个信号都是实值信号，公式里面的共轭就可以省略。
对于两个实数离散时间信号
x
[
n
]
,
y
[
n
]
,
n
∈
{
0
,
.
.
.
,
N
−
1
}
x\left[ n \right],y\left[ n \right],n \in \left\{ {0,...,N - 1} \right\}
x[n],y[n],n∈{0,...,N−1}，它们之间的普通相关运算定义为：
R
x
y
[
n
]
=
∑
m
,
m
−
n
∈
[
0
,
N
−
1
]
x
[
m
]
⋅
y
[
m
−
n
]
R_{xy} \left[ n \right] = \sum\limits_{m,m - n \in \left[ {0,N - 1} \right]}^{} {x\left[ m \right] \cdot y\left[ {m - n} \right]}
Rxy[n]=m,m−n∈[0,N−1]∑x[m]⋅y[m−n]
从相关运算定义来看，计算两个长度为 N N N的序列相关运算，乘法、加法的计算复杂度与 N 2 N^2 N2成正比。
（2）相关运算与卷积运算的关系
在信号运算中，还有一个应用更广泛的运算：卷积运算。
x
(
t
)
,
y
(
t
)
x\left( t \right),y\left( t \right)
x(t),y(t)之间的卷积运算定义为：
x
(
t
)
∗
y
(
t
)
=
∫
−
∞
∞
x
(
τ
)
⋅
y
(
t
−
τ
)
d
τ
x\left( t \right) * y\left( t \right) = \int_{ - \infty }^\infty {x\left( \tau \right) \cdot y\left( {t - \tau } \right)d\tau }
x(t)∗y(t)=∫−∞∞x(τ)⋅y(t−τ)dτ
对比一下信号的相关运算和卷积运算的定义，可以看出它们之间的关系：
R
x
y
(
t
)
=
x
(
t
)
∗
y
∗
(
−
t
)
R_{xy} \left( t \right) = x\left( t \right) * y^* \left( { - t} \right)
Rxy(t)=x(t)∗y∗(−t)
（3）快速卷积数值计算
之所以讨论相关运算与卷积运算之间的关系，是为了寻找相关运算的快速算法。
计算一个序列 x [ n ] , n ∈ { 0 , . . . , N − 1 } x\left[ n \right],n \in \left\{ {0,...,N - 1} \right\} x[n],n∈{0,...,N−1}的离散傅里叶变换 D F T { x [ n ] } DFT \left\{ {x\left[ n \right]} \right\} DFT{x[n]}有相应的快速算法-快速傅里叶变换 F F T { x [ n ] } FFT\left\{ {x\left[ n \right]} \right\} FFT{x[n]}，在 N N N为2的整数次幂的情况下，计算FFT的乘法，加法的复杂度都在 log 2 N \log _2 N log2N的数量级别。正变换和反变换的复杂度相同。
在根据傅里叶变换的卷积定理，序列的时域卷积（和）运算，在频域是乘积运算。基于此，再利用前面讨论的相关与卷积运算之间的关系，可以得到计算两个序列的相关运算的快速算法：
R x , y ( t ) = F F T − 1 { F F T [ x ( t ) ] ⋅ F F T [ y ( t ) ] ∗ } R_{x,y} \left( t \right) = FFT^{ - 1} \left\{ {FFT\left[ {x\left( t \right)} \right] \cdot FFT\left[ {y\left( t \right)} \right]^* } \right\} Rx,y(t)=FFT−1{FFT[x(t)]⋅FFT[y(t)]∗}
在实际工程中，往往参与卷积的两个信号实现已知，比如在利用声音定位的时候，发送声音信号往往是事先确定好的固定的Chirp信号，每次参与计算的新的信号是接收到的回声信号。所以在上面卷积的快速算法中，对于已知信号的FFT可以事先计算好并存储，实际运算中只需要完成对新采集到信号的FFT计算，以及对乘积结果的反FFT计算。
最后需要补充一下，在利用上面公式计算的时候，还需要将两个信号通过补0，变成长度等于两个序列长度之和减一。
2. 快速相关运算结果
使用快速计算，比直接在时域中进行卷积速度大大提高了。通过计算获得100组采样数据，每组2048个数据的相关峰值点，使用FFT需要大约：0.36秒钟；而使用普通的相关运算则需要200秒左右。
下图显示了通过两种方法所得到的相关峰值位置随着麦克风距离音箱位置变化而产生的延迟，这两种方法所得到的结果是一样的。
通过FFT得到的相关结果峰值位置
快速算法和直接计算两种方法的速度对比，可以通过下面两个动图所显示的计算过程能够体会出来。
下面是应用FFT快速计算机相关运算的过程。
使用FFT计算100组，每组2048点相关结果只需要0.37秒
下面则是直接计算卷积的过程：
直接计算100组，每组2048数据点相关结果则需要192秒左右
提高相关运算的空间精度
1. 为什么前面测量结果曲线中出现台阶
如果直接根据序列的相关结果峰值位置确定声音延迟，那么声音延迟的时间分辨率就是声音信号的采样时间 T s T_s Ts，再根据声音速度 v s o u n d v_{sound} vsound，可以计算出所对应的测量距离的空间分辨率 Δ d = T s × v s o u n d \Delta d = T_s \times v_{sound} Δd=Ts×vsound。
在采样时间 T s = 0.1 m s T_s = 0.1ms Ts=0.1ms，20℃空气速度 v 2 0 o c = 343 ( m / s ) v_{20^o c} = 343\,\,\left( {m/s} \right) v20oc=343(m/s)，对应的空间分辨率 Δ d = 3.4 ( c m ) \Delta d = 3.4\left( {cm} \right) Δd=3.4(cm)，也就是当收音麦克与音源之间的距离变化小于3.4厘米时，所测量得到的结果是一样的。这也就解释了前面100个位置点测量声音延迟曲线出现了很多台阶的原因。
2. 如何提高测量结果的空间分辨率
提高基于声音采样数据相关方法测量距离的空间分辨率，可以通过提高AD采样速率来解决。但这需要更高速的AD转换器，更多的数据存储内存以及更快速数据计算能力。
除此之外，还可以通过数据插值处理的方法来提高测量结果的空间分辨率。
数据插值可从离散时间采样数据 x [ n T s ] x\left[ {nT_s } \right] x[nTs]中获得时间更加密集的数据 x [ n T s 1 ] x\left[ {nT_{s1} } \right] x[nTs1]。插值分解成两个过程：第一个过程是将离散时间信号恢复成一个连续时间信号 x ( t ) x\left( t \right) x(t)；第二个过程就是在连续时间信号的基础上采用更加密集的是时间间隔采样 x [ n T s 1 ] x\left[ {nT_{s1} } \right] x[nTs1]。
恢复成连续时间信号可以有零阶保持、一阶保持、理想插值等不同方法，它们都可以看成是离散时间采样脉冲信号与一个插值函数进行卷积的结果： x ( t ) r e c o n s t r u c t e d = h ( t ) ∗ x s ( t ) x\left( t \right)_{reconstructed} = h\left( t \right) * x_s \left( t \right) x(t)reconstructed=h(t)∗xs(t)
零阶保持、一阶保持、理想插值分别对应的卷积信号是矩形信号、三角信号以及 sin c ( t ) \sin c\left( t \right) sinc(t)信号等。
不同的离散时间信号重建成连续时间信号的方法
使用理想插值所获得的结果更加平滑，但计算起来相对比较复杂。但如果是从离散时间信号的傅里叶变换结果中恢复插值信号的话，则有一个非常方便的方法，那就是通过对数据的DFT结果补零，获得更长的频谱数据，再通过反离散傅里叶变换，就可以得到原来数据的理想插值结果了。具体的 原理在信号与系统课程中会进行介绍的。
由于前面在快速计算相关结果的时候，就利用了快速傅里叶变换，所以可以在最后一步进行反傅里叶变换的时候，先进行补零，然后在进行。
3. 实验数据对比
下面给出了插值10倍之后所获得的相关峰值位置结果，对比原始计算方法，可以看到经过插值之后的结果明显平滑多了。通过插值后的结果所获得的空间分辨率就从原来的3.6厘米降低到3.6毫米了。
经过空间插值细化后10倍后的相关峰值位置计算结构
由于实验环境是在室内，存在着很多反射波的干扰，所以当距离远了之后，距离测量出现了很多的波动，它们反映了空间中的很多驻波干扰。下图给出了不同插值倍率下结果曲线。
随着插值倍数增加，所得到延迟曲线变化
通过简单的差值就可以轻松提高测距的空间分辨率，所需要的代价就是计算时间加长了。下面图给出了插值的倍数与结果计算消耗的时间之间的关系，整体上呈现线性比例关系。
细化倍数和计算时间
室内环境反射波对于测量结果的影响
由于是在室内进行实验，麦克风所接受到的声音信号除了直接来自于声源之外，可能还包括有四周墙壁的反射信号。如果声源距离比较近，四周的反射声波强度受到衰减，对于测量结果影响较小。
下面通过控制接受麦克的不同方向，考察一下测量结果是否受到影响。
使用舵机控制MIC的方向
将麦克放置距离音箱30厘米左右，方向从左到右旋转180°。对所获得的声音数据使用前面给出的插值细分方法，得到的声音延迟时间。下图给出了不同方向声音延迟结果曲线：
不同指向对应的相关延迟结果
延迟时间对应的最大值和最小值分别是：
t
max
=
4895.5
,
t
min
=
4896.4
t_{\max } = 4895.5,\,\,\,\,t_{\min } = 4896.4
tmax=4895.5,tmin=4896.4
在10kHz的采样率下，上述时间差所对应的距离变化为：
Δ
d
=
(
t
max
−
t
min
)
⋅
v
2
0
0
c
=
0.9
×
1
0
−
4
×
343
=
3
(
c
m
)
\Delta d = \left( {t_{\max } - t_{\min } } \right) \cdot v_{20^0 c} = 0.9 \times 10^{ - 4} \times 343 = 3\,\,\left( {cm} \right)
Δd=(tmax−tmin)⋅v200c=0.9×10−4×343=3(cm)
这个3厘米的变化距离和麦克风转动过程中所引起的距离变化大体相当，说明此事四周的反射波对于测量结果影响不大。同时所使用的MIC的方向指向性并不强，可以对来自于180°方向的声波都能够很好的探测。
下面是分别将麦克在距离音箱10厘米和50厘米处重新测量不同的指向对于测距结果的影响。
在10厘米距离下麦克不同方向对应的声音延迟
结论
根据前面的一些初步实验数据，对于普通的喇叭发出的Chirp声音，通过小型全向麦克所获得的声音信号可以获得空间分辨率在0.5厘米之内的检测精度。当然，这个精度会受到麦克风的与声源之间的指向的影响，上述的测量误差就会引入3个厘米的误差。
3个厘米误差对于在赛场上比赛的智能车模来说，问题不是很大。可以基于此设计出未来新型的室外声音引导的信标赛题组。
AD\DA采样板的设计说明：https://blog.csdn.net/zhuoqingjoking97298/article/details/104180878 ↩︎
扫一扫
6955
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
