CN104361894A - 一种基于输出的客观语音质量评估的方法 - Google Patents

Abstract

一种基于输出的客观语音质量评估的方法，其步骤如下：首先提取干净语音的非均匀线形预测倒谱系数用来对GMM-HMM模型进行训练，通过训练给干净语音建立参考模型；然后由所述参考模型和失真语音的非均匀线形预测倒谱系数向量可以得到它们之间的一致性测度；最后，通过多元非线性回归模型建立主观MOS分和一致性测度之间的映射关系，可以得到对MOS分的客观预测模型，通过所述客观预测模型进行语音质量的客观评价。本发明建立主观MOS分和客观测度之间的映射关系，得到对主观MOS分的预测模型，使得分更接近主观质量。

Description

一种基于输出的客观语音质量评估的方法

技术领域

本发明涉及语音信号处理技术领域，具体为一种基于输出的客观语音质量评估的方法。

背景技术

语音质量客观评价从评价方法上可分为侵入式(intrusive)和非侵入式(non-intrusive)。侵入式的评价通常也被称为双端或基于输入-输出的评价，它是以语音***的输入信号和输出信号之间的误差大小来判断语音质量的好坏，通过提取两端语音信号的特征参量来建立评价模型。非侵入式的评价也被称为单端或基于输出的评价，它仅根据语音***的输出信号来进行质量评价。基于输入-输出的评价，在评价语音质量时必须要有原始语音。而在实际应用中，某些场合中难以或者不便于采集到原始语音材料。因此，对该的语音材料的评估需要有一种科学客观的语音质量评估方法。

在客观语音质量评估中，选择特征参数是至关重要的，它影响着评估***的性能。由于人耳对声音的感知特性并不是以线性频率为尺度来变化的，而是用称之为巴克的临界频带的一系列有限频段来表达。故一般在特征参数提取时要通过频谱弯折来实现对人的听觉特性的模拟。如在MFCC参数时，将频谱转化为基于Mel评标的非线性频谱。带宽类似人耳听觉临界谱带的滤波器组来实现频谱弯折。而PLP参数将频域划分为22个频带，取中间20个频带。将每个频带内的能量谱与知觉驱动的加权重函数相乘得到临界听觉谱带，以此实现频谱弯折。显然在上面两种特征参数提取时，线性频率和弯曲频率间的变换计算量是比较大的。

发明内容

本申请的目的在于克服现有技术的不足，提供一种基于输出的客观语音质量评估的方法，这种方法在评估时不需要原始语音，通过结合非均匀线性预测倒谱系数和GMM-HMM的语音质量客观评估方法，建立主观MOS分和客观测度之间的映射关系，得到对主观MOS分的预测模型，从而可以用来客观评价编码或经信道传输后的语音质量。

为了达到上述目的，本发明的技术方案如下：

一种基于输出的客观语音质量评估的方法，其步骤如下：

首先提取干净语音的非均匀线形预测倒谱系数用来对GMM-HMM模型进行训练，通过训练给干净语音建立参考模型；

然后由所述参考模型和失真语音的非均匀线形预测倒谱系数向量可以得到它们之间的一致性测度；

最后，通过多元非线性回归模型建立主观MOS分和一致性测度之间的映射关系，可以得到对MOS分的客观预测模型，通过所述客观预测模型进行语音质量的客观评价。

进一步，所述非均匀线形预测倒谱系数为5阶的非均匀感知线性预测倒谱(NLPC)的特征参数，其处理步骤如下：

1、预处理，在提取NLPC参数前要进行预处理，预处理包括预加重，分帧和加窗：

(1)预加重，所述预加重用具有6dB/倍频程的提升高频特性的数字滤波器来实现，其传递函数为：

H(z)＝1-μz^-1   (1)

其中μ为预加重系数，可取为1或比1稍小的值，一般取μ＝0.95；

(2)分帧，将语音信号划分为一个一个的短时段，每一短时段称为一帧，每一帧的长度为10-30ms，所述分帧为交叠分段的方法，即每一帧的帧尾与下一帧的帧头是重叠的；

(3)加窗，设帧信号为x(n)，窗函数为w(n)，则加窗后的信号y(n)为:

y(n)＝x(n)w(n),0≤n≤N-1   (2)

其中，N为每帧的取样点数,所述窗函数是Hamming窗，即升余弦窗：

w(n)＝0.54-0.46cos[2πn/(N-1)],0≤n≤N-1   (3)

2、基于Bark双线性变换的频谱弯折，设线性频域为z平面，弯曲频域为l平面，则选择三个特殊的点(ζ₁,z₁)、(ζ₂,z₂)和(ζ₃,z₃)，通过一阶双线性变换可以确定z平面和l平面的映射关系：

$\frac{(\mathrm{\ζ}-{\mathrm{\ζ}}_1)({\mathrm{\ζ}}_2-{\mathrm{\ζ}}_3)}{({\mathrm{\ζ}}_2-{\mathrm{\ζ}}_1)(\mathrm{\ζ}-{\mathrm{\ζ}}_3)}=\frac{(z-z_1)(z_2-z_3)}{(z_2-z_1)(z-z_3)}---\left(4\right)$

考虑到线性频域和Bark谱域的特点，则取z₁＝ζ₁＝1，z₂＝ζ₂＝-1，可以得到下式：

$z=A_p\left(\mathrm{\ζ}\right)=\frac{\mathrm{\ζ}+\mathrm{\ρ}}{1+\mathrm{\ζ\ρ}}---\left(5\right)$

其中z和ω为复数，常数ρ为弯折系数，且0＜ρ＜1；显然式(7-4)为全通***，它将z平面的单位圆映射为ζ平面的单位圆；在进行频谱弯折时，z为线性频域，ζ为Bark谱域，则有z＝e^iω,ζ＝e^iθ；故弯曲频率θ可以由下式求出：

$\mathrm{\θ}=\arctan \left[\frac{(1-{\mathrm{\ρ}}^2)\sin \mathrm{\ω}}{-2\mathrm{\ρ}+(1+{\mathrm{\ρ}}^2)\cos \mathrm{\ω}}\right]---\left(\text{6}\right)$

若原始***具有截止频率为ω_p的分段恒定低通频率响应，经过双线性变换，***将同样具有类似的截止频率为θ_p的低通响应，ρ、ω_p和θ_p满足如下公式：

