紧接上一个博客文章，此为第三部分。上一部分见：麦克风阵列技术 二 （自动增益控制 自动噪声抑制 回声消除 语音活动检测）

麦克风阵列技术详解

声源定位

麦克风阵列可以自动检测声源位置，跟踪说话人，声源定位信息既可以用于智能交互，也可以用于后续的空域滤波，对目标方向进行语音增强。

利用麦克风阵列可以实现声源到达方向估计（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/pixc12 = zeros(2*N-1,1);
m = 0;
for i = -(N-1):N-1m = m+1;for t = 1:Nif 0<(i+t)&&(i+t)<=Nxc12(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(ω)|A(ω)=1|P(ω)|

基本思想就是求时间差只需要相位信息，舍弃不相关的幅度信息以提高健壮性，可以看到当A(ω)=1A(ω)=1的情况下就是经典互相关
P(ω)P(ω)为复数，可以表示为|P(ω)| * e - jωp|P(ω)| * e - jωp,去掉幅度信息后，就只剩相位信息e - jωpe - jωp了，要得到相位信息，可以用P(ω)abs(P(ω))P(ω)abs(P(ω))计算，也可以直接用matlab中的angle函数计算，即angle(P(ω))angle(P(ω))，

几行代码验证下：

x1 = [1,2,3,7,9,8,3,7]';
x2 = [4,5,6,5,4,3,8,2]';[tau,R,lag] = gccphat(x1,x2) N = length(x1)+length(x2)-1;
NFFT = 32;
P = (fft(x1,NFFT).*conj(fft(x2,NFFT)));
A = 1./abs(P);
R_est1 = fftshift(ifft(A.*P));
range = NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2;
R_est1 = R_est1(range);R_est2 = fftshift(ifft(exp(1i*angle(P))));
R_est2 = R_est2(range);

可以看到，三种不同写法得到的R_est1 、R_est2 与matlab自带函数gccphat计算得到的R相等。

那上面例子中的宽带语音信号，用GCC-PHAT方法得到具有陡峭峰值互相关函数，找到互相关最大时的点，结合采样频率Fs与与麦克风间距dFs与与麦克风间距d，就可以得到方向信息。

角度计算

上面的内容计算了两个麦克风的延时，实际中假设阵列中麦克风个数为NN，则所有麦克风间两两组合共有N(N-1)/2N(N-1)/2对，记第kk个麦克风坐标为(xk,yk,zk)(xk,yk,zk),声源单位平面波传播向量u? =(u,v,w)u→=(u,v,w),如果麦克风k,jk,j之间的延时为τkjτkj，则根据向量关系有下式，其中c为声速，
C * τkj = -(xk→-xj→)?u? c?Τkj = -(xk→-xj→) * u→

这样看起来不够直观，那就代入坐标写成标量形式如下：
Cτkj=u(xk-xj)+v*(yk-yj)+w*(zk-zj)cτkj=u(xk-xj)+v*(yk-yj)+w*(zk-zj)

当有多个麦克风时，每两个麦克风就可以得到一组上式，N个麦克风就会有N*(N-1)/2个等式N个麦克风就会有N*(N-1)/2个等式，声源单位传播向量u? =(u,v,w)u→=(u,v,w)有三个未知数，因此最少只需要三组等式，也就是三个麦克风就可以计算出声源方向，这里就先假定N=3N=3，可以得到方程组如下：

写成矩阵形式

求出u? =(u,v,w)u→=(u,v,w)后，由正余弦关系就有了角度值了
θ=acos(1w)θ=acos(1w)
α=acos(usin(acos(1w)))α=acos(usin(acos(1w)))

当麦克风数量N>3N>3时，其实所有组合信息对于角度值的计算是有冗余的，这个时候可以求出所有组合的角度值，然后利用最小二乘求出最优解，这样可以利用到所有的麦克风的信息来提高角度估计的稳定性。

波束形成

DBF是Digital Beam Forming的缩写，译为数字波束形成 或数字波束合成。数字波束形成技术是天线波束形成原理与数字信号处理技术相结合的产物，其广泛应用于阵列信号处理领域。
专业：电子、通信与自动控制技术阵列信号处理最主要的研究内容包括DOA估计和波束形成。较早的DOA估计方法又称为波束形成方法，而该波束形成方法利用了空域维纳滤波的匹配概念，由阵列流形在信号空间中的投影大小判定信号方向，后来随着研究的深入，高分辨谱估计技术的发展，才把DOA估计和波束形成加以区分，DOA估计是为了确定信号的方位，从接收数据中测出信号方向，不管信号是有用信号还是干扰信号，在DOA估计方向图中都表现为峰值，而此峰值并不是实际阵列输出功率；波束形成是传统滤波的空域拓展，其根本目的是有效提取有用信号并抑制噪声和干扰，在方向图中表现为有用信号方向形成峰值、干扰方向形成零陷，可以认为DOA估计为波束形成的前端处理，确定期望信号和干扰方向后，阵列对期望信号方向形成波束并在干扰方向形成零陷。

DBF体现的是声源信号的空域选择性，许多传统波束形成方法具有线性处理结构；波束形成需要考虑三个方面：
1)麦克风阵列个数；
2)性能；
3)鲁棒性。

