OMLSA&IMCRA学习笔记

前言-MMSE-LSA

这两个搭配好像这个是经典算法的集大成者。首先回忆一下MMSE-LSA的经典公式如下
$$\begin{array}{rl}\widehat{A_k}& =\frac{{\xi }_k}{1+{\xi }_k}exp\left\{\frac{1}{2}{\int }_{v_k}^{\infty }\frac{e^{-t}}{t}\mathrm{d}t\right\}{R}_k\end{array}$$

OMLSA

而所谓的最优改进（Optimal Modified）其实增加了话音不确定性分析，先果后因的看一下OMLSA的公式：
G ( k , l ) = { G H 1 ( k , l ) } p ( k , l ) × G m i n 1 − p ( k , l ) G(k,l)=\{G_{H_1}(k,l)\}^{p(k,l)}\times G_{min}^{1-p(k,l)} G(k,l)={GH1​​(k,l)}p(k,l)×Gmin1−p(k,l)​

OMLSA公式推导

上式中的$G_{H_1}(k,l)=\frac{\xi (k,l)}{1+\xi (k,l)}exp\left\{\frac{1}{2}{\int }_{v(k,l)}^{\infty }\frac{e^{-t}}{t}\mathrm{d}t\right\}$，很显然这个增益的算法比LSA复杂多了，而语音出现概率$p(k,l)$变成了指数运算。这一结果的来源是cohen提出了$H_0(k,l)$和$H_1(k,l)$作为语音缺席和语音出现两种假设。另外假设STFT的bin是复高斯（多元？）分布，所以被观测信号的条件概率密度被定义如下：
$$\begin{array}{rl}p(Y(k,l)|H_0(k,l))& =\frac{1}{\pi {\lambda }_d(k,l)}exp\{-\frac{|Y(k,l){|}^2}{{\lambda }_d(k,l)}\}\\ p(Y(k,l)|H_1(k,l))& =\frac{1}{\pi ({\lambda }_x(k,l)+\pi ({\lambda }_d(k,l))}exp\{-\frac{|Y(k,l){|}^2}{\pi ({\lambda }_x(k,l)+{\lambda }_d(k,l)}\}\end{array}$$
$k$表示频域的分量，$l$表示帧的编号。${\lambda }_x(k,l)=\mathbf{E}[|X(k,l){|}^2|H_1(k,l)]$，${\lambda }_d(k,l)=\mathbf{E}[|D(k,l){|}^2]$在此基础上，利用bayes规则，“以条件语音出现概率作为每一帧的语音出现概率”，这句话很重要，它表明每一帧的语音出现概率还是有据可循不是瞎猜的。定义的公式如下：
$$p(k,l)\triangleq P(H_1(k,l)|Y(k,l))$$
基于统计独立的假设，cohen定义了LSA估计器如下的形式：
$$\hat{A}(k,l)=exp\{\mathbf{E}[logA(k,l)|Y(k,l)]\}\triangleq G(k,l)|Y(k,l)|$$
在没有话音的时候，增益被限定到一个$G_{min}$，得到
$$exp\{\mathbf{E}[logA(k,l)|Y(k,l),H_0(k,l)](1-p(k,l))\}=G_{min}|Y(k,l)|$$
而在话音出现的时候，利用LSA的增益公式得到
$$exp\{\mathbf{E}[logA(k,l)|Y(k,l),H_1(k,l)](p(k,l))\}=G_{H_1}|Y(k,l)|$$
利用bayes概率公式：
$$\begin{array}{rl}\mathbf{E}[logA(k,l)|Y(k,l)]& =\mathbf{E}[logA(k,l)|Y(k,l),H_1(k,l)]p(k,l)\\ & +\mathbf{E}[logA(k,l)|Y(k,l),H_0(k,l)](1-p(k,l))\\ & =log(G_{H_1}|Y(k,l)|)p(k,l)+log(G_{min}|Y(k,l)|)(1-p(k,l))\\ & =log((G_{H_1}|Y(k,l)|{)}^{p(k,l)})+log((G_{min}|Y(k,l)|{)}^{(1-p(k,l))})\end{array}$$
那么最终的幅度谱估计公式如下。
$$\begin{array}{rl}\hat{A}(k,l)& =exp\{\mathbf{E}[logA(k,l)|Y(k,l)]\}\\ & =exp\{log((G_{H_1}|Y(k,l)|{)}^{p(k,l)})+log((G_{min}|Y(k,l)|{)}^{(1-p(k,l))})\}\\ & \triangleq G(k,l)|Y(k,l)|\end{array}$$
上式利用了log和exp的技巧，最终得出了新的公式。新公式相对于老公式多了一个$p(k,l)$语音存在概率，那么这个语音存在概率如何获得呢？

语音存在概率

语音出现概率是基于似然比的统计方法得出的，定义这两个假设的被观测信号频域表达：
$$\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}$$
令${\sigma }_d^2(k,l)=E[|D(k,l){|}^2]$表示噪声的方差，那么在最小均方误差意义上最优的噪声功率密度${\sigma }_d^2(k,l)$的估计表示为：

σ ^ d 2 ( k , l ) = E { σ d 2 ( k , l ) ∣ Y ( k , l ) } = E { σ d 2 ( k , l ) ∣ H 0 ( k , l ) } P ( H 0 ( k , l ) ∣ Y ( k , l ) ) + E { σ d 2 ( k , l ) ∣ H 1 ( k , l ) } P ( H 1 ( k , l ) ∣ Y ( k , l ) ) \begin{aligned} \hat\sigma^2_d(k,l)&=E\{\sigma^2_d(k,l)|Y(k,l)\}\\ &=E\{\sigma^2_d(k,l)|H_0(k,l)\}P(H_0(k,l)|Y(k,l)) + E\{\sigma^2_d(k,l)|H_1(k,l)\}P(H_1(k,l)|Y(k,l)) \end{aligned} σ^d2​(k,l)​=E{σd2​(k,l)∣Y(k,l)}=E{σd2​(k,l)∣H0​(k,l)}P(H0​(k,l)∣Y(k,l))+E{σd2​(k,l)∣H1​(k,l)}P(H1​(k,l)∣Y(k,l))​
其中$P(H_0(k,l)|Y(k,l))$利用贝叶斯条件概率公式可以写成如下形式：
$$\begin{array}{rl}P(H_0(k,l)|Y(k,l))& =\frac{P(Y(k,l)|H_0(k,l))P(H_0(k,l))}{P(Y(k,l))}\\ & =\frac{P(Y(k,l)|H_0(k,l))P(H_0(k,l))}{P(Y(k,l)|H_0(k,l))P(H_0(k,l))+P(Y(k,l)|H_1(k,l))P(H_1(k,l))}\\ & =\frac{1}{1+r\Lambda (k,l)}\end{array}$$
这里
$$r\triangleq \frac{P(H_1(k,l))}{P(H_0(k,l))},\Lambda (k,l)\triangleq \frac{P(Y(k,l)|H_1(k,l))}{P(Y(k,l)|H_0(k,l))}$$
同理可以求得：
$$P(H_1(k,l)|Y(k,l))=\frac{r\Lambda (k,l)}{1+r\Lambda (k,l)}$$

带回原来的公式
σ ^ d 2 ( k , l ) = E { σ d 2 ( k , l ) ∣ Y ( k , l ) } = 1 1 + r Λ ( k , l ) E { σ d 2 ( k , l ) ∣ H 0 ( k , l ) } + r Λ ( k , l ) 1 + r Λ ( k , l ) E { σ d 2 ( k , l ) ∣ H 1 ( k , l ) } ≈ 1 1 + r Λ ( k , l ) ∣ Y ( k , l ) ∣ 2 + r Λ ( k , l ) 1 + r Λ ( k , l ) σ d 2 ( k , l − 1 ) \begin{aligned} \hat\sigma^2_d(k,l)&=E\{\sigma^2_d(k,l)|Y(k,l)\}\\ &=\frac{1}{1+r\Lambda(k,l)}E\{\sigma^2_d(k,l)|H_0(k,l)\} + \frac{r\Lambda(k,l)}{1+r\Lambda(k,l)}E\{\sigma^2_d(k,l)|H_1(k,l)\}\\ &\approx\frac{1}{1+r\Lambda(k,l)}|Y(k,l)|^2 + \frac{r\Lambda(k,l)}{1+r\Lambda(k,l)}\sigma^2_d(k,l-1) \end{aligned} σ^d2​(k,l)​=E{σd2​(k,l)∣Y(k,l)}=1+rΛ(k,l)1​E{σd2​(k,l)∣H0​(k,l)}+1+rΛ(k,l)rΛ(k,l)​E{σd2​(k,l)∣H1​(k,l)}≈1+rΛ(k,l)1​∣Y(k,l)∣2+1+rΛ(k,l)rΛ(k,l)​σd2​(k,l−1)​
根据观察得出上式最后一项，这样就得出了时间平滑的噪声估计公式。这是语音出现概率
$$p(k,l)=\frac{r\Lambda (k,l)}{1+r\Lambda (k,l)}$$
假定噪声的频域估计是一个零均值方差为${\lambda }_d(k)$的复高斯分布，可以得到概率密度函数
p ( Y ( k , l ) ∣ H 0 ( k , l ) ) = 1 π λ d ( k ) e x p { − Y 2 ( k , l ) λ d ( k ) } p(Y(k,l)|H_0(k,l)) = \frac{1}{\pi \lambda_d(k)}exp\{-\frac{Y^2(k,l)}{\lambda_d(k)}\} p(Y(k,l)∣H0​(k,l))=πλd​(k)1​exp{−λd​(k)Y2(k,l)​}
同理假设噪声和语音是不相关的零均值方差为${\lambda }_x(k)+{\lambda }_d(k)$复高斯分布，可以得到概率密度函数
p ( Y ( k , l ) ∣ H 1 ( k , l ) ) = 1 π [ λ x ( k ) + λ d ( k ) ] e x p { − Y 2 ( k , l ) λ x ( k ) + λ d ( k ) } p(Y(k,l)|H_1(k,l)) = \frac{1}{\pi [\lambda_x(k)+\lambda_d(k)]}exp\{-\frac{Y^2(k,l)}{\lambda_x(k)+\lambda_d(k)}\} p(Y(k,l)∣H1​(k,l))=π[λx​(k)+λd​(k)]1​exp{−λx​(k)+λd​(k)Y2(k,l)​}
Λ ( k , l ) ≜ P ( Y ( k , l ) ∣ H 1 ( k , l ) ) P ( Y ( k , l ) ∣ H 0 ( k , l ) ) = p ( Y ( k , l ) ∣ H 1 ( k , l ) ) p ( Y ( k , l ) ∣ H 0 ( k , l ) ) = 1 π [ λ x ( k ) + λ d ( k ) ] e x p { − Y 2 ( k , l ) λ x ( k ) + λ d ( k ) } 1 π λ d ( k ) e x p { − Y 2 ( k , l ) λ d ( k ) } = λ d ( k ) ( λ x ( k ) + λ d ( k ) ) e x p { Y 2 ( k , l ) λ d ( k ) − Y 2 ( k , l ) λ x ( k ) + λ d ( k ) } = 1 1 + ξ k ( l ) e x p { λ x ( k ) Y 2 ( k , l ) λ d ( k ) ( λ x ( k ) + λ d ( k ) ) } = 1 1 + ξ k ( l ) e x p { ξ k ( l ) 1 + ξ k ( l ) γ k ( l ) } \begin{aligned} \Lambda(k,l)&\triangleq\frac{P(Y(k,l)|H_1(k,l))}{P(Y(k,l)|H_0(k,l))}\\ &= \frac{p(Y(k,l)|H_1(k,l))}{p(Y(k,l)|H_0(k,l))}\\ &= \frac{ \frac{1}{\pi [\lambda_x(k)+\lambda_d(k)]}exp\{-\frac{Y^2(k,l)}{\lambda_x(k)+\lambda_d(k)}\}}{ \frac{1}{\pi \lambda_d(k)}exp\{-\frac{Y^2(k,l)}{\lambda_d(k)}\}} \\ &= \frac{\lambda_d(k) }{ (\lambda_x(k)+\lambda_d(k))} exp\{\frac{Y^2(k,l)}{\lambda_d(k)}-\frac{Y^2(k,l)}{\lambda_x(k)+\lambda_d(k)}\}\\ &= \frac{1 }{ 1+\xi_k(l)} exp\{\frac{\lambda_x(k)Y^2(k,l) }{ \lambda_d(k)(\lambda_x(k)+\lambda_d(k))} \}\\ &= \frac{1 }{ 1+\xi_k(l)} exp\{\frac{\xi_k(l) }{ 1+\xi_k(l)}\gamma_k(l) \}\\ \end{aligned} Λ(k,l)​≜P(Y(k,l)∣H0​(k,l))P(Y(k,l)∣H1​(k,l))​=p(Y(k,l)∣H0​(k,l))p(Y(k,l)∣H1​(k,l))​=πλd​(k)1​exp{−λd​(k)Y2(k,l)​}π[λx​(k)+λd​(k)]1​exp{−λx​(k)+λd​(k)Y2(k,l)​}​=(λx​(k)+λd​(k))λd​(k)​exp{λd​(k)Y2(k,l)​−λx​(k)+λd​(k)Y2(k,l)​}=1+ξk​(l)1​exp{λd​(k)(λx​(k)+λd​(k))λx​(k)Y2(k,l)​}=1+ξk​(l)1​exp{1+ξk​(l)ξk​(l)​γk​(l)}​ 此时得到
p ( k , l ) = r 1 1 + ξ k ( l ) e x p { ξ k ( l ) 1 + ξ k ( l ) γ k ( l ) } 1 + r 1 1 + ξ k ( l ) e x p { ξ k ( l ) 1 + ξ k ( l ) γ k ( l ) } = { 1 + q ( k , l ) 1 − q ( k , l ) ( 1 + ξ ( k , l ) ) e x p ( − v ( k , l ) ) } − 1 \begin{aligned} p(k,l) &= \frac{r \frac{1 }{ 1+\xi_k(l)} exp\{\frac{\xi_k(l) }{ 1+\xi_k(l)}\gamma_k(l) \}}{1+r \frac{1 }{ 1+\xi_k(l)} exp\{\frac{\xi_k(l) }{ 1+\xi_k(l)}\gamma_k(l) \}}\\ &= \big \{1+\frac{q(k,l)}{1-q(k,l)}(1+\xi(k,l))exp(-v(k,l)) \big \}^{-1} \end{aligned} p(k,l)​=1+r1+ξk​(l)1​exp{1+ξk​(l)ξk​(l)​γk​(l)}r1+ξk​(l)1​exp{1+ξk​(l)ξk​(l)​γk​(l)}​={1+1−q(k,l)q(k,l)​(1+ξ(k,l))exp(−v(k,l))}−1​

其中$q(k,l)\triangleq P(H_0(k,l))$被定义为语音不存在的先验概率。$\xi (k,l)\triangleq {\lambda }_x(k,l)/{\lambda }_d(k,l)$为先验信噪比，$\gamma (k,l)\triangleq |Y(k,l){|}^2/{\lambda }_d(k,l)$，$v(k,l)\triangleq \gamma (k,l)\xi (k,l)/(1+\xi (k,l))$。众所周知，噪声抑制算法除了一个有效的增益消除公式，剩下最难的就是噪声的估计和跟踪了。上面公式内容来解析，这是所谓的先验信噪比估计器和语音缺席概率SAP联合 估计器。作者推导的先验信噪比估计算法为：
ξ ^ ( k , l ) = α G H 1 2 ( k , l − 1 ) γ ( k , l − 1 ) + ( 1 − α ) m a x { γ ( k , l ) − 1 , 0 } \hat{\xi}(k,l)=\alpha G_{H_1}^2(k,l-1)\gamma(k,l-1)+(1-\alpha)max\{\gamma(k,l)-1,0\} ξ^​(k,l)=αGH1​2​(k,l−1)γ(k,l−1)+(1−α)max{γ(k,l)−1,0}
SAP算法论文也给出，偷懒不再敲键盘了。

MCRA是怎么一回事

递归平均

后来cohen又发表Noise Estimation by Minima Controlled Recursive Averaging for Robust Speech Enhancement和Noise Spectrum Estimation in Adverse Environments:Improved Minima Controlled Recursive Averaging了两篇论文，提出并改进了MCRA算法，令噪声估计在低SNR的情况下误差更小。回顾这个算法，首先要了解就是他其实是将最小值统计跟踪和递归平均两种噪声估计的理论进行了有机的整合。继续以上的种种假设，又引入了一个新的假设公式：
$$Y(k,l)=\sum _{n=0}^{N-1}y(n+lM)h(n)e^{-j\frac{2\pi }{N}nk}$$
前文给出了$H_0(k,l)$和$H_1(k,l)$假设下被观测信号的条件概率密度公式，另外根据加性噪声的假设，定义这两个假设的被观测信号频域表达：
$$\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}$$
令${\lambda }_d(k,l)=E[|D(k,l){|}^2]$表示噪声的方差，这样一个非常常用的时间循环迭代平滑方法可以用于估计噪声的变化，设定一个新的条件：
$$\begin{array}{rl}& H_0^{\prime }(k,l):{\hat{\lambda }}_d(k+1,l)={\alpha }_d{\hat{\lambda }}_d(k,l)+(1-{\alpha }_d)|Y(k,l){|}^2\\ & H_1^{\prime }(k,l):{\hat{\lambda }}_d(k+1,l)={\hat{\lambda }}_d(k,l)\end{array}$$
如果假设$p^{\prime }(k,l)\triangleq P(H^{\prime }(k,l)|Y(k,l))$是在当前$Y(k,l)$下的信号出现条件概率，结合上边的假设可以设计出继续信号出现概率的平滑公式如下：
$$\begin{array}{rl}{\hat{\lambda }}_d(k+1,l)& ={\hat{\lambda }}_d(k,l)p^{\prime }(k,l)+({\alpha }_d{\hat{\lambda }}_d(k,l)+(1-{\alpha }_d)|Y(k,l){|}^2)(1-p^{\prime }(k,l))\\ & ={\hat{\alpha }}_d{\hat{\lambda }}_d(k,l)+(1-{\hat{\alpha }}_d)|Y(k,l){|}^2\end{array}$$
新的平滑因子加入了信号出现概率的影响，而且是一个时频变化的变量：
$${\hat{\alpha }}_d(k,l)={\alpha }_d+(1-{\alpha }_d)p^{\prime }(k,l)$$
以上表达式是一个非常典型的基于语音信号出现概率的一阶时间递归方程。这里MCRA的Recursive Average算是有了，那么minima controlled从何谈起呢？答案就在语音信号概率的计算上面。OMLSA的语音出现概率是基于似然比检验得到的，MCRA的语音出现概率是基于带噪语音功率谱与其局部最小值的比来得到的。局部最小值的统计算法

