单通道语音增强之统计信号模型

1. 信号估计理论简述
2. 最大似然估计ML
3.贝叶斯估计
- 3.1 最小均方估计（MMSE）
- 3.2 最大后验估计 MAP
  - 3.2.1 幅度和相位混合最大后验估计
  - 3.2.2 纯幅度最大后验估计

本文代码位于GitHub。

1. 信号估计理论简述

信号估计理论是现代统计处理的基础课题[@ZhangXianDa2002ModernSP]，在通信、语音、图像领域均有广泛应用。语音增强，就是从带噪的语音测量信号中估计原始的无噪语音，这是典型的信号估计问题。《语音增强–理论与实践》[@loizou2007speech]一书中列举了用于语音增强的一系列统计模型。

假设麦克风采集到的带噪语音序列为 $y[n]$ ，并且噪声都是加性噪声。则带噪语音序列为无噪语音序列与噪声序列的和。原始语音信号与噪声均可视为随机信号。

$y[n] = x[n] + d[n]$

在时域对$x[n]$进行估计是非常困难的，通过傅立叶变换，我们可以将信号分解为频域上互相独立的系数。信号估计模型转变为对每一个频点的系数进行估计的模型，不同频点之间的参数是相互独立的。

$Y(\omega_{k}) = X(\omega_{k}) + D(\omega_{k})$

这个方法就叫做统计信号谱分析（Statistical Spectral Analysis)。显然地，纯净信号谱$X(\omega_{k})$带有幅度与相位两参数，我们实际上是对幅度$X_{k}$和相位$\theta_{y}(k)$进行参数估计。

$Y_{k}e^{j\theta_{y}(k)} = X_{k}e^{j\theta_{x}(k)} + D_{k}e^{j\theta_{d}(k)}$

重组谱幅度和谱相位估计值即可恢复纯净语音谱，估计值用上标来表示：$\hat{X}(\omega_{k})=\hat{X}(k)e^{j\hat{\theta}_{x}(k)}$。

实际信号的幅度和相位是不方便直接用在运算过程中的，因为信号取值范围不定，且瞬时变化。在噪声抑制领域，更常用的语音谱估计方法是对抑制增益(Suppression Gain)进行估计，不同的估计准则称为抑制准则(Supression Rule)。

$\hat{X}(\omega_{k}) = H_{k}Y(\omega_{k})$

通常会根据先验信噪比、后验信噪比来估计抑制增益$H_{k}$。并且可以在只有噪声出现的时刻更新$H_{k}$，在语音存在的时刻进行抑制，无须每帧去调用噪声抑制算法，计算过程比直接估计信号谱灵活。

综上，语音增强的典型流程就是：

对带噪语音y[n]分帧，每一帧进行DFT得到$Y(\omega_{k})$。
估计或者沿用上一帧的抑制增益$H_{k}$，得到纯净语音谱$\hat{X}(\omega_{k})$
对$\hat{X}(\omega_{k})$进行IDFT,得到纯净语音序列的估计$x[n]$。

为了估计模型的建模，对测量信号、估计信号、噪声信号都需要作一些数学上的假设和简化。其中对噪声一般会作以下假设：

噪声是与语音独立的加性噪声；
每一帧噪声的统计分布是稳态的；
噪声的傅立叶级数是零均值复高斯分布。

$\left[D(\omega_{k})\right] = \frac{1}{\pi\lambda_{d}(k)} \exp\left[-\frac{|D(\omega_{k})|^2}{\lambda_{d}(k)}\right]$

2. 最大似然估计ML

如果不考虑信号的先验分布，即认为信号值是确定信号，而不是随机信号，我们只需要分析含有信号$x$为参数的带噪信号$y$的概率分布$p(y;x)$，并使之最大。这种估计方法叫最大似然估计器（Maximum Likelihood Estimation）。将$y$的带参概率分布$p(y;x)$称为似然函数（Likelihood Function）。对纯净信号$x$的估计，表达为求解合适的$x$值，使得似然函数$p(y;x)$最大。

