OMLSA算法推导

OM-LSA(Optimally-modified log-spectral amplitude)是经典的降噪算法，这里做一个学习总结，把中间的一些公式推导一下，写一些自己的想法，水平有限，欢迎指正。

LSA估计增益

OM-LSA(Optimally-modified log-spectral amplitude)从名字上就可以看出来估计的是对数幅度谱，LSA做为最小优化目标，即，
e = E { ( l o g ( A ^ ( k , l ) ) − l o g ( A ( k , l ) ) ) 2 } (1) e=E\{(log(\hat{A}(k,l))- log(A(k,l)))^2\} \tag{1} e=E{(log(A^(k,l))−log(A(k,l)))2}(1)

这里使用的其实是贝叶斯MSE准则，通过对$log(\hat{A}(k,l))$求导得到，
$$\begin{gathered} log(\hat{A}(k,l))=E(log(A(k,l)|Y(k,l))) \\ \hat{A}(k,l)=exp(E(logA(k,l)|Y(k,l)))\text{\hspace{1em}(2)} \end{gathered}$$

每一帧可以分为存在语音和不存在语音两种情况，即，
$$\begin{array}{rl}{H}_0(k,l):Y(k,l)& =D(k,l)\\ H_1(k,l):Y(k,l)& =X(k,l)+D(k,l)\end{array}\text{\hspace{1em}(3)}$$

式(2)中需要求解E(log A(k,l)|Y(k,l)),在存在语音和不存在语音两种假设情况下可以写成，
$$\begin{array}{rl}E(logA(k,l)|Y(k,l))& =E(logA(k,l)|Y(k,l),H_1(k,l))p(k,l)\\ & +E(logA(k,l)|Y(k,l),H_0(k,l))(1-p(k,l))\end{array}\text{\hspace{1em}(4)}$$

将式(4)带人(2),并且利用$e^{xy}=(e^x{)}^y,e^{x+y}=e^x{e}^y$的性质，可以得到，
$$\begin{array}{rl}\hat{A}(k,l)& =exp(E(logA(k,l)|Y(k,l),H_1(k,l)){)}^{p(k,l)}\\ & +exp(E(logA(k,l)|Y(k,l),H_0(k,l)){)}^{(1-p(k,l))}\end{array}\text{\hspace{1em}(5)}$$

当语音不存在时，使用一个最小值超参进行约束，即，
$$exp(E(logA(k,l)|Y(k,l),H_0(k,l)))=G_{min}|Y(k,l)|\text{\hspace{1em}(6)}$$

当语音存在时，可以推导出，
$$\begin{gathered} exp(E(logA(k,l)|Y(k,l),H_1(k,l)))=G_{H_1}|Y(k,l)| \\ {G}_{H_1}=\frac{\zeta (k,l)}{1+\zeta (k,l)}{\int }_{v(k,l)}^{\infty }\frac{1}{2}\frac{e^{-t}}{t}dt\text{\hspace{1em}(7)} \end{gathered}$$

其中$p(k,l)$为条件语音存在概率，$\zeta (k,l)$为先验信噪比，$\gamma (k,l)$为后验信噪比，$v(k,l)=\frac{\gamma (k,l)\zeta (k,l)}{1+\zeta (k,l)}$是先验后验信噪比的函数。将式(7),(6)带入式(5)得，
A ^ ( k , l ) = { G H 1 ( k , l ) } p ( k , l ) G m i n 1 − p ( k , l ) ∣ Y ( k , l ) ∣ = G ( k , l ) ∣ Y ( k , l ) ∣ (8) \begin{aligned} \hat{A}(k,l)&= \{G_{H_1}(k,l)\}^{p(k,l)}G_{min}^{1-p(k,l)} |Y(k,l)| \\ &= G(k,l) |Y(k,l)| \end{aligned} \tag{8} A^(k,l)​={GH1​​(k,l)}p(k,l)Gmin1−p(k,l)​∣Y(k,l)∣=G(k,l)∣Y(k,l)∣​(8)

要想得到$G(k,l))$,需要估计出条件语音存在概率和先后验信噪比（间接需要估计出底噪）。

贝叶斯定理估计条件语音存在概率

假设干净语音和噪声的短时傅里叶变换系数满足复高斯分布，且不相干，则可以得到概率密度函数为，
$$\begin{array}{rl}p(Y(k,l)|H_0(k,l))& =\frac{1}{\pi {\lambda }_x(k,l)}exp(-\frac{|Y(k,l){|}^2}{{\lambda }_d(k)})，\\ p(Y(k,l)|H_1(k,l))& =\frac{1}{\pi ({\lambda }_x(k,l)+{\lambda }_d(k,l))}exp(-\frac{|Y(k,l){|}^2}{{\lambda }_x(k)+{\lambda }_d(k)})\end{array}\text{\hspace{1em}(9)}$$

