关于 IMCRA+OMLSA 语音降噪算法的详细解释
关于 IMCRA+OMLSA 语音降噪算法的详细解释
- 概述
- OMLSA 算法
- IMCRA 算法
概述
IMCRA+OMLSA 算法及其一些变体是目前语音降噪中常用的算法。很多文献在解释这两种算法的时候，条理有一些混乱，使得读者难以弄清楚问题的本质。一些文献从 IMCRA 算法开始介绍，我个人觉得这样反而容易引起理解上的混乱。就好比画画，总是要先画出轮廓来，才能逐步细化。围棋也是一样，一定要先布局，才能进入到中盘。我觉得先介绍 OMLSA 算法，才能起到纲举目张的作用。所以本文打算先讲解 OMLSA 算法，再讲解 IMCRA 算法，以给读者一个明确的思路。
OMLSA 算法
OMLSA 算法的英文名称是 optimally modified log-spectral amplitude estimator，文献1中将其翻译为 “最优改进对数谱幅度估计”。从这个名字可以看出，该方法一定是对信号的幅度谱进行了一些操作。OMLSA 算法的计算准则实际上就是最小化实际干净语音和估计出来的干净语音的差异2，这一差异如果使用公式可以表示为：
Delta=E(∣logA(k,l)−logA^(k,l)∣2) {Delta} = E\left( \left| logA( k,l)-log\hat{A}(k,l) \right|^2 \right)\, Delta=E(∣∣∣logA(k,l)−logA^(k,l)∣∣∣2)
其中 A(k,l)A(k,l)A(k,l) 是干净语音信号 xxx 的频谱幅值，A^(k,l)\hat{A}(k,l)A^(k,l) 是估计出来的频谱幅值。如果使用传统的方法，那 A^(k,l)\hat{A}(k,l)A^(k,l) 可以表示为：
A^(k,l)=G(k,l)⋅∣Y(k,l)∣,\hat{A}(k,l) = G(k,l)\cdot |Y(k,l)|, A^(k,l)=G(k,l)⋅∣Y(k,l)∣,
其中 G(k,l)G(k,l)G(k,l) 为带噪语音的估计增益，Y(k,l)Y(k,l)Y(k,l) 就是含噪语音的幅度谱。但是现在OMLSA肯定不会这么简单，它会变化成一种复杂一点的形式。为了避免推导过程中使读者的关注点发生混乱，我们直接给出此时 A^(k,l)\hat{A}(k,l)A^(k,l) 的计算公式
A^(k,l)=(Gmin⋅∣Y(k,l)∣)(1−p(k,l))×(GH1(k,l)⋅∣Y(k,l)∣)p(k,l)\hat{A}(k,l) = (G_{\rm min}\cdot |Y(k,l)|)^{(1-p(k,l))}\times (G_{\rm H1}(k,l)\cdot |Y(k,l)|)^{p(k,l)} A^(k,l)=(Gmin⋅∣Y(k,l)∣)(1−p(k,l))×(GH1(k,l)⋅∣Y(k,l)∣)p(k,l)
可以看出，这个公式里有三个变量：GminG_{\rm min}Gmin，p(k,l)p(k,l)p(k,l) 和 GH1(k,l)G_{\rm H1}(k,l)GH1(k,l)。GminG_{\rm min}Gmin 表示的是语音不存在时的增益函数，一般将 GminG_{\rm min}Gmin 的值设为带噪信号的最小信噪比1，比如假设最小信噪比为 −18dB-18 \rm dB−18dB 时，GminG_{\rm min}Gmin 取值为 0.1257。p(k,l)p(k,l)p(k,l) 是语音存在的后验概率。那么 p(k,l)p(k,l)p(k,l) 和 GH1(k,l)G_{\rm H1}(k,l)GH1(k,l) 又分别怎样求呢？我们这里再次直接给出公式23：
p(k,l)={1+q(k,l)1−q(k,l)(1+ξ(k,l))exp(−γ(k,l)⋅ξ(k,l)1+ξ(k,l))}−1p(k,l) = \left \{ 1+\frac{q(k,l)}{1-q(k,l)} \left( 1+\xi(k,l) \right){\rm exp}(-\frac{\gamma(k,l) \cdot \xi(k,l)}{1+\xi(k,l)}) \right \}^{-1} p(k,l)={1+1−q(k,l)q(k,l)(1+ξ(k,l))exp(−1+ξ(k,l)γ(k,l)⋅ξ(k,l))}−1
GH1(k,l)=ξ(k,l)1+ξ(k,l)exp(12∫v(k,l)∞e−xxdx),withv(k,l)=ξ(k,l)γ(k,l)1+ξ(k,l)G_{\rm H1}(k,l) = \frac{\xi(k,l)}{1+\xi(k,l)}{\rm exp} \left ( \frac{1}{2} \int_{v(k,l)}^\infty \frac{e^{-x}}{x}dx \right), {\rm with}~~v(k,l)=\frac{\xi(k,l)\gamma(k,l)}{1+\xi(k,l)} GH1(k,l)=1+ξ(k,l)ξ(k,l)exp(21∫v(k,l)∞xe−xdx),with v(k,l)=1+ξ(k,l)ξ(k,l)γ(k,l)
擦，这下突然多出来了这么多变量，抓狂。很多文献在这一步的时候，对于这些变量都没有解释清楚，这就导致了理解起来的混乱。现在来一个个梳理。q(k,l)q(k,l)q(k,l) 指的是这一帧的语音不存在的先验概率4，ξ(k,l)\xi(k,l)ξ(k,l) 是一帧语音的先验信噪比，γ(k,l)\gamma(k,l)γ(k,l) 是一帧语音的后验信噪比。也就是说，只要了这三个值，我们就可以求出 p(k,l)p(k,l)p(k,l) 和 GH1(k,l)G_{\rm H1}(k,l)GH1(k,l) 了。我们先来看一下 γ(k,l)\gamma(k,l)γ(k,l)，它的表达式为：
γ(k,l)=∣Y(k,l)∣2λd(k,l)\gamma(k,l)= \frac{|Y(k,l)|^2}{\lambda_{d}(k,l)} γ(k,l)=λd(k,l)∣Y(k,l)∣2
注意到这里出现了一个非常重要的变量 λd(k,l)\lambda_{d}(k,l)λd(k,l)，这个变量是噪声谱估计，它才是融合OMLSA和IMCRA方法的精华所在。后面我们会讲到，它是上一个时刻推导出来的。此时我们先假定它是已知的。而 ξ(k,l)\xi(k,l)ξ(k,l) 的表达式为：
ξ(k,l)=αGH12(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α 是权重因子，用来控制降噪和语音之间的平衡。关于这个因子的选择，有人进行过研究，在这里我们不进行探讨，只需要知道它是一个常数就可以了。 而 GH1(k,l−1)G_{H1}(k,l-1)GH1(k,l−1) 是上一帧的对数谱增益函数，在求解本帧（lll 帧）的时候，它是已经知道的。γ(k,l−1)\gamma(k,l-1)γ(k,l−1) 也是同样的道理。所以现在根本问题就在于求解 γ(k,l)\gamma(k,l)γ(k,l)，而求解 γ(k,l)\gamma(k,l)γ(k,l) 的根本问题就在于求解 λd(k,l)\lambda_{d}(k,l)λd(k,l)。那么在这里就要和IMCRA方法建立一种联系了。
IMCRA 算法
现在继续来讲解如何求解 λd(k,l)\lambda_{d}(k,l)λd(k,l)。实际上，λd(k,l)\lambda_{d}(k,l)λd(k,l) 的求解要用到上一帧的 qqq，也就是 q(k,l−1)q(k,l-1)q(k,l−1)。很多文献里面，关于第 lll 帧和第 l−1l-1l−1 帧的表述非常混乱。为了避免引起表述过程中的混乱，我们假定当前正处于第 lll 帧，即将推导第 l+1l+1l+1 帧的情况。首先，我们要定义一个新的变量 Sf(λ,k)S_{f}(\lambda,k)Sf(λ,k)，它表示的是在频域对每一帧含噪语音谱平滑后的功率谱值5，其表达式为：
Sf(k,l)=∑i=−ωωb(i)∣Y(k−i,l)∣2,S_{f}(k,l)=\sum\limits_{i=-\omega}^{\omega}b (i)|Y(k-i,l)|^2, Sf(k,l)=i=−ω∑ωb(i)∣Y(k−i,l)∣2,
其中 b(i)b(i)b(i) 表示标准化窗函数，窗函数的长度为 2ω+12\omega+12ω+1，Y(k−i,l)Y(k-i,l)Y(k−i,l) 表示含噪语音在时频域作短时傅里叶变换的幅度值。接下来，我们可以通过平滑求出当前帧含噪语音的功率谱：
S(k,l)=αsS(k,l−1)+(1−αs)Sf(k,l),S(k,l)=\alpha_{s}S(k,l-1)+(1-\alpha_{s})S_{f}(k,l), S(k,l)=αsS(k,l−1)+(1−αs)Sf(k,l),
其中 αs\alpha_{s}αs 为平滑参数，我们可以将其设置为 0.8，S(k,l−1)S(k,l-1)S(k,l−1) 表示前一帧含噪语音的功率谱，此时它是已知的（因为是前一帧的信息，之前肯定已经求出来了）。接下来跟踪平滑功率谱 S(k,l)S(k,l)S(k,l) 的最小值：
Smin(k,l)=min{Smin(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)I(k,l)I(k,l) 用于对语音存在进行粗略估计。在此之前，我们先定义两个变量：
γmin(k,l)=∣Y(k,l)∣2BminSmin(k,l),\gamma_{\rm min}(k,l)=\frac{|Y(k,l)|^2}{B_{\min}S_{\min}(k,l)}, γmin(k,l)=BminSmin(k,l)∣Y(k,l)∣2,
ζ(k,l)=S(k,l)BminSmin(k,l),\zeta(k,l)=\frac{S(k,l)}{B_{\min}S_{\min}(k,l)}, ζ(k,l)=BminSmin(k,l)S(k,l),
其中 Bmin=1.66B_{\min}=1.66Bmin=1.66 是噪声谱最小值估计偏置补偿因子。接下来 I(k,l)I(k,l)I(k,l) 的定义为：
I(k,l)={1,ifγmin(k,l)<γ0andζ(k,l)<ζ0(speechisabsent)0,otherwise(speechispresent)I(k,l)=\left\{ \begin{aligned} &1,~~{\rm if}~\gamma_{\min}(k,l)<\gamma_{0}~{\rm and}~\zeta(k,l)<\zeta_0~~~(\rm speech~is ~absent) \\ &0,~~{\rm otherwise}~(\rm speech~is ~present) \end{aligned} \right. I(k,l)={1, if γmin(k,l)<γ0 and ζ(k,l)<ζ0 (speech is absent)0, otherwise (speech is present)
其中 γ0=4.6\gamma_0=4.6γ0=4.6，ζ0=1.67\zeta_0=1.67ζ0=1.67。接下来，就要用上这个 I(k,l)I(k,l)I(k,l)。对功率谱频点间进行二次平滑得到平滑功率谱值：
S~f(k,l)={∑i=−ωωb(i)I(k−i,l)∣Y(k−i,l)∣2∑i=−ωωb(i)I(k−i,l),if∑i=−ωωI(k−i,l)≠0S~(k,l−1),otherwise\tilde{S}_f(k,l)=\left\{ \begin{aligned} &\frac{\sum\limits_{i=-\omega}^{\omega}b(i)I(k-i,l)|Y(k-i,l)|^2}{\sum\limits_{i=-\omega}^{\omega}b(i)I(k-i,l)}, ~~{\rm if}~\sum\limits_{i=-\omega}^{\omega}I(k-i,l)\neq 0\\ &\tilde{S}(k,l-1), {\rm otherwise} \end{aligned} \right. S~f(k,l)=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧i=−ω∑ωb(i)I(k−i,l)i=−ω∑ωb(i)I(k−i,l)∣Y(k−i,l)∣2, if i=−ω∑ωI(k−i,l)̸=0S~(k,l−1),otherwise
其中 S~(k,l)\tilde{S}(k,l)S~(k,l) 表示二次平滑时，在时域采用一阶递归平滑方法得到的平滑功率谱值，其表达式为：
S~(k,l)=αsS~(k,l−1)+(1−αs)S~f(k,l)\tilde{S}(k,l)=\alpha_{s}\tilde{S}(k,l-1)+(1-\alpha_{s})\tilde{S}_{f}(k,l) S~(k,l)=αsS~(k,l−1)+(1−αs)S~f(k,l)
跟踪平滑功率谱 S~min(k,l)\tilde{S}_{\min}(k,l)S~min(k,l) 的最小值：
S~min(k,l)=min{Smin(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)}
在此之后，我们还要定义两个变量：
γ~min(k,l)=∣Y(k,l)∣2BminS~min(k,l)\tilde{\gamma}_{\min}(k,l)=\frac{|Y(k,l)|^2}{B_{\min}\tilde{S}_{\min}(k,l)} γ~min(k,l)=BminS~min(k,l)∣Y(k,l)∣2
ζ~(k,l)=S(k,l)BminS~min(k,l)\tilde{\zeta}(k,l)=\frac{S(k,l)}{B_{\min}\tilde{S}_{\min}(k,l)} ζ~(k,l)=BminS~min(k,l)S(k,l)
接下来，终于到激动人心的一步，可以求先验语音不存在的概率 q(k,l)q(k,l)q(k,l) 了。它的计算公式如下：
q~(k,l)={1,ifγ~min(k,l)≤1andζ~(k,l)<ζ~0γ1−γ~min(k,l)γ1−1,if1<γ~min(k,l)<γ1andζ~(k,l)<ζ00,otherwise\tilde{q}(k,l)=\left\{ \begin{aligned} &1,~~{\rm if}~\tilde{\gamma}_{\min}(k,l)\leq 1~{\rm and}~\tilde{\zeta}(k,l)<\tilde{\zeta}_0 \\ &\frac{\gamma_{1}-\tilde{\gamma}_{\min}(k,l)}{\gamma_{1}-1},~~{\rm if}~1<\tilde{\gamma}_{\min}(k,l)<\gamma_{1}~{\rm and}~\tilde{\zeta}(k,l)<\zeta_{0}\\ &0, {\rm otherwise} \end{aligned} \right. q~(k,l)=⎩⎪⎪⎪⎨⎪⎪⎪⎧1, if γ~min(k,l)≤1 and ζ~(k,l)<ζ~0γ1−1γ1−γ~min(k,l), if 1<γ~min(k,l)<γ1 and ζ~(k,l)<ζ00,otherwise
其中 ζ0=1.67\zeta_{0}=1.67ζ0=1.67，γ1=3\gamma_{1}=3γ1=3。截止到目前为止，q(k,l)q(k,l)q(k,l)，ξ(k,l)\xi(k,l)ξ(k,l) 和 γ(k,l)\gamma(k,l)γ(k,l) 这三个变量都已经求出来了，那么对于当前时刻的降噪来说已经足够了。但是为了可持续发展，我们还需要递推下一个时刻 l+1l+1l+1 的噪声谱估计 λd(k,l+1)\lambda_{d}(k,l+1)λd(k,l+1)。首先使用 q(k,l)q(k,l)q(k,l) 来求 p(k,l)p(k,l)p(k,l)，这个公式前面已经说过：
p(k,l)={1+q(k,l)1−q(k,l)(1+ξ(k,l))exp(−γ(k,l)⋅ξ(k,l)1+ξ(k,l))}−1p(k,l) = \left \{ 1+\frac{q(k,l)}{1-q(k,l)} \left( 1+\xi(k,l) \right){\rm exp}(-\frac{\gamma(k,l) \cdot \xi(k,l)}{1+\xi(k,l)}) \right \}^{-1} p(k,l)={1+1−q(k,l)q(k,l)(1+ξ(k,l))exp(−1+ξ(k,l)γ(k,l)⋅ξ(k,l))}−1
然后定义语音存在条件概率计算时变平滑参数 α~d(k,l)\tilde{\alpha}_{d}(k,l)α~d(k,l)，公式如下：
α~d(k,l)=αd+(1−αd)p(k,l)\tilde{\alpha}_{d}(k,l)=\alpha_{d}+(1-\alpha_{d})p(k,l) α~d(k,l)=αd+(1−αd)p(k,l)
其中 αd=0.85\alpha_{d}=0.85αd=0.85。根据这里计算出来的 α~d(k,l)\tilde{\alpha}_{d}(k,l)α~d(k,l)，就可以求解下一个时刻的 λˉd(k,l+1)\bar{\lambda}_{d}(k,l+1)λˉd(k,l+1)
λˉd(k,l+1)=α~d(k,l)λˉd(k,l)+[1−α~d(k,l)]∣Y(k,l)∣2\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 λˉd(k,l+1)=α~d(k,l)λˉd(k,l)+[1−α~d(k,l)]∣Y(k,l)∣2
为了避免噪声谱估计的过低，需要采用一个偏置因子对噪声谱估计进行补偿，公式如下：
λ^d(k,l+1)=β⋅λˉd(k,l+1)\hat{\lambda}_{d}(k,l+1)=\beta\cdot \bar{\lambda}_{d}(k,l+1) λ^d(k,l+1)=β⋅λˉd(k,l+1)
到这里，就已经求出来了下一个周期的 λ^d(k,l+1)\hat{\lambda}_{d}(k,l+1)λ^d(k,l+1)，就可以进入下一个周期的计算。
本文的参考文件仅仅列举出来了一些中文的文献。关于 IMCRA+OMLSA 算法的英文经典文献也有很多，在此不一一列举，大家可以在网上查找。
OM-LSA和小波阈值去噪结合的语音增强，刘凤增，李国辉，李博，计算机科学与探索，2011, 5(6), 541-552. ↩︎ ↩︎
针对于家居环境的语音增强系统的研究与开发，陈成斌，华南理工大学工程硕士学位论文，2014. ↩︎ ↩︎
特定人语音增强算法的研究，郭栗，上海交通大学硕士学位论文，2015. ↩︎
基于改进谱平滑策略的IMCRA算法及其语音增强，张建伟，陶亮，周健，王华彬，计算机工程与应用，2017, 53(1): 153-157. ↩︎
基于维纳过滤的IMCRA算法，吴进，赵隽，李乔深，西安邮电大学学报，2017, 22(5): 73-77. ↩︎
关于 IMCRA+OMLSA 语音降噪算法的详细解释相关推荐
- 传统语音增强——基于先验信噪比的维纳滤波语音降噪算法
一.基于先验信噪比的维纳滤波语音降噪算法的基本概念 改进的维纳滤波器为基于先验信噪比的维纳滤波器,其原理框图下图所示. 对于第m帧带噪语音信号ym(n)=sm(n)+nm(n) 式中,sm(n)是第m ...
- 传统语音增强——基于小波分解的语音降噪算法
一.小波分析的意义 在传统的傅里叶分析中,信号完全是在频域展开的,不包含任何时频的信息.因为丢弃的时域信息对某些应用同样重要,所以出现很多能表征时域和频域信息的信号分析方法,如短时傅里叶变换.Gabo ...
- 传统语音增强——基本的维纳滤波语音降噪算法
一.维纳滤波的基本原理 基本维纳滤波就是用来解决从噪声中提取信号问题的一种过滤(或滤波)方法.它基于平稳随机过程模型,且假设退化模型为线性空间不变系统的.实际上这种线性滤波问题,可以看成是一种估计问题 ...
- [投稿]一个频域语音降噪算法实现及改进方法
本文是音频处理的朋友icoolmedia(QQ:314138065)的投稿.对音频处理有兴趣的朋友可以通过下面的方式与他交流: 作者:icoolmedia QQ:314138065 音视频算法讨论 ...
- KMP算法的详细解释及实现
这是我自己学习算法时有关KMP的学习笔记,代码注释的十分的详细,分享给大家,希望对大家有所帮助 在介绍KMP算法之前, 先来介绍一下朴素模式匹配算法: 朴素模式匹配算法: 假设要从主串S=" ...
- Gauss-Newton算法代码详细解释(转载+自己注释)
这篇博客是对[1]中不详细的地方进行细节上的阐述, 并且每句代码都加了注释,使得更加容易理解 下面的论述(包括伪代码和算法)特指被最小化的目标函数是MSE的时候 需要注意,如果不是MSE为目标函数,那 ...
- 反向传播算法的详细解释（下）
上一篇文章算是用"链式法则"给我们开了个头,下一篇文章则是将反向传播算法应用到神经网络. 原文出处:知乎 https://zhuanlan.zhihu.com/p/25081671 ...
- SORT跟踪算法的详细解释，不容错过
转载自:https://blog.csdn.net/HaoBBNuanMM/article/details/85555547 SORT - SIMPLE ONLINE AND REALTIME TRA ...
- 狄克斯特拉（Dijkstra）算法原理详细解释与实现(python)
目录 写在前面 1. 简介 2. 原理 2.1 找出最便宜的节点 2.2 计算前往该节点的各个邻居的开销 2.3 重复上面的步骤 实现 总结 写在前面 本文原理摘自<算法图解>这本书. 其 ...
最新文章
- nsTimer的简单用法
- [置顶] 2014年八大最热门IT技能
- asp ed什么意思 j_这部洗脑ED动画是如何创作出来的？
- 绍兴市一男子醉酒驾车还冲上公交车暴打司机
- 第五章：创建自定义绑定
- mysql-表记录之增删改操作
- java g1的并行_Java 11好用吗
- vi 之行号操作---显示行号、跳到指定行
- HTTP、Asp.net管道与IIS
- ArcGIS操作系列（一）之地理配准
- 计算机毕业设计（附源码）python学科竞赛赛场安排系统
- codevs2069 油画 — 动态维护优先队列
- 解决黑苹果和Windows双系统时，时钟不同步的问题
- SP4354 TWINSNOW - Snowflakes
- AM335x DDR3配置
- 群翔ShopNum1分销系统V8.1升级版，更优更全更盈利
- 【技术美术图形部分】图形渲染管线2.0-GPU管线概述几何阶段
- kafka维护工具使用指南
- CVE-2017-8464震网3.0漏洞分析与复现
- vue element-ui 中走马灯自适应图片高度