13.4 非特定人语音识别算法——HMM

来源：互联网发布：天则经济研究知乎编辑：程序博客网时间：2024/05/21 05:22

与DTW相比，HMM一方面用隐含的状态对应于声学层各相对稳定的发音单位，并通过状态转移和状态驻留来描述发音的变化；另一方面，它引入了概率统计模型，不再用动态时间对齐的方法求匹配距离，而是用概率密度函数计算语音参数对HMM模型的输出概率，通过搜索最佳状态序列，以最大后验概率为准则找到识别结果。
HMM模型较为完整地表达了语音的声学模型，并且采用统计的训练方法将底层的声学模型和上层的语言模型融入统一的语音识别搜索算法中，可以获得较好的效果。
由于HMM所涉及的内容较多，本章只针对孤立词识别，以连续高斯混合HMM模型为例，对HMM模型结构、模型参数的训练及识别算法进行介绍。
首先以一天的气温变化为例，简单介绍什么是HMM模型。
1. 对于空气的温度，一人的感觉为标准，可以比较主观地划分为冷、凉爽、温暖、热、酷热等几类，这里称其为五个感觉状态，温度是可以用温度计定量地测定的，然而温度的数值与这五个状态却不是一一对应的因为人的感觉是不能准确地量化的，比如不能简单地将20℃强行称为凉爽或是温暖，但是20℃肯定不是冷或酷热。对于这种温度到感觉之间的对应关系，可以用模糊数学中的隶属度来描述。但是在概率模型中，则是用统计的方法，用概率密度函数来描述。比如，20℃对“凉爽”的输出概率可能为0.7，而对“温暖”的输出概率可能为0.3。
2. 一天中的气温是在随时变化的，这一变化是随机的，但是也有一个基本的规律。比如一个晴天夏日为例，凌晨的温度比较低，到了上午达到20℃左右，而到了正午则可能高达35℃，傍晚温度又降到20℃左右，深夜则又降得更低。假设每隔10分钟测量一次温度，那么一天24小时将得到60/10*24=144个温度采样，把它们按时间排起来，就可以看到一天温度变化的规律。
3. 温度的变化是随机的，它的数值随着时间的变化本身就是一个随机过程。温度给人的感觉所对应的五种状态随时间的变化则是另一个随机过程。可以注意到，温度的变化有时是稳定的，可以认为它驻留于某个状态中；也可能忽冷忽热，这时说它发生了从一个状态到另一个状态的转移。由于当前状态只与前面一时刻所处的状态有关，所以这一随机过程为马尔科夫过程。在温度变化的随机过程中，可以用温度计观测到温度的数值，而在状态变化的随机过程中，却无法准确地指出某一时刻到底是什么状态，也就是说，状态的变化是隐含的，这就是所谓“隐含”的来历。
4. 总之，HMM的精髓就在于：观察可测，状态隐含。这里观察就是温度，状态就是人的感觉。
在语音识别中，所谓的观察序列就是通过计算得到的一帧帧的语音参数，如MFCC参数。而状态则是在训练阶段事先规定好的不同语音单元。对于汉语普通话来说，语音单元可以是一个完整的音节，也可以是声母或韵母，还可以是更为精细的音素。
一个HMM模型由若干个状态组成，随着时间的变化，各个状态之间可以发生转移，也可以在一个状态内驻留。每个观察向量对不同的状态都有相应的输出概率。如上图所示的HMM，包含有四个状态S1~S4，状态之间或状态自身的转移概率用aij表示，输入观察序列为o1，o2，...，oT。每个观察序列是一帧MFCC参数。在这个模型中，序列o1，o2，...，oT是可观测的输入序列，称为观察序列，而每一时刻所处的状态却是隐含的。
一个连续混合高斯HMM的基本元素组合：
N——模型的状态数；
A={aij}——状态转移概率矩阵；
π={πi}——各状态的起始概率分布
B={bj(o)}——输出概率密度函数
输出概率密度函数中参数描述如下表所示

