一种基于高斯均值超矢量与谱聚类的录音设备聚类方法与流程

本发明涉及智能语音信号处理、模式识别与音频取证技术领域，特别是涉及一种基于高斯均值超矢量与谱聚类的录音设备聚类方法。

背景技术：

随着语音取证技术的发展，基于语音样本的录音设备辨识已取得较好效果，在司法取证中具有重要意义。录音设备采集的语音证据已成为常见证据形式之一，被大量递交到法庭或者其他执法机构，对破案具有重要作用。

在实际案例中，可能由于录音设备标签信息丢失、设备损坏、设备识别结果不确定等因素，导致录音设备的识别结果受到影响，从而降低其在法庭取证中的认可度；另一方面，当法官面对大量被递交上来的语音样本时，首要关心的可能并不是录音设备所属的类别，而是想知道哪些语音样本来自于相同的录音设备。此时，法官所要面对的问题变为：在未知任何录音设备先验信息的情况下，如何估计用于采集语音样本的录音设备个数并将相同录音设备的语音样本合并在一起。

技术实现要素：

本发明的目的在于克服现有技术的不足和缺点，提供了一种基于高斯均值超矢量与谱聚类的录音设备聚类方法：对读入的语音样本进行预加重、分帧和加窗等预处理；从每帧语音中提取梅尔频率倒谱系数mfcc特征；采用期望最大化em算法训练一个通用背景模型ubm；通过最大后验概率map算法对ubm进行自适应调整，得到每个语音样本的高斯混合模型gmm；将每个gmm的所有高斯分量的均值矢量依次拼接构成高斯均值超矢量特征；利用谱聚类算法对所有语音样本的高斯均值超矢量进行聚类，得到语音样本所包含的录音设备个数并将相同录音设备的语音样本合并在一起。

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

一种基于高斯均值超矢量与谱聚类的录音设备聚类方法，包括以下步骤：

s1、读入语音：读入由不同录音设备采集的语音样本；

s2、预处理：对读入的语音样本进行预加重、分帧和加窗处理；

s3、提取特征：从每帧语音中提取mfcc特征；

s4、训练通用背景模型：以全部语音样本的mfcc特征作为输入，采用em算法训练生成一个ubm；

s5、训练高斯混合模型：以某个语音样本的mfcc特征作为输入，采用map算法更新ubm参数，得到该语音样本对应的gmm；

s6、构造高斯均值超矢量：将某个gmm的所有高斯分量的均值矢量依次拼接，构成该gmm的高斯均值超矢量；

s7、录音设备聚类：采用谱聚类算法对全部语音样本的高斯均值超矢量进行聚类，估计录音设备个数并将相同录音设备的语音样本合并。

优选的，所述步骤s2中语音预处理包括以下步骤：

s2.1、预加重：设置数字滤波器的传递函数为h(z)＝1-αz^-1，其中α为一个系数且取值为：0.9≤α≤1，读入的语音通过该数字滤波器后实现预加重；

s2.2、分帧：设置语音帧的帧长为25毫秒、帧移为10毫秒，帧长和帧移所对应的采样点个数分别为n＝0.025×fs和s＝0.01×fs，其中fs为语音采样频率，将读入的语音切分成语音帧xt'(n)，1≤t≤t，0≤n≤n-1，其中t和n分别表示语音帧数和每帧语音的采样点数；

s2.3、加窗：窗函数ω(n)为汉明窗：

将每帧语音xt'(n)与汉明窗ω(n)相乘得到加窗后的语音xt(n)：

xt(n)＝ω(n)×xt'(n)n＝0,1,...,n-1；t＝1,2,...,t。

优选的，所述步骤s3中提取语音信号特征的步骤如下：

s3.1、对第t帧语音xt(n)做离散傅立叶变换(discretefouriertransformation,dft)得到线性频谱xt(k)：

s3.2、将上述线性频谱xt(k)通过梅尔频率滤波器组得到梅尔频谱，再进行对数运算得到对数频谱st(m)，其中梅尔频率滤波器组为若干个带通滤波器hm(k)，0≤m＜d，d为滤波器的个数，每个滤波器具有三角形滤波特性，其中心频率为f(m)，当m值较小时相邻f(m)的间隔较小，随着m的增加相邻f(m)的间隔逐渐变大，每个带通滤波器的传递函数为：

其中，0≤m＜d，f(m)定义如下：

其中，fl、fh为滤波器的最低频率和最高频率，b^-1为b的逆函数：

b^-1(b)＝700(e^b/1125-1)

