weixin_40725706

这个屌丝很懒，什么也没留下！

热门标签

【SLP·Python】基于 DTW GMM HMM 三种方法实现的语音分类识别系统_语音分类模型

作者：weixin_40725706 | 2024-08-02 21:46:25

踩

语音分类模型

前言

Github 项目地址：https://github.com/Sherry-XLL/Digital-Recognition-DTW_HMM_GMM

在孤立词语音识别(Isolated Word Speech Recognition) 中，DTW，GMM 和 HMM 是三种典型的方法：

动态时间规整(DTW, Dyanmic Time Warping)
高斯混合模型(GMM, Gaussian Mixed Model)
隐马尔可夫模型(HMM, Hidden Markov Model)

本文并不介绍这三种方法的基本原理，而是侧重于 Python 版代码的实现，针对一个具体的语音识别任务——10 digits recognition system，分别使用 DTW、GMM 和 HMM 建立对 0～9 十个数字的孤立词语音分类识别模型。

一、音频数据处理

1. 数据集

本文中的语音识别任务是对 0 ~ 9 这十个建立孤立词语言识别模型，所使用的数据集类似于语音版的 MNIST 手写数字数据库。原始数据集是 20 个人录制的数字 0 ~ 9 的音频文件，共 200 条 .wav 格式的录音，同样上传到了 Github 项目中，分支情况如下：

records 文件夹下包括 digit_0 ~ digit_9 十个子文件夹。每个子文件夹代表一个数字，里面有 20 个人对该数字的录音音频，1_0.wav 就代表第 1 个人录数字 0 的音频文件。

为了比较不同方法的性能，我们将数据集按照 4:1 的比例分为训练集和测试集。我们更想了解的是模型在未知数据上的表现，于是前16个人的数据都被划分为了训练集，共 160 条音频文件；17-20 这四个人的音频为测试集，共 40 条 .wav 格式的录音。

2. 音频预处理

preprocess.py referencing
https://github.com/rocketeerli/Computer-VisionandAudio-Lab/tree/master/lab1

当录数字的音频时，录音前后会有微小时间的空隙，并且每段音频的空白片段长度不一。如果忽略这个差异，直接对原音频进行建模运算，结果难免会受到影响，因此本文选择基于双阈值的语音端点检测对音频进行预处理。

对于每段 .wav 音频，首先用 wave 读取，并通过 np.frombuffer 和 readframes 得到二进制数组。然后，编写计算能量 (energy) 的函数 calEnergy 和计算过零率 (zero-crossing rate,ZCR) 的函数 calZeroCrossingRate，帧长 framesize 设置为常用的 256：

def calEnergy(wave_data):
    """
    :param wave_data: binary data of audio file
    :return: energy
    """
    energy = []
    sum = 0
    for i in range(len(wave_data)):
        sum = sum + (int(wave_data[i]) * int(wave_data[i]))
        if (i + 1) % 256 == 0:
            energy.append(sum)
            sum = 0
        elif i == len(wave_data) - 1:
            energy.append(sum)
    return energy


def calZeroCrossingRate(wave_data):
    """
    :param wave_data: binary data of audio file
    :return: ZeroCrossingRate
    """
    zeroCrossingRate = []
    sum = 0
    for i in range(len(wave_data)):
        sum = sum + np.abs(int(wave_data[i] >= 0) - int(wave_data[i - 1] >= 0))
        if (i + 1) % 256 == 0:
            zeroCrossingRate.append(float(sum) / 255)
            sum = 0
        elif i == len(wave_data) - 1:
            zeroCrossingRate.append(float(sum) / 255)
    return zeroCrossingRate
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

代入音频的二进制数据得到能量和过零率数组之后，通过语音端点检测得到处理之后的音频数据，端点检测函数 endPointDetect 主要为三个步骤：

先根据 energy 和 zeroCrossingRate 计算高能阈值 (high energy threshold, MH)
第一步是计算较高的短时能量作为高能阈值 (high energy threshold, MH)，可用于识别语音的浊音部分，得到初步检测数据 A。
第二步是计算较低的能量阈值 (low energy threshold, ML)，使用该阈值搜索两端，第二次检测从 A 得到 B。
第三步是计算过零率阈值 (zero crossing rate threshold, Zs)，在考虑过零率的基础上进一步处理第二步的结果 B，返回处理好的数据 C。

