关于 IMCRA+OMLSA 语音降噪算法的详细解释

概述

IMCRA+OMLSA 算法及其一些变体是目前语音降噪中常用的算法。很多文献在解释这两种算法的时候，条理有一些混乱，使得读者难以弄清楚问题的本质。一些文献从 IMCRA 算法开始介绍，我个人觉得这样反而容易引起理解上的混乱。就好比画画，总是要先画出轮廓来，才能逐步细化。围棋也是一样，一定要先布局，才能进入到中盘。我觉得先介绍 OMLSA 算法，才能起到纲举目张的作用。所以本文打算先讲解 OMLSA 算法，再讲解 IMCRA 算法，以给读者一个明确的思路。

OMLSA 算法

OMLSA 算法的英文名称是 optimally modified log-spectral amplitude estimator，文献¹中将其翻译为 “最优改进对数谱幅度估计”。从这个名字可以看出，该方法一定是对信号的幅度谱进行了一些操作。OMLSA 算法的计算准则实际上就是最小化实际干净语音和估计出来的干净语音的差异²，这一差异如果使用公式可以表示为：

$$Delta=E\left({\left|logA(k,l)-log\hat{A}(k,l)\right|}^2\right)$$

其中 $A(k,l)$ 是干净语音信号 $x$ 的频谱幅值，$\hat{A}(k,l)$ 是估计出来的频谱幅值。如果使用传统的方法，那 $\hat{A}(k,l)$ 可以表示为：

$$\hat{A}(k,l)=G(k,l)\cdot |Y(k,l)|,$$

其中 $G(k,l)$ 为带噪语音的估计增益，$Y(k,l)$ 就是含噪语音的幅度谱。但是现在OMLSA肯定不会这么简单，它会变化成一种复杂一点的形式。为了避免推导过程中使读者的关注点发生混乱，我们直接给出此时 $\hat{A}(k,l)$ 的计算公式

$$\hat{A}(k,l)=(G_{\mathrm{m}\mathrm{i}\mathrm{n}}\cdot |Y(k,l)|{)}^{(1-p(k,l))}\times (G_{\mathrm{H}1}(k,l)\cdot |Y(k,l)|{)}^{p(k,l)}$$

可以看出，这个公式里有三个变量：$G_{\mathrm{m}\mathrm{i}\mathrm{n}}$，$p(k,l)$ 和 $G_{\mathrm{H}1}(k,l)$。$G_{\mathrm{m}\mathrm{i}\mathrm{n}}$ 表示的是语音不存在时的增益函数，一般将 $G_{\mathrm{m}\mathrm{i}\mathrm{n}}$ 的值设为带噪信号的最小信噪比¹，比如假设最小信噪比为 $-18\mathrm{d}\mathrm{B}$ 时，$G_{\mathrm{m}\mathrm{i}\mathrm{n}}$ 取值为 0.1257。$p(k,l)$ 是语音存在的后验概率。那么 $p(k,l)$ 和 $G_{\mathrm{H}1}(k,l)$ 又分别怎样求呢？我们这里再次直接给出公式²³：

$$p(k,l)={\left\{1+\frac{q(k,l)}{1-q(k,l)}\left(1+\xi (k,l)\right)\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{\gamma (k,l)\cdot \xi (k,l)}{1+\xi (k,l)})\right\}}^{-1}$$

$$G_{\mathrm{H}1}(k,l)=\frac{\xi (k,l)}{1+\xi (k,l)}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{1}{2}{\int }_{v(k,l)}^{\infty }\frac{e^{-x}}{x}dx\right),\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}v(k,l)=\frac{\xi (k,l)\gamma (k,l)}{1+\xi (k,l)}$$

擦，这下突然多出来了这么多变量，抓狂。很多文献在这一步的时候，对于这些变量都没有解释清楚，这就导致了理解起来的混乱。现在来一个个梳理。$q(k,l)$ 指的是这一帧的语音不存在的先验概率⁴，$\xi (k,l)$ 是一帧语音的先验信噪比，$\gamma (k,l)$ 是一帧语音的后验信噪比。也就是说，只要了这三个值，我们就可以求出 $p(k,l)$ 和 $G_{\mathrm{H}1}(k,l)$ 了。我们先来看一下 $\gamma (k,l)$，它的表达式为：

$$\gamma (k,l)=\frac{|Y(k,l){|}^2}{{\lambda }_d(k,l)}$$

注意到这里出现了一个非常重要的变量 ${\lambda }_d(k,l)$，这个变量是噪声谱估计，它才是融合OMLSA和IMCRA方法的精华所在。后面我们会讲到，它是上一个时刻推导出来的。此时我们先假定它是已知的。而 $\xi (k,l)$ 的表达式为：