最小值控制的语音出现概率估计

在给定一帧的子带中，这个语音出现概率决定于：附近受噪话音能量（local energy）与世间窗口中最小能量的比值。这个local energy通过平滑时频窗口的STFT幅度平方来获取，具体方法是先定义一个窗函数获取当前频率的附近能量：
$$S_f(k,l)=\sum _{i=-w}^{w}b(i)|Y(k-i,l){|}^2$$
而时域的平滑采用一阶递归方程：$$S(k,l)={\alpha }_{s}S(k,l-1)+(1-{\alpha }_s)S_f(k,l)$$
论文里还有用到贝叶斯判决规则来证明这个比例是单调的，但似乎多此一举，估计和跟踪基本方法在【4】中很好地归纳，先定义最小值搜索的伪代码：

if mod(l/L)=0
	S_min(k,l) = min(S_tmp(k,l-1),S(k,l))
	S_tmp(k,l) = S_f(k,l)
else
	S_min(k,l) = min(S_min(k,l-1),S(k,l))
	S_tmp(k,l) = min(S_tmp(k,l-1),S(k,l))
end	

这个算法还有一些优化的改进，一般窗口取0.8-1.4s左右。
有了最小值搜索算法，MCRA的整个流程大致梳理如下：

1. 先计算$S_f(k,l)$和平滑后的$S(k,l)$。
2. 启用 最小值搜索算法来得到$S_{min}(k,l)$
3. 计算$S_r(k,l)=\frac{S(k,l)}{S_{min}(k,l)}$，并根据阈值$\delta$进行比较
   $$\begin{gathered} if{S}_r(k,l)>\delta \\ p(k,l)=1\text{语音存在} \\ else \\ p(k,l)=0\text{语音存在} \end{gathered}$$