$\mathrm{\ρ}=\frac{\sin [({\mathrm{\θ}}_p-{\mathrm{\ω}}_p)/2]}{\sin [({\mathrm{\θ}}_p+{\mathrm{\ω}}_p)/2]}---\left(7\right)$

对ρ的优化取决于语音信号的采样频率f_s，本发明采用的ρ的优化形式：

${\mathrm{\ρ}}_Y=1.067{\left[\frac{2}{\mathrm{\π}}\arctan \left(0.6583f_s\right)\right]}^{\frac{1}{2}}-0.191---\left(8\right)$

3、非均匀线性预测倒谱的特征参数提取，非均匀线性预测倒谱特征参数的计算步骤如下：

1)通过每帧语音信号，采用p阶线性预测分析计算出声道全极点传输函数H(z)

$H\left(z\right)=\frac{G}{A\left(z\right)}=\frac{G}{1-\underset{k=1}{\overset{P}{\mathrm{\Σ}}}{a}_k{z}^{-k}}---\left(9\right)$

式中G为增益常数，a_k为第k阶线性预测参数，取p＝12；

2)对向量a(n)＝[1,-a₁,…-a_p]进行DFT变换得到A(e^iω)；在A(e^iω)上非均匀取M个点，则有

$\tilde{A}\left(k\right)=\underset{n=0}{\overset{p}{\mathrm{\Σ}}}{a}_n{e}^{-j{\mathrm{\θ}}_{k}n},k=\mathrm{0,1},...,M-1---\left(10\right)$

其中θ_k为利用双线性变换逼近的巴克频率刻度，它与原始线性频率之间的关系为：

${\mathrm{\θ}}_k=\arctan \left[\frac{e^{j2\mathrm{\πk}/M}+\mathrm{\ρ}}{1+{\mathrm{\ρe}}^{j2\mathrm{\πk}/M}}\right]---\left(11\right)$

其中ρ是由采样频率决定的频率弯折因子，考虑Bark谱域有22个频带，本发明中M取为88；

3)弯折后的LPC谱为

$\tilde{P}\left(k\right)=\frac{G^2}{{\left|\tilde{A}\left(k\right)\right|}^2}---\left(12\right)$

用q阶AR模型来逼近得到声道时域相应的自相关函数：

$r\left(n\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}\tilde{P}\left(k\right)\cos (2\mathrm{\πkn}/M),n=\mathrm{0,1},...,q---\left(\text{13}\right)$

4)根据莱文逊-德宾(Levinson-Durbin)算法，从r(n)可以求得一组新的线性预测系数n＝0,1,…,q，根据式(7-14)可以由线性预测系数得出倒谱系数c_n，该参数即为NLPC参数：

$c_1={\tilde{a}}_1$

$c_n={\tilde{a}}_n+\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},1<n\&le;p$

$c_n=\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},n>p---\left(14\right)\text{.}$

进一步，所述参考语音模型是通过采用HMM模型对特征向量进行分类，然后通过GMM对每一类特征向量建模，其方法如下：

1、基于GMM-HMM的语音建模和训练，设观测的特征向量序列为O＝o₁,o₂,…,o_T，该序列的状态模型序列为S＝s₁,s₂,…,s_N，则该序列的HMM模型可表示为：

λ＝(π,A,B)   (15)

其中，π＝{π_i＝P(s₁＝i),i＝1,2,…,N}为初始状态概率矢量；A＝{a_ij}为状态间跳转的转移概率矩阵，a_ij为从状态i跳转到状态j的概率；B＝{b_i(o_t)＝P(o_t|s_t＝i),2≤i≤N-1}为状态的输出概率分布集；

对连续HMM模型，观测序列为连续信号，其与状态j对应的信号空间用M个混合高斯密度函数的和来表示：

$b_j\left(o_t\right)=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}}),1\&le;j\&le;N---\left(16\right)$

$N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})={\left(2\mathrm{\&pi;}\right)}^{-D/2}{\left|C_{\mathrm{jk}}\right|}^{-1/2}\exp (-\frac{1}{2}{(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T{C}_{\mathrm{jk}}^{-1}(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}))---\left(17\right)$

其中，c_jk表示状态j的第k个高斯混合密度函数的系数；μ_jk是高斯密度函数的均值向量；C_jk为协方差矩阵，D为观测序列O的维数；HMM参数可以观测序列O＝o₁,o₂,…,o_T估计得到，估计的目标是使模型与训练数据的似然函数P(O|λ)最大化来估计最新的λ，即使这可以采用EM算法来实现，所述EM算法包含两部分：前向后向概率计算和HMM参数及高斯混合参数的再估计；

上述似然函数p(O|λ)的前向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_T\left(i\right)---\left(18\right)$

其中

a₁(i)＝πb_i(o₁),1≤i≤N

${\mathrm{\&alpha;}}_{t+1}\left(i\right)=\left[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_{t+1}\right),1\&le;t\&le;T-1,1\&le;j\&le;N$

似然函数p(O|λ)的后向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{b}_1\left(o_i\right){\mathrm{\&beta;}}_1\left(i\right)---\left(19\right)$

其中

β_i(T)＝1,1≤i≤N

${\mathrm{\&beta;}}_t\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}{b}_j\left(o_{i+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right),t=T-1,T-2,...,\mathrm{1,1}\&le;j\&le;N$

对给定观测序列O＝o₁,o₂,…,o_T通过重估计得到最新的λ，在此定义ξ_t(i,j)为t时刻状态为s_i且t+1时刻状态为s_j的概率，可由下式得到

$\begin{array}{c}{\mathrm{\&xi;}}_t(i,j)=P_{S|O,\mathrm{\&lambda;}}(s_t=s_i,s_{t+1}=s_j|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}---\left(20\right)$

在给定模型λ和观测序列O的条件下，状态s_i在时刻t的后验概率为

$\begin{array}{c}{\mathrm{\&gamma;}}_t\left(i\right)=P_{S|O,\mathrm{\&lambda;}}(s_t=i|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\\ =\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)\end{array}---\left(21\right)$

由此，HMM参数λ重估计如下

${\stackrel{\&OverBar;}{\mathrm{\&pi;}}}_i={\mathrm{\&gamma;}}_1\left(i\right)---\left(22\right)$

${\stackrel{\&OverBar;}{a}}_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t\left(i\right)}---\left(23\right)$

