基于线性预测的实时去混响以及降噪——利用交替的卡尔曼滤波

基于线性预测的实时去混响以及降噪——利用交替的卡尔曼滤波

2. 信号模型建立以及问题表述

​ 假设混响环境下有未知数量的声源，使用固定在任意位置的$M$个子麦进行采集。给出采集信号的$stft$域的表达：
$$\mathbf{y}(k,n)=[Y_1(k,n),\cdots ,Y_M(k,n){]}^T$$
​ 其中，$Y_m(k,n)$是第$m$路信号第$k$个子频带，第$n$帧的频域表达。我们假设多通道麦克风信号有两部分组成：
$$y(k,n)=x(k,n)+v(k,n)$$
​ 其中，向量$x(k,n)$以及$v(k,n)$分别表示阵列上麦克风采集的混响语音信号以及加性噪声。

A. 多通道自回归混响模型

​ 假设混响信号的产生能够在各个子频带上独立，即就是采集信号中的混响部分的某一帧的某个子频带可以通过且仅需要该子频带的过去几帧来预测，给出表达：
$$x(k,n)=\underset{r(k,n)}{\underbrace{\sum _{l=D}^L{C}_l(k,n)x(k,n-l)}}+s(k,n)$$
​ 其中，向量$s(k,n)=[S_1(k,n),\cdots ,S_M(k,n){]}^T$表示采集信号中希望获取的直达声和早期混响信号的$stft$域，$S_m(k,n)$表示第$m$个子麦的第$n$帧第$k$个子频带频域表达。矩阵$C_l(k,n)\in C^{M\times M}$表示对于第$n-l,l\in [D,D+1,\cdots ,L]$帧的采集信号$stft$域$x(k,n-l)\in C^{M,1}$的滤波参数。该自回归模型中的早期混响语音信号$s(k,n)$（同时也可以认为是线性预测算法中的预测误差）是卡尔曼滤波过程中的新息。延迟参数$D$的选择决定了期望信号中我们希望保留的早期混响的量，并且需要基于$STFT$的帧选择的重叠率来调整来确保期望信号$s(k,n)$中的直达声成分和晚期混响信号$r(k,n)$没有相关性。

​ 我们假设$s(k,n)\sim N(0_{M\times 1},{\Phi }_s(k,n))$且噪声向量$v(k,n)\sim N(0_{M\times 1},{\Phi }_v(k,n))$服从零均值多维复高斯分布的随机向量，它们的协方差矩阵分别为 Φ s ( k , n ) = E { s ( k , n ) s H ( k , n ) } \Phi_s(k,n)=E\{s(k,n)s^H(k,n)\} Φs​(k,n)=E{s(k,n)sH(k,n)}以及 Φ v ( k , n ) = E { v ( k , n ) v H ( k , n ) } \Phi_v(k,n)=E\{v(k,n)v^H(k,n)\} Φv​(k,n)=E{v(k,n)vH(k,n)}；除此之外，我们假设$s(k,n)，v(k,n)$在不同时刻是不相关的。

B. 建立信号模型的两种紧凑形式表达

​ 为了建立一个损失函数（下文会介绍，实际是分解成两个子损失函数，分别使用两个交替卡尔曼滤波来求解）。我们首先引入两种结果相同，写法不同的矩阵表达式来描述公式2的麦克风采集信号向量。为了描述更加紧凑，论文后续省略掉频段索引$k$；我们首先来看第一种表达方式：
$$X(n)=I_M\otimes [x^T(n-L+D)\dots x^T(n)]$$

$$c(n)=Vec\{[C_L(n)\dots C_D(n){]}^T\}$$

其中，$I_M$是$M$维的单位矩阵，$\otimes$表示Kronecker乘积，定义如下：

V e c { . } Vec\{.\} Vec{.}是矩阵拉直操作，表示将大括号内的矩阵的列从左到右按顺序首位拼接起来，得到的一个新的列向量。因此，$c(n)$是一个长度为$L_c=M^2(L-D+1)$的列向量，$X(n)$是一个形状为$M\times L_c$的稀疏矩阵，基于定义式4以及式5，结合信号模型式2和式3。给出麦克风观测信号：
$$y(n)=\underset{r(n)}{\underbrace{X(n-D)c(n)}}+\underset{u(n)}{\underbrace{s(n)+v(n)}}$$
其中，$u(n)$包含了早期语音混响加上噪声信号。得到$u(n)$的方差矩阵为 Φ u ( k , n ) = E { u ( k , n ) u H ( k , n ) } ～ N ( 0 M × 1 , Φ u ( k , n ) ) \Phi_u(k,n)=E\{u(k,n)u^H(k,n)\}～N(0_{M×1},\Phi_u(k,n)) Φu​(k,n)=E{u(k,n)uH(k,n)}～N(0M×1​,Φu​(k,n))

