基于最小均方误差短时谱估计的语音增强

#1 基于最小值控制的噪声估计

Cohen 和 Berdugo 提出了一种基于最小值控制的递归平均（MCRA: Minima Controlled Recursive Averageing）方法来估计噪声。该方法由带噪语音的局部能量值与一定时间范围内的最小值的比较判断某个子带（频点）是否存在语音。如果该子带存在语音则不更新噪声谱，如果不存在则跟带噪语音的功率谱进行加权更新。

具体实现流程：
1、 对每一帧的信号的各点功率谱进行加权平滑，对平滑的结果，然后通过一个时间范围内（range）的最小值来粗略估计噪声功率, 这个最小值为上一个range到本次range当前帧的中的最小值； range大小理论保证说话人语音存在间歇期，以期望在非语音段估计到噪声最小值。
2、 用当前带噪功率谱与最小值估计的噪声功率比较，判断是否该频点是否为语音的可能性。

static void update_noise_prob(NsState *ns)
{
	int i;
	int min_range;
	int N = ns->ps_size;

	for (i=1;i<N-1;i++)
		ns->S[i] = 0.8f * ns->S[i] + 0.05f * ns->ps[i-1] + 0.1f * ns->ps[i] + 0.5f * ns->ps[i +1];
	ns->S[0] = 0.8f * ns->S[0] + 0.2f * ns->ps[0];
	ns->S[N-1] = 0.8f * ns->S[N-1] + 0.2f *  ns->ps[N-1];


	if (ns->nb_adapt==1)
	{
		//for (i=0;i<N;i++)
		//	ns->Smin[i] = ns->Stmp[i] = 0; // stmp[] 若初始化为0， 那么其实下面比较大小在第一个range毫无意义
		// lucas change
		for (i =0; i < N ; i++)
		{
			ns->Smin[i] = 0;
			ns->Stmp[i] = ns->S[i];
		}
		
	}

	if (ns->nb_adapt < 100)
		min_range = 15;
	else if (ns->nb_adapt < 1000)
		min_range = 50;
	else if (ns->nb_adapt < 10000)
		min_range = 150;
	else
		min_range = 300;
	///  Stmp 是一个range内的最小，Smin是上一个range和当前的最小
	if (ns->min_count > min_range)
	{
		ns->min_count = 0; 
		for (i=0;i<N;i++)
		{
			ns->Smin[i] = (ns->Stmp[i] < ns->S[i]) ? ns->Stmp[i] : ns->S[i];
			ns->Stmp[i] = ns->S[i];
		}
	} else {
		for (i=0;i<N;i++)
		{
			ns->Smin[i] = (ns->Smin[i] < ns->S[i]) ? ns->Smin[i] : ns->S[i];
			ns->Stmp[i] = (ns->Stmp[i] < ns->S[i]) ? ns->Stmp[i] : ns->S[i];    
		}
	}
	for (i=0;i<N;i++)
	{
		if(0.4f * ns->S[i] > ns->Smin[i])
			ns->update_prob[i] = 1;
		else
			ns->update_prob[i] = 0;

		/*fprintf (stderr, "%f ", st->S[i]/st->Smin[i]);*/
		/*fprintf (stderr, "%f ", st->update_prob[i]);*/
	}

}

3、更新噪声功率谱noise。当判断该频点被判断为非语音声，则更新噪声功率，由于一段音频一开始基本是非语音段而且大部分背景噪声其实分布相对稳定，更新参数beta 在一开始设计比较大，随着时间变小到一定程度。

	for (i=0;i<N;i++)
	{ //update_prob 声音概率，又么声音概率为零，又么此时能量小于噪声能量则更新噪声
		if (!ns->update_prob[i] || ns->ps[i] < ns->noise[i])
		{
			ns->noise[i] = beta_1* ns->noise[i] + beta * ns->ps[i];
			if (ns->noise[i] < 0)
				ns->noise[i] = 0.0f;
		}
	}

#2 临界带（Bark scale）
声学研究中，使用滤波器组模拟临界带对语音信号进行处理更符合人耳听觉结构模型，Eberhard Zwicker 进行大量的听音实验，对听觉临界带带宽进行合理评估，提出声频范围按临界带划分的方法，即Bark 域。

将线性频域非线性转换到bark 域，如下方程：