ξ ( k , l ) = α G H 1 2 ( k , l − 1 ) γ ( k , l − 1 ) + ( 1 − α ) max ⁡ { γ ( k , l ) − 1 , 0 } , \xi(k,l)= \alpha G_{H1}^{2}(k,l-1)\gamma(k,l-1)+(1-\alpha)\max \{ \gamma(k,l)-1,0 \}, ξ(k,l)=αGH12​(k,l−1)γ(k,l−1)+(1−α)max{γ(k,l)−1,0},

其中 $\alpha$ 是权重因子，用来控制降噪和语音之间的平衡。关于这个因子的选择，有人进行过研究，在这里我们不进行探讨，只需要知道它是一个常数就可以了。 而 $G_{H1}(k,l-1)$ 是上一帧的对数谱增益函数，在求解本帧（$l$ 帧）的时候，它是已经知道的。$\gamma (k,l-1)$ 也是同样的道理。所以现在根本问题就在于求解 $\gamma (k,l)$，而求解 $\gamma (k,l)$ 的根本问题就在于求解 ${\lambda }_d(k,l)$。那么在这里就要和IMCRA方法建立一种联系了。

IMCRA 算法

现在继续来讲解如何求解 ${\lambda }_d(k,l)$。实际上，${\lambda }_d(k,l)$ 的求解要用到上一帧的 $q$，也就是 $q(k,l-1)$。很多文献里面，关于第 $l$ 帧和第 $l-1$ 帧的表述非常混乱。为了避免引起表述过程中的混乱，我们假定当前正处于第 $l$ 帧，即将推导第 $l+1$ 帧的情况。首先，我们要定义一个新的变量 $S_f(\lambda ,k)$，它表示的是在频域对每一帧含噪语音谱平滑后的功率谱值⁵，其表达式为：

$$S_f(k,l)=\sum _{i=-\omega }^{\omega }b(i)|Y(k-i,l){|}^2,$$

其中 $b(i)$ 表示标准化窗函数，窗函数的长度为 $2\omega +1$，$Y(k-i,l)$ 表示含噪语音在时频域作短时傅里叶变换的幅度值。接下来，我们可以通过平滑求出当前帧含噪语音的功率谱：

$$S(k,l)={\alpha }_{s}S(k,l-1)+(1-{\alpha }_s)S_f(k,l),$$

其中 ${\alpha }_s$ 为平滑参数，我们可以将其设置为 0.8，$S(k,l-1)$ 表示前一帧含噪语音的功率谱，此时它是已知的（因为是前一帧的信息，之前肯定已经求出来了）。接下来跟踪平滑功率谱 $S(k,l)$ 的最小值：

S m i n ( k , l ) = min ⁡ { S m i n ( k , l − 1 ) , S ( k , l ) } . S_{\rm min}(k,l)=\min \{ S_{\rm min}(k,l-1),S(k,l) \}. Smin​(k,l)=min{Smin​(k,l−1),S(k,l)}.

接下来，我们需要计算标识函数 $I(k,l)$ 用于对语音存在进行粗略估计。在此之前，我们先定义两个变量：

$${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}(k,l)=\frac{|Y(k,l){|}^2}{B_{\min }{S}_{\min }(k,l)},$$

$$\zeta (k,l)=\frac{S(k,l)}{B_{\min }{S}_{\min }(k,l)},$$

其中 $B_{\min }=1.66$ 是噪声谱最小值估计偏置补偿因子。接下来 $I(k,l)$ 的定义为：

$$I(k,l)=\{\begin{array}{rl}& 1,\mathrm{i}\mathrm{f}{\gamma }_{\min }(k,l)<{\gamma }_0\mathrm{a}\mathrm{n}\mathrm{d}\zeta (k,l)<{\zeta }_0(\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{e}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t})\\ & 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}(\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{e}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t})\end{array}$$

其中 ${\gamma }_0=4.6$，${\zeta }_0=1.67$。接下来，就要用上这个 $I(k,l)$。对功率谱频点间进行二次平滑得到平滑功率谱值：