因此由线性谱xt(k)到对数谱st(m)的变换为：

s3.3、采用离散余弦变换(discretecosinetransformation,dct)将上述对数频谱st(m)变换到倒谱域，得到第t帧的mfcc特征，即ct(p)：

s3.4、对每帧语音重复步骤s3.1-s3.3，得到所有t帧语音的mfcc特征，将它们按帧的顺序组合成一个mfcc特征矩阵。

优选的，所述步骤s4采用em算法训练ubm的步骤包括：

s4.1、包含g个高斯分量的ubm表示为：

θ＝{wi,μi,σi}1≤i≤g

其中wi表示第i个高斯分量的权重系数，μi表示第i个高斯分量的均值矢量，∑i表示第i个高斯分量的协方差矩阵，第t帧语音的特征矩阵ft输入θ的输出概率为：

其中bi(ft)表示第i个高斯分量的概率分布，表示为：

其中d和ft分别表示特征维数和第t帧语音的特征矩阵，所有t帧语音的特征矩阵f＝[f1,…,ft,…,ft]^t，t表示矩阵转置运算，f输入θ的输出概率为：

s4.2、求期望：计算函数q(θ,θ^*)：

q(θ,θ^*)＝e{lnp[(f,i)|θ]}

整理得：

其中θ^*表示迭代更新后的通用背景模型，ln(·)表示自然对数函数，p(i|ft,θ^*)表示更新后的ubm的第i个高斯对特征矩阵ft的输出概率，p(i|ft,θ)表示更新前的ubm的第i个高斯对特征矩阵ft的输出概率，p(i|ft,θ^*)与p(i|ft,θ)的表达式相同，只需用当前迭代后的ubm参数代替当前迭代前的ubm参数即可，根据贝叶斯公式，特征矩阵ft在第i个高斯的概率为：

其中pi(ft)和pj(ft)分别为：

s4.3、最大化：根据q(θ,θ^*)函数估计θ＝{wi,ui,σi}：

s4.4、ubm参数迭代：当似然函数值达到最大时停止迭代，即当前迭代的输出概率p(f|θ^*)相对上次迭代时的输出概率p(f|θ)增幅小于设定的阈值10^-4，则停止迭代，得到最终的模型参数：

混合权重系数：

均值矢量：

协方差矩阵：

优选的，步骤s5采用map算法更新ubm参数，得到gmm的步骤包括：

s5.1、计算gmm参数：已知ubm和某个语音样本的特征矩阵f＝[f1,…,ft,…,ft]^t，第t帧语音的特征矩阵ft在第i个高斯的概率为：

s5.2、由于录音设备之间的差别主要由各高斯分量的均值矢量描述，因此只更新ubm的均值矢量，ubm的第i个高斯分量的均值矢量更新如下：

其中和μi分别表示更新后和更新前的第i个高斯分量的均值矢量，而且

其中系数ai与特征矩阵有关，控制着均值矢量的原值与新值之间的均衡。τ是一个常数相关因子，控制着每个高斯分量的更新程度，通常取值为14；

s5.3、参数迭代：设map的迭代次数为m，则重复m次步骤s5.1和步骤s5.2，得到更新的均值矢量，m通常取值为5～10；

s5.4、将ubm的权重系数、协方差矩阵及更新后的均值矢量分别作为gmm的权重系数、协方差矩阵及均值矢量，得到gmm的参数。

优选的，步骤s6的高斯均值超矢量的提取步骤如下：

s6.1、表示第j个gmm的第i个高斯分量的均值矢量，1≤j≤j，1≤i≤g，其中j和g分别表示语音样本总个数和gmm的高斯分量个数，表示如下：

s6.2、将第j个gmm的g个高斯分量的均值矢量依次拼接，得到第j个gmm的高斯均值超矢量：

优选的，步骤s7中的谱聚类算法步骤如下：

s7.1、所有待聚类语音样本的高斯均值矢量集合u＝{u¹,u²,...,u^j}，j为语音样本总个数，根据u^j构造亲和矩阵a∈r^j×j,a的第(j，k)个元素ajk定义如下：

其中d(u^j,u^k)是高斯均值超矢量u^j与u^k之间的欧氏距离，σj是一个尺度参数，定义为第j个高斯均值超矢量u^j与其他j-1个高斯均值超矢量之间的欧氏距离矢量的方差；

s7.2、构造对角矩阵h，其第(j，k)个元素代表亲和矩阵中第j行所有元素之和，再根据对角矩阵h和a构造归一化的亲和矩阵：

