懷舊篇, 單通道降噪, MMSE-STSA, MMSE-LSA 方法


記錄一下單通道降噪的一個經典方法, MMSE-STSA, MMSE-LSA, 已經是 1984 左右的文章了. 單通道降噪 OMLSA 也從這衍生出來的. 我們先從 MMSE-STSA 說起, 全名是 minimum mean-square error short time spectral amplitude.
$y(t)=x(t)+d(t),0\leq t\leq T$
$x$, $d$, $y$ 分別是 speech, noise, 和收到的 noisy signal, 其中 $x$, $d$ 相互獨立. 相對應的第 $k$ 個 frequency bin 如下:
$$X_k=A_k\exp(j\alpha_k) \\ D_k \\ Y_k=R_k\exp(j\theta_k)$$


MMSE-STSA $^{[1]}$

目標函式為
$$\begin{align} \arg\min_{\hat{A}_k}{\mathbb{E}\left[\left(A_k-\hat{A}_k\right)^2\vert y(t),0\leq t\leq T\right]} \end{align}$$ 最佳解為
$$\begin{align} \hat{A}_k=\mathbb{E}\left[A_k\vert y(t),0\leq t \leq T\right] \end{align}$$ 但關鍵是我們不知道 clean speech 的 amplitude $A_k$, 那該怎麼估呢?

首先我們對每個 frequency bin 的分布假設為 Gaussian distribution (complex).

引用原文 “Since the Fourier coefficient is, after all, a weighted sum (or integral) of random variables resulting from the random process samples”, 在一個短時的 frame 中大致上是 stationary, 因此可以看作是一個 WSS 的 ramdom process, 再加上 cental limit theorem, 就當作高斯分布吧.

套用 Guassian distribution 假設, 做如下推導
$$\begin{align} \hat{A}_k=\mathbb{E}\left[A_k\vert y(t),0\leq t \leq T\right]=\mathbb{E}\left[A_k\vert Y_0,Y_1,...\right] \\ =\mathbb{E}\left[A_k\vert Y_k\right] \\ =\int_0^{\infty}\int_0^{2\pi}a_k p(a_k,\alpha_k\vert Y_k)d\alpha_k d a_k = \int_0^{\infty}\int_0^{2\pi}a_k \frac{p(a_k,\alpha_k,Y_k)}{p(Y_k)}d\alpha_k d a_k \\ =\frac{ \int_0^{\infty}\int_0^{2\pi}a_k p(Y_k\vert a_k,\alpha_k) p(a_k,\alpha_k) d\alpha_k d a_k }{ \int_0^{\infty}\int_0^{2\pi} p(Y_k\vert a_k,\alpha_k) p(a_k,\alpha_k) d\alpha_k d a_k } \end{align}$$ 其中 (3) 到 (4) 我們假設每個 frequency bin 是獨立的
由於我們假設每個 frequency bin 都是 complex Gaussian distribution, 因此 (6) 的機率分佈如下定義:
$$\begin{align} p(Y_k\vert a_k,\alpha_k)=\frac{1}{\pi\lambda_d (k)}\exp\left[ -\frac{1}{\lambda_d (k)}\vert Y_k - a_k e^{j\alpha_k} \vert^2 \right] \\ p(a_k,\alpha_k)=\frac{1}{\pi\lambda_x (k)}\exp\left[-\frac{a_k^2}{\lambda_x (k)}\right] \end{align}$$ 注意到 (7) 能這麼寫是因為我們知道 $x$ and $d$ 互相獨立, 因此在給定 $x$ 的情形下, 只是改變 mean 的位置, 其 variance 仍由 $d$ 來決定. 另外:
$$\begin{align} \lambda_x (k)=\mathbb{E}\left[\vert X_k \vert ^2\right]=A_k^2 \\ \lambda_d (k)=\mathbb{E}\left[\vert D_k \vert ^2\right] \end{align}$$ 表示第 $k$ 個 bin 的 speech and noise 的 variance
將 (7) and (8) 帶入 (6) 並感謝偉大的作者推導得到:
$$\begin{align} \hat{A}_k=\Gamma(1.5)\frac{\sqrt{\upsilon_k}}{\gamma_k}M(-0.5;1;-\upsilon_k)R_k \\ \hat{A}_k=\Gamma(1.5)\frac{\sqrt{\upsilon_k}}{\gamma_k}\exp\left(-\frac{\upsilon_k}{2}\right)\left[(1+\upsilon_k)I_0(\frac{\upsilon_k}{2})+\upsilon_k I_1(\frac{\upsilon_k}{2})\right]R_k \end{align}$$ 其中 $\Gamma$ 表示 gamma function, $\Gamma(1.5)=\sqrt{\pi}/2$; $M(a;c;x)$ 是 confluent hypergeometric function (這是外星符號吧), $I_0$ and $I_1$ 是 modified Bessel funciton of zero and first order. 總之就是能帶入計算的東西, 最重要, 也是需要我們估計的變數如下:
$$\begin{align} \upsilon_k\triangleq \frac{\xi_k}{1+\xi_k}\gamma_k \\ \color{orange}{ \xi_k\triangleq\frac{\lambda_x (k)}{\lambda_d (k)} } \\ \color{orange}{ \gamma_k\triangleq\frac{R_k^2}{\lambda_d (k)} } \\ \end{align}$$ $\xi_k$ 和 $\gamma_k$ 分別稱為 prior SNR 和 posterior SNR. 總之如能估出 $\xi_k$ 和 $\gamma_k$, 我們就能計算出 gain 值, 之後的方法如 LSA, OMLSA 也都如此. 文章後面會使用 MCRA 來估算這兩個 SNR.