最后，将编号 1-16 被试者的音频文件存入 processed_train_records ，编号 17-20 存入 processed_test_records。

3. MFCC 特征提取

音频文件无法直接进行语音识别，需要对其进行特征提取。在语音识别（Speech Recognition）领域，最常用到的语音特征就是梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients，MFCC）。

在 Python 中，可以用封装好的工具包 librosa 或 python_speech_features 实现对 MFCC 特征的提取。本文使用 librosa 提取给定音频的 MFCC 特征，每帧提取 39 维 MFCC 特征：

import librosa
def mfcc(wav_path, delta = 2):
    """
    Read .wav files and calculate MFCC
    :param wav_path: path of audio file
    :param delta: derivative order, default order is 2
    :return: mfcc
    """
    y, sr = librosa.load(wav_path)
    # MEL frequency cepstrum coefficient
    mfcc_feat = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
    ans = [mfcc_feat]
    
    # Calculate the 1st derivative
    if delta >= 1:
        mfcc_delta1 = librosa.feature.delta(mfcc_feat, order=1, mode ='nearest')
        ans.append(mfcc_delta1)
        
    # Calculate the 2nd derivative
    if delta >= 2:
        mfcc_delta2 = librosa.feature.delta(mfcc_feat, order=2, mode ='nearest')
        ans.append(mfcc_delta2)
        
    return np.transpose(np.concatenate(ans, axis=0), [1,0])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

可以参考官网 librosa 和 python-speech-features 了解更多音频特征提取工具。

二、Isolated Speech Recognition - DTW

1. DTW 算法步骤

动态时间规整(DTW, Dyanmic Time Warping) 本质上是一种简单的动态规划算法，它可以分为两个步骤，一个是计算两种模式的帧之间的距离，即得到帧匹配距离矩阵；另一个是在帧匹配距离矩阵中找到最优路径。算法步骤如下：

首先，使用欧几里德距离 (Euclidean distance) 初始化成本矩阵 (cost matrix)。
其次，使用动态规划计算成本矩阵，将每个数据的向量与模板向量进行比较，动态转移方程为
$c o s t [i, j] + = m i n ([c o s t [i, j], c o s t [i + 1, j], c o s t [i, j + 1]]) .$
第三，计算规整路径 (warp path)，最小距离为标准化距离。将第一个训练序列设置为标准 (template)，合并所有训练序列并求其平均值，返回 0-9 位数字的模型。
测试过程中，使用 DTW 衡量每条音频到 0-9 这十个模板的”距离“，并选择模板距离最小的数字作为音频的预测值。

2. DTW 函数编写

编写 dtw 函数的关键在于两段 MFCC 序列的动态规划矩阵和规整路径的计算。令输入的两段 MFCC 序列分别为 M1 和 M2，M1_len 和 M2_len 表示各自的长度，则 cost matrix 的大小为 M1_len * M2_len，先用欧式距离对其进行初始化，再根据转移式计算 cost matrix：

def dtw(M1, M2):
    """
    Compute Dynamic Time Warping(DTW) of two mfcc sequences.
    :param M1, M2: two mfcc sequences
    :return: the minimum distance and the wrap path
    """
    # length of two sequences
    M1_len = len(M1)
    M2_len = len(M2)
    cost_0 = np.zeros((M1_len + 1, M2_len + 1))
    cost_0[0, 1:] = np.inf
    cost_0[1:, 0] = np.inf
    
    # Initialize the array size to M1_len * M2_len
    cost = cost_0[1:, 1:]
    for i in range(M1_len):
        for j in range(M2_len):
            cost[i, j] = calEuclidDist(M1[i], M2[j])
            
    # dynamic programming to calculate cost matrix
    for i in range(M1_len):
        for j in range(M2_len):
            cost[i, j] += min([cost_0[i, j], \
                               cost_0[min(i + 1, M1_len - 1), j], \
                               cost_0[i, min(j + 1, M2_len - 1)]])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

欧式距离的计算函数 calEuclidDist 可以用一行代码完成：

def calEuclidDist(A, B):
    """
    :param A, B: two vectors
    :return: the Euclidean distance of A and B
    """
    return sqrt(sum([(a - b) ** 2 for (a, b) in zip(A, B)]))
1
2
3
4
5
6

cost matrix 计算完成后还需计算 warp path，MFCC 序列长度为 1 时的路径可以单独分情况讨论，最小代价的 path_1 和 path_2 即为所求，最后返回数组 path：

