紧接上一个博客文章,此为第三部分。上一部分见:麦克风阵列技术 二 (自动增益控制 自动噪声抑制 回声消除 语音活动检测)
麦克风阵列技术详解
声源定位
麦克风阵列可以自动检测声源位置,跟踪说话人,声源定位信息既可以用于智能交互,也可以用于后续的空域滤波,对目标方向进行语音增强。
利用麦克风阵列可以实现声源到达方向估计(direction-of-arrival (DOA) estimation),DOA估计的其中一种方法是计算到达不同阵元间的时间差,另外一种可以看这里,这篇主要介绍经典的GCC-PHAT方法
简单说明问题背景,信号模型如下图,远场平面波,二元阵列
要计算得到θθ,其实就是要求两个阵元接收到的信号时间差,现在问题变成到达时间差估计(Time-Difference-of-Arrival Estimation),因此,基于延时估计的DOA方法,其实也可以看做是分两步进行的,第一步是估计延时,第二步是计算角度,与之相对应的基于空间谱估计的DOA方法就是一步完成的。下面就分两步进行介绍
延时估计
一、 互相关函数(cross-correlation function)
计算y1(k)y1(k)与y2(k)y2(k)的时间差,可以计算两个信号的互相关函数,找到使互相关函数最大的值即是这两个信号的时间差
离散信号的互相关函数
R(τ)=E[x1(m)x2(m+τ)]R(τ)=E[x1(m)x2(m+τ)]
求时间差就是找到互相关函数最大时的点
D=argmaxR(n)D=argmaxR(n)
MATLAB代码如下:
%%
% Load the chirp signal.
load chirp;
c = 340.0;
Fs = 44100;
%%
d = 0.25;
N = 2;
mic = phased.OmnidirectionalMicrophoneElement;
% array = phased.URA([N,N],[0.0724,0.0418],'Element',mic);
array = phased.ULA(N,d,'Element',mic);
%%
% Simulate the incoming signal using the |WidebandCollector| System
% object(TM).
arrivalAng = 42;
collector = phased.WidebandCollector('Sensor',array,'PropagationSpeed',c,...
'SampleRate',Fs,'ModulatedInput',false);
signal = collector(y,arrivalAng);
x1 = signal(:,1);
x2 = signal(:,2);
N =length(x2);
xc = xcorr(x1,x2,'biased');
[k,ind] = max(xc);
an = acos((ind-N)/Fs*340/d)*180/pi
xc12 = zeros(2*N-1,1);
m = 0;
for i = -(N-1):N-1
m = m+1;
for t = 1:N
if 0<(i+t)&&(i+t)<=N
xc12(m) = xc12(m) + x2(t)*x1(t+i);
end
end
end
xc12 = xc12/N;
以上程序中的循环就是上面的定义公式,运行程序可以看到循环部分计算的互相关与直接调用matlab的xcorr结果相同(注意matlab中互相关默认没做归一化),找到互相关函数的最大值就可以得到时间差 。
二、广义互相关(generalized cross-correlation)
理论上使用上面个介绍的CCF方法就可以得到时间差,但是实际的信号会有噪声,低信噪比会导致互相关函数的峰值不够明显,这会在找极值的时候造成误差。
为了得到具有更陡峭极值的互相关函数,一般在频域使用一个加权函数来白化输入信号,这就是经典的广义互相关方法。
由维纳-辛钦定理可知,随机信号的自相关函数和功率谱密度函数服从一对傅里叶变换的关系,即x1、x2x1、x2的互功率谱可由下式计算
P(ω)=∫+∞?∞R(τ)e?jωτdτP(ω)=∫?∞+∞R(τ)e?jωτdτ
R(τ)=∫+∞?∞P(ω)ejωτdωR(τ)=∫?∞+∞P(ω)ejωτdω
这一步是把互相关函数变换到了频域,哦对,上面说到是想白化互相关函数,那就把上面第二式添加一个系数
R(τ)=∫+∞?∞A(ω)P(ω)ejωτdωR(τ)=∫?∞+∞A(ω)P(ω)ejωτdω
设计不同的频域系数A(ω)A(ω)对应着不同方法,这里只介绍 PHAT(phase transform)方法,即取系数如下:
A(ω)=1|P