4. 再次建立一个一阶递归方程$$\hat{p}(k,l)={\alpha }_p\hat{p}(k,l-1)+(1-{\alpha }_p)p(k,l)$$得出的估计值作为$p^{\prime }(k,l)$带入$${\hat{\alpha }}_d(k,l)={\alpha }_d+(1-{\alpha }_d)p^{\prime }(k,l)$$得出${\hat{\alpha }}_d(k,l)$
5. 利用${\hat{\alpha }}_d(k,l)$求得新的噪声估计$${\hat{\lambda }}_d(k+1,l)={\hat{\alpha }}_d{\hat{\lambda }}_d(k,l)+(1-{\hat{\alpha }}_d)|Y(k,l){|}^2$$

从MCRA到IMCRA

MCRA-2

这时罗爱洲提出的对MCRA的改进，主要是替换了最小值搜索算法，其次是对$\delta$门限做出了频率相关的改进，最小值搜索采用了连续频谱最小值跟踪，伪代码如下：
$$\begin{array}{rl}& if{S}_{min}(k,l-1)<S(k,l)\\ & S_{min}(k,l)=\gamma S_{min}(k,l-1)+\frac{1-\gamma }{1-\beta }(S(k,l)-\beta S(k,l-1))\\ & else\\ & S_{min}(k,l)=S(k,l)\\ & end\end{array}$$
这里对${\alpha }_d$提出了0.7~0.9的取值范围，$beta=0.96,\gamma =0.998$作为经验值。而对$\delta$按照频率调整的公式如下：
$$\delta (k)=\{\begin{array}{l}2,1\le k\le LF\\ 2,LF\le k\le MF\\ 5,MF\le k\le Fs/2\end{array}$$这里$LF=1kHz，MF=3kHz，Fs=NyqustFreq$，具体的取值也可以根据应用场景噪声调整。
下图是一个频率分量上的噪声跟踪图，从实际对比来看，其实两者差距并不大，各有所唱吧

IMCRA

水到渠成，懒得写了，最后贴一张图，比较一下三种跟踪方法在omlsa加持下的对比。从语音存在概率来看，imcra确实很优秀，但最终结果其实mcra好像反而更稳。

参考：

1.Optimal Speech Enhancement Under Signal Presence Uncertainty Using Log-Spectral Amplitude Estimator 和
2.Noise Estimation by Minima Controlled Recursive Averaging for Robust Speech Enhancement
3.关于 IMCRA+OMLSA 语音降噪算法的详细解释
4.语音增强理论与实践， （美）罗爱洲， （译）高毅等