在t时刻状态j第k个高斯混合分量的参数c_jk，μ_jk和C_jk可以由下列式子重新估计：

${\stackrel{\&OverBar;}{c}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(24\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)o_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(25\right)$

${\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(26\right)$

其中，γ_t(j,k)表示在t时刻状态j第k个高斯混合分量的概率，可由下式得到

${\mathrm{\&gamma;}}_t(j,k)=\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\left[\frac{c_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jm}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}})}\right]---\left(27\right)$

2、基于多元非线性自回归分析的MOS分预测，在对参考语音建模及训练后，对失真语音的特征向量序列可以由Viterbi译码算法计算其状态模型序列 $S_t=\{s_{t_1},s_{t_2},...,s_{t_T}\},$ 使

$P\{O_t,S_t|\mathrm{\&lambda;}\}=\underset{S}{\max }P\{O_t,S\left|\mathrm{\&lambda;}\right\}---\left(28\right)$

同时也可得到第i帧的失真语音特征向量序列处于对应状态模型的概率其中S＝s₁,s₂,…,s_N，当所有失真语音特征向量序列的状态模型确定后，此段语音中所有处于状态模型s_i的语音帧数由可以计算出来，状态模型s_i下参考语音参数模型与失真语音参数向量的一致性测度的计算公式如下：

$c_i=\frac{1}{N_{s_i}}\underset{j=1}{\overset{N_{s_i}}{\mathrm{\&Sigma;}}}\log p\left(o_j\right|s_i)---\left(\text{29}\right)$

其中o_j是处于状态模型s_i的第j个失真语音参数向量，c_i表示处于状态模型s_i语音帧的一致性测度，其值越大表示失真语音与参考语音一致性越好，对每一类状态模型的失真语音帧都计算出其一致性测度，然后通过MOS映射来得到最终对失真语音MOS分的预测值；

求出每一类语音帧的一致性测度后可以用一种函数映射关系来表示一致性测度与MOS分的映射关系：

其中，f(·)为预测函数，它可以是线性或非线性回归关系，也可以是多项式拟合关系；ε为预测方差；是通过f(·)预测出的MOS分，假定MOS分映射关系f(·)为非线性关系，通过多元非线性回归模型来建立；

主观MOS分与一致性测度间的映射关系f(·)的非线性关系可用如下多元非线性回归模型描述：

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{ik}}{c}_i^k+{\mathrm{\&beta;}}_0---\left(31\right)$

其中，表示处于状态模型s_i语音帧的一致性测度c_i的k次方；令则上述多元非线性回归模型可以转化为多元线性回归模型，即

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{\&beta;}}_{i1}{c}_{i1}+{\mathrm{\&beta;}}_{i2}{c}_{i2}+...+{\mathrm{\&beta;}}_{\mathrm{iK}}{c}_{\mathrm{iK}})+{\mathrm{\&beta;}}_0---\left(32\right)$

令

F＝[f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)]   (33)

β＝[β₀,β₁₁,…,β_1K,…,β_N1,…,β_NK]^T   (34)

$C=\left[\begin{array}{cccccccc}1& c_{11}^{\left(1\right)}& ...& c_{1K}^{\left(1\right)}& ...& c_{N1}^{\left(1\right)}& ...& c_{\mathrm{NK}}^{\left(1\right)}\\ 1& c_{11}^{\left(2\right)}& ...& c_{1K}^{\left(2\right)}& ...& c_{N1}^{\left(2\right)}& ...& c_{\mathrm{NK}}^{\left(2\right)}\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ 1& c_{11}^{\left(n\right)}& ...& c_{1K}^{\left(n\right)}& ...& c_{N1}^{\left(n\right)}& ...& c_{\mathrm{NK}}^{\left(n\right)}\end{array}\right]---\left(\text{35}\right)$

其中，f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)表示1,2,…,n个观测值；矩阵C中各元素的上标代表各模型一致性测度的1,2,…,n个观测值；

根据最小二乘法原理，回归系数的估计值应满足

$\left(C^{T}C\right)\widehat{\mathrm{\&beta;}}=C^{T}F---\left(36\right)$

式中，C是n行(N×K+1)列的矩阵；N为HMM参考模型的个数；K为多项式的阶数；(N×K+1)表示多元非线性回归模型转化为多元线性回归模型的变量数；所述观测值的个数n要大于变量数，所述矩阵C是行满秩矩阵，C^T是列满秩矩阵，根据矩阵的秩的性质有

R(C)＝R(C^T)＝R(C^TC)   (37)

矩阵C^TC是满秩矩阵，即对称矩阵C^TC可逆，从而(C^TC)^-1存在，β的最小二乘估计为

$\widehat{\mathrm{\&beta;}}={\left(C^{T}C\right)}^{-1}{C}^{T}F---\left(38\right)$

对式(38)根据最小二乘法进行求解，便可求得回归系数。

本发明具有以下优点：

1)用Bark双线性变换来逼近Bark频率刻度，从而拉伸语音的低频信息和压制高频信息，得到一种称为非均匀线性预测倒谱的特征参数。它可用于语音鲁棒性分析和语音识别，抑制依赖于说话人的特征，保留语音段的语言学的质量；

2)本发明建立主观MOS分和客观测度之间的映射关系，得到对主观MOS分的预测模型，使得分更接近主观质量。

附图说明

图1为基于GMM-HMM模型的语音质量评估算法的原理框图。

具体实施方式

1.NLPC倒谱特征参数提取

为了减少计算量，本发明采用Bark双线性变换来逼近Bark频率刻度，从而拉伸语音的低频信息和压制高频信息，得到一种称为非均匀线性预测倒谱的特征参数。它可用于语音鲁棒性分析和语音识别，抑制依赖于说话人的特征，保留语音段的语言学的质量。考虑到高阶参数中包含了较多依赖于说话人的特征信息，故本发明采用5阶的非均匀感知线性预测倒谱(NLPC,Non-uniform Linear Prediction Cepstrum)的特征参数来进行语音质量分析。

1.1.预处理

在提取NLPC参数前要进行预处理，预处理包括预加重，分帧和加窗：

(1)预加重

