部分beamforming知识汇总

文章目录

• 1. 基础知识以及概念
• 1.1. 通用模型假设：
• 1.2. 频率-波数响应
• 1.3. 均匀线性阵列（Uniform Linear Array）
• 1.4. 主波束调向(beam steering)
• 1.5. 方向图合成
• 1.6. 零点调向（null steering/ sidelobe canceller）
• 小总结
• 2. 常见的Beamformer（DSB，MVDR，GSC，GEV）
• 2.1. Delay-Sum Beamforming
• 2.2. Minimum Variance Distortionless Response (MVDR)
• 2.3. Generalized Sidelobe Canceller (GSC)
• 2.4. Generalized Eigen Value(GEV)
• 3 声源定位算法
• 4.

这段时间做了一些和增强相关的任务，虽然我用的方法是NN-based方法，但是这套传统理论我看了之后还是很感兴趣。学了一些深蓝学院的课程之后，结合之前做的，大概梳理了一些东西出来。

1. 基础知识以及概念

1.1. 通用模型假设：

这里可以假设声音属于平面波，方向为$\textbf a$a;假设有M个麦克风组成麦克风阵列$[\textbf p_0,...,\textbf p_{M-1}]$[p0​,...,pM−1​]
则接收到的信号为:
$\textbf x(n,\textbf p) = \left[ \begin{array}{c} x(n,\textbf p_0) \\ x(n,\textbf p_1) \\ ... \\ x(n, \textbf p_{M-1}) \end{array} \right]$x(n,p)=⎣⎢⎢⎡​x(n,p0​)x(n,p1​)...x(n,pM−1​)​⎦⎥⎥⎤​
由于入射方向$\textbf a$a是已知的，就可以通过入射方向估计出信号到达每个麦克风的时延$\tau$τ
$\tau = \frac{\textbf a^T\textbf p_i}{c}, i = 0,...,N-1$τ=caTpi​​,i=0,...,N−1
意为每个麦克风的坐标在入射方向上的投影距离，处以声速$c$c，如果规定$\theta$θ是入射角与z轴的夹角，$\phi$ϕ为向量在x0y的投影与x的夹角，那么$\textbf a$a可以表示为
$\textbf a = \left[\begin{array}{c} -sin\theta cos \phi \\ -sin\theta sin \phi \\ -cos\theta \end{array} \right]$a=⎣⎡​−sinθcosϕ−sinθsinϕ−cosθ​⎦⎤​

由于实际接收到的信号有延时差异。那么，接收到的信号在时域上就可以表示为
$\textbf x(n,\textbf p) = \left[ \begin{array}{c} x(n-\tau_0) \\ x(n-\tau_1) \\ ... \\ x(n-\tau_{M-1}) \end{array} \right]$x(n,p)=⎣⎢⎢⎡​x(n−τ0​)x(n−τ1​)...x(n−τM−1​)​⎦⎥⎥⎤​
时域上的延时可以直接在频域上的相位中体现出来
$\textbf X(\omega) = \left[ \begin{array}{c} e^{-j\omega\tau_0}\\ e^{-j\omega\tau_1}\\ ...\\ e^{-j\omega\tau_{M-1}} \end{array} \right]X(\omega)$X(ω)=⎣⎢⎢⎡​e−jωτ0​e−jωτ1​...e−jωτM−1​​⎦⎥⎥⎤​X(ω)
$\omega\tau_i = \omega \frac{\textbf a^T \textbf p_i}{c}= (\frac{w}{c}\textbf a)^T\textbf p_i = \textbf k^T \textbf p_i = \frac {2\pi}{\lambda}\textbf a^T\textbf p_i$ωτi​=ωcaTpi​​=(cw​a)Tpi​=kTpi​=λ2π​aTpi​
其中，$\lambda$λ是波长，$\textbf k$k就被称为波数，是后面会经常用到的一个概念，因为它同时包含了声源方向和频率两个信息。则，上面的频域接收信号表达式可以表示为
$\textbf X(\omega) = \left[ \begin{array}{c} e^{-j\textbf k^T\textbf p_0}\\ e^{-j\textbf k^T\textbf p_1}\\ ...\\ e^{-j\textbf k^T\textbf p_{M-1}} \end{array} \right]X(\omega) = \textbf V_k(\textbf k)X(\omega)$X(ω)=⎣⎢⎢⎡​e−jkTp0​e−jkTp1​...e−jkTpM−1​​⎦⎥⎥⎤​X(ω)=Vk​(k)X(ω)
其中$\textbf V_k(\textbf k)$Vk​(k)被称为阵列流形矢量（array mainifold vector）