权系数满足下面条件：
这种连续混合高斯HMM通常简称为CHMM。对于每一个状态，都用若干个正态高斯概率密度函数（简称为pdf）的线性组合来表示，每个pdf有各自的均值矢量和协方差矩阵，这些都是通过对大量的MFCC参数进行统计得到的。
对于HMM模型，有三个基本问题需要解决：
1）输出概率的计算问题：给定观察序列O=(o1,o2,...,oT)和HMM模型λ=(Α,Β,π)，如果已知状态转移序列q=(q1,q2,...,qT)，则HMM模型λ以状态序列q输出观察序列O概率为：

HMM模型λ输出序列q的概率为：

而这里需要的是对所有可能的状态转移序列q，模型输出观察序列O的概率P(O|λ)。
由全概率公式可得：

该式需要进行2TN^T次计算，这在实际中是无法承受的。为了降低计算复杂度，可以采用前向算法和后向算法。
首先定义HMM的前向概率为：

表示给定HMM模型参数λ，部分观察序列{o1o2...ot}在t时刻处于状态i的概率。
前向概率αt(i)可以用下面的递推公式计算：
（1）初始化

（2）迭代计算

（3）终止计算
与前向概率相对应，还有后向概率，定义后向概率为：

表示给定HMM模型参数λ，观察序列在t时刻处于状态i，系统输出部分观察序列{ot+1ot+2...oT}的概率。
后向概率βt(i)也有类似的递推公式计算：
（1）初始化

（2）迭代计算
前向概率和后向概率的递推关系可用下图说明：
利用前向概率和后向概率计算输出概率
前向概率公式和后向概率公式巧妙地将整个观察序列对HMM模型的输出概率分成两个部分观察序列的输出概率的乘积，而且它们各自都有相应的递推公式，可以大大简化计算。经过分析，可以得到下面的输出概率计算公式：

实际上，这就是HMM三个基本问题中第一个问题的解答，它的另一种常用的形式是：

实际计算中首先计算出对于每个t和每个状态i的前向概率和后向概率，然后套用上面的公式，计算出该观察序列对模型的输出概率。这两个公式也成为全概率公式。
Viterbi算法是一种广泛用于通信领域中的动态规划算法，该算法在语音识别中也得到了应用。利用全概率公式，可以计算出系统的输出概率，但是无法找到一条最佳的状态转移路径。而用Viterbi算法，不仅可以找到一条足够好的状态转移路径，还可以得到该路径所对应的输出概率。同时，用Viterbi算法计算输出概率所需要的计算量要比全概率公式的计算量小很多。应该注意，这里是说”足够好“，而不是说”最优“，因为动态规划算法得到的结果通常是满意的，但并不保证它是最优的。
Viterbi算法也有递推的形式：
（1）初始化

（2）迭代计算

（3）终止计算

（4）回溯最佳路径

这里，δt(i)为t时刻第i状态的累计输出概率，ψt(i)为t时刻第i状态的前续状态号，为最优状态序列中t时刻所处的状态，为最终的输出概率。
实际上，对于最常使用的状态转移只能局限于自身和下一个状态的HMM模型来说，在第二步计算中每次可能的路径选择只有两个，与全概率算法相比并没有减少多少计算量。实际使用中，通常用对数形式的Viterbi算法，这样将避免进行大量的乘法计算，真正地减少了计算量，同时还可以保证有很高的动态范围，而不会由于过多的连乘而导致溢出问题。
对数形式的Viterbi算法如下：
（1）预处理

（2）初始化

（3）递推

（4）终止

（5）回溯最佳路径

在识别阶段，如果HMM模型为整词模型，就没有必要保存前续节点矩阵和状态转移路径，可以进一步减少计算量。
Baum-Welch算法实际上是极大似然（ML）准则的一个应用，它采用了一种多次迭代的优化算法。用Lagrange数乘法构造一个目标优化函数Q，其中包含了所有HMM参数作为变量，然后令Q对各变量的偏导数为0，推导出Q达到极点时新的HMM参数对应于旧的模型参数之间的关系，从而得到HMM各参数的估计。用新旧HMM模型参数之间的函数关系反复迭代运算，知道HMM模型参数不再发生明显的变化为止。
这里不再详细讨论具体的公式推导过程，只考虑具体实现中的一个具体问题，也就是在计算过程中经常遇到的溢出问题。还记得前向概率的定义：