由于语音信号的平均功率谱受声门激励和口鼻辐射的影响，高频端大约在800Hz以上按6dB/倍频程跌落，为此要进行预加重。预加重的目的是将更为有用的高频部分的频谱进行提升，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱，以便于进行频谱分析或声道参数分析。预加重一般用具有6dB/倍频程的提升高频特性的数字滤波器来实现，其传递函数为：

H(z)＝1-μz^-1   (1)

其中μ为预加重系数，可取为1或比1稍小的值，一般取μ＝0.95。

(2)分帧

由于语音信号的准平稳特性，使得只在短时段上才可视为是一个平稳过程，可以用对平稳过程的分析方法进行分析，因此需要将语音信号划分为一个一个的短时段，每一短时段称为一帧，每一帧的长度大概为10-30ms。分帧可以采用连续分段的方法，但为了使帧与帧之间平滑过渡，保持其连贯性，一般采用交叠分段的方法，即每一帧的帧尾与下一帧的帧头是重叠的。

(3)加窗

为了减小语音帧的截断效应，降低帧两端的坡度，使语音帧的两端不引起急剧变化而平滑过渡到零，就要让语音帧乘以一个窗函数。设帧信号为x(n)，窗函数为w(n)，则加窗后的信号y(n)为:

y(n)＝x(n)w(n),0≤n≤N-1   (2)

其中，N为每帧的取样点数。

常用的窗函数是Hamming窗，即升余弦窗：

w(n)＝0.54-0.46cos[2πn/(N-1)],0≤n≤N-1   (3)

1.2.基于Bark双线性变换的频谱弯折

频谱弯折的目的在于使得低频部分被扩展，而高频部分被压缩，完成了线性频率到弯曲频率的变换。若能确立线性频域到弯曲频域的映射关系，就能很方便地实现频谱弯折。设线性频域为z平面，弯曲频域为l平面，则选择三个特殊的点(ζ₁，z₁)、(ζ₂,z₂)和(ζ₃,z₃)，通过一阶双线性变换可以确定z平面和l平面的映射关系：

$\frac{(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_1)({\mathrm{\&zeta;}}_2-{\mathrm{\&zeta;}}_3)}{({\mathrm{\&zeta;}}_2-{\mathrm{\&zeta;}}_1)(\mathrm{\&zeta;}-{\mathrm{\&zeta;}}_3)}=\frac{(z-z_1)(z_2-z_3)}{(z_2-z_1)(z-z_3)}---\left(4\right)$

考虑到线性频域和Bark谱域的特点，则取z₁＝ζ₁＝1，z₂＝ζ₂＝-1。可以得到下式：

$z=A_p\left(\mathrm{\&zeta;}\right)=\frac{\mathrm{\&zeta;}+\mathrm{\&rho;}}{1+\mathrm{\&zeta;\&rho;}}---\left(5\right)$

其中z和ω为复数，常数ρ为弯折系数，且0＜ρ＜1。显然式(7-4)为全通***，它将z平面的单位圆映射为ζ平面的单位圆。在进行频谱弯折时，z为线性频域，ζ为Bark谱域，则有z＝e^iω,ζ＝e^iθ。故弯曲频率θ可以由下式求出：

$\mathrm{\&theta;}=\arctan \left[\frac{(1-{\mathrm{\&rho;}}^2)\sin \mathrm{\&omega;}}{-2\mathrm{\&rho;}+(1+{\mathrm{\&rho;}}^2)\cos \mathrm{\&omega;}}\right]---\left(\text{6}\right)$

若原始***具有截止频率为ω_p的分段恒定低通频率响应，经过双线性变换，***将同样具有类似的截止频率为θ_p的低通响应，ρ、ω_p和θ_p满足如下公式：

$\mathrm{\&rho;}=\frac{\sin [({\mathrm{\&theta;}}_p-{\mathrm{\&omega;}}_p)/2]}{\sin [({\mathrm{\&theta;}}_p+{\mathrm{\&omega;}}_p)/2]}---\left(7\right)$

对ρ的优化取决于语音信号的采样频率f_s，本发明采用的ρ的优化形式：

${\mathrm{\&rho;}}_Y=1.067{\left[\frac{2}{\mathrm{\&pi;}}\arctan \left(0.6583f_s\right)\right]}^{\frac{1}{2}}-0.191---\left(8\right)$

1.3.非均匀线性预测倒谱的特征参数提取

非均匀线性预测倒谱特征参数的计算步骤如下：

1)通过每帧语音信号，采用p阶线性预测分析计算出声道全极点传输函数H(z)

$H\left(z\right)=\frac{G}{A\left(z\right)}=\frac{G}{1-\underset{k=1}{\overset{P}{\mathrm{\&Sigma;}}}{a}_k{z}^{-k}}---\left(9\right)$

式中G为增益常数，a_k为第k阶线性预测参数。本发明中取p＝12。

2)对向量a(n)＝[1,-a₁,…-a_p]进行DFT变换得到A(e^iω)。在A(e^iω)上非均匀取M个点，则有

$\tilde{A}\left(k\right)=\underset{n=0}{\overset{p}{\mathrm{\&Sigma;}}}{a}_n{e}^{-j{\mathrm{\&theta;}}_{k}n},k=\mathrm{0,1},...,M-1---\left(10\right)$

其中θ_k为利用双线性变换逼近的巴克频率刻度，它与原始线性频率之间的关系为：

${\mathrm{\&theta;}}_k=\arctan \left[\frac{e^{j2\mathrm{\&pi;k}/M}+\mathrm{\&rho;}}{1+{\mathrm{\&rho;e}}^{j2\mathrm{\&pi;k}/M}}\right]---\left(11\right)$

其中ρ是由采样频率决定的频率弯折因子，考虑Bark谱域有22个频带，本发明中M取为88。

3)弯折后的LPC谱为

$\tilde{P}\left(k\right)=\frac{G^2}{{\left|\tilde{A}\left(k\right)\right|}^2}---\left(12\right)$

用q阶AR模型来逼近得到声道时域相应的自相关函数：

$r\left(n\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\tilde{P}\left(k\right)\cos (2\mathrm{\&pi;kn}/M),n=\mathrm{0,1},...,q---\left(\text{13}\right)$

4)根据莱文逊-德宾(Levinson-Durbin)算法，从r(n)可以求得一组新的线性预测系数n＝0,1,…,q。根据式(7-14)可以由线性预测系数得出倒谱系数c_n，该参数即为NLPC参数。