$\hat{x} = \argmax_{x}p(y;x)$

文献[@mcaulay1980speech]最早将最大似然估计法用在语音增强领域。对于纯净语音，可以假设纯净语音幅度$X_{k}$和相位$\theta_{y}(k)$是未知但确定，无需考虑其先验概率分布。最大似然语音增强模型表达为：

$\hat{X}_{k},\hat{\theta}_{k} = \argmax_{X_{k},\theta_{k}} p\left[Y(\omega_{k});X_{k},\theta_{k}\right]$

把$D_{k}e^{\theta_{d}(k)}=Y_{k}e^{\theta_{y}(k)} - X_{k}e^{\theta_{x}(k)} $代入噪声零均值复高斯分布公式中，得到：

$p\left[Y(\omega_{k});X_{k},\theta_{k}\right] = \frac{1}{\pi \lambda_{d}(k)}\exp\left[-\frac{|Y(\omega_{k}) - X_{k}e^{j\theta_{x}(k)}|^2}{\lambda_{d}(k)}\right]$ $H_{k} = \frac{1}{2}+\frac{1}{2}\sqrt{\frac{\xi_{k}}{1+\xi_{k}}}$

3.贝叶斯估计

如果比最大似然估计更进一步，考虑待估计量$x$也是随机变量，且$x$的先验分布为$p(x)$，这种假设下的估计方法叫做贝叶斯估计[@ZhangXianDa2002ModernSP]。定义估计值$\hat{x}$与实际值$x$之间的误差函数为$c(\hat{x},x)$，贝叶斯估计器的目标即为找出是平均误差$E[c(\hat{x},x)]$最小的估计值$x$。

$\hat{x} = \argmin_{\hat{x}} E[c(\hat{x}, x)]$

对于待估计的纯净语音谱，贝叶斯估计器可以表达为：

$\hat{X}(\omega_{k}) = \argmin_{\hat{X_{k}}} E[c(\hat{X}(\omega_{k}), X(\omega_{k}))]$

误差$c(\hat{x},x)$的期望值取决于待测信号与测量信号的联合概率分布。

$E[c(\hat{x}, x)] = \int_{x}\int_{y}c(\hat{x}, x)p(x, y)dxdy$

根据条件概率密度公式对联合概率密度进行分解：

$E[c(\hat{x}, x)] = \int_{x}\int_{y}c(\hat{x}, x)p(x|y)p(y)dxdy = \int_{y} \left[ \int_{x}c(\hat{x}, x)p(x|y)dx\right] p(y)dy$

因为估计值$\hat{x}$与$p(y)$相互独立，可以把外层对$y$的积分消除，也就是：

$\hat{x} = \argmin_{\hat{x}} E[c(\hat{x}, x)|y]$

误差函数$c(\hat{x}, x)$的计算模型，会引出不同种类的估计器。典型的误差函数有几种类型[@ZhangXianDa2002ModernSP]：

平方误差函数，对应最小均方估计。
$c(\hat{x}, x) = (\hat{x} - x)^2$
绝对值误差函数，对应条件中位数估计。
$c(\hat{x}, x) = |\hat{x} - x|$
均匀误差函数，对应最大后验估计。