$${\tilde{S}}_f(k,l)=\{\begin{array}{rl}& \frac{\sum _{i=-\omega }^{\omega }b(i)I(k-i,l)|Y(k-i,l){|}^2}{\sum _{i=-\omega }^{\omega }b(i)I(k-i,l)},\mathrm{i}\mathrm{f}\sum _{i=-\omega }^{\omega }I(k-i,l)̸=0\\ & \tilde{S}(k,l-1),\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$$

其中 $\tilde{S}(k,l)$ 表示二次平滑时，在时域采用一阶递归平滑方法得到的平滑功率谱值，其表达式为：

$$\tilde{S}(k,l)={\alpha }_s\tilde{S}(k,l-1)+(1-{\alpha }_s){\tilde{S}}_f(k,l)$$

跟踪平滑功率谱 ${\tilde{S}}_{\min }(k,l)$ 的最小值：

S ~ min ⁡ ( k , l ) = min ⁡ { S min ⁡ ( k , l − 1 ) , S ~ ( k , l ) } \tilde{S}_{\min}(k,l)=\min\{ S_{\min}(k,l-1),\tilde{S}(k,l)\} S~min​(k,l)=min{Smin​(k,l−1),S~(k,l)}

在此之后，我们还要定义两个变量：

$${\tilde{\gamma }}_{\min }(k,l)=\frac{|Y(k,l){|}^2}{B_{\min }{\tilde{S}}_{\min }(k,l)}$$

$$\tilde{\zeta }(k,l)=\frac{S(k,l)}{B_{\min }{\tilde{S}}_{\min }(k,l)}$$

接下来，终于到激动人心的一步，可以求先验语音不存在的概率 $q(k,l)$ 了。它的计算公式如下：

$$\tilde{q}(k,l)=\{\begin{array}{rl}& 1,\mathrm{i}\mathrm{f}{\tilde{\gamma }}_{\min }(k,l)\le 1\mathrm{a}\mathrm{n}\mathrm{d}\tilde{\zeta }(k,l)<{\tilde{\zeta }}_0\\ & \frac{{\gamma }_1-{\tilde{\gamma }}_{\min }(k,l)}{{\gamma }_1-1},\mathrm{i}\mathrm{f}1<{\tilde{\gamma }}_{\min }(k,l)<{\gamma }_1\mathrm{a}\mathrm{n}\mathrm{d}\tilde{\zeta }(k,l)<{\zeta }_0\\ & 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$$

其中 ${\zeta }_0=1.67$，${\gamma }_1=3$。截止到目前为止，$q(k,l)$，$\xi (k,l)$ 和 $\gamma (k,l)$ 这三个变量都已经求出来了，那么对于当前时刻的降噪来说已经足够了。但是为了可持续发展，我们还需要递推下一个时刻 $l+1$ 的噪声谱估计 ${\lambda }_d(k,l+1)$。首先使用 $q(k,l)$ 来求 $p(k,l)$，这个公式前面已经说过：

$$p(k,l)={\left\{1+\frac{q(k,l)}{1-q(k,l)}\left(1+\xi (k,l)\right)\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{\gamma (k,l)\cdot \xi (k,l)}{1+\xi (k,l)})\right\}}^{-1}$$

然后定义语音存在条件概率计算时变平滑参数 ${\tilde{\alpha }}_d(k,l)$，公式如下：

$${\tilde{\alpha }}_d(k,l)={\alpha }_d+(1-{\alpha }_d)p(k,l)$$

其中 ${\alpha }_d=0.85$。根据这里计算出来的 ${\tilde{\alpha }}_d(k,l)$，就可以求解下一个时刻的 ${\bar{\lambda }}_d(k,l+1)$

$${\bar{\lambda }}_d(k,l+1)={\tilde{\alpha }}_d(k,l){\bar{\lambda }}_d(k,l)+[1-{\tilde{\alpha }}_d(k,l)]|Y(k,l){|}^2$$

为了避免噪声谱估计的过低，需要采用一个偏置因子对噪声谱估计进行补偿，公式如下：

$${\hat{\lambda }}_d(k,l+1)=\beta \cdot {\bar{\lambda }}_d(k,l+1)$$

到这里，就已经求出来了下一个周期的 ${\hat{\lambda }}_d(k,l+1)$，就可以进入下一个周期的计算。

本文的参考文件仅仅列举出来了一些中文的文献。关于 IMCRA+OMLSA 算法的英文经典文献也有很多，在此不一一列举，大家可以在网上查找。

---

1. OM-LSA和小波阈值去噪结合的语音增强，刘凤增，李国辉，李博，计算机科学与探索，2011, 5(6), 541-552. ↩︎ ↩︎

2. 针对于家居环境的语音增强系统的研究与开发，陈成斌，华南理工大学工程硕士学位论文，2014. ↩︎ ↩︎

3. 特定人语音增强算法的研究，郭栗，上海交通大学硕士学位论文，2015. ↩︎

4. 基于改进谱平滑策略的IMCRA算法及其语音增强，张建伟，陶亮，周健，王华彬，计算机工程与应用，2017, 53(1): 153-157. ↩︎

5. 基于维纳过滤的IMCRA算法，吴进，赵隽，李乔深，西安邮电大学学报，2017, 22(5): 73-77. ↩︎