$c_1={\tilde{a}}_1$

$c_n={\tilde{a}}_n+\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},1<n\&le;p$

$c_n=\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},n>p---\left(14\right)$

2.基于GMM-HMM非侵入式语音质量评估

在基于非侵入式的语音质量评估方法中，因为没有原始语音，所以必须建立参考语音模型，通过对比失真语音特征向量与参考语音模型的距离来预测失真语音的质量。为了建立参考语音模型，本发明中采用HMM模型对特征向量进行分类，然后通过GMM对每一类特征向量建模。

2.1.基于GMM-HMM的语音建模和训练

设观测的特征向量序列为O＝o₁,o₂,…,o_T，该序列的状态模型序列为S＝s₁,s₂,…,s_N，则该序列的HMM模型可表示为：

λ＝(π,A,B)   (15)

其中，π＝{π_i＝P(s₁＝i),i＝1,2,…,N}为初始状态概率矢量；A＝{a_ij}为状态间跳转的转移概率矩阵，a_ij为从状态i跳转到状态j的概率；B＝{b_i(o_t)＝P(o_t|s_t＝i),2≤i≤N-1}为状态的输出概率分布集。

对连续HMM模型，观测序列为连续信号，其与状态j对应的信号空间用M个混合高斯密度函数的和来表示：

$b_j\left(o_t\right)=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}}),1\&le;j\&le;N---\left(16\right)$

$N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})={\left(2\mathrm{\&pi;}\right)}^{-D/2}{\left|C_{\mathrm{jk}}\right|}^{-1/2}\exp (-\frac{1}{2}{(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T{C}_{\mathrm{jk}}^{-1}(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}))---\left(17\right)$

其中，c_jk表示状态j的第k个高斯混合密度函数的系数；μ_jk是高斯密度函数的均值向量；C_jk为协方差矩阵，D为观测序列O的维数。HMM参数可以观测序列O＝o₁,o₂,…,o_T估计得到，估计的目标是使模型与训练数据的似然函数P(O|λ)最大化来估计最新的λ，即使这可以采用EM算法(又称为Baum-Welch算法)来实现。这一算法包含两部分：前向后向概率计算和HMM参数及高斯混合参数的再估计。

似然函数p(O|λ)的前向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_T\left(i\right)---\left(18\right)$

其中

a₁(i)＝πb_i(o₁),1≤i≤N

${\mathrm{\&alpha;}}_{t+1}\left(i\right)=\left[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_{t+1}\right),1\&le;t\&le;T-1,1\&le;j\&le;N$

似然函数p(O|λ)的后向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{b}_1\left(o_i\right){\mathrm{\&beta;}}_1\left(i\right)---\left(19\right)$

其中

β_i(T)＝1,1≤i≤N

${\mathrm{\&beta;}}_t\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}{b}_j\left(o_{i+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right),t=T-1,T-2,...,\mathrm{1,1}\&le;j\&le;N$

对给定观测序列O＝o₁,o₂,…,o_T通过重估计得到最新的λ。在此定义ξ_t(i,j)为t时刻状态为s_i且t+1时刻状态为s_j的概率，可由下式得到

$\begin{array}{c}{\mathrm{\&xi;}}_t(i,j)=P_{S|O,\mathrm{\&lambda;}}(s_t=s_i,s_{t+1}=s_j|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}---\left(20\right)$

在给定模型λ和观测序列O的条件下，状态s_i在时刻t的后验概率为

$\begin{array}{c}{\mathrm{\&gamma;}}_t\left(i\right)=P_{S|O,\mathrm{\&lambda;}}(s_t=i|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\\ =\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)\end{array}---\left(21\right)$

由此，HMM参数λ重估计如下

${\stackrel{\&OverBar;}{\mathrm{\&pi;}}}_i={\mathrm{\&gamma;}}_1\left(i\right)---\left(22\right)$

${\stackrel{\&OverBar;}{a}}_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t\left(i\right)}---\left(23\right)$

在t时刻状态j第k个高斯混合分量的参数c_jk，μ_jk和C_jk可以由下列式子重新估计：

${\stackrel{\&OverBar;}{c}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(24\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)o_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(25\right)$

${\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(26\right)$

其中，γ_t(j,k)表示在t时刻状态j第k个高斯混合分量的概率，可由下式得到

${\mathrm{\&gamma;}}_t(j,k)=\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\left[\frac{c_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jm}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}})}\right]---\left(27\right)$

2.2.基于多元非线性自回归分析的MOS分预测

在对参考语音建模及训练后，对失真语音的特征向量序列可以由Viterbi(维特比)译码算法计算其状态模型序列使

$P\{O_t,S_t|\mathrm{\&lambda;}\}=\underset{S}{\max }P\{O_t,S\left|\mathrm{\&lambda;}\right\}---\left(28\right)$

同时也可得到第i帧的失真语音特征向量序列处于对应状态模型的概率其中S＝s₁,s₂,…,s_N。当所有失真语音特征向量序列的状态模型确定后，此段语音中所有处于状态模型s_i的语音帧数由可以计算出来。状态模型s_i下参考语音参数模型与失真语音参数向量的一致性测度的计算公式如下：

$c_i=\frac{1}{N_{s_i}}\underset{j=1}{\overset{N_{s_i}}{\mathrm{\&Sigma;}}}\log p\left(o_j\right|s_i)---\left(\text{29}\right)$

其中o_j是处于状态模型s_i的第j个失真语音参数向量。c_i表示处于状态模型s_i语音帧的一致性测度，其值越大表示失真语音与参考语音一致性越好。对每一类状态模型的失真语音帧都计算出其一致性测度，然后通过MOS映射来得到最终对失真语音MOS分的预测值。

求出每一类语音帧的一致性测度后可以用一种函数映射关系来表示一致性测度与MOS分的映射关系：

其中，f(·)为预测函数，它可以是线性或非线性回归关系，也可以是多项式拟合关系；ε为预测方差；是通过f(·)预测出的MOS分。本发明中假定MOS分映射关系f(·)为非线性关系，通过多元非线性回归模型来建立。

