記錄一下單通道降噪的一個經典方法, 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
- Speech Enhancement Using a Minimum Mean-Square Error Short-Time Spectral Amplitude Estimator by Yariv Ephraim and David Malah
- Speech Enhancement Using a Minimum Mean-Square Error Log-Spectral Amplitude Estimator by Yariv Ephraim and David Malah
- [A Noise Reduction Pre-processor for Mobile Voice Communication] by R. Martin …
- Speech enhancement for non-stationary noise environments by Israel Cohen and Baruch Berdugo