$c(\hat{x}, x) = \left\{ \begin{array}{lr} 1 \quad (|\hat{x}-x|\geq\Delta/2) & \\ &\Delta>0\\ 0 \quad (|\hat{x}-x|<\Delta/2) & \end{array} \right.$

3.1 最小均方估计（MMSE）

最小均方估计(Minimum Mean Square Error Estimation)[@ephraim1985speech]使用平方误差函数，平均误差为：

$E[c(\hat{x}, x)|y]=\int_{x}(\hat{x} - x)^2 p(x|y)dx$

为求平均误差的极大值，可对估计量$\hat{x}$求导，并求极值点。

$\frac{\partial E[c(\hat{x}, x)|y]}{\partial \hat{x}}=\int_{x}2(\hat{x} - x) p(x|y)dx =0$

显然极值点的$\hat{x}$为被估计量$x$的条件均值：

$\hat{x} =\int_{x}x p(x|y)dx = E[ x| y]$

亦可以使用贝叶斯公式展开，得到：

$\begin{aligned} \hat{x} = \int_{x}x \frac{p(y|x)p(x)}{p(y)} dx \\ = \frac{\int_{x}xp(y|x)p(x)dx}{\int_{x}p(y|x)p(x)dx} \end{aligned}$

在语音增强的模型里，纯净语音谱的估计为其在带噪语音谱下的条件均值。

$\hat{X}(\omega_{k}) = E[X(\omega_{k})| Y(\omega_{k})]$

上式中，为得到纯净语音谱需要分别估计谱幅度和谱相位: $\hat{X}, \hat{\theta}_{x}$。

$\hat{X}_{k}e^{j\hat{\theta}_{x}(k)} = E[ X_{k}e^{j\theta_{x}(k)}| Y(\omega_{k})]$

同时估计谱幅度和谱相位是很难的，研究者提出了许多分别估计谱幅度和谱相位的方法，估计完成后再用两者重组复语音信号。

$\hat{X}_{k} = E[ X_{k}| Y(\omega_{k})]$ $e^{j\hat{\theta}_{x}(k)} = E[ e^{j\theta_{x}(k)}| Y(\omega_{k})]$

3.1.1 MMSE谱幅度估计

最小均方根估计器MMSE short-time spectral amplitude

$\hat{X_{k}} = E[ X_{k}| Y(\omega_{k})]$ $v_{k} = \frac{\xi_{k}}{1+\xi_{k}}\gamma_{k}$ $H_{k} =\frac{\sqrt{\pi v_{k}}}{ 2\gamma_{k}}[(1+v_{k})I_{0}(\frac{v_k}{2})+v_{1}I_{1}(\frac{v_{k}}{2})]\exp(\frac{-v_{k}}{2})$

	def mmse_stsa_gain(parameters=None):
		gamma = parameters['gamma']
		ksi = parameters['ksi']
			
		vk = ksi * gamma / (1 + ksi)
		j0 = np.i0(vk / 2)
		j1 = np.i1(vk / 2)
			
		A = sqrt(pi * vk) / 2 / gamma
		B = (1 + vk) * j0 + vk * j1
		C = np.exp(-0.5 * vk)
		gain = A * B * C
		return gain

3.1.2 MMSE对数谱幅度估计

对数最小均方根估计器The MMSE log spectral amplitude (MMSE-LSA)，或者缩写为LogMMSE估计器。

$c(\hat{X_{k}}, X_{k}) = (\log{X_{k}}- \log{X_{k}})^2$ $\log{\hat{X_{k}}} = E[ \log{X_{k}}| Y(\omega_{k})]$ $H_{k} = \frac{\xi_{k}}{1+\xi_{k}}\exp(\frac{1}{2}\int_{v_{k}}^{\infty}\frac{e^{-t}}{t}dt)$

	def logmmse_gain(parameters=None):
		gamma = parameters['gamma']
		ksi = parameters['ksi']
		A = ksi / (1 + ksi)
		vk = A * gamma
		ei_vk = 0.5 * expn(1, vk)
		gain = A * np.exp(ei_vk)
		return gain

3.1.3 MMSE平方谱幅度估计

频谱幅度平方估计器MMSE magnitude squared[@wolfe2003efficient]

$\hat{X_{k}^2} = E[ X_{k}^2| Y_{k}]$ $\hat{X_{k}^2} = \frac{\xi_{k}}{1+\xi_{k}} (\frac{1+v_{k}}{\gamma_{k}})Y_{k}^2$ $H_{k} = \sqrt{\frac{\xi_{k}}{1+\xi_{k}} (\frac{1+v_{k}}{\gamma_{k}})}$

	def mmse_sqr_gain(parameters=None):
		gamma = parameters['gamma']
		ksi = parameters['ksi']
			
		vk = ksi * gamma / (1 + ksi)
		j0 = np.i0(vk / 2)
		j1 = np.i1(vk / 2)
			
		A = ksi / (1 + ksi)
		B = (1 + vk) / gamma
		gain = sqrt(A * B)
		return gain

3.2 最大后验估计 MAP

当贝叶斯估计器采用均匀误差函数时，平均误差为：

$\begin{aligned} E[c(\hat{x}, x)|y] &= \int_{-\infty}^{\hat{x}-\Delta/2} p(x|y)dx + \int_{\hat{x}+\Delta/2}^{\infty} p(x|y)dx \\ &= 1-\int_{\hat{x}\Delta/2}^{\hat{x}+\Delta/2} p(x|y)dx \end{aligned}$

显然要使得平均误差最小，就是要求目标估计$\hat{x}$，使得 $p(x|y)$最大。这种估计模型称作最大后验估计(Maximum A Posteriori Estimation)。这个模型的意思是只有估计值$\hat{x}$与原始值$x$相等，误差才为0，其他时候误差均匀为1。估计值可以表达为：

$\hat{x} =\argmin_{x} p(x|y)$

在语音增强的模型里，纯净语音谱的估计为其在带噪语音谱下的条件均值。

$\hat{X}(\omega_{k}) = \argmin_{ X(\omega_{k})}\quad p\left[X(\omega_{k})| Y(\omega_{k})\right]$

文献[@wolfe2003efficient]提出了两种基于最大后验估计(Maximum A Posteriori Estimation)的语音增强算法。一种是同时求解幅度和相位的混合最大后验估计：另外一种是单纯估计幅度的方法，两种估计器最后的噪声抑制增益略有不同。

3.2.1 幅度和相位混合最大后验估计

[@wolfe2003efficient]

$H_{k} = \frac{\xi_{k}+\sqrt{\xi_{k}^2+2(1+\xi_{k})\frac{\xi_{k}}{\gamma_{k}}}}{2(1+\xi_{k})}$

def map_joint_gain(parameters=None):
	gamma = parameters['gamma']
	ksi = parameters['ksi']
	
	eps = 1e-6
	gain = (ksi + sqrt(ksi^ 2 + 2 * (1.0 + ksi)* ksi/ (gamma + eps))) / 2.0/ (1.0 + ksi)
	return gain

3.2.2 纯幅度最大后验估计

[@wolfe2003efficient]

$H_{k} = \frac{\xi_{k}+\sqrt{\xi_{k}^2+(1+\xi_{k})\frac{\xi_{k}}{\gamma_{k}}}}{2(1+\xi_{k})}$

def map_sa_gain(parameters=None):
	gamma = parameters['gamma']
	ksi = parameters['ksi']
	
	eps = 1e-6
	gain = (ksi + sqrt(ksi^ 2 + (1.0 + ksi)* ksi/ (gamma + eps))) / 2.0/ (1.0 + ksi)
	return gain

参考文献

[1] 张贤达, 现代信号处理. 清华大学出版社有限公司, 2002.

[2] P. C. Loizou, Speech enhancement: theory and practice. CRC press, 2007.

[3] R. McAulay and M. Malpass, Speech enhancement using a soft-decision noise suppression filter,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 28, no. 2, pp. 137-145,1980.

[4] Y. Ephraim and D. Malah, Speech enhancement using a minimum mean-square error log-spectral amplitude estimator,” IEEE transactions on acoustics, speech, and signal processing,vol. 33, no. 2, pp. 443-445, 1985.

[5] P. J. Wolfe and S. J. Godsill, Efficient alternatives to the ephraim and malah suppression rule for audio signal enhancement,” EURASIP Journal on Applied Signal Processing, vol. 2003, pp.1043-1051, 2003.

wjchen