在麦克风较少时，波束形成的空域选择性差，当麦克风数量较多时，其波束3dB带宽较为窄，如果估计的目标声源方向有稍有偏差，带来的影响也更大，鲁棒性不好。通常鲁棒性和性能是对矛盾体，需要均衡来看。

WebRTC使用了如下几个点以提高鲁棒性和性能（其算法性能优先）：
1.可以使用多个后置滤波器而非一个，2.每个后置滤波使用新的结构。

每个后置滤波器为每个声学场景的时频域bin在均方误意义上提供了最优实增益。在WebRTC中后置滤波器根据声源的空域协方差矩阵，干扰源协方差矩阵，绕射场（零阶贝塞尔函数计算）以及临近麦克风的时频信号信息求得。

这样的话就可以为每个声源和干扰场景计算出不同的最优后置滤波器，也可以使用级联的方式灵活使用多个不同选择性的后置滤波器。
当前现存的波束形成算法的鲁棒性成为它们使用的一道门槛，如MVDR和多通道维纳滤波。

WebRTC为了增强鲁棒性，在求最优矩阵时，对声源信号添加了限制条件，使用Gabor frame将声源变成时频bin的系数，对这些bin按照目标声源和干扰声源附加了条件，如果满足条件，则门操作让目标声源通过，而让干扰源乘以零以实现选择最优目标信号。

在WebRTC中这些增益系数称为自适应标量（上面的实）乘法增益，均方误差准则被用来做为计算的准则。由于阵列方向响应随频率是变换的，而语音信号又是宽带信号，所以WebRTC中使用了gabor变换来表示声音信号。增益源于目标信号和干扰的比例。

波束形成模型

以均匀线阵为例：

按窄带模型分析：

可以写成矩阵形式：

其中a(θ)为方向矢量或导向矢量（Steering Vector），波束形成主要是针对各个接收信号X进行权重相加。

波束形成基本理论

1. A-波束形成
   权重相加：
   不同的波束形成，就是不同的权重W。

2. B-瑞利限
   以均匀直角窗为例：
   
   得出方向图：
   
   主瓣宽度正比于孔径宽度的倒数：
   
   因为孔径的限制，造成波束宽度存在限制（不会无限制小），近而落在主瓣波束内部的两个信号便会混在一起而分不清，这就存在瑞利限的问题。

   直角窗主瓣宽度为：
   
   其中λ为入射波长，theta1为入射角，Md为阵列孔径。

3. C-常见窗函数

   对于空间不同的阵列信号，类似采样分析（空域采样），自然可以加窗进行处理，不加窗可以认为是直角窗，另外也可以选择汉明窗、hanning窗等等。

   加窗可以改变波束宽度以及主瓣、副瓣等特性，可以借助MATLAB 的wvtool观察不同窗函数特性。

N = 192;
w = window(@blackmanharris,N);
wvtool(w)

4. D-DFT实现
   阵列的采样间隔是相位信息：
   
   这就类似于频域变换，只不过这里的相位信息：对应的不是频率，而是不同位置，可以看作空域的变换。

   分别对阵列信号进行直接加权、加窗、DFT实现：