# calculate the warp path
   if len(M1) == 1:
       path = np.zeros(len(M2)), range(len(M2))
   elif len(M2) == 1:
       path = range(len(M1)), np.zeros(len(M1))
   else:
       i, j = np.array(cost_0.shape) - 2
       path_1, path_2 = [i], [j]
       # path_1, path_2 with the minimum cost is what we want
       while (i > 0) or (j > 0):
           arg_min = np.argmin((cost_0[i, j], cost_0[i, j + 1], cost_0[i + 1, j]))
           if arg_min == 0:
               i -= 1
               j -= 1
           elif arg_min == 1:
               i -= 1
           else:
               j -= 1
           path_1.insert(0, i)
           path_2.insert(0, j)
       # convert to array
       path = np.array(path_1), np.array(path_2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

3. 模型训练与预测

在模型训练阶段，对于每个数字的 16 条音频文件，现计算 MFCC 序列到 mfcc_list 列表中，将第一个音频文件的 MFCC 序列设置为标准，计算标准模板和每条模板之间的动态时间规整路径，再对其求平均进行归一化，将结果 append 到 model 列表中。返回的 model 包含 0-9 这 10 个数字的 DTW 模型：

 # set the first sequence as standard, merge all training sequences
 mfcc_count = np.zeros(len(mfcc_list[0]))
 mfcc_all = np.zeros(mfcc_list[0].shape)
 for i in range(len(mfcc_list)):
     # calculate the wrap path between standard and each template
     _, path = dtw(mfcc_list[0], mfcc_list[i])
     for j in range(len(path[0])):
         mfcc_count[int(path[0][j])] += 1
         mfcc_all[int(path[0][j])] += mfcc_list[i][path[1][j]]

 # Generalization by averaging the templates
 model_digit = np.zeros(mfcc_all.shape)
 for i in range(len(mfcc_count)):
     for j in range(len(mfcc_all[i])):
         model_digit[i][j] = mfcc_all[i][j] / mfcc_count[i]
         
 # return model with models of 0-9 digits
 model.append(model_digit)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

测试阶段比较简单，将训练得到的 model 和需要预测音频的路径输入到函数 predict_dtw 中，分别计算 model 中每位数字到音频 MFCC 序列的最短距离 min_dist，min_dist 所对应的数字即为该音频的预测标签。

def predict_dtw(model, file_path):
    """
    :param model: trained model
    :param file_path: path of .wav file
    :return: digit
    """
    # Iterate, find the digit with the minimum distance
    digit = 0
    min_dist, _ = dtw(model[0], mfcc_feat)
    for i in range(1, len(model)):
        dist, _ = dtw(model[i], mfcc_feat)
        if dist < min_dist:
            digit = i
            min_dist = dist
    return str(digit)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

三、Isolated Speech Recognition - GMM

1. GMM 模型训练流程

参考 sklearn 高斯混合模型中文文档：https://www.scikitlearn.com.cn/0.21.3/20/#211

高斯混合模型(GMM, Gaussian Mixed Model) 是一个假设所有的数据点都是生成于有限个带有未知参数的高斯分布所混合的概率模型，包含了数据的协方差结构以及隐高斯模型的中心信息，同样可以用于解决本文的数字识别任务。

scikit-learn 是基于 Python 语言的机器学习工具，sklearn.mixture 是其中应用高斯混合模型进行非监督学习的包，GaussianMixture 对象实现了用来拟合高斯混合模型的期望最大化 (EM, Expectation-Maximum) 算法。

如果想要直接调用封装好的库，只需在代码中加入 from sklearn.mixture import GaussianMixture，GaussianMixture.fit 便可以从训练数据中拟合出一个高斯混合模型，并用 predict 进行预测，详情参见官方文档。

GMM 模型的训练和预测过程与 DTW 类似，此处不再赘述，流程图如下所示：
GMM 训练流程图

2. EM 算法步骤

参考文档：https://blog.csdn.net/nsh119/article/details/79584629?spm=1001.2014.3001.5501

GMM 模型的核心在于期望最大化 (EM, Expectation-Maximum) 算法，模型参数的训练步骤主要为 Expectation 的 E 步和 Maximum 的 M 步：

(1) 取参数的初始值开始迭代。
(2) E 步：依据当前模型参数，计算分模型 $k$ 对观测数据 $y_j$ 的相应度
$\hat{\gamma}_{jk}=\dfrac{\alpha_j\phi(y_j|\theta_k)}{\Sigma^K_{k=1}\alpha_k\phi(y_j|\theta_k)}，j=1,2,...,N; k=1,2,...,K$
(3) M 步：计算新一轮迭代的模型参数

\begin{aligned} {\hat{μ}}_{k} & = \frac{Σ_{j = 1}^{N} {\hat{γ}}_{j k} y_{j}}{Σ_{j = 1}^{N} {\hat{γ}}_{j k}} ， k = 1, 2, . . ., K \\ {\hat{σ}}_{k}^{2} & = \frac{Σ_{j = 1}^{N} {\hat{γ}}_{j k} (y_{j} - μ_{k})^{2}}{Σ_{j = 1}^{N} {\hat{γ}}_{j k}} ， k = 1, 2, . . ., K \\ {\hat{α}}_{k} & = \frac{Σ_{j = 1}^{N} {\hat{γ}}_{j k}}{N} ， k = 1, 2, . . ., K \end{aligned}

$\begin{aligned} \hat{\mu}_k&=\dfrac{\Sigma^N_{j=1}\hat{\gamma}_{jk}y_j}{\Sigma^N_{j=1}\hat{\gamma}_{jk}}，k=1,2,...,K \\ \hat{\sigma}^2_k&=\dfrac{\Sigma^N_{j=1}\hat{\gamma}_{jk}(y_j-\mu_k)^2}{\Sigma^N_{j=1}\hat{\gamma}_{jk}}，k=1,2,...,K \\ \hat{\alpha}_k&=\dfrac{\Sigma^N_{j=1}\hat{\gamma}_{jk}}{N}，k=1,2,...,K \end{aligned}$

\overset{μ}{^}_{k} \overset{σ}{^}_{k}^{2} \overset{α}{^}_{k} = \frac{Σ _{j = 1}^{N} γ ^ _{j k} y _{j}}{Σ _{j = 1}^{N} γ ^ _{j k}} ， k = 1, 2, . . ., K = \frac{Σ _{j = 1}^{N} γ ^ _{j k} ( y _{j} - μ _{k} ) ^{2}}{Σ _{j = 1}^{N} γ ^ _{j k}} ， k = 1, 2, . . ., K = \frac{Σ _{j = 1}^{N} γ ^ _{j k}}{N} ， k = 1, 2, . . ., K

3. GMM 实现代码

sklearn.mixture/_gaussian_mixture.py：https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/mixture/_gaussian_mixture.py

GMM 模型实现可以参考 GaussianMixture 的源码，本文对 sklearn 的实现方式进行了简化，关键在于 GMM 这个类的函数实现。

我们实现的 GMM 类包含三个参数 (parameter) 和三个属性 (attribute)：n_components 表示模型中状态的数目；mfcc_data 为音频提取的 MFCC 数据；random_state 是控制随机数的参数；means 表示在每个状态中 mixture component 的平均值；var 为方差；weights 是每个 component 的参数矩阵。

class GMM:
    """Gaussian Mixture Model.
    Parameters
    ----------
    n_components : int
        Number of states in the model.
    mfcc_data : array, shape (, 39)
        Mfcc data of training wavs.
    random_state: int
        RandomState instance or None, optional (default=0)
    Attributes
    ----------
    means : array, shape (n_components, 1)
        Mean parameters for each mixture component in each state.
    var : array, shape (n_components, 1)
        Variance parameters for each mixture component in each state.
    weights : array, shape (1, n_components)
        Weights matrix of each component.
    """
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

根据均值和方差可以计算对数高斯概率：

def log_gaussian_prob(x, means, var):
    """
    Estimate the log Gaussian probability
    :param x: input array
    :param means: The mean of each mixture component.
    :param var: The variance of each mixture component.
    :return: the log Gaussian probability
    """
    return (-0.5 * np.log(var) - np.divide(np.square(x - means), 2 * var) - 0.5 * np.log(2 * np.pi)).sum()
1
2
3
4
5
6
7
8
9

GMM 类中包含五个部分：__init__ 用于初始化；e_step 是 Expectation 步骤；m_step 是 Maximization 步骤；train 用于调用 E 步和 M 步进行参数训练；log_prob 计算每个高斯模型的对数高斯概率。

def __init__(self, mfcc_data, n_components, random_state=0):
    # Initialization
    self.mfcc_data = mfcc_data
    self.means = np.tile(np.mean(self.mfcc_data, axis=0), (n_components, 1))
    # randomization
    np.random.seed(random_state)
    for k in range(n_components):
        randn_k = np.random.randn()
        self.means[k] += 0.01 * randn_k * np.sqrt(np.var(self.mfcc_data, axis=0))
    self.var = np.tile(np.var(self.mfcc_data, axis=0), (n_components, 1))
    self.weights = np.ones(n_components) / n_components
    self.n_components = n_components

def e_step(self, x):
    """
    Expectation-step.
    :param x: input array-like data
    :return: Logarithm of the posterior probabilities (or responsibilities) of
             the point of each sample in x.
    """
    log_resp = np.zeros((x.shape[0], self.n_components))
    for i in range(x.shape[0]):
        log_resp_i = np.log(self.weights)
        for j in range(self.n_components):
            log_resp_i[j] += log_gaussian_prob(x[i], self.means[j], self.var[j])
        y = np.exp(log_resp_i - log_resp_i.max())
        log_resp[i] = y / y.sum()
    return log_resp

def m_step(self, x, log_resp):
    """
    Maximization step.
    :param x: input array-like data
    :param log_resp: Logarithm of the posterior probabilities (or responsibilities) of
                  the point of each sample in data
    """
    self.weights = np.sum(log_resp, axis=0) / np.sum(log_resp)
    denominator = np.sum(log_resp, axis=0, keepdims=True).T
    means_num = np.zeros_like(self.means)
    for k in range(self.n_components):
        means_num[k] = np.sum(np.multiply(x, np.expand_dims(log_resp[:, k], axis=1)), axis=0)
    self.means = np.divide(means_num, denominator)
    var_num = np.zeros_like(self.var)
    for k in range(self.n_components):
        var_num[k] = np.sum(np.multiply(np.square(np.subtract(x, self.means[k])),
                                        np.expand_dims(log_resp[:, k], axis=1)), axis=0)
    self.var = np.divide(var_num, denominator)

def train(self, x):
    """
    :param x: input array-like data
    """
    log_resp = self.e_step(x)
    self.m_step(x, log_resp)

def log_prob(self, x):
    """
    :param x: input array-like data
    :return: calculated log gaussian probability of single modal Gaussian
    """
    sum_prob = 0
    for i in range(x.shape[0]):
        prob_i = np.array([np.log(self.weights[j]) + log_gaussian_prob(x[i], self.means[j], self.var[j])
                           for j in range(self.n_components)])
        sum_prob += np.max(prob_i)
    return sum_prob
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

四、Isolated Speech Recognition - HMM

1. 隐马尔可夫模型

参考 hmmlearn 文档：https://hmmlearn.readthedocs.io/en/stable

隐马尔可夫模型(HMM, Hidden Markov Model) 是一种生成概率模型，其中可观测的 $X$ 变量序列由内部隐状态序列 $Z$ 生成。隐状态 (hidden states) 之间的转换假定为（一阶）马尔可夫链 (Markov chain) 的形式，可以由起始概率向量 $π$ 和转移概率矩阵 $A$ 指定。可观测变量的发射概率 (emission probability) 可以是以当前隐藏状态为条件的任何分布，参数为 $\theta$ 。HMM 完全由 $π$ 、 $A$ 和 $\theta$ 决定。

HMM 有三个基本问题：

给定模型参数和观测数据，估计最优隐状态序列 (estimate the optimal sequence of hidden states)
给定模型参数和观测数据，计算模型似然 (calculate the model likelihood)
仅给出观测数据，估计模型参数 (estimate the model parameters)

第一个和第二个问题可以通过动态规划中的维特比算法 (Viterbi algorithm) 和前向-后向算法 (Forward-Backward algorithm) 来解决，最后一个问题可以通过迭代期望最大化 (EM, Expectation-Maximization) 算法求解，也称为 Baum-Welch 算法。

如果想要直接调用封装好的库，只需在代码中加入 from hmmlearn import hmm，详情参见官方文档和源码。

2. 代码实现

对于本文的数字语音识别任务，HMM 模型的训练预测过程与 DTW 和 GMM 类似，不再重复说明。HMM 模型实现可以参考 hmmlearn 的源码，本文的实现关键在于 HMM 这个类中的函数。与 GMM 的 EM 算法不同的是，在 HMM 的实现中我们使用的是维特比算法 (Viterbi algorithm) 和前向-后向算法 (Forward-Backward algorithm)。

HMM 类包含两个参数 (parameter) 和四个属性 (attribute)：n_components 表示模型中状态的数目；mfcc_data 为音频提取的 MFCC 数据；means 表示在每个状态中 mixture component 的平均值；var 为方差；trans_mat 是状态之间的转移矩阵。

class HMM():
    """Hidden Markov Model.
    Parameters
    ----------
    n_components : int
        Number of states in the model.
    mfcc_data : array, shape (, 39)
        Mfcc data of training wavs.
    Attributes
    ----------
    init_prob : array, shape (n_components, )
        Initial probability over states.
    means : array, shape (n_components, 1)
        Mean parameters for each mixture component in each state.
    var : array, shape (n_components, 1)
        Variance parameters for each mixture component in each state.
    trans_mat : array, shape (n_components, n_components)
        Matrix of transition probabilities between states.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

HMM 类中包含五个部分：

__init__ 对参数进行初始化 (Initialization)；
viterbi Viterbi 算法用于解决 HMM 的解码问题，即找到观测序列中最可能的标记序列，本质上是一种动态规划算法 (dynamic programming)；
forward 计算所有时间步中所有状态的对数概率 (log probabilities)；
forward_backward 是参数估计的前向-后向算法 (Forward-Backward algorithm)；
train 用于调用 forward_backward 和 viterbi 进行参数训练；
log_prob 计算每个模型的对数概率 (log likelihood)。

def viterbi(self, data):
    # Viterbi algorithm is used to solve decoding problem of HMM
    # That is to find the most possible labeled sequence of the observation sequence
    for index, single_wav in enumerate(data):
        n = single_wav.shape[0]
        labeled_seq = np.zeros(n, dtype=int)
        log_delta = np.zeros((n, self.n_components))
        psi = np.zeros((n, self.n_components))
        log_delta[0] = np.log(self.init_prob)
        for i in range(self.n_components):
            log_delta[0, i] += log_gaussian_prob(single_wav[0], self.means[i], self.var[i])
        log_trans_mat = np.log(self.trans_mat)
        for t in range(1, n):
            for j in range(self.n_components):
                temp = np.zeros(self.n_components)
                for i in range(self.n_components):
                    temp[i] = log_delta[t - 1, i] + log_trans_mat[i, j] + log_gaussian_prob(single_wav[t], self.means[j], self.var[j])
                log_delta[t, j] = np.max(temp)
                psi[t, j] = np.argmax(log_delta[t - 1] + log_trans_mat[:, j])
        labeled_seq[n - 1] = np.argmax(log_delta[n - 1])
        for i in reversed(range(n - 1)):
            labeled_seq[i] = psi[i + 1, labeled_seq[i + 1]]
        self.states[index] = labeled_seq

def forward(self, data):
    # Computes forward log-probabilities of all states at all time steps.
    n = data.shape[0]
    m = self.means.shape[0]
    alpha = np.zeros((n, m))
    alpha[0] = np.log(self.init_prob) + np.array([log_gaussian_prob(data[0], self.means[j], self.var[j]) for j in range(m)])
    for i in range(1, n):
        for j in range(m):
            alpha[i, j] = log_gaussian_prob(data[i], self.means[j], self.var[j]) + np.max(
                np.log(self.trans_mat[:, j].T) + alpha[i - 1])
    return alpha

def forward_backward(self, data):
    # forward backward algorithm to estimate parameters
    gamma_0 = np.zeros(self.n_components)
    gamma_1 = np.zeros((self.n_components, data[0].shape[1]))
    gamma_2 = np.zeros((self.n_components, data[0].shape[1]))
    for index, single_wav in enumerate(data):
        n = single_wav.shape[0]
        labeled_seq = self.states[index]
        gamma = np.zeros((n, self.n_components))
        for t, j in enumerate(labeled_seq[:-1]):
            self.trans_mat[j, labeled_seq[t + 1]] += 1
            gamma[t, j] = 1
        gamma[n - 1, self.n_components - 1] = 1
        gamma_0 += np.sum(gamma, axis=0)
        for t in range(n):
            gamma_1[labeled_seq[t]] += single_wav[t]
            gamma_2[labeled_seq[t]] += np.square(single_wav[t])
    gamma_0 = np.expand_dims(gamma_0, axis=1)
    self.means = gamma_1 / gamma_0
    self.var = (gamma_2 - np.multiply(gamma_0, self.means ** 2)) / gamma_0
    for j in range(self.n_components):
        self.trans_mat[j] /= np.sum(self.trans_mat[j])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

五、实验结果

Github 项目地址：https://github.com/Sherry-XLL/Digital-Recognition-DTW_HMM_GMM

1. 实验环境

  macOS Catalina Version 10.15.6, Python 3.8, PyCharm 2020.2 (Professional Edition).
1

2. 文件简介

  dtw.py: Implementation of Dynamic Time Warping (DTW)
  gmm.py: Implementation of Gaussian Mixture Model (GMM)
  hmm.py: Implementation of Hidden Markov Model (HMM)

  gmm_from_sklearn.py: Train gmm model with GaussianMixture from sklearn
  hmm_from_hmmlearn.py: Train hmm model with hmm from hmmlearn

  preprocess.py: preprocess audios and split data
  processed_test_records: records with test audios
  processed_train_records: records with train audios
  records: original audios
  utils.py: utils function
1
2
3
4
5
6
7
8
9
10
11
12

若想运行代码，只需在命令行输入如下指令，以 DTW 为例：

  eg:
  python preprocess.py (mkdir processed records)
  python dtw.py 
1
2
3

3. 结果对比

各个方法的数据集和预处理部分完全相同，下面是运行不同文件的结果：

python dtw.py
----------Dynamic Time Warping (DTW)----------
Train num: 160, Test num: 40, Predict true num: 31
Accuracy: 0.78
1
2
3
4

python gmm_from_sklearn.py:
---------- GMM (GaussianMixture) ----------
Train num: 160, Test num: 40, Predict true num: 34
Accuracy: 0.85
1
2
3
4

python gmm.py
---------- Gaussian Mixture Model (GMM) ----------
confusion_matrix: 
 [[4 0 0 0 0 0 0 0 0 0]
 [0 4 0 0 0 0 0 0 0 0]
 [0 0 4 0 0 0 0 0 0 0]
 [0 0 0 4 0 0 0 0 0 0]
 [0 0 0 0 4 0 0 0 0 0]
 [0 0 0 0 0 4 0 0 0 0]
 [0 0 0 0 0 0 4 0 0 0]
 [0 0 0 0 0 0 0 4 0 0]
 [0 0 0 0 0 0 0 0 4 0]
 [0 0 0 0 0 0 0 0 0 4]]
Train num: 160, Test num: 40, Predict true num: 40
Accuracy: 1.00
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

python hmm_from_hmmlearn.py
---------- HMM (GaussianHMM) ----------
Train num: 160, Test num: 40, Predict true num: 36
Accuracy: 0.90
1
2
3
4

python hmm.py
---------- HMM (Hidden Markov Model) ----------
confusion_matrix: 
 [[4 0 0 0 0 0 0 0 0 0]
 [0 4 0 0 0 0 0 0 0 0]
 [0 0 3 0 1 0 0 0 0 0]
 [0 0 0 4 0 0 0 0 0 0]
 [0 0 0 0 4 0 0 0 0 0]
 [0 0 0 0 1 3 0 0 0 0]
 [0 0 0 0 0 0 3 0 1 0]
 [0 0 0 0 0 0 0 4 0 0]
 [0 0 0 1 0 0 1 0 2 0]
 [0 0 0 0 0 0 0 0 0 4]]
Train num: 160, Test num: 40, Predict true num: 35
Accuracy: 0.875
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

Method	DTW	GMM from sklearn	Our GMM	HMM from hmmlearn	Our HMM
Accuracy	0.78	0.85	1.00	0.90	0.875

值得注意的是，上面所得正确率仅供参考。我们阅读源码就会发现，不同文件中 n_components 的数目并不相同，最大迭代次数 max_iter 也会影响结果。设置不同的超参数 (hyper parameter) 可以得到不同的正确率，上表并不是各种方法的客观对比。事实上 scikit-learn 的实现更完整详细，效果也更好，文中三种方法的实现仅为基础版本。

完整代码详见笔者的 Github 项目，如有问题欢迎在评论区留言指出。

六、参考链接

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/weixin_40725706/article/detail/920074