多通道混响抑制：GWPE(Generalized Weighted Prediction Error)

1.模型构建

考虑一个M个说话人，N个麦克风的场景。
$s^m[k](1\le m\le M)$：第m个说话人k时刻的信号。
$y^n[k]$和$v^n[k](1\le n\le N)$：第n个麦克风k时刻的接收信号和噪声。
则有：
$$\mathbf{y}[k]=\sum _{\tau =0}^{J-1}{\mathbf{H}}^T[\tau ]\mathbf{s}[k-\tau ]+\mathbf{v}[k]（1）$$
其中
$\mathbf{y}[k]={\left[y^1[k],\dots ,y^N[k]\right]}^T$,$\mathbf{s}[k]={\left[s^1[k],\dots ,s^M[k]\right]}^T$,$\mathbf{v}[k]={\left[v^1[k],\dots ,v^N[k]\right]}^T$

$\mathbf{H}[\tau ]$为房间的脉冲响应矩阵，J表示响应的阶数。
$$\mathbf{H}[\tau ]={\left[\begin{array}{ccc}{h}^{1,1}[\tau ]& \cdots & h^{N,1}[\tau ]\\ ⋮& \ddots & ⋮\\ h^{1,M}[\tau ]& \cdots & h^{N,M}[\tau ]\end{array}\right]}_{M\times N}$$
{ h n , m [ τ ] } 0 ≤ τ ≤ J − 1 \left\{h^{n, m}[\tau]\right\} 0 \leq \tau \leq J-1 {hn,m[τ]}0≤τ≤J−1为说话人m与麦克风n之间的冲击响应。
式（1）其实就是卷积的矩阵形式。
为了减少滤波器的阶数，原论文中还对$\mathbf{y}[k]$进行了子带分解，本文不做考虑。

我们希望能够减少混响响应的长度来对混响效果进行削弱。
将$\mathbf{y}[k]$写为：
$$\mathbf{y}[k]=\sum _{\tau =0}^{\Delta }{\mathbf{H}}^T[\tau ]\mathbf{s}[k-\tau ]+\sum _{\tau =\Delta }^{J-1}{\mathbf{H}}^T[\tau ]\mathbf{s}[k-\tau ]+\mathbf{v}[k]$$
希望削减混响响应的长度到$\Delta$，即只保留第一项。

首先采用多通道线性预测技术用$\Delta -(J-1)$的麦克风接收信号去估计$\Delta -(J-1)$部分的混响：
$$\tilde{\mathbf{y}}(k)=\sum _{\tau =\Delta }^{\Delta +K-1}{G}^T(\tau )\mathbf{y}(k-\tau )$$
其中 { G ( τ ) } Δ ≤ τ ≤ Δ + K − 1 \left\{G(\tau)\right\}_{\Delta \leq \tau \leq \Delta+K-1} {G(τ)}Δ≤τ≤Δ+K−1​为$N\times N$的预测矩阵，由于实际中混响的效果随着时间逐渐减弱，所以只需要估计长度为K的混响信号。

然后与$\mathbf{y}[k]$作差：
$$\mathbf{x}(k)=\mathbf{y}(k)-\tilde{\mathbf{y}}(k)(2)$$
如果最小化（2），由于$\tilde{y}(k)$只与$y(k)$的$\Delta \le \tau \le \Delta +K-1$的部分相关，所以会对消掉。
理想的预测误差$\mathbf{x}(k)$等于
$${\mathbf{x}}_{}(k)=\sum _{\tau =0}^{\Delta -1}{\mathbf{H}}^T(\tau )\mathbf{s}(k-\tau )+\text{ noise }(3)$$

2.GWPE损失函数定义

假设我们收集T点的预测误差数据${\left(\mathbf{x}(k)\right)}_{k\in T}$，干净的语音信号在延迟时间大于几十毫秒后，其自相关系数基本为0，但是混响信号在延迟时间较长时自相关系数仍然很大，所以GWPE损失函数就是通过最小化T点麦克风阵列接收信号的相关性来达到去混响的目的。

• GWPE使用Hadamard-Fischer mutual correlation作为损失函数，其定义如下：
  $$C_{HF}\left({\mathbf{U}}_1,\dots ,{\mathbf{U}}_N\right)=\frac{1}{N}\sum _n\log \left(\det E\left({\mathbf{U}}_n{\mathbf{U}}_n^H\right)\right)-\log \left(\det E\left(\mathbf{U}{\mathbf{U}}^H\right)\right)$$
  其中${\mathbf{U}}_1,\dots ,{\mathbf{U}}_N$表示复随机变量，当${\mathbf{U}}_1,\dots ,{\mathbf{U}}_N$互不相关时：
  $$C_{\mathrm{HF}}\left({\mathbf{U}}_1,\dots ,{\mathbf{U}}_N\right)=0$$

