Beam Pattern

在看前端方面的论文的时候，一定会遇到beam pattern这个概念，本篇文章主要讲述beam pattern的意义，应用以及可视化方面的东西。

beamformer很重要的一个性质就是它的指向性，指向性可以反映波束形成器对每个方向的信号，哪些频率进行增强，哪些频率进行抑制。如果可以准确的估计出DoA的话，增强声源方向的信号，抑制其他方向，就成为了beamformer的主要任务。

根据上述表述，beamformer要想表现的完美，准确的DoA不可或缺。如果可以实时的根据信号环境调整DoA，那么在工作中，beamformer的指向性便会发生改变，这种波束形成器被称为adaptive beamformer。反之，如果固定其指向性，在工作中只增强事先设定的方向信号，便被称为fixed beamformer。反映到beam weights上，如果beam weights被计算出并保持不变，对应后者，如果有声学环境的参数掺杂在其中，或者不断的update DoA，则对应前者，比如MVDR，PMWF，Max-SNR等等。

上篇文章中也提到了，用NN做前端，一般用mask模型估计adaptive beamformer中target的协相关/功率谱密度矩阵，并用其最大特征向量来估计steer vector。我们更加在意系统最终对WER或者PESQ，SDR，SNR等等（根据任务不同）的贡献，而不会刻意的去分析其指向性。这种不涉及麦克风拓扑结构，不涉及声源信息的方法称为Blind方法。

Beam Pattern

现在做如下定义：存在$M$个麦克风的线阵中，第$i$个麦克风相对第一个的距离为$d_i(d_0 = 0, d_i < d_j \;\text{if} \; i < j)$，如果DoA为$\theta(\theta \in [0, \pi])$，声速为$c$，那么，第$i$个麦克风相对第一个的延时$\tau_i$定义为：

$$\tau_{\theta i} = \frac{\cos(\theta)d_i}{c}$$

以MVDR为例，已知相对延迟$\tau$，beam weights由下式得出：

$$\mathbf{w}(\omega) = \frac{\mathbf{R}_{vv}^{-1}(\omega) \mathbf{d}_\tau(\omega)}{\mathbf{d}_\tau^H(\omega)\mathbf{R}_{vv}^{-1} (\omega)\mathbf{d}_\tau(\omega)}$$

其中：

$$\mathbf{d}_\tau = [e^{-j\pi \omega\tau_0}, e^{-j\pi \omega\tau_1}, \cdots, e^{-j\pi \omega\tau_{M - 1}}]^T$$

如果我们想要分析在该$\mathbf{w}, \mathbf{w} \in \mathbf{C}^{F \times M}​$下，beamformer对任意方向信号的作用，就可以借助beam pattern进行。

令$Z \in \mathbf{C}^{F \times T}$为beamformer输出的信号，那么有（详细见文章“Overview of beamformer”）：

$$Z(\omega) = \mathbf{w}(\omega)^H \mathbf{y}(\omega)$$

其中$\mathbf{y} = [Y_0, Y_1, \cdots, Y_{M-1}] \in \mathbf{C}^{F \times M \times T}$。

把假设DoA为$\theta$下计算的beam weights表示为$\mathbf{w}(\omega, \theta)$：

$$\mathbf{w}_\theta(\omega) = \frac{\mathbf{R}_{vv}^{-1}(\omega) \mathbf{d}_\theta(\omega)}{\mathbf{d}_\theta^H(\omega)\mathbf{R}_{vv}^{-1}(\omega) \mathbf{d}_\theta(\omega)}$$

$\mathbf{d}_\phi​$表示DoA为$\phi​$时的steer vector：

$$\mathbf{d}_\phi = [e^{-j\pi \omega\tau_{\phi 0}}, e^{-j\pi \omega\tau_{\phi 1}}, \cdots, e^{-j\pi \omega\tau_{\phi{M - 1}}}]^T$$

考虑到

$$\mathbf{y}(\omega) = \mathbf{d}_\phi(\omega) X(\omega) + \mathbf{v}(\omega)$$

那么在beam weights为$\mathbf{w}(\omega, \theta)$下，有

$$Z(\omega) = \mathbf{w}_\theta(\omega)^H \mathbf{y}(\omega) = \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega) X(\omega) + \mathbf{w}_\theta(\theta)^H\mathbf{v}(\omega)$$

把上式中$X(\omega)$的系数称为beam pattern $\mathcal{B}(\omega, \phi)$：

$$\mathcal{B}(\omega, \phi) = \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega)$$