第二个表达式利用了带下划线表示的拼接向量：
$$\underline{x}(n)=[x^T(n-L+1)\dots x^T(n){]}^T$$

$$\underline{s}(n)=[0_{1\times M(L-1)}{s}^T(n){]}^T$$

式7和式8表示的列向量长度均为$ML$，给出对应卡尔曼滤波中的状态传播矩阵以及观测矩阵（状态传播矩阵：状态量经过传播矩阵左乘后得到对下一时刻状态量的预估；观测矩阵：状态量通常是隐变量，不能直接被观测到，将状态量左乘观测矩阵得到对观测量的估计），分别表示为：
$$F(n)=\left[\begin{array}{ccc}{0}_{M(L-1)\times M}& & I_{M(L-1)}\\ C_L(n)& & C_D(n)0_{M\times M(D-1)}\end{array}\right]$$

$$H=[0_{M\times M(L-1)}{I}_M]$$

其中，形状为$ML\times ML$的状态传播矩阵$F(n)$的底部$M$行包含了自回归参数$C_l(n)$，从$F(n)$的分块部分可以看出，上面的$M(L-1)$行的目的是将$\underline{x}(n)$左移一个时间单位，同时下面$M$行包含了自回归参数，利用了$(n-L+1)$时刻到$(n-D)$时刻的无噪声混响信号$x(.)$以及自回归参数来预测当前时刻的无噪混响信号并放入$\underline{x}(n)$的最右侧。观测矩阵$M$是一个选择矩阵，它将$\underline{x}$的最右侧对应$n$时刻的信号$x(n)$取出。结合式9和式10，我们可以将式3和式2用另一种方式描述：
$$\underline{x}(n)=F(n)\underline{x}(n-1)+\underline{s}(n)$$

$$y(n)=H\underline{x}(n)+v(n)$$

注意到式6和式12表达含义相同但是表达方式不同。

C. 自回归参数的随机状态空间模型

考虑到自回归参数有时变特性以及$STFT$模型引入的误差，我们使用一阶Markov模型对自回归参数建模：
$$c(n)=Ac(n-1)+w(n)$$
我们假设传递矩阵$A=I_{L_C}$是一个单位阵。同时，过程噪声$w(n)$用来表征参数$c(n)$在不同时刻的不确定性。我们假设$w(n)～N(0_{M\times 1},{\Phi }_w(n))$服从零均值，方差矩阵为$\Phi w(n)$的复高斯分布。同时，$w(n)$在各个时刻独立且和$u(n)$不相关。

图1展示了麦克风观测信号的产生过程以及潜在的混响信号$x(n)$和自回归参数$c(n)$的交互过程。

D. 问题模型建立

​ 我们的目标是预估多通道早期混响信号$s(n)$。我们提出的方法是首先预估无噪声的混响信号$x(n)$以及多通道自回归参数$c(n)$而不是直接去预估早期混响语音$s(n)$。我们在符号上面加个帽子符号$\hat{X}$来表示对参数的预估，例如$\hat{x}(n),\hat{c}(n)$分别表示对$x(n),c(n)$的预估；接下来，我们可以通过使用预估的自回归参数以及无噪声混相信号FIR滤波得到晚期混响信号$\hat{r}(n)$，将预估的晚期混相信号从预估的无噪声混响信号中减去得到我们期望预估的早期混响信号$\hat{s}(n)$，即就是：
$$\hat{s}(n)=\hat{x}(n)-\hat{X}(n-D)\hat{c}(n)$$
​ 其中$\hat{X}(n)$由式4定义，$\hat{r}(n)$认为是预估的晚期混响信号。下一节我们将介绍如何同时预估$x(n)$以及$c(n)$。

3. 通过交替最小化进行最小均方差估计