x = 0:16000;
y = 13.1*atan(0.00074*x)+2.24*atan( x.^2*1.85e-8)+1e-4*x;
plot(x,y);

则在bark域的滤波器组可由下面函数实现：

function Filter = newFilterBank(banks, samplerate,len)
df = samplerate/(2* len);
max_mel = toBark(samplerate/2);
mel_interval = max_mel / (banks-1);
Filter.nb_banks = banks;
Filter.len = len;
for i=1:Filter.len
    curr_freq = (i-1)*df;
    mel = toBark(curr_freq);
    if(mel > max_mel)
        break;
    end
    
    id1 = floor(mel/mel_interval);
    if id1>(banks-2)
        id1 = banks-2;
        val = 1;
    else
        val = (mel - id1*(mel_interval)) / mel_interval;
    end
    
    id2 = id1 +1;
    
    Filter.bank_left(i) = id1;
    Filter.filter_left(i) = 1 - val;
    Filter.bank_right(i) = id2;
    Filter.filter_right(i) = val;
end
end

function mel = toBark(freq)
mel = 13.1*atan(0.00074*freq)+2.24*atan( freq.^2*1.85e-8)+1e-4*freq;
end

上述函数建立的滤波组如下图：

len = 8000*10;
filterBank = newFilterBank(24, 16000, len);

x=0:len-1;
figure
plot(x,filterBank.filter_left,'.r', x,filterBank.filter_right,'.b')

下文将在Bark域计算降噪增益，用对来线性频域计算得出的增益进行约束，使之不偏离声学模型。

#3 最小均方误差法
最小均方误差法由Ephrain 和 Malah 提出，根据其论文，带噪语音信号模型，若$Y(k)$ 、$S(k)$和 $D(k)$ 分别表示 $y(n)$、 $s(n)$ 和 $d(n)$ 的傅里叶变换后的频谱，$Y(k)$ 的幅值是 $R(k)$, 相位是 $ \varphi_k$，$S(k)$ 的幅值是 $A(k)$ , 并假设噪声时平稳加性高斯白噪声，由贝叶斯公式可得 $A(k)$ 的估计值$\hat{A}(k)$为：
A ^ ( k ) = E { A ( k ) ∣ Y ( k ) } \hat{A}(k) = E\{ A(k) | Y(k) \} A^(k)=E{A(k)∣Y(k)} $$={\int }_0^{\infty }p[a_k|Y(k)]a_{k}d{a}_k$$ $$={\int }_0^{\infty }\frac{p[a_k,Y(k)]}{p[Y(k)]}{a}_{k}d{a}_k$$ $$=\frac{{\int }_0^{2\pi }{\int }_0^{\infty }{a}_{k}p[Y(k),(a_k,{\phi }_k)]d{\phi }_{k}d{a}_k}{p[Y(k)]}$$ $$=\frac{{\int }_0^{2\pi }{\int }_0^{\infty }{a}_{k}p[Y(k)|(a_k,{\phi }_k)]p(a_k,{\phi }_k)d{\phi }_{k}d{a}_k}{{\int }_0^{2\pi }{\int }_0^{\infty }p[Y(k)|(a_k,{\phi }_k)]p(a_k,{\phi }_k)d{\phi }_{k}d{a}_k}$$
式中， E { . } E\{.\} E{.} 代表参数的期望； p { . } p\{.\} p{.} 代表概率密度函数；$p(a_k)$ 是幅值$A(k)$的概率密度函数； $ p(a_k, \varphi_k)$ 为幅相联合概率分布概率。