将其分解就可以看到，它包含很多项形如：

的乘积项的和，其中每个乘积项都是小于1的。由于多项连乘，导致随着t的增加而越来越小，当t达到一定程度，的数值动态范围将超过计算机所能表示的范围，造成下溢。类似的，后向概率也有这个问题。前面的迭代计算公式只是能够保证快速地计算出和，但是并不能保证计算结果不溢出。为了解决这个问题，在实际运算中通常都要采用动态标定的办法。
引入标定系数c用表示标定后的α，用作为一个中间变量，可以得到下面的标定过程。
（1）当t=1时，计算，令。然后计算。

（2）对于2<=t<=T，计算，以及标定后的前向概率：

标定前后的前向概率之间的关系为：

后向概率标定过程类似，所不同的是它采用了与前向概率相同的标定系数，引入两个心的参数，使得：

也就是：

在后面将看到，标定后的前向概率和后向概率巧妙地通过分式中分子和分母的约分，即保持了与未标定前计算公式结果的一致性，又能够避免计算过程中出现溢出问题。
经典的Baum-Welch算法中，参数重估公式是在假设只有一个观察序列的条件下推导出来的。而在实际应用中，都是有大量观察序列参与训练的，也就是对于每个HMM模型，都会收集大量的语音数据，分别计算出各自的MFCC参数序列，再用来对该HMM的参数进行重估。例如，对数字”0“建立HMM模型，就要找很多人，每人录制多个”0“的wav文件，进行端点检测之后，计算出MFCC参数序列，也就是所谓的观察序列，就可以对模型的参数进行训练了。
设有K个观察序列构成的集合为：

其中为第k个观察序列，其时间长度为Tk：

定义过渡概率为观察序列在时刻t处于状态i，在时刻t+1处于状态j的概率，即：

利用标定后的前向概率和后向概率，可以为每个观察序列分别计算过渡概率，而不会有溢出问题：

利用这一公式，可以很快计算出HMM模型的转移概率：

而对于常用的从左到右单向结构的HMM来说，起始概率应该有下面默认的数值：

因此往往在实际计算中不对它进行重估。
再定义为某个观察序列在t时刻处于状态j的对于第l个混合高斯元输出概率，姑且成为混合输出概率，即：

在但观察序列的情况下，由于上式的分母是恒定的，在后续的计算中，将在分式的分母和分子中被约去，可以省略。但是再多观察序列中，由于不同观察序列对应的分母是不同的，因此不能简单地略去。因此可以得到下面的高斯概率密度函数的重估公式：
HMM参数的存储格式：将HMM的所有参数都保存在一个结构中，如：

这里hmm.N为状态数；hmm.M为一N*1数组，对应于每个状态所包含的高斯混合的个数；hmm.init为各状态的初始概率，一般来说，只有第一个状态的初始概率为1，其他为0：

hmm.trans为转移概率矩阵，通常只有上的数据有效，如：

最后，hmm.mix是一个结构，其中保存了各个HMM状态对应的高斯概率密度函数的均值和方差：

如，该HMM第一个状态为：

表示该状态有3个混合高斯概率密度函数，均值mean为3*24的矩阵，每一行是一个pdf的均值。方差var也是一个3*24的矩阵，每一行是一个pdf的方差。注意，这里用的是对角阵形式，没有用协方差矩阵。权重weight包含了三个pdf的权系数，它们的和为1，如：

有了以上这些信息，就足以描述一个HMM模型了。
高斯混合的输出概率的计算
对于单个的高斯密度函数pdf，假设其协方差矩阵为对角阵，函数pdf.m计算给定观察向量对该pdf的输出概率：
*****************************************************************
function p = pdf(m, v, x)
%计算多元高斯密度函数
%输入:
% m -- 均值向量, SIZE*1
% v -- 方差向量, SIZE*1
% x -- 输入向量, SIZE*1
%输出:
% p -- 输出概率