主观MOS分与一致性测度间的映射关系f(·)的非线性关系可用如下多元非线性回归模型描述：

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{ik}}{c}_i^k+{\mathrm{\&beta;}}_0---\left(31\right)$

其中，表示处于状态模型s_i语音帧的一致性测度c_i的k次方。令则上述多元非线性回归模型可以转化为多元线性回归模型，即

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{\&beta;}}_{i1}{c}_{i1}+{\mathrm{\&beta;}}_{i2}{c}_{i2}+...+{\mathrm{\&beta;}}_{\mathrm{iK}}{c}_{\mathrm{iK}})+{\mathrm{\&beta;}}_0---\left(32\right)$

令

F＝[f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)]   (33)

β＝[β₀,β₁₁,…,β_1K,…,β_N1,…,β_NK]^T   (34)

$C=\left[\begin{array}{cccccccc}1& c_{11}^{\left(1\right)}& ...& c_{1K}^{\left(1\right)}& ...& c_{N1}^{\left(1\right)}& ...& c_{\mathrm{NK}}^{\left(1\right)}\\ 1& c_{11}^{\left(2\right)}& ...& c_{1K}^{\left(2\right)}& ...& c_{N1}^{\left(2\right)}& ...& c_{\mathrm{NK}}^{\left(2\right)}\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ 1& c_{11}^{\left(n\right)}& ...& c_{1K}^{\left(n\right)}& ...& c_{N1}^{\left(n\right)}& ...& c_{\mathrm{NK}}^{\left(n\right)}\end{array}\right]---\left(\text{35}\right)$

其中，f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)表示1,2,…,n个观测值；矩阵C中各元素的上标代表各模型一致性测度的1,2,…,n个观测值。

根据最小二乘法原理，回归系数的估计值应满足

$\left(C^{T}C\right)\widehat{\mathrm{\&beta;}}=C^{T}F---\left(36\right)$

式中，C是n行(N×K+1)列的矩阵；N为HMM参考模型的个数；K为多项式的阶数；(N×K+1)表示多元非线性回归模型转化为多元线性回归模型的变量数。为了获得比较准确的结果，通常观测值的个数n要大于变量数。因此，矩阵C是行满秩矩阵，C^T是列满秩矩阵。根据矩阵的秩的性质有

R(C)＝R(C^T)＝R(C^TC)   (37)

因此，矩阵C^TC是满秩矩阵，即对称矩阵C^TC可逆，从而(C^TC)^-1存在。因此β的最小二乘估计为

$\widehat{\mathrm{\&beta;}}={\left(C^{T}C\right)}^{-1}{C}^{T}F---\left(38\right)$

对式(38)根据最小二乘法进行求解，便可求得回归系数。

3.基于NLPC系数和GMM-HMM模型的非侵入式语音质量评估方法

基于NLPC系数和GMM-HMM模型的非侵入式语音质量评估方法原理框图如图1所示。具体步骤为：

1)原始语音经过通信***后得到失真语音，为了得到失真语音的MOS分估计，采用非均匀线性预测倒谱进行特征参数提取。

2)采用HMM模型对参考语音信号的特征向量进行分类，然后通过GMM对每一类特征向量可能分布建模。

3)提取参考语音的特征矢量送入参考模型，计算与每一类GMM模型的一致性测度。

以一致性量度值作为语音质量的指标，通过多元非线性函数建立主观MOS分和一致性测度之间的映射关系，与MOS分的映射得到MOS分估计值。

虽然上面的举例了一些特定实施例来说明和描述本发明，但并不意味着本发明仅局限于其中的各种细节。相反地，在等价于权利要求书的范畴和范围内可以不偏离本发明精神地在各种细节上做出各种修改。

Claims (3)

1.一种基于输出的客观语音质量评估的方法，其步骤如下：

首先提取干净语音的非均匀线形预测倒谱系数用来对GMM-HMM模型进行训练，通过训练给干净语音建立参考模型；

然后由所述参考模型和失真语音的非均匀线形预测倒谱系数向量可以得到它们之间的一致性测度；

最后，通过多元非线性回归模型建立主观MOS分和一致性测度之间的映射关系，可以得到对MOS分的客观预测模型，通过所述客观预测模型进行语音质量的客观评价。

2.根据权利要求1所述的一种基于输出的客观语音质量评估的方法，其特征在于：所述非均匀线形预测倒谱系数为5阶的非均匀感知线性预测倒谱(NLPC)的特征参数，其处理步骤如下：

一)预处理，在提取NLPC参数前要进行预处理，预处理包括预加重，分帧和加窗：

(1)预加重，所述预加重用具有6dB/倍频程的提升高频特性的数字滤波器来实现，其传递函数为：

H(z)＝1-μz^-1    (1)

其中μ为预加重系数，可取为1或比1稍小的值，一般取μ＝0.95；

(2)分帧，将语音信号划分为一个一个的短时段，每一短时段称为一帧，每一帧的长度为10-30ms，所述分帧为交叠分段的方法，即每一帧的帧尾与下一帧的帧头是重叠的；

(3)加窗，设帧信号为x(n)，窗函数为w(n)，则加窗后的信号y(n)为:

y(n)＝x(n)w(n),0≤n≤N-1    (2)

其中，N为每帧的取样点数,所述窗函数是Hamming窗，即升余弦窗：

w(n)＝0.54-0.46cos[2πn/(N-1)],0≤n≤N-1    (3)

二)基于Bark双线性变换的频谱弯折，设线性频域为z平面，弯曲频域为l平面，则选择三个特殊的点和通过一阶双线性变换可以确定z平面和l平面的映射关系：

考虑到线性频域和Bark谱域的特点，则取可以得到下式：

其中z和ω为复数，常数ρ为弯折系数，且0＜ρ＜1；显然式(7-4)为全通***，它将z平面的单位圆映射为平面的单位圆；在进行频谱弯折时，z为线性频域，为Bark谱域，则有z＝e^iω,故弯曲频率θ可以由下式求出：

$\mathrm{\&theta;}=\arctan \left[\frac{(1-{\mathrm{\&rho;}}^2)\sin \mathrm{\&omega;}}{-2\mathrm{\&rho;}+(1+{\mathrm{\&rho;}}^2)\cos \mathrm{\&omega;}}\right]---\left(6\right)$