s7.3、计算矩阵l的前kmax个最大的特征值及其特征矢量其中1≤k≤kmax，vk为列矢量，根据相邻特征值之间的差值估计最优类别数(即录音设备个数)k：

根据估计出来的录音设备个数k，构造矩阵v＝[v1,v2,...,vk]∈r^j×k；

s7.4、归一化矩阵v的每一行，得到矩阵y∈r^j×k,y的第(j,k)个元素yjk：

s7.5、将矩阵中的每一行当作空间r^k中的一个点，利用k均值算法(k-meansalgorithm)将这j行聚类成k类。将高斯均值超矢量u^j所对应的语音样本判为第k类(即第k个录音设备)，当且仅当矩阵y的第j行被聚类在第k类中；

s7.6、根据上述聚类结果，得到全部待聚类语音样本所对应的录音设备个数，并将相同录音设备的语音样本合并。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明采用高斯均值超矢量特征表征各语音样本对应的录音设备，比传统mfcc特征更有效刻画录音设备的特性差异。

2、本发明采用谱聚类的录音设备聚类方法是一种无监督方法，与传统的有监督分类方法相比，不需要训练复杂分类器且无须知道待聚类录音设备的任何先验知识，适用范围更广。

附图说明

图1是本发明实施例一种基于高斯均值超矢量和谱聚类的录音设备聚类方法的流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

图1是本发明的结构流程图。包括如下步骤：

1、首先读入记录有录音设备信息的语音样本。

2、对读入的语音样本进行预处理，预处理包括对语音信号的预加重、分帧和加窗等步骤。预处理具体包括以下步骤：

2.1、预加重：设置数字滤波器的传递函数为h(z)＝1-αz^-1，其中α为一个系数且取值为：0.9≤α≤1，读入的语音通过该数字滤波器后实现预加重；

2.2、分帧：设置语音帧的帧长为25毫秒、帧移为10毫秒，帧长和帧移所对应的采样点个数分别为n＝0.025×fs和s＝0.01×fs，其中fs为语音采样频率，将读入的语音切分成语音帧xt'(n)，1≤t≤t，0≤n≤n-1，其中t和n分别表示语音帧数和每帧语音的采样点数；

2.3、窗函数ω(n)为汉明窗：

将每帧语音xt'(n)与汉明窗ω(n)相乘得到加窗后的语音xt(n)：

xt(n)＝ω(n)×xt'(n)n＝0,1,...,n-1；t＝1,2,...,t。

3、提取mfcc特征，具体步骤如下：

3.1、对第t帧语音xt(n)做离散傅立叶变换(discretefouriertransformation,dft)得到线性频谱xt(k)：

3.2、将上述线性频谱xt(k)通过梅尔频率滤波器组得到梅尔频谱，再进行对数运算得到对数频谱st(m)，其中梅尔频率滤波器组为若干个带通滤波器hm(k)，0≤m＜d，d为滤波器的个数，每个滤波器具有三角形滤波特性，其中心频率为f(m)，当m值较小时相邻f(m)的间隔较小，随着m的增加相邻f(m)的间隔逐渐变大，每个带通滤波器的传递函数为：

其中，0≤m＜d，f(m)定义如下：

其中，fl、fh为滤波器的最低频率和最高频率，b^-1为b的逆函数：

b^-1(b)＝700(e^b/1125-1)，

因此由线性谱xt(k)到对数谱st(m)的变换为：

3.3、采用离散余弦变换(discretecosinetransformation,dct)将上述对数频谱st(m)变换到倒谱域，得到第t帧的mfcc特征，即ct(p)：

3.4、对每帧语音重复步骤3.1-3.3，得到所有t帧语音的mfcc，将它们按帧的顺序组合成一个mfcc矩阵。

4、训练通用背景模型(ubm)，以全部语音样本的mfcc特征作为输入，采用em算法训练生成一个ubm，具体步骤如下：

4.1、包含g个高斯分量的ubm表示为：

θ＝{wi,μi,σi}1≤i≤g

其中wi表示第i个高斯分量的权重系数，μi表示第i个高斯分量的均值矢量，∑i表示第i个高斯分量的协方差矩阵，第t帧语音的特征矩阵ft输入θ的输出概率为：

其中bi(ft)表示第i个高斯分量的概率分布，表示为：

其中d和ft分别表示特征维数和第t帧语音的特征矩阵，所有t帧语音的特征矩阵f＝[f1,…,ft,…,ft]^t，t表示矩阵转置运算，f输入θ的输出概率为：

4.2、求期望：计算函数q(θ,θ^*)：