1.2. 频率-波数响应

一个频率响应函数表示的是一个滤波器的频率相应，但是在阵列系统中，最关键的就是引入了空间的概念，也就是有频率-空间，两个维度。那么在求响应函数的时候也需要把空间的信息考虑进来。
$\textbf Y(\omega) = \textbf H^T(\omega)\textbf X(\omega) = \textbf H^T(\omega)\textbf V_k(\textbf k)X(\omega)$Y(ω)=HT(ω)X(ω)=HT(ω)Vk​(k)X(ω)
其中$\textbf H(\omega)$H(ω)就是频率响应函数。那么就可以规定频率-波数响应（frequency-wavnumber response）函数为
$\gamma(\omega,\textbf k) = \textbf H^T(\omega)\textbf V_k(\textbf k)$γ(ω,k)=HT(ω)Vk​(k)
那我们经常听到另一个概念是beam pattern，意味波束方向图，实际上是频率-波数响应在球面$k = 2\pi / \lambda$k=2π/λ上进行取值
$B(\omega:\theta, \phi) = \gamma(\omega,\textbf k)|_{\textbf k = \frac{2\pi}{\lambda}\textbf a}$B(ω:θ,ϕ)=γ(ω,k)∣k=λ2π​a​

1.3. 均匀线性阵列（Uniform Linear Array）

ULA适用于麦克风阵列是一排的这种情况，这种就被称为线性阵列。假设所有的麦克风阵列均匀分布在z轴上，相邻两麦克风的间距为d，也就是每个麦克风的坐标就是
$p_{z_i} =(i - \frac {M-1}2)d, i=0,1,…,M-1 $
$p_{x_i} = p_{y_i} = 0$pxi​​=pyi​​=0
把坐标带入阵列流形矢量可以得到
$\textbf V_k(\textbf k) =[e^{-j(0-\frac {M-1}2)k_zd},e^{-j(1-\frac {M-1}2)k_zd},...,e^{-j(M-\frac {M-1}2)k_zd}]^T$Vk​(k)=[e−j(0−2M−1​)kz​d,e−j(1−2M−1​)kz​d,...,e−j(M−2M−1​)kz​d]T
$k_z = -\frac {2\pi}\lambda cos\theta$kz​=−λ2π​cosθ
则它的波束方向图是
$B_\theta(\theta) = \textbf H^T\textbf V_\theta(\theta) = e^{-j(\frac {M-1}2)\frac{2\pi d}{\lambda}cos\theta}\sum_{i=0}^{N-1}H_i e^{ji\frac{2\pi d}\lambda cos\theta}, 0\le\theta\le\pi$Bθ​(θ)=HTVθ​(θ)=e−j(2M−1​)λ2πd​cosθi=0∑N−1​Hi​ejiλ2πd​cosθ,0≤θ≤π
若$N=11,H_i=\frac 1 N,d = \frac \lambda 2$N=11,Hi​=N1​,d=2λ​

可以看出，频率-波数响应函数和权重向量之间是傅立叶变化的关系
$B_\psi(\psi) = \textbf H^T\textbf V_\psi(\psi) = e^{-j(\frac {M-1}2)\psi}\sum_{i=0}^{N-1}H_i e^{ji\psi}= e^{-j(\frac {M-1}2)\psi}(\sum_{i=0}^{N-1}H_i e^{-ji\psi})^*$Bψ​(ψ)=HTVψ​(ψ)=e−j(2M−1​)ψi=0∑N−1​Hi​ejiψ=e−j(2M−1​)ψ(i=0∑N−1​Hi​e−jiψ)∗

1.4. 主波束调向(beam steering)

如果我们希望$\textbf k = \textbf k_T$k=kT​,那么频率-波数响应函数为$\gamma(\omega,\textbf k -\textbf k_T)$γ(ω,k−kT​)，阵列流形矢量
$\left[ \begin{array}{c} e^{-j(\textbf k-\textbf k_T)^T\textbf p_0}\\ e^{-j(\textbf k-\textbf k_T)^T\textbf p_1}\\ ...\\ e^{-j(\textbf k-\textbf k_T)^T\textbf p_{M-1}} \end{array} \right]$⎣⎢⎢⎡​e−j(k−kT​)Tp0​e−j(k−kT​)Tp1​...e−j(k−kT​)TpM−1​​⎦⎥⎥⎤​
实际上也就可以视为，对权重信息$H$H(或者可以称为$w$w)，进行修改
$\textbf w_T = \left[ \begin{array}{c} e^{-j\textbf k^T\textbf p_0}\\ e^{-j\textbf k^T\textbf p_1}\\ ...\\ e^{-j\textbf k^T\textbf p_{M-1}} \end{array} \right]\bigodot \textbf w$wT​=⎣⎢⎢⎡​e−jkTp0​e−jkTp1​...e−jkTpM−1​​⎦⎥⎥⎤​⨀w