現在就算傳統方法一般也很少使用 MMSE-STSA, 至少會使用 LSA 取代. LSA 有近似的計算方式, 因此我們也不糾結 (12) 到底怎麼算出來.


MMSE-LSA $^{[2]}$

大致想法跟流程跟上面一樣(只是我算不出來), 只是目標函數針對 log 值來計算
$$\begin{align} \arg\min_{\hat{A}_k}{\mathbb{E}\left[\left(\log A_k-\log\hat{A}_k\right)^2\vert y(t),0\leq t\leq T\right]} \end{align}$$ 同樣經過不是人類的計算後得到:
$$\begin{align} \hat{A}_k=\frac{\xi_k}{1+\xi_k}\exp\left[\frac{1}{2}\int_{\upsilon_k}^{\infty}\frac{e^{-t}}{t}dt\right]R_k \end{align}$$ [3] 給出了一個好算的近似結果
$$\begin{align} \int_{\upsilon_k}^{\infty}\frac{e^{-t}}{t}dt\approx \left\{ \begin{array}{rcl} -2.31\log_{10}(\upsilon_k)-0.6\mbox{ for }\upsilon_k<0.1 \\ -1.544\log_{10}(\upsilon_k)+0.166\mbox{ for }0.1\leq\upsilon_k\leq 1 \\ 10^{-(0.52\upsilon_k+0.26)}\mbox{ for }\upsilon_k>1 \\ \end{array}\right. \end{align}$$ 另外還有 optimally-modified log-spectral amplitude (OMLSA) [4] 方法, 作者有提供 MATLAB codes. 這算單通道降噪標配了, 但實驗結果對聽覺有幫助, 對 WER 不一定降低. 總之不管哪一種方法, 都必須很好的估出 prior and posterior SNR.


MCRA Prior/Posterior SNR 估計

針對 STFT 時間 $l$, frequency bin $k$ 來說, 假設我們已估出來 speech presence probability $p(k,l)$, 我們可以這麼 update noise 的 variance:
$$\begin{align} \hat{\lambda}_d(k,l+1)=\hat{\lambda}_d(k,l)p(k,l)+\left[\alpha_d\hat{\lambda}_d(k,l)+(1-\alpha_d)|Y(k,l)|^2\right](1-p(k,l)) \end{align}$$ 這很好理解, 如果有 speech 的話, noise variance 就沿用原來舊的, 而如果沒有 speech, nosie vaiance 就要用當前 frame 透過 $\alpha_d$ 平滑地更新一下 (就稱這樣的平滑為 $\alpha$ 平滑).

估計 $p(k,l)$ 之前, 文章的做法是都先針對 time and frequency 做平滑. frequency 可選用一個 window (可用類似 Gaussian window), 而時間上的平滑可使用 $\alpha$ 平滑. 令 $S(k,l)$ 為我們平滑後的 spectrum power, 然後對每個 bin 都 tracking 一小段時間的最小值, 令為 $S’(k,l)$. 則很明顯如果 $S(k,l)>\delta S’(k,l)$, 我們就可以認為有 speech, 機率為 1, 否則為 0. 這樣的 speech 機率過了 $\alpha$ 平滑的結果就是 $p(k,l)$. 明確一點寫下為:
$$\begin{align} p(k,l)=\alpha_p p(k,l-1)+(1-\alpha_p)\mathbf{I}[S(k,l)>\delta S'(k,l)] \end{align}$$ 其中 $\mathbf{I}[.]$ 為 indicator function


MCRA 有哪些調整的參數

實際情形有一些需要調整的參數, 列在下面

  • $\alpha_d$: noise variance smoothing
  • $\alpha_p$: speech probability smoothing
  • STFT 的 time and frequency smoothing 參數
  • $\delta$: 判斷當前 frame and bin 是否為 speech 的 threshold
  • tracking minimal power $S’(k,l)$ 的參數, 譬如要用多少個 frame 來找 minimum

待做些實驗才會知道效果…


Reference

  1. Speech Enhancement Using a Minimum Mean-Square Error Short-Time Spectral Amplitude Estimator by Yariv Ephraim and David Malah
  2. Speech Enhancement Using a Minimum Mean-Square Error Log-Spectral Amplitude Estimator by Yariv Ephraim and David Malah
  3. [A Noise Reduction Pre-processor for Mobile Voice Communication] by R. Martin …
  4. Speech enhancement for non-stationary noise environments by Israel Cohen and Baruch Berdugo