从这里可以看出，beam pattern是频率$\omega​$和DoA $\phi​$的函数，$\mathcal{B}(\omega, \phi)​$越大，表示对应方向上和频率处的信号增益越大。也就是说，通过beam pattern我们可以看出一组空间滤波器在方位轴和频率轴上的表现。这是分析滤波器性能的一个很方便的手段。

Fixed Beamformer Design

根据beam pattern的定义，我们可以由$\mathcal{B}(\omega, \phi)$反推$\mathbf{w}_\theta(\omega)$。

假设我们想要增强$\theta$方向的信号，对应的在beam pattern中反映的应该是$\delta \in [\theta - \epsilon,\theta +\epsilon ]$处的beam pattern较大。将其作为beam weights的target，在频率$\omega$处，令：

$$\mathcal{B}_{\text{target}}(\omega, \delta) = 1$$

现在，想要求出$\mathbf{w}_\theta(\omega)​$，使得 $ \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega)​$的结果接近$\mathcal{B}_{\text{target}}​$，转化为下面的优化问题：

$$\mathbf{w}_\theta(\omega) = \arg \min_{\mathbf{w}_\theta(\omega)} \Vert \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega) - \mathcal{B}_{\text{target}}(\omega) \Vert_2^2$$

$\mathbf{d}_\phi(\omega) $的定义和上一部分相同，表示假设DoA为$\phi$时，算出的steer vector矩阵。和MVDR一样，需要限制DoA方向的语音不失真，即：

$$\mathbf{w}_\theta(\omega)^H\mathbf{d}_\theta(\omega) = 1$$

求解上述问题即可得到指向$\theta$方向的beam weights。不过这样求出的beam weights在低频处权值很大，往往上百上千，beamforming之后，低频信号放大的很厉害。通过加入WNG（白噪声增益，White Noise Gain）约束可以缓解这种情况的发生。

WNG定义如下：

$$A(\omega) = \frac{|\mathbf{w}_\theta(\omega)^H\mathbf{d}_\theta(\omega)|^2}{\mathbf{w}_\theta(\omega)^H\mathbf{w}_\theta(\omega)}$$

约束$A(\omega) \geqslant \gamma$可以改写为：

$$\mathbf{w}_\theta(\omega)^H\mathbf{w}_\theta(\omega) = \Vert\mathbf{w}_\theta(\omega) \Vert_2^2 \leqslant 1/\gamma$$

通常设置$20 \log10(\gamma) = 0, -10 \text{dB}$。

因此，Fixed Beamformer的设计流程如下：

1. 根据DoA，定义$\mathcal{B}_{\text{target}}(\omega, \delta)$。

2. 求解凸优化问题：

   $$\begin{align} \mathbf{w}_\theta(\omega) &= \arg \min_{\mathbf{w}_\theta(\omega)} \Vert \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega) - \mathcal{B}_{\text{target}}(\omega) \Vert_2^2 \\ & \text{s.t} \quad \mathbf{w}_\theta(\omega)^H\mathbf{d}_\theta(\omega) = 1, \; \Vert\mathbf{w}_\theta(\omega) \Vert_2^2 \leqslant 1/\gamma \end{align}$$

3. 分析$\mathcal{B}(\omega, \phi) = \mathbf{w}_\theta(\omega)^H\mathbf{d}_\phi(\omega)$验证

Visualize Beam Pattern

将$\theta = 60^\circ$时的$\mathcal{B}(\omega, \phi)$整体绘制，纵轴表示DoA，横轴表示频率

黄色区域表示响应越大，即这个部分的语音gain越多，可以看出，60度方向，全频带范围内gain值都很高，和设计初衷相符。如果固定频率比如5kHz，画出极坐标图，可以得到下面的图（图的幅值方向单位为dB）：

在信号处理中，$\mathbf{w}_\theta(\omega)$可以准确的求出，因此，beam pattern是可以求出解析的表达式的，那么往往喜欢用极坐标图的形式显示。

总结下来，这种套路应该是信号那边常见的方法，将beam weights的设计转为凸优化问题进行求解。实际应用中，不同DoA下设计的weights作用在相同的输入上，结果应该会有明显的不同。

另外，需要注意$\omega \in [0, 2 \pi]$这个量，角频率，和频率$f$的关系是：

$$\omega = 2\pi f$$

DFT中，时域，频域采样数相同，设为$N$，$f$为信号的采样频率，那么

$$\omega_k = 2 \pi \frac{kN}{T} = 2\pi f_k，f_k = \frac{kN}{T} = 2\pi f \quad 0 \leqslant k < N$$

在短时域中，考虑DFT的共轭对称性，只保留前$\frac{N}{2} +1$个点，即$0 \leqslant k \leqslant \frac{N}{2} +1$。

具体编程实现的时候，需要注意它们之间的转换关系。