q(θ,θ^*)＝e{lnp[(f,i)|θ]}

整理得：

其中θ^*表示迭代更新后的通用背景模型，ln(·)表示自然对数函数，p(i|ft,θ^*)表示更新后的ubm的第i个高斯对特征矩阵ft的输出概率，p(i|ft,θ)表示更新前的ubm的第i个高斯对特征矩阵ft的输出概率，p(i|ft,θ^*)与p(i|ft,θ)的表达式相同，只需用当前迭代后的ubm参数代替当前迭代前的ubm参数即可，根据贝叶斯公式，特征矩阵ft在第i个高斯的概率为：

其中pi(ft)和pj(ft)分别为：

4.3、最大化：根据q(θ,θ^*)函数估计θ＝{wi,ui,σi}：

4.4、ubm参数迭代：当似然函数值达到最大时停止迭代，即当前迭代的输出概率p(f|θ^*)相对上次迭代时的输出概率p(f|θ)增幅小于设定的阈值10^-4，则停止迭代，得到最终的模型参数：

混合权重系数：

均值矢量：

协方差矩阵：

5、训练高斯混合模型，以某个语音样本的mfcc特征作为输入，采用map算法更新ubm参数，得到该语音样本对应的gmm。得到gmm的步骤具体包括：

5.1、计算gmm参数：已知ubm和某个语音样本的特征矩阵f＝[f1,…,ft,…,ft]^t，第t帧语音的特征矩阵ft在第i个高斯的概率为：

5.2由于录音设备之间的差别主要由各高斯分量的均值矢量描述，因此只更新ubm的均值矢量，ubm的第i个高斯分量的均值矢量更新如下：

其中和μi分别表示更新后和更新前的第i个高斯分量的均值矢量，

系数ai与特征数据有关，并且控制着均值矢量在原值与新值之间的均衡。τ是一个常数相关因子，控制着每个高斯分量的更新程度，通常取值为14；

5.3、设map的迭代次数为m，则重复m次步骤1)和步骤2)，得到更新的均值矢量，m通常取值为5～10；

5.4、将ubm的权重系数、协方差矩阵及更新后的均值矢量分别作为gmm的权重系数、协方差矩阵及均值矢量，得到gmm的参数。

6、构造高斯均值超矢量，将某个gmm的所有高斯分量的均值矢量依次拼接，构成该gmm的高斯均值超矢量。高斯均值超矢量具体提取步骤如下：

6.1、表示第j个gmm的第i个高斯分量的均值矢量，1≤j≤j，1≤i≤g，其中j和g分别表示语音样本总个数和gmm的高斯分量个数。表示如下：

6.2、将第j个gmm的g个高斯分量的均值矢量依次拼接，得到第j个gmm的高斯均值超矢量：

7、录音设备聚类：采用谱聚类算法对全部语音样本的高斯均值超矢量进行聚类，估计录音设备个数并将相同录音设备的语音样本合并。谱聚类算法步骤如下：

7.1、所有待聚类语音样本的高斯均值矢量集合u＝{u¹,u²,...,u^j}，j为语音样本总个数，根据u^j构造亲和矩阵a∈r^j×j,a的第(j，k)个元素ajk定义如下：

其中d(u^j,u^k)是高斯均值超矢量u^j与u^k之间的欧氏距离，σj是一个尺度参数，定义为第j个高斯均值超矢量u^j与其他j-1个高斯均值超矢量之间的欧氏距离矢量的方差；

7.2、构造对角矩阵h，其第(j，k)个元素代表亲和矩阵中第j行所有元素之和，再根据对角矩阵h和a构造归一化的亲和矩阵：

7.3、计算矩阵l的前kmax个最大的特征值及其特征矢量其中1≤k≤kmax，vk为列矢量，根据相邻特征值之间的差值估计最优类别数(即录音设备个数)k：

根据估计出来的录音设备个数k，构造矩阵v＝[v1,v2,...,vk]∈r^j×k；

7.4、归一化矩阵v的每一行，得到矩阵y∈r^j×k,y的第(j,k)个元素yjk：

7.5、将矩阵中的每一行当作空间r^k中的一个点，利用k均值算法(k-meansalgorithm)将这j行聚类成k类。将高斯均值超矢量u^j所对应的语音样本判为第k类(即第k个录音设备)，当且仅当矩阵y的第j行被聚类在第k类中；

7.6、根据上述聚类结果，得到全部待聚类语音样本所对应的录音设备个数，并将相同录音设备的语音样本合并。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。