拼接的混响语音信号向量$\underline{x}(n)$以及自回归参数向量$c(n)$（同时也是$F(n)$内部的一部分）可以通过最小均方差算法来预估；损失函数定义为：
J ( x ‾ , c ) = E { ∣ ∣ x ‾ ( n ) − F ^ ( n ) x ‾ ^ ( n − 1 ) + s ‾ ^ ( n ) ∣ ∣ 2 2 } J(\underline{x}, c)=E\{||\underline{x}(n)-\hat{F}(n)\hat{\underline{x}}(n-1)+\hat{\underline{s}}(n)||^2_2\} J(x​,c)=E{∣∣x​(n)−F^(n)x​^​(n−1)+s​^​(n)∣∣22​}
注意到损失函数有两个因变量分别为晚期混响信号$\underline{x}$以及自回归参数$c$，同时调整两个因变量来最小化该损失函数较为困难。为了简化预估问题来得到问题的闭解，我们使用一种交替最小化的技术[26]来重新定义两个子损失函数。对于两个因变量$a,b$的损失函数，该技术首先使得某个因变量$a$固定，然后寻找另一个因变量$b$使得损失函数最小；紧接着固定因变量$b$，同样寻找另一个因变量$a$来使得损失函数最小；以此往复，交替进行。因此，我们定义两个子损失函数：
J c ( c ( n ) ∣ x ‾ ( n ) ) = E { ∣ ∣ c ( n ) − c ^ ( n ) ∣ ∣ 2 2 } J_c(c(n)|\underline{x}(n)) = E\{||c(n)-\hat{c}(n)||^2_2\} Jc​(c(n)∣x​(n))=E{∣∣c(n)−c^(n)∣∣22​}

J x ( x ‾ ( n ) ∣ c ( n ) ) = E { ∣ ∣ x ‾ ( n ) − x ‾ ^ ( n ) ∣ ∣ 2 2 } J_x(\underline{x}(n)|c(n))=E\{||\underline{x}(n)-\hat{\underline{x}}(n)||^2_2\} Jx​(x​(n)∣c(n))=E{∣∣x​(n)−x​^​(n)∣∣22​}

注意到第$n$帧时的式16以及信号模型式6，由于$n$时刻的麦克风采集信号$y(n)$ 仅依赖于使用过去延迟大于$D$的晚期混响信号$x(n)$来构建的$X(n-D)$；因此我们可以将式16改为$J_c(c(n)|\underline{x}(n))=J_c(c(n)|\underline{x}(n-D)) $；

使用$n$时刻可用的预估值替换损失函数16, 17的不变量依赖$\underline{x}(n)$以及$c(n)$，我们就能得到每一个时刻$n$的交替最小化过程：
$$1)\hat{c}(n)=\underset{c}{\mathrm{argmin}}{J}_c(c(n)|\underline{x}(n-D))$$

$$2)\hat{\underline{x}}(n)=\underset{\underline{x}}{\mathrm{argmin}}{J}_x(\underline{x}(n)|\hat{c}(n))$$

这里需要特别注意的是求解顺序，必须先求解$\hat{c}$再求解$\hat{\underline{x}}(n)$；这是因为自回归参数的时变的，如果先使用上一时刻的自回归参数$\hat{c}(n-1)$来预估$\hat{\underline{x}}(n)$，近似误差很大并且会影响到$\hat{\underline{x}}(n)$的估计。所以先利用过去的混响信号来估计$\hat{c}(n)$，再利用更新过的$\hat{c}(n)$来预估$\hat{\underline{x}}(n)$。尽管全局损失函数15很难进行收敛求解，但是如果式16, 17可以独自交替地最小化，全局损失函数可以收敛于局部最小点。对于给定模型，式16，17可以通过卡尔曼滤波来求解。

最终用来预估期望信号向量$s(n)$遵循以下三个步骤，标识在图2上，说明如下：