利用贝叶斯定理，
$$\begin{array}{rl}p(H_1(k,l)|Y(k,l))& =\frac{p(Y(k,l)|H_1(k,l))p(H_1)}{p(Y(k,l))}\\ & =\frac{p(Y(k,l)|H_1(k,l))p(H_1)}{p(Y(k,l)|H_1(k,l))p(H_1)+p(Y(k,l)|H_0(k,l))p(H_0)}\end{array}\text{\hspace{1em}(10)}$$

将式(9)带入式(10),得到，
p ( k , l ) = { 1 + q ( k , l ) 1 − q ( k , l ) ( 1 + ζ ( k , l ) ) e x p ( − v ( k , l ) ) } − 1 (11) p(k,l)=\{1+\frac{q(k,l)}{1-q(k,l)}(1+\zeta(k,l))exp(-v(k,l)) \}^{-1} \tag{11} p(k,l)={1+1−q(k,l)q(k,l)​(1+ζ(k,l))exp(−v(k,l))}−1(11)

其中$q(k,l)$为先验语音缺失概率，$\zeta (k,l)=\frac{{\lambda }_x(k,l)}{{\lambda }_d(k,l)}$为先验信噪比，$\gamma (k,l)=\frac{|Y(k,l){|}^2}{{\lambda }_d(k,l)}$为后验信噪比，$v(k,l)=\frac{\gamma (k,l)\zeta (k,l)}{1+\zeta (k,l)}$是先验后验信噪比的函数(和上面的定义一样)。

根据式(11),要想得到条件语音存在概率，需要估计出先后验信噪比（需要估计底噪）和先验语音缺失概率。

DD准则估计先验信噪比

在假设噪声估计出来的情况下，后验信噪比的估计只需要利用当前帧的能量除以噪声能量就可以了，而先验信噪比的估计需要利用DD准则进行估计，即，
ζ ^ ( k , l ) = α G H 1 2 ( k , l − 1 ) γ ( k , l − 1 ) + ( 1 − α ) m a x { γ ( k , l ) − 1 , 0 } (12) \hat{\zeta}(k,l)=\alpha G_{H_1}^2(k,l-1)\gamma(k,l-1) + (1-\alpha)max\{\gamma(k,l)-1,0\} \tag{12} ζ^​(k,l)=αGH1​2​(k,l−1)γ(k,l−1)+(1−α)max{γ(k,l)−1,0}(12)

式(12)的第一部分是上一帧的先验信噪比估计,
G H 1 2 ( k , l − 1 ) γ ( k , l − 1 ) = { G H 1 ( k , l − 1 ) ∣ Y ( k , l − 1 ) ∣ } 2 λ d ( k , l − 1 ) = λ x ( k , l − 1 ) λ d ( k , l − 1 ) = ζ ( k , l − 1 ) \begin{aligned} G_{H_1}^2(k,l-1)\gamma(k,l-1)&=\frac{\{G_{H_1}(k,l-1)|Y(k,l-1)|\}^2}{\lambda_d(k,l-1)}\\ &= \frac{\lambda_x(k,l-1)}{\lambda_d(k,l-1)} \\ &= \zeta(k,l-1) \end{aligned} GH1​2​(k,l−1)γ(k,l−1)​=λd​(k,l−1){GH1​​(k,l−1)∣Y(k,l−1)∣}2​=λd​(k,l−1)λx​(k,l−1)​=ζ(k,l−1)​

式(12)的第二部分是当前帧的先验信噪比估计,
$$\gamma (k,l)-1=\frac{|Y(k,l){|}^2}{{\lambda }_d(k,l)}-1=\frac{|Y(k,l){|}^2-{\lambda }_d(k,l)}{{\lambda }_d(k,l)}=\frac{{\lambda }_x(k,l)}{{\lambda }_d(k,l)}=\zeta (k,l)$$

将两部分做一个平滑得到最终的先验信噪比估计。

通过先验信噪比的软决策估计先验语音缺失概率

首先将估计出来的先验信噪比进行平滑，即，
$$\zeta (k,l)=\beta \zeta (k,l-1)+(1-\beta )\hat{\zeta }(k,l-1)\text{\hspace{1em}(13)}$$

接着在频域进行综合，考虑频带间的影响，定义
$${\zeta }_{\lambda }(k,l)=\sum _{i=-w_{\lambda }}^{w_{\lambda }}{h}_{\lambda }(i)\zeta (k-i,l)$$