若原始***具有截止频率为ω_p的分段恒定低通频率响应，经过双线性变换，***将同样具有类似的截止频率为θ_p的低通响应，ρ、ω_p和θ_p满足如下公式：

$\mathrm{\&rho;}=\frac{\sin [({\mathrm{\&theta;}}_p-{\mathrm{\&omega;}}_p)/2]}{\sin [({\mathrm{\&theta;}}_p+{\mathrm{\&omega;}}_p)/2]}---\left(7\right)$

对ρ的优化取决于语音信号的采样频率f_s，本方法采用的ρ的优化形式：

${\mathrm{\&rho;}}_Y=1.067{\left[\frac{2}{\mathrm{\&pi;}}\arctan \left(0.6583f_s\right)\right]}^{\frac{1}{2}}-0.191---\left(8\right)$

三)非均匀线性预测倒谱的特征参数提取，非均匀线性预测倒谱特征参数的计算步骤如下：

1)通过每帧语音信号，采用p阶线性预测分析计算出声道全极点传输函数H(z)

$H\left(z\right)=\frac{G}{A\left(z\right)}=\frac{G}{1-\underset{k=1}{\overset{P}{\mathrm{\&Sigma;}}}{a}_k{z}^{-k}}---\left(9\right)$

式中G为增益常数，a_k为第k阶线性预测参数，取p＝12；

2)对向量a(n)＝[1,-a₁,…-a_p]进行DFT变换得到A(e^iω)；在A(e^iω)上非均匀取M个点，则有

$\tilde{A}\left(k\right)=\underset{n=0}{\overset{p}{\mathrm{\&Sigma;}}}{a}_n{e}^{-j{\mathrm{\&theta;}}_{k}n}k=\mathrm{0,1},...,M-1---\left(10\right)$

其中θ_k为利用双线性变换逼近的巴克频率刻度，它与原始线性频率之间的关系为：

${\mathrm{\&theta;}}_k=\arctan \left[\frac{e^{j2\mathrm{\&pi;k}/M}+\mathrm{\&rho;}}{1+{\mathrm{\&rho;e}}^{j2\mathrm{\&pi;k}/M}}\right]---\left(11\right)$

其中ρ是由采样频率决定的频率弯折因子，考虑Bark谱域有22个频带，本发明中M取为88；

3)弯折后的LPC谱为

$\tilde{P}\left(k\right)=\frac{G^2}{{\left|\tilde{A}\left(k\right)\right|}^2}---\left(12\right)$

用q阶AR模型来逼近得到声道时域相应的自相关函数：

$r\left(n\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\tilde{P}\left(k\right)\cos (2\mathrm{\&pi;kn}/M)n=\mathrm{0,1},...,q---\left(13\right)$

4)根据莱文逊-德宾(Levinson-Durbin)算法，从r(n)可以求得一组新的线性预测系数n＝0,1,…,q，根据式(7-14)可以由线性预测系数得出倒谱系数c_n，该参数即为NLPC参数：

$c_1={\tilde{a}}_1$

$c_n={\tilde{a}}_n+\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},1<n\&le;p$

$c_n=\underset{k=1}{\overset{n-1}{\mathrm{\&Sigma;}}}(1-\frac{k}{n}){\tilde{a}}_k{c}_{n-k},n>p---\left(14\right).$

3.根据权利要求2所述的一种基于输出的客观语音质量评估的方法，其特征在于：所述参考语音模型是通过采用HMM模型对特征向量进行分类，然后通过GMM对每一类特征向量建模，其方法如下：

一)基于GMM-HMM的语音建模和训练，设观测的特征向量序列为O＝o₁,o₂,…,o_T，该序列的状态模型序列为S＝s₁,s₂,…,s_N，则该序列的HMM模型可表示为：

λ＝(π,A,B)    (15)

其中，π＝{π_i＝P(s₁＝i),i＝1,2,…,N}为初始状态概率矢量；A＝{a_ij}为状态间跳转的转移概率矩阵，a_ij为从状态i跳转到状态j的概率；B＝{b_i(o_t)＝P(o_t|s_t＝i),2≤i≤N-1}为状态的输出概率分布集；

对连续HMM模型，观测序列为连续信号，其与状态j对应的信号空间用M个混合高斯密度函数的和来表示：

$b_j\left(o_t\right)=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}}),1\&le;j\&le;N---\left(16\right)$

$N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})={\left(2\mathrm{\&pi;}\right)}^{-D/2}{\left|C_{\mathrm{jk}}\right|}^{-1/2}\exp (-\frac{1}{2}{(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T{C}_{\mathrm{jk}}^{-1}(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}))---\left(17\right)$

其中，c_jk表示状态j的第k个高斯混合密度函数的系数；μ_jk是高斯密度函数的均值向量；C_jk为协方差矩阵，D为观测序列O的维数；HMM参数可以观测序列O＝o₁,o₂,…,o_T估计得到，估计的目标是使模型与训练数据的似然函数P(O|λ)最大化来估计最新的λ，即使这可以采用EM算法来实现，所述EM算法包含两部分：前向后向概率计算和HMM参数及高斯混合参数的再估计；

上述似然函数p(O|λ)的前向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_T\left(i\right)---\left(18\right)$

其中

a₁(i)＝πb_i(o₁),1≤i≤N

${\mathrm{\&alpha;}}_{t+1}\left(i\right)=\left[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_{t+1}\right),1\&le;t\&le;T-\mathrm{1,1}\&le;j\&le;N$

似然函数p(O|λ)的后向概率计算公式如下：

$P\left(O\right|\mathrm{\&lambda;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{b}_1\left(o_i\right){\mathrm{\&beta;}}_1\left(i\right)---\left(19\right)$

其中

β_i(T)＝1,1≤i≤N

${\mathrm{\&beta;}}_t\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}{b}_j\left(o_{i+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right),t=T-1,T-2,...,\mathrm{1,1}\&le;j\&le;N$

对给定观测序列O＝o₁,o₂,…,o_T通过重估计得到最新的λ，在此定义ξ_t(i,j)为t时刻状态为s_i且t+1时刻状态为s_j的概率，可由下式得到

$\begin{array}{c}{\mathrm{\&xi;}}_t(i,j)=P_{S|O,\mathrm{\&lambda;}}(s_t=s_i,s_{t+1}=s_j|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}---\left(20\right)$

在给定模型λ和观测序列O的条件下，状态s_i在时刻t的后验概率为

$\begin{array}{c}{\mathrm{\&gamma;}}_t\left(i\right)=P_{S|O,\mathrm{\&lambda;}}(s_t=i|O,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\\ =\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)\end{array}---\left(21\right)$

由此，HMM参数λ重估计如下

${\stackrel{\&OverBar;}{\mathrm{\&pi;}}}_i={\mathrm{\&gamma;}}_1\left(i\right)---\left(22\right)$

${\stackrel{\&OverBar;}{a}}_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t\left(i\right)}---\left(23\right)$

在t时刻状态j第k个高斯混合分量的参数c_jk，μ_jk和C_jk可以由下列式子重新估计：

${\stackrel{\&OverBar;}{c}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(24\right)$

${\stackrel{\&OverBar;}{\mathrm{\&mu;}}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)o_t}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(25\right)$

${\stackrel{\&OverBar;}{C}}_{\mathrm{jk}}=\frac{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}}){(o_t-{\mathrm{\&mu;}}_{\mathrm{jk}})}^T}{\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_t(j,k)}---\left(26\right)$