在假设噪声是平稳加性高斯白噪声的统计模型下，$ p[Y(k) | (a_k, \varphi_k)]$ 和 $p(a_k,{\phi }_k)$ 分别满足：
p [ Y ( k ) ∣ ( a k , φ k ) ] = 1 π p d ( k ) e x p { − 1 p d ( k ) ∣ Y ( k ) − a k e j φ k ∣ 2 } p[Y(k) | (a_k, \varphi_k)] = \frac{1}{\pi p_d(k)}exp\{-\frac{1}{p_d(k)}|Y(k) -a_ke^{j\varphi_k}|^2\} p[Y(k)∣(ak​,φk​)]=πpd​(k)1​exp{−pd​(k)1​∣Y(k)−ak​ejφk​∣2} p ( a k , φ k ) = a k π p s ( k ) e x p { − a k 2 p s ( k ) } p(a_k, \varphi_k) = \frac{a_k}{\pi p_s(k)}exp\{-\frac{a_k^2}{p_s(k)}\} p(ak​,φk​)=πps​(k)ak​​exp{−ps​(k)ak2​​}
式中，$p_s(k)$、$p_d(k)$ 分别表示纯语音信号和噪声的第k个频谱分量的能量，计算如下：
p s ( k ) = E { ∣ S ( K ) ∣ 2 } p_s(k) = E\{|S(K)|^2\} ps​(k)=E{∣S(K)∣2} p d ( k ) = E { ∣ D ( K ) ∣ 2 } p_d(k) = E\{|D(K)|^2\} pd​(k)=E{∣D(K)∣2}
综合上式可得$S(K)$的幅值的估计值$\hat{A}(k)$， 如下所示：
$$\hat{A}(k)=\Gamma (1.5)\frac{\sqrt{{\nu }_k}}{{\gamma }_k}M(-0.5;1;-{\nu }_k)R(k)$$ $$\hat{A}(k)=G_{MMSE}R(k)$$
式中，$ \Gamma(.)$ 是伽马函数， $M(a;c;x)$ 是合流超几何函数。${\nu }_k$定义为：
$${\nu }_k=\frac{{\xi }_k}{1+{\xi }_k}{\gamma }_k$$
式中，${\xi }_k$ 和 ${\gamma }_k$ 分别为是先验和后验信噪比，他们分别定义为：
$${\xi }_k=\frac{p_s(k)}{p_d(k)}$$ $${\gamma }_k=\frac{|y_k{|}^2}{p_d(k)}=\frac{|R(k){|}^2}{p_d(k)}$$
先验信噪比通过Ephrain 和 Malah 提出的使用前一帧信号信息估计当前帧先验信噪比的反馈方法：
ξ ^ ( i , k ) = η A 2 ^ ( i − 1 , k ) p d ( k ) + ( 1 − η ) m a x { γ ( i , k ) − 1 , 0 } \hat{\xi}(i, k) = \eta\frac{\hat{A^2}(i-1,k)}{p_d(k)}+(1-\eta )max\{\gamma(i,k)-1,0\} ξ^​(i,k)=ηpd​(k)A2^(i−1,k)​+(1−η)max{γ(i,k)−1,0}
综合上述的式子，可计算得到增益$G_{MMSE}$，其中$ \Gamma(1.5)M(-0.5; 1; -\nu_k)$可如下近似计算：

x = 0:0.5:10
H = hypergeom(-.5,1,-x)
y = gamma(1.5)* H

通过上面使用matlab 对$ \Gamma(1.5)M(-0.5; 1; -\nu_k)$采样21个点生成表：

static const float table[21] = {
   0.88623f, 1.09501f, 1.28192f, 1.45141f, 1.60682f, 1.75066f, 1.88487f,
   2.01098f, 2.13015f, 2.24335f, 2.35134f, 2.45474f, 2.55407f, 2.64977f,
   2.74218f, 2.83163f, 2.91837f, 3.00263f, 3.08461f, 3.16449f, 3.24241f};
   
static inline float hypergeom_gain(float xx)
{
   int ind;
   float integer, frac;
   float x;
   x = xx;
   integer = floor(2*x);
   ind = (int)integer;
   if (ind<0)
      return 1.f;
   if (ind>19)
      return 1.f*(1+.1296/x);
   frac = 2*x-integer;
   return 1.f*((1-frac)*table[ind] + frac*table[ind+1])/sqrt(x+.0001f);
}

实际算法的实现，分别计算了线性频域的增益和非线性频域（Bark 域）的增益，并将通过Bark 域的增益对线性频域增益进行约束，最终得到增强的语音信号。

#参考文献
【1】I. Cohen and B. Berdugo, “Speech enhancement for non-stationary noise environments”. Signal Processing, vol. 81, no. 2, pp. 2403-2418, 2001.
【2】Y. Ephraim and D. Malah, “Speech enhancement using minimum mean-square error short-time spectral amplitude estimator”. IEEE Transactions on Acoustics, Speech and Signal Processing, vol. ASSP-32, no. 6, pp. 1109-1121, 1984.