其中$h_{\lambda }$为hanning窗，下标$\lambda$可以等于"local",也可以等于"global",当$w_{\lambda }$=1时，为"local"，代表局部信噪比的平均，当$w_{\lambda }$=15时，为"global"，代表全局信噪比的平均，根据${\zeta }_{\lambda }(k,l)$的数值和超参${\zeta }_{min}$,${\zeta }_{max}$定义，
$$p_{\lambda }(k,l)=\{\begin{array}{ll}& 0,\text{if}{\zeta }_{\lambda }(k,l)\le {\zeta }_{min}\\ & 1,{\zeta }_{\lambda }(k,l)\ge {\zeta }_{max}\\ & \frac{log({\zeta }_{\lambda }(k,l)/{\zeta }_{min})}{log({\zeta }_{max}/{\zeta }_{min})},otherwise\end{array}\text{\hspace{1em}(14)}$$

为了进一步在噪声帧加强噪声抑制，定义，
ζ f r a m e ( l ) = m e a n 1 ≤ k ≤ M / 2 + 1 { ζ ( k , l ) } \zeta_{frame}(l)=\underset{1\leq k \leq M/2+1}{mean}\{\zeta(k,l)\} ζframe​(l)=1≤k≤M/2+1mean​{ζ(k,l)}

$$u(l)=\{\begin{array}{ll}& 0，\text{if},{\zeta }_{frame}(l)\le {\zeta }_{peak}{\zeta }_{min}\\ & 1,\text{if},{\zeta }_{frame}(l)\ge {\zeta }_{peak}{\zeta }_{max}\\ & \frac{log({\zeta }_{frame}(l)/{\zeta }_{peak}(l)/{\zeta }_{min}(l))}{log({\zeta }_{max}/{\zeta }_{min})},otherwise\end{array}\text{\hspace{1em}(15)}$$

$p_{frame}$的计算如下图所示，

先验语音缺失概率根据下式求得，
$$\hat{q}(k,l)=1-p_{local}(k,l)p_{global}(k,l)p_{frame}(l)\text{\hspace{1em}(16)}$$

MCRA方法估计底噪

常见的噪声估计方法有递归平滑，最小值跟踪，直方图统计和分位数噪声估计等方法，这里的MCRA(Minima controlled recursive averaging)算法使用了递归平均和最小值跟踪相结合的方法，这里的最小值跟踪体现在语音存在概率是由最小值跟踪确定的。

递归平滑

在噪声段进行语音平滑处理，在语音段不更新噪声，则，
$$\begin{array}{rl}{H}_0^{{}^{\prime }}(k,l):\widehat{{\lambda }_d}(k,l+1)& ={\alpha }_d\widehat{{\lambda }_d}(k,l)+(1-{\alpha }_d)|Y(k,l){|}^2\\ H_1^{{}^{\prime }}(k,l):\widehat{{\lambda }_d}(k,l+1)& =\widehat{{\lambda }_d}(k,l)\end{array}\text{\hspace{1em}(17)}$$

${\alpha }_d$是平滑参数，数值过大会导致跟踪较慢，数值过小，容易产生音乐噪声。利用条件语音存在概率$p^{{}^{\prime }}(k,l)$进行综合，可得，
$$\begin{array}{rl}\widehat{{\lambda }_d}(k,l+1)& =({\alpha }_d\widehat{{\lambda }_d}(k,l)+(1-{\alpha }_d)|Y(k,l){|}^2)(1-p^{{}^{\prime }}(k,l))+\widehat{{\lambda }_d}(k,l)p^{{}^{\prime }}(k,l)\\ & ={\tilde{\alpha }}_d\widehat{{\lambda }_d}(k,l)+(1-{\tilde{\alpha }}_d)|Y(k,l){|}^2\end{array}\text{\hspace{1em}(18)}$$

其中，$\widetilde{{\alpha }_d}={\alpha }_d+(1-{\alpha }_d)p^{{}^{\prime }}(k,l)$

最小值控制语音存在概率

首先对幅度谱在频域进行平滑，得，
$$S_f(k,l)=\sum _{i=-w}^{w}b(i)|Y(k-i,l){|}^2\text{\hspace{1em}(19)}$$

接着在时域进行平滑，
$$S(k,l)={\alpha }_{s}S(k,l-1)+(1-{\alpha }_s)S_f(k,l)\text{\hspace{1em}(20)}$$

利用最小值跟踪法求出$S_{min}(k,l)$,并且定义比值$S_r(k,l)=S(k,l)/S_{min}(k,l)$,如果$S_r(k,l)$大于阈值，则令$I(k,l)=1$,否则为0，条件语音存在概率根据I(k,l)平滑得到，即，
$$p^{{}^{\prime }}(k,l)={\alpha }_p{p}^{{}^{\prime }}(k,l-1)+(1-{\alpha }_{)}I(k,l)\text{\hspace{1em}(21)}$$

将式(21)带入式(18)得到最终的噪声估计值。

总结

算法流程为：首先使用MCRA方法估计出噪声，接着估计出后验信噪比和先验信噪比(DD准则)，利用先验信噪比的软决策估计出先验语音缺失概率，之后利用贝叶斯准则求出条件语音存在概率，将前面求出的值带入增益函数表达式得到增益值。