1）基于麦克风采集的带噪混响语音$y(n)$以及延迟的无噪声信号 x ( n ′ ) , n ′ ∈ { 1 , … , n − D } x(n'), n'\in\{1, \dots,n-D\} x(n′),n′∈{1,…,n−D}预估自回归参数$c(n)$，（带噪混响语音$y(n)$以及延迟的无噪声混响信号$x(n^{\prime })$都是已知且固定的）。具体实现中，将这些依赖量使用通过步骤2中的第二个卡尔曼预估值$\hat{x}(n^{\prime })$来代替。

2）利用自回归模型来预估无噪声混响信号$\underline{x}(n)$. 该步骤可以认为是消噪阶段。注意到这里的自回归参数$c(n)$被假设为不变的以及已知的。具体实现中，自回归参数通过步骤1的预估来获得。

3）利用预估的自回归参数$\hat{c}(n)$以及延迟的无噪声混响信号$\hat{x}(n)$，我们可以得到对晚期混响信号$r(n)$的预估。将预估的无噪声语音信号减去预估的晚期混响信号，就可以得到期望的无噪声语音早期混响信号。

步骤2中消噪阶段要求二阶噪声统计量${\Phi }_v(n)$已知，表示在图2中的灰色方块。由于目前存在很多优秀的二阶噪声统计量估计方法，例如文献[28]-[30]，因此对噪声的统计估计不在本论文的讨论范围内，并且接下来都假设噪声的统计特征${\Phi }_v(n)$是已知的。

本文提出的结构解决了顺序结构中通常面临的因果问题：每一个预测步骤需要其他变量的当前帧预估。图三展示了该模型，其中的噪声消除阶段利用的是延迟的自回归参数。该结构适用与自回归参数稳定或者随时间变化缓慢的情形，但是在时变的自回归参数模型中通常存在很大的误差。

A. 最优顺序结构——预估自回归参数

给定了图二展示的延迟晚期混响信号$x(n)$，我们本节来推导用来预估自回归参数的卡尔曼滤波器。

1）预估自回归参数的卡尔曼滤波器：假设$n$时刻下，我们已经有了包含在矩阵$X(n-D)$中的过去一段时刻的无噪声混响信号。接下来，我们考虑式13和式6分别作为状态方程和观测方程。假设过程噪声$w(n)$以及观测噪声$u(n)$符合零均值复高斯分布，并且两者不相关。我们可以通过最小化状态向量估计误差矩阵来得到对自回归参数向量的最优顺序估计，给出状态向量误差方差矩阵（实际上就是卡尔曼滤波最小化的损失函数）：
$${\Phi }_{{\Delta }_c}(n)=E\{[c(n)-\hat{c}(n)][c(n)-\hat{c}(n){]}^H\}$$
得到的解即就是卡尔曼滤波黄金五条：

$${\hat{\Phi }}_{{\Delta }_c}(n|n-1)=A{\hat{\Phi }}_{{\Delta }_c}(n-1)A^H+{\Phi }_w(n)$$

$$\hat{c}(n|n-1)=A\hat{c}(n-1)$$

$$e(n)=y(n)-X(n-D)\hat{c}(n|n-1)$$

$$\begin{gathered} K(n)={\hat{\Phi }}_{{\Delta }_c}(n|n-1)X^H(n-D) \\ [X(n-D){\hat{\Phi }}_{{\Delta }_c}(n|n-1)X^H(n-D)+{\Phi }_u(n){]}^{-1} \end{gathered}$$

$${\hat{\Phi }}_{{\Delta }_c}(n)=[I_{L_c}-K(n)X(n-D)]{\hat{\Phi }}_{{\Delta }_c}(n|n-1)$$

$$\hat{c}(n)=\hat{c}(n|n-1)+K(n)e(n)$$

其中$K(n)$为卡尔曼增益，$e(n)$为预测误差。注意到预测误差即就是利用预估自回归参数$\hat{c}(n|n-1)$得到的语音的早期混响加上噪声，因此预测误差$e(n)$实际上就是利用上一时刻的$u(n-1)$对此时刻$u(n)$的估计，即$e(n)=u(n|n-1)$.

2）参数估计：仅包含延迟的混响信号$x(n)$的矩阵$X(n-D)$是通过$III-B$节中的第二个卡尔曼滤波器来估计的。

我们假设状态传递矩阵为单位阵$A=I_{L_c}$，并且过程噪声的方差矩阵${\Phi }_w(n)={\varphi }_w(n)I_{L_c}$，其中的${\varphi }_w(n)$通过下式计算：
$${\hat{\varphi }}_w(n)=\frac{1}{L_c}||\hat{c}(n)-\hat{c}(n-1)|{|}_2^2+\eta$$
其中$\eta$是一个极小的正数，当连续两帧预估的状态向量完全相同时，增加一个极小数来表征自回归参数的连续变化特性。

观测噪声方差矩阵${\Phi }_u(n)$可以通过给定概率密度函数$f(y(n)|\hat{\Theta }(n))$后，利用极大似然法来预估。其中先验参数 Θ ^ ( n ) = { x ^ ( n − L , … , x ^ ( n − 1 ) , c ^ ( n ) } \hat{\Theta}(n)=\{\hat{x}(n-L, \dots, \hat{x}(n-1), \hat{c}(n)\} Θ^(n)={x^(n−L,…,x^(n−1),c^(n)}为当前已知的预估参数。假设$N$帧内，${\Phi }_u(n)$具有稳定性，极大似然法通过当前可用的信息得到的估计为：
$${\hat{\Phi }}_u^{ML}(n)=\frac{1}{N}(\sum _{l=N-1}^1\hat{u}(n-l){\hat{u}}^H(n-l)+e(n)e^H(n))$$
其中$\hat{u}(n)=y(n)-\hat{X}(n-D)\hat{c}(n)$；$e(n)=u(n|n-1)$是预测的噪声加语音信号。由于$n$时刻$\hat{c}(n)$还没有计算出来，因此我们无法得到后验预估值$\hat{u}(n)$，因此式28求和项上限为1，同时使用$e(n)=u(n|n-1)$来近似$\hat{u}(n)$.

在具体实现中，式28的算数均值可以用迭代平均来代替，由此得到迭代版本的极大似然估计
$${\hat{\Phi }}_u(n)=\alpha {\hat{\Phi }}_u^{pos}(n-1)+(1-\alpha )e(n)e^H(n)$$
其中，后验方差估计${\hat{\Phi }}_u^{pos}(n)$依赖先前帧的迭代平滑得到，有：
$${\hat{\Phi }}_u^{pos}(n)=\alpha {\hat{\Phi }}_u^{pos}(n-1)+(1-\alpha )\hat{u}(n){\hat{u}}^H(n)$$
其中的迭代均值参数$\alpha =e^{-\frac{\Delta t}{\tau }}$依赖于以秒为单位指数平滑常数$\tau$以及以秒为单位帧偏移参数$\Delta t$.由于$u(n)$可以仅仅可以在一小段时间内（几帧内）是稳态的，式29的近似估计对极大似然的估计的影响微乎其微。

B. 最优顺序结构——降噪

给定了当前时刻自回归参数$c(n)$;本节我们来推导第二个卡尔曼滤波来预估无噪声的混响信号$\underline{x}(n)$.

1）利用卡尔曼滤波降噪：假定自回归参数$c(n)$以及矩阵$F(n)$（由$c(n)$组成）是已知的，由$x(n)$的最新L帧拼接构成的混响信号向量$\underline{x}(n)$认为是状态向量，我们将式11和式12认为是状态和观测方程。由于假设语音信号$s(n)$以及（8）的定义，$\underline{s}(n)$也是零均值的高斯随机变量并且它的协方差矩阵 Φ s ‾ ( n ) = E { s ‾ ( n ) s ‾ H ( n ) } \Phi_{\underline{s}}(n)=E\{\underline{s}(n)\underline{s}^H(n)\} Φs​​(n)=E{s​(n)s​H(n)}，其中右下角包含了${\Phi }_s(n)$，其他位置则全是零。

已知了$\underline{s}(n)$以及$v(n)$是零均值高斯噪声且相互没有相关性，我们通过最小化状态向量误差矩阵来得到对$\underline{x}(n)$的最优估计，状态向量误差矩阵有：
$${\Phi }_{{\Delta }_x}(n)=E\{[\underline{x}(n)-\widehat{\underline{x}(n)}][\underline{x}(n)-\hat{\underline{x}}(n){]}^H\}$$
给出卡尔曼滤波的预测方程：
$${\hat{\Phi }}_{{\Delta }_x}(n|n-1)=F(n){\hat{\Phi }}_{{\Delta }_x}(n-1)F^H(n)+{\Phi }_{\underline{s}}(n)$$

$$\hat{\underline{x}}(n|n-1)=F(n)\hat{\underline{x}}(n-1)$$

更新方程：
$$\begin{array}{c}\begin{array}{rl}{K}_x(n)=& {\hat{\Phi }}_{{\Delta }_x}(n|n-1)H^H\\ & \times [H{\hat{\Phi }}_{{\Delta }_x}(n|n-1)H^H+{\Phi }_v(n)]\end{array}\end{array}$$

$$e_x(n)=y(n)-H\hat{\underline{x}}(n|n-1)$$

4. 抑制能力控制

在实际应用中，对音频中不期望的成份（例如混响、环境噪声）的抑制能力进行控制通常是有好处的。在很多情况下，通过控制抑制的量来掩盖音频处理痕迹以及缓和语音失真效果，从而显著提升主观音频听感。在通话场景中，通常倾向于保留一定量的残留噪声，否则的话听者会有一种通话断开的感觉。对于去混响而言，如果早期反射声保留了同时晚期混响完全被抑制，主观听感上就会很不自然。因此考虑到主观感受，通常倾向于保留一定的晚期混响。在本节中，我们来推导一个具有控制噪声和混响抑制能力的输出方案。

算法简化记录