1.5. 方向图合成

主要是可以通过调整权重向量，使得旁瓣尽可能小
比如汉明窗

实际上我们所想要得到的一种情况，就是目标方向的响应很大，同时旁瓣很小，再有就是主瓣的宽度越小越好。这都是通过调整权重矩阵w来实现的。包括beamforming的一些方法。

1.6. 零点调向（null steering/ sidelobe canceller）

假设知道了监听的方向，同时知道了干扰噪声的方向，那么我就会希望干扰噪声出现一个零点，也就是对干扰噪声方向不进行任何响应。着就用到的是零点调向技术。数学符号可以表达为$\textbf w^T\textbf V_k(\textbf k_T) =1$wTVk​(kT​)=1同时满足$\textbf w^T\textbf V_k(\textbf k_i) =0, i=1,2,...,M_0$wTVk​(ki​)=0,i=1,2,...,M0​也就是有$M_0$M0​个方向我们希望没有响应。 一般来说这种条件数大于未知变量数的问题，我们只有一个最小二乘解，也就是对所有的条件都不一定满足，但是各个条件的误差平方和是最小的。然而这里我们明确希望的一点，一定要满足零点约束的话，问题就会随之改变。

于是这个问题可以近似为，已知一个期望权值向量$w_d$wd​，这个向量一般来说我们可以就假设$w_T$wT​为期望的权值向量，它满足了无畸变约束。现在要求一个向量$w$w，使得两者尽可能相近，同时满足零点约束的限制。于是公式可以写成
$\left \{ \begin{aligned} &min ||\textbf w_d - \textbf w||^2 \\ &s.t. \textbf w^H\textbf C = \textbf 0 \end{aligned}\right .${​min∣∣wd​−w∣∣2s.t.wHC=0​
其中$\textbf C$C是零点向量长成的空间。也就是零点列向量组成的矩阵。

这种问题的解就很经典，使用拉格朗日乘子法即可，这个方法后面也会经常用到。最后得到的结果是
$\textbf w^H = \textbf w_d^H(\textbf I_M - \textbf C(\textbf C^H \textbf C)^{-1}\textbf C^H)$wH=wdH​(IM​−C(CHC)−1CH)
实际上可以很明显地看出来，最终的结果是$\textbf w_d$wd​减去一个东西。I后面的一坨是一个投影矩阵，所以最终的结果就是$\textbf w_d$wd​减去了$\textbf w_d$wd​在零点约束空间$\textbf C$C上的投影。

投影矩阵：
满足“幂等矩阵”同时又是“共轭对称阵”，那么它就是投影矩阵
$\textbf C(\textbf C^H\textbf C)^{-1}\textbf C^H$C(CHC)−1CH和$\textbf I_M - \textbf C(\textbf C^H\textbf C)^{-1}\textbf C^H$IM​−C(CHC)−1CH都是投影矩阵，且两者相互正交

小总结

阵列麦克风信号处理实际上最关键的就是引入了空间这个概念，如何体现空间，就是通过信号到达麦克风的时间差体现的，时域的时间差可以直接等价地表现在频域的相移上。然后我们可以通过频域的窗函数去控制它的beam pattern，最终的目的是表现出一种期望方向的信号很大，其他方向或者干扰方向信号为0或者尽可能小的一种情况。

2. 常见的Beamformer（DSB，MVDR，GSC，GEV）

2.1. Delay-Sum Beamforming

我们会经常听到这个名字，它的原理是最简单的，只要知道入射方向，就能得到信号到达每个麦克风的时间差，只要将这些麦克风的信号在时间上对齐即可。
所以它的做法就是，首先对每一个麦克风的信号进行延时操作。
$y'(n) = \frac 1 N \sum_{i=0}^{M-1}x(n+\tau_i)$y′(n)=N1​i=0∑M−1​x(n+τi​)
那么相应地，它的窗函数就是$h_i(n) = \frac 1N \delta(n+\tau_i)$hi​(n)=N1​δ(n+τi​)，频域表示为$H(\omega) = \frac 1Ne^{-j \omega\tau_i}$H(ω)=N1​e−jωτi​

D-S 在很多地方都会使用，它虽然没有限制任何噪声方向，但是它很稳定。

2.2. Minimum Variance Distortionless Response (MVDR)

这个方法的核心思想就是，找到权值向量$\textbf w$w，保持住目标方向响应不变，同时让让输出信号功率最小，公式表达为
$\left \{ \begin{aligned} &min \quad \textbf w \textbf R_{xx} \textbf w \\ &s.t. \textbf w^H\textbf V_k(\textbf k_s) = 1 \end{aligned}\right .${​minwRxx​ws.t.wHVk​(ks​)=1​
这个方程的解法就是拉格朗日乘子法，得到的解为
$\textbf w = \frac {\textbf R_{xx}^{-1}\textbf V_k(\textbf k_s)}{\textbf V_k^H(\textbf k_s)\textbf R_{xx}^{-1}\textbf V_k(\textbf k_s)}$w=VkH​(ks​)Rxx−1​Vk​(ks​)Rxx−1​Vk​(ks​)​
保持目标方向不变，最小化输出功率就相当于最小化噪声功率，因此，也可以写成

$\textbf w = \frac {\textbf R_{nn}^{-1}\textbf V_k(\textbf k_s)}{\textbf V_k^H(\textbf k_s)\textbf R_{nn}^{-1}\textbf V_k(\textbf k_s)}$w=VkH​(ks​)Rnn−1​Vk​(ks​)Rnn−1​Vk​(ks​)​
这里面的$\textbf R_{nn}$Rnn​就是噪声自相关矩阵。
这个方法高度依赖于声源方向的估计，方向预测的准，效果就好。

2.3. Generalized Sidelobe Canceller (GSC)

我对GSC的理解就是，在DSB上减去了噪声项。给我的感觉有点像AEC的一些做法，不同的是，GSC的噪声是估计出来的。

可以直接来看一个GSC的实例，我觉得比理论还好理解一些。整个框架氛围三个模块：FBF、ABM、AIC。

1. 对于FBF可以使用DSB，这个实际上就要求我们去预先知道或者估计出时延信息或者入射方向信息。$y_f(n)$yf​(n)就是DS之后的结果。
2. 然后在通过自适应滤波的技术去更新ABM的各个滤波器，这些滤波器的作用就是屏蔽掉语音信号成分，保留噪声成分。由于滤波器希望对语音信号进行一个屏蔽，那么滤波器的更新条件就是语音信号存在。
3. 最后是AIC模块，这个模块就是通过前面每个麦克风的噪声估计，对噪声进行一个滤波，然后和DS的信号$y_f(n)$yf​(n)相减。就得到了最后干净的语音。那么这个模块为了准确估计，就需要在语音不存在的时候进行更新。

其中会提及到语音存在概率这个东西，传统做法会有Minima Controlled Recursive Averaging (MCRA)方法。神经网络方法会有估计mask，对mask进行sigmoid去估计语音存在概率。

现在再从理论上来说说GSC，它的思想是把权重$\textbf w$w投影在两个相互正交的平面上，且其中一个平面是我们定义的（称作“约束子空间”），另一个平面是正交补空间，被称为“自适应子空间”。$\textbf w = \textbf w_q - \textbf w_p$w=wq​−wp​，其中$\textbf w_q$wq​是约束子空间的投影，$\textbf w_p$wp​是自适应子空间的投影。

在上面的例子中，$\textbf w_q$wq​就是DSB的窗函数，$w_p$wp​就是ABM和AIC组合而成的。$w_p = \textbf B (\textbf B^H\textbf B)^{-1}\textbf B^H\textbf w = \textbf B \textbf w_a$wp​=B(BHB)−1BHw=Bwa​具体来说，ABM对应了$\textbf B$B,AIC对应了$\textbf w_a$wa​。两个变量都是通过自适应滤波器的方法求解得到的。

2.4. Generalized Eigen Value(GEV)

这个beamformer的目标函数为，最大化每个频点的信噪比
$\textbf w_{GEV} = arg \max_{w(k)}\frac{\textbf w^H(k)\textbf R_{zz}\textbf w(k)}{\textbf w^H(k)\textbf R_{nn}\textbf w(k)}$wGEV​=argw(k)max​wH(k)Rnn​w(k)wH(k)Rzz​w(k)​
这里我们用$z$z表示干净语音信号。

这个问题实际上是一个广义瑞丽商问题。求解方法依旧是拉格朗日成子法
我们可以定义$R(x) = \frac {x^TMx}{x^TQx}$R(x)=xTQxxTMx​，如果我们限制$x^TQx = 1$xTQx=1，那么就可以用拉格朗日成子法来求极值。
结果是极点为$x = Q^{-1}M$x=Q−1M，对应的极值为$Q^-1M$Q−1M的特征值

最后解出来的结果为$\textbf w_{GEV}(\omega)= \textbf R_{nn}^{-1}(\omega)\textbf R_{zz}(\omega)$wGEV​(ω)=Rnn−1​(ω)Rzz​(ω)。

3 声源定位算法

前面提到，阵列信号处理可以看成是空间、频率两个维度共同的响应。那么对于声源定位算法来说，就是查找适合的或者说是响应最大的空间方位。
所以很多空间定位算法，就是遍历一圈方向，找到最大的那个值，只不过是定义的值是不同的。

1. 如果我们对DS的信号功率，在方向上进行遍历，选择一个最大的方向，这个算法就叫SRP算法。
   $\textbf k_s = arg \max_{k} \sum_\omega|Y(\omega,\textbf k)|^2$ks​=argkmax​ω∑​∣Y(ω,k)∣2

2. 但是人们发现，更多时候，方向更多的时候是与相位相关的，与幅度无关，于是更倾向于使用SRP-PHANT的方法。由于两个麦克风的信号是高度相似的，所以可以对一个信号进行时间延迟（频域相移），计算相关函数，如果某个方向所有麦克风两两的相关响应最大，那么这个方向就是入射角。PHANT的方法使用的是广义互相关函数，也就是不考虑幅度大小。计算公式可以表示为：
   $p(\textbf k) = \sum_{i=0}^{M-1}\sum_{l=0}^{M-1}\sum_\omega\frac{X_i(\omega)X^*_l(\omega)}{|X_i(\omega)X^*_l(\omega)|}e^{j\omega\tau_{il}}$p(k)=i=0∑M−1​l=0∑M−1​ω∑​∣Xi​(ω)Xl∗​(ω)∣Xi​(ω)Xl∗​(ω)​ejωτil​
   前面的$\textbf k$k体现在$\tau_{il}$τil​上。
   这个$\tau_{il}$τil​就是TDOA。
   
   $\tau_{ij} = \frac{d cos\theta}{c}$τij​=cdcosθ​
   还有一个有意思的理解，可以分享一下，公式前面是$\frac{XX^*}{|XX^*|}$∣XX∗∣XX∗​，实际上就是相位差（IPD），我们谈绝对相位是没有意义的，就像我们提绝对时间也是没有意义的是一样的，总是拿一个麦克风作为参考麦克风再去谈相对时延或者相位差，这才是有意义的。也就是前面是信号算出来的相位差，后面的$\tau_{il}$τil​是我们估计的两个麦克风的到达时间差在相位的体现，去算没有归一化的cosine距离。

3. 我们如果不用DS去算，而是使用MVDR去进行计算，那么得到的就是Capon谱
   最优权重矩阵$\textbf w_c$wc​下的阵列输出功率为
   $P(\textbf k_s) = P_{capon} (\textbf k)= \textbf w_c^H\textbf R_{xx}(\omega)\textbf w_c = \frac 1{\textbf V_{\textbf k}^H(\textbf k_s)\textbf R_{xx}^{-1}\textbf V_{\textbf k}(\textbf k_s)}$P(ks​)=Pcapon​(k)=wcH​Rxx​(ω)wc​=VkH​(ks​)Rxx−1​Vk​(ks​)1​
   $\textbf k_s = arg \max_k P_{capon}(\textbf k)$ks​=argkmax​Pcapon​(k)

4.

现在在回过头来看看当时的google和amazon的做法，的确是很相似。google是让神经网络自己学一组steer vector，每个表示不同方向，在频域直接点乘。亚马逊是预设好几个方向的steer vector，这个是根据拓扑结构算出来的。后面就是方向选择问题，可以是拼一起，可以根据能量进行选择，或是attention。
在神经网络化的beamformer中，感觉大部分的神经网络的作用是用来预测speech或者是noise，或是预测T-F mask来得到speech和（或）noise。我们也看到了上面的很多公式中需要大量的speech或者是noise的自相关矩阵，这样估计出了mask之后就可以送给后面继续进行beamformer的估计。