所以我们就可以通过上述损失函数来最小化${\left(\mathbf{x}(k)\right)}_{k\in T}$的相关性。
令${\mathbf{x}}_T={\left[{\mathbf{x}}^T(T),\dots ,{\mathbf{x}}^T(1)\right]}^T$损失函数可以写为：
$\begin{array}{rl}F\left(\mathcal{G}\right)=& C_{\mathrm{HF}}\left(\mathbf{x}(1),\dots ,\mathbf{x}(T)\right)\\ =& \frac{1}{|T|}\sum _{t\in T}\log \left(\det E\left({\mathbf{x}}_T(t){\mathbf{x}}_T^T(t)\right)\right)-\log \left(\det E\left({\mathbf{x}}_T{\mathbf{x}}_T^T\right)\right)\end{array}（4）$
其中 G = { G ( τ ) } Δ ≤ τ ≤ Δ + K − 1 \mathcal{G}=\left\{\boldsymbol{G}(\tau)\right\}_{\Delta \leq \tau \leq \Delta+K-1} G={G(τ)}Δ≤τ≤Δ+K−1​.

损失函数第二项实际上是常数，证明如下：
我们将式（2）写为矩阵形式：
$$x_T=G\cdot y_T$$
其中${\mathbf{y}}_T={\left[{\mathbf{y}}^T(T),\dots ,{\mathbf{y}}^T(1)\right]}^T$,
G = { I , O , … , O , − G ( Δ ) , … , − G ( Δ + K l − 1 ) } G=\left\{\boldsymbol{I}, \boldsymbol{O}, \ldots, \boldsymbol{O},-\boldsymbol{G}(\Delta), \ldots,-\boldsymbol{G}\left(\Delta+K_{l}-1\right)\right\} G={I,O,…,O,−G(Δ),…,−G(Δ+Kl​−1)}因此有：$$\begin{array}{rl}\det E\left({\mathbf{x}}_T{\mathbf{x}}_T^T\right)& ={\left|\det \mathbf{G}\right|}^2\det E\left({\mathbf{y}}_T{\mathbf{y}}_T^T\right)\\ & =\det E\left({\mathbf{y}}_T{\mathbf{y}}_T^T\right)\\ & =\text{ constant. }\end{array}$$
（上标T表示转置/共轭转置）所以式（4）的第二项可以去掉。
$$\begin{array}{rl}F\left(\mathcal{G}\right)=& \frac{1}{|T|}\sum _{t\in T}\log \left(\det E\left({\mathbf{x}}_T(t){\mathbf{x}}_T^T(t)\right)\right)\end{array}（5）$$
上式就是GWPE的损失函数，现在混响抑制的目标可以写为：
$${\mathcal{G}}_{}=\arg \underset{{\mathcal{G}}_{}}{\min }F\left({\mathcal{G}}_{}\right)$$

3.优化算法

最小化（5）来获得$\mathcal{G}$的估计，但是由于（5）没有解析解，所以论文中给出辅助函数的方法，构造辅助函数：
$$\begin{array}{r}\tilde{F}\left(\mathcal{G},\mathcal{L}\right)=\frac{1}{|T|}{\sum }_{t\in T}(E\left(x^T(t)\mathbf{\Lambda }(t{)}^{-1}{\mathbf{x}}_T(t)\right)-N+\log \left(\det \mathbf{\Lambda }(\mathbf{t})\right))\end{array}$$
其中$$\mathbf{\Lambda }=E\left({\mathbf{x}}_T{\mathbf{x}}_T^T\right)$$
为空间相关矩阵，$\mathbf{\Lambda }$通常通过时间平均的方法得到
$$\mathbf{\Lambda }=\sum _{k=t-\delta }^{t+\delta }\frac{1}{2\delta +1}{\mathbf{x}}_T{\mathbf{x}}_T^T$$
然后估计$\mathcal{G}$：
$$\begin{array}{rl}\hat{\mathbf{R}}& =\sum _{t\in \mathcal{T}}\overline{\mathbf{\psi }(t-\Delta )}\hat{\mathbf{\Lambda }}(t{)}^{-1}\overline{{\mathbf{\psi }}^{\ast }(t-\Delta )}\\ \hat{\mathbf{r}}& =\sum _{t\in \mathcal{T}}\overline{\mathbf{\psi }(t-\Delta )}\hat{\mathbf{\Lambda }}(t{)}^{-1}\mathbf{y}(t)\end{array}$$
$$\hat{\mathbf{g}}={\hat{\mathbf{R}}}^{-1}\hat{\mathbf{r}}$$
$g$为预测矩阵 G = { G ( τ ) } Δ ≤ τ ≤ Δ + K − 1 \mathcal{G}=\left\{\boldsymbol{G}(\tau)\right\}_{\Delta \leq \tau \leq \Delta+K-1} G={G(τ)}Δ≤τ≤Δ+K−1​的列向量拼接起来的向量：
$$\mathbf{g}=\left[\begin{array}{c}{\mathbf{g}}^1(\Delta )\\ ⋮\\ {\mathbf{g}}^N(\Delta )\\ ⋮\\ {\mathbf{g}}^1\left(\Delta +K-1\right)\\ ⋮\\ g^N\left(\Delta +K-1\right)\end{array}\right]$$
其中${\mathbf{g}}^n(\tau )$代表$\mathbf{G}(\tau )$的第n列。