記錄一下單通道降噪的一個經典方法, MMSE-STSA, MMSE-LSA, 已經是 1984 左右的文章了. 單通道降噪 OMLSA 也從這衍生出來的. 我們先從 MMSE-STSA 說起, 全名是 minimum mean-square error short time spectral amplitude.
$y(t)=x(t)+d(t),0leq tleq T$
$x$, $d$, $y$ 分別是 speech, noise, 和收到的 noisy signal, 其中 $x$, $d$ 相互獨立. 相對應的第 $k$ 個 frequency bin 如下:
$$X_k=A_kexp(jalpha_k) \
D_k \
Y_k=R_kexp(jtheta_k)$$
MMSE-STSA $^{[1]}$
目標函式為
$$begin{align}
argmin_{hat{A}_k}{mathbb{E}left[left(A_k-hat{A}_kright)^2vert y(t),0leq tleq Tright]}
end{align}$$ 最佳解為
$$begin{align}
hat{A}_k=mathbb{E}left[A_kvert y(t),0leq t leq Tright]
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_kvert y(t),0leq t leq Tright]=mathbb{E}left[A_kvert Y_0,Y_1,...right] \
=mathbb{E}left[A_kvert Y_kright] \
=int_0^{infty}int_0^{2pi}a_k p(a_k,alpha_kvert Y_k)dalpha_k d a_k = int_0^{infty}int_0^{2pi}a_k frac{p(a_k,alpha_k,Y_k)}{p(Y_k)}dalpha_k d a_k \
=frac{ int_0^{infty}int_0^{2pi}a_k p(Y_kvert a_k,alpha_k) p(a_k,alpha_k) dalpha_k d a_k }{ int_0^{infty}int_0^{2pi} p(Y_kvert a_k,alpha_k) p(a_k,alpha_k) dalpha_k d a_k }
end{align}$$ 其中 (3) 到 (4) 我們假設每個 frequency bin 是獨立的
由於我們假設每個 frequency bin 都是 complex Gaussian distribution, 因此 (6) 的機率分佈如下定義:
$$begin{align}
p(Y_kvert a_k,alpha_k)=frac{1}{pilambda_d (k)}expleft[ -frac{1}{lambda_d (k)}vert Y_k - a_k e^{jalpha_k} vert^2 right] \
p(a_k,alpha_k)=frac{1}{pilambda_x (k)}expleft[-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 ^2right]=A_k^2 \
lambda_d (k)=mathbb{E}left[vert D_k vert ^2right]
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}expleft(-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_ktriangleq frac{xi_k}{1+xi_k}gamma_k \
color{orange}{
xi_ktriangleqfrac{lambda_x (k)}{lambda_d (k)}
} \
color{orange}{
gamma_ktriangleqfrac{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}
argmin_{hat{A}_k}{mathbb{E}left[left(log A_k-loghat{A}_kright)^2vert y(t),0leq tleq Tright]}
end{align}$$ 同樣經過不是人類的計算後得到:
$$begin{align}
hat{A}_k=frac{xi_k}{1+xi_k}expleft[frac{1}{2}int_{upsilon_k}^{infty}frac{e^{-t}}{t}dtright]R_k
end{align}$$ [3] 給出了一個好算的近似結果
$$begin{align}
int_{upsilon_k}^{infty}frac{e^{-t}}{t}dtapprox left{
begin{array}{rcl}
-2.31log_{10}(upsilon_k)-0.6mbox{ for }upsilon_k<0.1 \
-1.544log_{10}(upsilon_k)+0.166mbox{ for }0.1lequpsilon_kleq 1 \
10^{-(0.52upsilon_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_dhat{lambda}_d(k,l)+(1-alpha_d)|Y(k,l)|^2right](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