p = (2 * pi * prod(v)) ^ -0.5 * exp(-0.5 * (x-m) ./ v * (x-m)');
*****************************************************************
而计算观察向量x对于某个HMM状态的输出概率，也就是x对该状态的若干高斯混合元的输出概率的线性组合，用下面的代码mixture.m来实现。
***********************************************
function prob = mixture(mix, x)
%计算输出概率
%输入:
% mix -- 混合高斯结构
% x    -- 输入向量, SIZE*1
%输出:
% prob -- 输出概率

prob = 0;
for j = 1:mix.M
    m = mix.mean(j,:);
    v = mix.var (j,:);
    w = mix.weight(j);
    prob = prob + w * pdf(m, v, x);
end

% 加上realmin, 以防止viterbi.m中计算log(prob)时溢出
if prob==0, prob=realmin; end
***********************************************
前向概率、后向概率及其他参数的计算
首先，给定了一个或多个观察序列，也就是语音参数的序列，要计算每个观察序列的前向概率、后向概率，以及各种中间参数。计算前向概率和后向概率要用标定后的公式。她们都是T*N的矩阵，分别保存在alpha和beta。另外还要保留计算过程中用到的标定系数数组c，这在Baum-Welch训练算法中会用到。
接下来计算过渡概率和混合输出概率，它们都是三维数组，分别保存在ksai和gama。
计算完毕，将c、alpha、beta、ksai和gama组装到结构param中，作为输出变量返回。下面就是计算各种参数的MATLAB程序getparam.m。
*********************************************************************
function param = getparam(hmm, O)
%给定输出序列O, 计算前向概率alpha, 后向概率beta, 标定系数c, 及ksai,gama
%输入:
% hmm -- HMM模型参数
% O   -- n*d 观察序列
%输出:
% param -- 包含各种参数的结构

T = size(O,1);  %序列的长度

init = hmm.init;   %初始概率
trans = hmm.trans;  %转移概率
mix   = hmm.mix;    %高斯混合
N     = hmm.N;      %HMM状态数

% 给定观察序列O, 计算前向概率alpha
alpha = zeros(T,N);

% t=1的前向概率
x = O(1,:);
for i = 1:N
    alpha(1,i) = init(i) * mixture(mix(i),x);
end

% 标定t=1的前向概率
c    = zeros(T,1);
c(1) = 1/sum(alpha(1,:));
alpha(1,:) = c(1) * alpha(1,:);

% t=2:T的前向概率和标定
for t = 2:T
    for i = 1:N
        temp = 0;
        for j = 1:N
            temp = temp + alpha(t-1,j) * trans(j,i);
        end
        alpha(t,i) = temp * mixture(mix(i),O(t,:));
    end
    c(t) = 1/sum(alpha(t,:));
    alpha(t,:) = c(t)*alpha(t,:);
end

% 给定观察序列O, 计算后向概率beta
beta = zeros(T,N);

% t=T的后向概率及标定
for l = 1:N
    beta(T,l) = c(T);
end

% t=T-1:1的后向概率和标定
for t = T-1:-1:1
    x = O(t+1,:);
    for i = 1:N
    for j = 1:N
        beta(t,i) = beta(t,i) + beta(t+1,j) * mixture(mix(j),x) * trans(i,j);
    end
    end
    beta(t,:) = c(t) * beta(t,:);
end

%过渡概率ksai
ksai = zeros(T-1,N,N);
for t = 1:T-1
    denom = sum(alpha(t,:).*beta(t,:));
    for i = 1:N-1
    for j = i:i+1
        nom = alpha(t,i) * trans(i,j) * mixture(mix(j),O(t+1,:)) * beta(t+1,j);
        ksai(t,i,j) = c(t) * nom/denom;
    end
    end
end

%混合输出概率:gama
gama = zeros(T,N,max(hmm.M));
for t = 1:T
    pab = zeros(N,1);
    for l = 1:N
        pab(l) = alpha(t,l) * beta(t,l);
    end
    x = O(t,:);
    for l = 1:N
        prob = zeros(mix(l).M,1);
        for j = 1:mix(l).M
            m = mix(l).mean(j,:);
            v = mix(l).var (j,:);
            prob(j) = mix(l).weight(j) * pdf(m, v, x);
        end
        tmp = pab(l)/sum(pab);
        for j = 1:mix(l).M
            gama(t,l,j) = tmp * prob(j)/sum(prob);
        end
    end
end

param.c     = c;
param.alpha = alpha;
param.beta = beta;
param.ksai = ksai;
param.gama = gama;
*********************************************************************
识别算法的实现
HMM的识别一般都是用Viterbi算法，这里用到了对数的形式。识别程序要求输入一个HMM模型的参数，以及一个测试用的语音观察序列，然后计算出它对该模型的输出概率，并给出最佳的状态路径。
在Viterbi算法中，对起始概率和转移概率有一个取对数的操作。为了防止出现对0取对数而造成下溢，程序中采用了下面的办法：
*******************************
ind1 = find(trans>0);
ind0 = find(trans<=0);
trans(ind0) = -inf;
trans(ind1) = log(trans(ind1));
*******************************
这段代码用来计算转移概率的对数。首先find函数分别找到大于0和小于等于0的元素的下标，然后用下标分别给trans矩阵的不同元素赋值，这样就不会出现对0取对数的警告错误了。
完整的Viterbi算法代码如下：
*********************************************************************
function [prob,q] = viterbi(hmm, O)
%Viterbi算法
%输入:
% hmm -- hmm模型
% O   -- 输入观察序列, N*D, N为帧数,D为向量维数
%输出:
% prob -- 输出概率
% q    -- 状态序列

init = hmm.init;   %初始概率
trans = hmm.trans;  %转移概率
mix   = hmm.mix;    %高斯混合
N     = hmm.N;      %HMM状态数
T     = size(O,1);  %语音帧数

% 计算log(init);
ind1 = find(init>0);
ind0 = find(init<=0);
init(ind0) = -inf;
init(ind1) = log(init(ind1));

% 计算log(trans);
ind1 = find(trans>0);
ind0 = find(trans<=0);
trans(ind0) = -inf;
trans(ind1) = log(trans(ind1));

% 初始化
delta = zeros(T,N);
fai   = zeros(T,N);
q     = zeros(T,1);

% t=1
x = O(1,:);
for i = 1:N
    delta(1,i) = init(i) + log(mixture(mix(i),x));
end

% t=2:T
for t = 2:T
for j = 1:N
    [delta(t,j) fai(t,j)] = max(delta(t-1,:) + trans(:,j)');
    x = O(t,:);
    delta(t,j) = delta(t,j) + log(mixture(mix(j),x));
end
end

% 最终概率和最后节点
[prob q(T)] = max(delta(T,:));

% 回溯最佳状态路径
for t=T-1:-1:1
    q(t) = fai(t+1,q(t+1));
end
*********************************************************************
函数返回为两个输出变量，prob为输出概率，q为状态路径。
训练过程的Baum-Welch算法
利用Baum-Welch算法，按照公式计算转移概率，用公式计算各pdf的均值、方差和权系数。函数baum.m就是进行一次迭代的Baum-Welch算法的MATLAB实现，首先是调用getparam函数，为每个输入观察序列计算各种参数，然后再对HMM的各种参数进行重估，将新的参数作为输出变量返回。
*********************************************************************
function hmm = baum(hmm, samples)

mix = hmm.mix;         %高斯混合
N    = length(mix);     %HMM状态数
K    = length(samples); %语音样本数
SIZE = size(samples(1).data,2); %参数阶数

% 计算前向、后向概率矩阵，考虑多观察序列和下溢问题
disp( '计算样本参数...');
for k = 1:K
    fprintf( '%d ',k)
    param(k) = getparam(hmm, samples(k).data);
end
fprintf( '\n')

% 重估转移概率矩阵A: trans
disp( '重估转移概率矩阵A...')
for i = 1:N-1
    denom = 0;
    for k = 1:K
        tmp   = param(k).ksai(:,i,:);
        denom = denom + sum(tmp(:));
    end
    for j = i:i+1
        nom = 0;
        for k = 1:K
            tmp = param(k).ksai(:,i,j);
            nom = nom   + sum(tmp(:));
        end
        hmm.trans(i,j) = nom / denom;
    end
end

% 重估混合高斯的参数
disp( '重估混合高斯的参数...')
for l = 1:N
for j = 1:hmm.M(l)
    fprintf( '%d,%d ',l,j)
    % 计算各pdf的均值和方差
    nommean = zeros(1,SIZE);
    nomvar = zeros(1,SIZE);
    denom   = 0;
    for k = 1:K
        T = size(samples(k).data,1);
        for t = 1:T
            x       = samples(k).data(t,:);
            nommean = nommean + param(k).gama(t,l,j) * x;
            nomvar = nomvar + param(k).gama(t,l,j) * (x-mix(l).mean(j,:)).^2;
            denom   = denom   + param(k).gama(t,l,j);
        end
    end
    hmm.mix(l).mean(j,:) = nommean / denom;
    hmm.mix(l).var (j,:) = nomvar / denom;

    % 计算各pdf的权
    nom   = 0;
    denom = 0;
    for k = 1:K
        tmp = param(k).gama(:,l,j);    nom   = nom   + sum(tmp(:));
        tmp = param(k).gama(:,l,:);    denom = denom + sum(tmp(:));
    end
    hmm.mix(l).weight(j) = nom/denom;
end
fprintf( '\n')
end
*********************************************************************
HMM参数的初始化
在进行Baum-Welch迭代前，首先要对HMM的参数进行初始化。函数inithmm.m就是完成这一工作的。初始化的过程中，首先将初始概率hmm.init赋值为只有第一个元素为1，其他元素为0的数组，其长度为该HMM的状态数N。接着将转移概率赋值为如下形式的矩阵：

最后是对各状态的混合高斯函数的均值、方差和权系数进行初始化。这里采用了按状态数对观察序列的参数平均分段，然后将所有观察序列中属于一个段的参数组成一个大的矩阵，调用voicebox工具箱中的聚类函数kmeans，由聚类结果，计算得到各pdf的均值、方差和权系数。
*********************************************************************
function hmm = inithmm(samples, M)

K = length(samples);    %语音样本数
N = length(M);          %状态数
hmm.N = N;
hmm.M = M;

% 初始概率矩阵
hmm.init    = zeros(N,1);
hmm.init(1) = 1;

% 转移概率矩阵
hmm.trans=zeros(N,N);
for i=1:N-1
    hmm.trans(i,i)   = 0.5;
    hmm.trans(i,i+1) = 0.5;
end
hmm.trans(N,N) = 1;

% 概率密度函数的初始聚类
% 平均分段
for k = 1:K
    T = size(samples(k).data,1);
    samples(k).segment=floor([1:T/N:T T+1]);
end

%对属于每个状态的向量进行K均值聚类，得到连续混合正态分布
for i = 1:N
    %把相同聚类和相同状态的向量组合到一个向量中
    vector = [];
    for k = 1:K
        seg1 = samples(k).segment(i);
        seg2 = samples(k).segment(i+1)-1;
        vector = [vector ; samples(k).data(seg1:seg2,:)];
    end
    mix(i) = getmix(vector, M(i));
end

hmm.mix = mix;

function mix = getmix(vector, M)

[mean esq nn] = kmeans(vector,M);

% 计算每个聚类的标准差，对角阵，只保存在对角线上的元素
for j = 1:M
    ind = find(j==nn);
    tmp = vector(ind,:);
    var(j,:) = std(tmp);
end

% 计算每个聚类的元素数，归一化为各pdf的权重
weight = zeros(M,1);
for j = 1:M
    weight(j) = sum(find(j==nn));
end
weight = weight/sum(weight);

% 保存结果
mix.M      = M;
mix.mean   = mean;      % M*SIZE
mix.var    = var.^2;    % M*SIZE
mix.weight = weight;    % M*1
*********************************************************************
训练主程序
上面的函数baum.m只是实现了一次迭代，实际应用中应该是多次迭代才能得到结果，同时还应该给出一个结束迭代的条件。本例中是计算所有观察序列的输出概率，对其进行累加，得到总和输出概率。当此概率的相对变化小到一定数值的时候，结束迭代。另外，设定最大迭代次数为常数，当迭代次数超过此常数时，也停止迭代。
下面是完整的训练程序train.m:
*********************************************************************
function [hmm, pout] = train(samples, M)
%输入:
% samples -- 样本结构
% M       -- 为每个状态指定pdf个数，如:[3 3 3 3]
%输出:
% hmm      -- 训练完成后的hmm

K   = length(samples);

% 计算语音参数
disp('正在计算语音参数');
for k = 1:K
    if isfield(samples(k), 'data') & ~isempty(samples(k).data)
        continue;
    else
        samples(k).data = mfcc(samples(k).wave);
    end
end

hmm = inithmm(samples, M);

for loop = 1:40
    fprintf('\n第%d遍训练\n\n',loop)
    hmm = baum(hmm, samples);

    %计算总输出概率
    pout(loop)=0;
    for k = 1:K
        pout(loop) = pout(loop) + viterbi(hmm, samples(k).data);
    end

    fprintf('总和输出概率(log)=%d\n', pout(loop))

    %比较两个HMM的距离
    if loop>1
        if abs((pout(loop)-pout(loop-1))/pout(loop)) < 5e-6
            fprintf( '收敛!\n');
            return
        end
    end
end

disp('迭代40次仍不收敛, 退出');
*********************************************************************
程序中，输入参数有两个，结构数组samples内包含了观察序列的信息，每个samples(k)都包含两个成员samples(k).wave和samples(k).data，分别为该观察序列的原始语音和参数。其中成员data可以在调用之前计算，也可以由train程序内计算。数组M包含了个状态对应的高斯混合数。
训练程序
利用训练函数train.m和识别函数viterbi.m，就可以对观察序列进行训练和识别了。下面是一段选连程序的脚本文件mian.m，假设原始语音信号已经保存在cell数组samples中，其中，samples{i}{1:K}表示第i个词的K个原始语音信号。循环程序中，将它们导入结构数组samples数组中，再交由训练函数train进行训练。
下面是训练程序main.m的代码：
*********************************************************************
for i=1:length(samples)
    sample=[];
    for k=1:length(samples{i})
        sample(k).wave=samples{i}{k};
        sample(k).data=[];
    end
    hmm{i}=train(sample,[3 3 3 3]);
end
*********************************************************************
识别结果hmm为一个cell数组，每个元素为一个hmm结构。下面是一段识别程序，识别所用的数据是wav文件，首先用函数wavread将其读入，再用函数vad进行端点检测，计算出MFCC参数之后，交由识别函数viterbi.m计算得到其对数形式的输出概率，最后用max函数找到识别结果。
下面是识别程序recog.m的代码：
*********************************************************************
for i=1:10
    fname = sprintf('..\\..\\ch6\\%db.wav',i-1);
    %fname = sprintf('tiger2.wav',i-1);
    x = wavread(fname);
    [x1 x2] = vad(x);
    m = mfcc(x);
    m = m(x1-2:x2-2,:);
    for j=1:10
        pout(j) = viterbi(hmm{j}, m);
    end
    [d,n] = max(pout);
    fprintf('第%d个词, 识别为%d\n', i-1,n)
end
*********************************************************************
利用录音工具，作者为汉语数码0~9录制了语音，每个数码录制了10个样本，保存在cell数组samples中，并保存为文件samples.mat，以备训练程序使用。

0 0