clc;clear all;close all
M = 32;
DOA = [-30 30];
SNR = 10;
theta = -90:.1:90;
len = length(theta);
SignalMode = 'Independent';
fc = 1e9;
c = 3e8;
lambda = c/fc;
d = lambda/2;
N = 100;%snap points
x = StatSigGenerate(M,N,DOA,SNR,SignalMode,lambda,d);
R_hat = 1/N*x*x';
output = zeros(3,len);
for i = 1:lena = exp(1j*2*pi*[0:M-1]'*sind(theta(i))*d/lambda);W = (inv(R_hat)*a)*(1./(a'*inv(R_hat)*a));output(1,i) = mean(abs(W'*x),2);output(2,i) = 1./(a'*inv(R_hat)*a);output(3,i) = a'*x*ones(N,1);
end
output = abs(output);
output = output - repmat(mean(output.')',1,size(output,2));
output = output./repmat(max(output.')',1,size(output,2));
%plot
plot(theta,output(1,:),'k',theta,output(2,:),'r--',theta,output(3,:),'b');
legend('MVDR 波束','MVDR 谱','固定权重 波束');

5. E-自适应波束形成
   直接相加也好、加窗也好，都是固定的权重系数，没有考虑到信号本身的特性，所以如果结合信号本身去考虑就形成了一系列算法：自适应波束形成。

   这类步骤通常是：
   1）给定准则函数；
   2）对准则函数进行求解。

   准则常用的有：信噪比（snr）最大准则、均方误差最小准则（MSE）、线性约束最小方差准则（LCMV）、最大似然准则（ML）等等；

   求解的思路大体分两类：1）直接求解，例如MVDR中的求解；2）也可以利用梯度下降的思想，如随机梯度下降、批量梯度下降、Newton-raphson等方法，不再详细说明。

   以MVDR举例：
   
   这里采用直接求解的思路：
   
   将求解的W带入
   
   即可得到波束形成。

6. F-栅瓣现象
   栅瓣是一类现象，对应干涉仪就是相位模糊（相位超过2*pi），对应到Beamforming就是栅瓣问题，具体不再论述，给出现象（同样的波束，在不同的位置分别出现）：
   

7. G-波束形成与空间谱

   之前分析过MVDR的方法，得到的输出（含有约束的最小均方误差准则）为：
   
   有时候也称这个输出为空间谱，其实就是|y2(t)|，但这个与MUSIC等算法的谱还不是一回事，只是有时候也被称作空间谱，所以这里多啰嗦几句，分析这个说法的来源。

   已知N个采样点的信号，对其进行傅里叶变换：
   
   进一步得到功率谱密度：
   
   根据上文的分析：y(t)其实对应的就是空域变换（可借助DFT实现），类比于时频处理中的频域变换。而这里又可以看到频域变换的平方/长度，对应就是功率谱，这是频域的分析。

   对应到空域，自然就是|y2(t)|/长度，对应空间谱，长度只影响比例关系，所以MVDR的最小方差输出被称作：空间谱也是合适的。

   给出一个测试（这里如果），对比MVDR的y(t)、MVDR功率谱以及普通Beamforming的结果：

clc;clear all;close all
M = 32;
DOA = [-30 30];
SNR = 10;
theta = -90:.1:90;
len = length(theta);
SignalMode = 'Independent';
fc = 1e9;
c = 3e8;
lambda = c/fc;
d = lambda/2;
N = 100;%snap points
x = StatSigGenerate(M,N,DOA,SNR,SignalMode,lambda,d);
R_hat = 1/N*x*x';
output = zeros(3,len);
for i = 1:lena = exp(1j*2*pi*[0:M-1]'*sind(theta(i))*d/lambda);W = (inv(R_hat)*a)*(1./(a'*inv(R_hat)*a));output(1,i) = mean(abs(W'*x),2);output(2,i) = 1./(a'*inv(R_hat)*a);output(3,i) = a'*x*ones(N,1);
end
output = abs(output);
output = output - repmat(mean(output.')',1,size(output,2));
output = output./repmat(max(output.')',1,size(output,2));
%plot
plot(theta,output(1,:),'k',theta,output(2,:),'r--',theta,output(3,:),'b');
legend('MVDR 波束','MVDR 谱','固定权重 波束');

如果将d = lambda/2;改为d = lambda/0.5;，自然就有了栅瓣：

去混响

混响是指声波在室内传播时，要被墙壁、天花板、地板等障碍物反射，当声源停止发声后，声波在室内要经过多次反射和吸收，最后才消失。这种现象称为混响。因此，当声源和麦克风之间的距离越远，反射声占的比例就越高，混响就严重。
经典的去混响方法包括形成拾音波束来减少反射声和基于反卷积的去混响方法。

麦克风阵列结构设计

MUC 孔的孔深孔径比越小越好，即开孔越大越好，深度越小越好，尽量向1:1靠近。孔深与孔径比值越大，麦克频响的震点越像低频靠近，要求震点在12KHz以上。最少也要在8KHz以上。喇叭腔体不能漏气。这是因为，喇叭正反两面的声波相位相差180度，当音腔有漏气时，声波会发生抵消，尤其是低频频段。
麦克和喇叭的失真都要小。麦克失真小于4%，喇叭失真小于10%，由于喇叭低频失真严重些，会超过10%，可以考虑增加滤波器滤掉低频成分。
喇叭腔体四周与其他机构件保留1mm的距离，防止腔壳与机构接触产生异音。
喇叭鼓膜上方与机构件保留1.5mm的距离，以防鼓膜振动碰到机构件产生异音。
喇叭与机构件有接触的地方，要增加泡面，以起到缓冲、减振的效果，防止喇叭振动时与机构件碰撞产生异音。

声学结构确认流程

1)远程会议或现场结构设计评估
确认麦克阵列构型，确认声腔及安装结构设计，确认进声孔深度、直径大小等；
2)声学实验室录音效果评估-第一阶段
计算裸麦和带声腔结构的麦克风之间的录音之间谐波程度，根据分析结果确定是否通过。
3)声学实验室录音效果评估-第二阶段
分别利用裸麦和带声腔结构的麦克风信号做基于相位的声源定位，如果两者定位误差小于5°，则认为通过该项测试。
4)声学实验室录音效果评估-第三阶段
分别利用裸麦和带声腔结构的麦克风录音数据进行识别，效果差距在2%以内，则认为远场识别方面无问题。