其中，γ_t(j,k)表示在t时刻状态j第k个高斯混合分量的概率，可由下式得到

${\mathrm{\&gamma;}}_t(j,k)=\frac{{\mathrm{\&alpha;}}_t\left(j\right){\mathrm{\&beta;}}_t\left(j\right)}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(i\right)}\left[\frac{c_{\mathrm{jk}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jk}},C_{\mathrm{jk}})}{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{\mathrm{jm}}N(o_t,{\mathrm{\&mu;}}_{\mathrm{jm}},C_{\mathrm{jm}})}\right]---\left(27\right)$

二)基于多元非线性自回归分析的MOS分预测，在对参考语音建模及训练后，对失真语音的特征向量序列可以由Viterbi译码算法计算其状态模型序列 $S_t=\{s_{t_1},s_{t_2},...,s_{t_T}\},$ 使

$P\{O_t,S_t|\mathrm{\&lambda;}\}=\underset{S}{\max }P\{O_t,S\left|\mathrm{\&lambda;}\right\}---\left(28\right)$

同时也可得到第i帧的失真语音特征向量序列处于对应状态模型的概率其中S＝s₁,s₂,…,s_N，当所有失真语音特征向量序列的状态模型确定后，此段语音中所有处于状态模型s_i的语音帧数由可以计算出来，状态模型s_i下参考语音参数模型与失真语音参数向量的一致性测度的计算公式如下：

$c_i=\frac{1}{N_{s_i}}\underset{j=1}{\overset{N_{s_i}}{\mathrm{\&Sigma;}}}\log p\left(o_j\right|s_i)---\left(29\right)$

其中o_j是处于状态模型s_i的第j个失真语音参数向量，c_i表示处于状态模型s_i语音帧的一致性测度，其值越大表示失真语音与参考语音一致性越好，对每一类状态模型的失真语音帧都计算出其一致性测度，然后通过MOS映射来得到最终对失真语音MOS分的预测值；

求出每一类语音帧的一致性测度后可以用一种函数映射关系来表示一致性测度与MOS分的映射关系：

其中，f(·)为预测函数，它可以是线性或非线性回归关系，也可以是多项式拟合关系；ε为预测方差；是通过f(·)预测出的MOS分，假定MOS分映射关系f(·)为非线性关系，通过多元非线性回归模型来建立；

主观MOS分与一致性测度间的映射关系f(·)的非线性关系可用如下多元非线性回归模型描述：

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{\mathrm{ik}}{c}_i^k+{\mathrm{\&beta;}}_0---\left(31\right)$

其中，表示处于状态模型s_i语音帧的一致性测度c_i的k次方；令则上述多元非线性回归模型可以转化为多元线性回归模型，即

$f\left(c_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{\&beta;}}_{i1}{c}_{i1}+{\mathrm{\&beta;}}_{i2}{c}_{i2}+...+{\mathrm{\&beta;}}_{\mathrm{iK}}{c}_{\mathrm{iK}})+{\mathrm{\&beta;}}_0---\left(32\right)$

令

F＝[f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)]    (33)

β＝[β₀,β₁₁,…,β_1K,…,β_N1,…,β_NK]^T    (34)

$C=\left[\begin{array}{cccccccc}1& c_{11}^{\left(1\right)}& ...& c_{1K}^{\left(1\right)}& ...& c_{N1}^{\left(1\right)}& ...& c_{\mathrm{NK}}^{\left(1\right)}\\ 1& c_{11}^{\left(2\right)}& ...& c_{1K}^{\left(2\right)}& ...& c_{N1}^{\left(2\right)}& ...& c_{\mathrm{Nk}}^{\left(2\right)}\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ 1& c_{11}^{\left(n\right)}& ...& c_{1K}^{\left(n\right)}& ...& c_{N1}^{\left(n\right)}& ...& c_{\mathrm{Nk}}^{\left(n\right)}\end{array}\right]---\left(35\right)$

其中，f⁽¹⁾(c_i),f⁽²⁾(c_i),…,f⁽ⁿ⁾(c_i)表示1,2,…,n个观测值；矩阵C中各元素的上标代表各模型一致性测度的1,2,…,n个观测值；

根据最小二乘法原理，回归系数的估计值应满足

$\left(C^{T}C\right)\widehat{\mathrm{\&beta;}}=C^{T}F---\left(36\right)$

式中，C是n行(N×K+1)列的矩阵；N为HMM参考模型的个数；K为多项式的阶数；(N×K+1)表示多元非线性回归模型转化为多元线性回归模型的变量数；所述观测值的个数n要大于变量数，所述矩阵C是行满秩矩阵，C^T是列满秩矩阵，根据矩阵的秩的性质有

R(C)＝R(C^T)＝R(C^TC)    (37)

矩阵C^TC是满秩矩阵，即对称矩阵C^TC可逆，从而(C^TC)^-1存在，β的最小二乘估计为

$\widehat{\mathrm{\&beta;}}={\left(C^{T}C\right)}^{-1}{C}^{T}F---\left(38\right)$

对式(38)根据最小二乘法进行求解，便可求得回归系数。