2018-11-27 23:31:21 entoon 阅读数 6592
• ###### 自然语言处理——实战分词

742 人正在学习 去看看 魏俊

matlab程序设计实例——语音识别

%设置模板
ncoeff=12; %mfcc系数的个数
fMatrix1 = cell(1,5);
fMatrix2 = cell(1,5);
fMatrix3 = cell(1,5);
fMatrix4 = cell(1,5);
fMatrix5 = cell(1,5);
fMatrix6 = cell(1,5);
fMatrix7 = cell(1,5);

for i = 1:5
q = [‘SpeechData\1’ num2str(i) ‘.wav’];
z=speechIn1(:,1);

fMatrix1(1,i) = {mfccf(ncoeff,speechIn1,FS1)};


end

for j = 1:5
q = [‘SpeechData\2’ num2str(j) ‘.wav’];
z=speechIn2(:,1);
fMatrix2(1,j) = {mfccf(ncoeff,speechIn2,FS2)};
end

for k = 1:5
q = [‘SpeechData\3’ num2str(k) ‘.wav’];
z=speechIn3(:,1);
fMatrix3(1,k) = {mfccf(ncoeff,speechIn3,FS3)};
end

for k = 1:5
q = [‘SpeechData\4’ num2str(k) ‘.wav’];
z=speechIn4(:,1);
fMatrix4(1,k) = {mfccf(ncoeff,speechIn4,FS4)};
end
for k = 1:5
q = [‘SpeechData\5’ num2str(k) ‘.wav’];
z=speechIn5(:,1);
fMatrix5(1,k) = {mfccf(ncoeff,speechIn5,FS5)};
end
for k = 1:5
q = [‘SpeechData\6’ num2str(k) ‘.wav’];
z=speechIn6(:,1);
fMatrix6(1,k) = {mfccf(ncoeff,speechIn6,FS6)};
end
for k = 1:5
q = [‘SpeechData\7’ num2str(k) ‘.wav’];
z=speechIn7(:,1);
fMatrix7(1,k) = {mfccf(ncoeff,speechIn7,FS7)};
end

%将数据保存为mat文件
fields = {‘One’,‘Two’,‘Three’,‘Four’,‘Five’};
s1 = cell2struct(fMatrix1, fields, 2); %fields项作为行
save Vectors1.mat -struct s1;
s2 = cell2struct(fMatrix2, fields, 2);
save Vectors2.mat -struct s2;
s3 = cell2struct(fMatrix3, fields, 2);
save Vectors3.mat -struct s3;
s4 = cell2struct(fMatrix4, fields, 2);
save Vectors4.mat -struct s4;
s5 = cell2struct(fMatrix5, fields, 2);
save Vectors5.mat -struct s5;
s6 = cell2struct(fMatrix6, fields, 2);
save Vectors6.mat -struct s6;
s7 = cell2struct(fMatrix7, fields, 2);
save Vectors7.mat -struct s7;

%%%%%匹配模板%%%
ncoeff = 12; %MFCC参数阶数
N = 5; %10个数字
fs=100000; %采样频率
duration2 = 5; %录音时长
k = 7; %训练样本的人数

speech = audiorecorder(fs,16,1);
disp(‘Press any key to start 5 seconds of speech recording…’);
pause
disp(‘Recording speech…’);
recordblocking(speech,duration2) % duration*fs 为采样点数
speechIn=getaudiodata(speech);
disp(‘Finished recording.’);
disp(‘System is trying to recognize what you have spoken…’);
z=speechIn(:,1);
rMatrix1 = mfccf(ncoeff,speechIn,fs); %采用MFCC系数作为特征矢量
rMatrix = CMN(rMatrix1); %归一化处理

Sco = DTWScores(rMatrix,N); %计算DTW值
[SortedScores,EIndex] = sort(Sco,2); %按行递增排序，并返回对应的原始次序
Nbr = EIndex(:,1:1) %得到每个模板匹配的2个最低值对应的次序

[Modal,Freq] = mode(Nbr(?); %返回出现频率最高的数Modal及其出现频率Freq

Word = char(‘hello’,‘matlab代码’,‘工作室’,‘欢迎’,‘你’);
if mean(abs(speechIn)) < 0.01
fprintf(‘No microphone connected or you have not said anything.\n’);
elseif (Freq <2) %频率太低不确定
fprintf(‘The word you have said could not be properly recognised.\n’);
else
fprintf(‘You have just said %s.\n’,Word(Modal,:));
end

%%%%来自PeP 工作室%%%%邮箱matlabcc2018@126.com
https://item.taobao.com/item.htm?spm=a1z10.1-c.w4004-21346770580.4.7c8172216LvrrQ&id=587158208204

2015-04-03 12:31:22 fandaoerji 阅读数 11232
• ###### 自然语言处理——实战分词

742 人正在学习 去看看 魏俊

# HMM+GMM语音识别技术详解级PMTK3中的实例

• 语音信号基础：语音信号的表示形式、分帧、特征(MFCC)、音素等等
• HMM模型：离散隐马尔科夫模型级3个问题的求解方法
• GMM：混合高斯模型，用于连续隐马尔科夫模型。

## 语音数据处理

• fs：采样率 8000Hz 115200Hz 等等，代表每1秒保存的语音数据点数
• bits:每个采样点用几个二进制保存
• 通道：很多音频都有左右2个通道，在语音识别中通常有一个通道的数据就够了。

[x fs bit]=wavread('apple.wave');
plot(x);
--------------
fs =8000
bits =16
--------------

    Fs = 8000;
framesize = 80;%真长
overlap = 20;%帧课重叠部分长度
frames = enframe(x, hamming(framesize), overlap);
>> size(frames)
131    80

80个数据点对于我们来说也比较多的，能不能再进一步减少呢？当然可以这就是特征提取了。现在语音识别中常用的特征是MFCC（可以参照：http://blog.csdn.net/xiaoding133/article/details/8106672）。我们需要知道80个数据点最后提取了多长的特征呢？答案是12个点。

## HMM模型

HMM模型的基本知识一定要读A tutorial on Hidden Markov Models and selected applications in speech recognition, L. Rabiner, 1989, Proc. IEEE 77(2):257–286. 相信我这个英文写的HMM说明比很多中文网站里写的要好懂很多。这里假设大家都懂了HMM模型是什么，我主要说一下，HMM模型怎么用到语音识别的吧。

## 少词汇量语音识别

• 状态个数N：整数
• 转移概率A：N×N矩阵
• 初始概率π：N×1矩阵
• 观察序列概率B：N×T矩阵 T是观察符号集个数
• 观察符号：模型输出的东西，可以是数值也可以是向量

N=2分别代表目前用的是红球还是绿球
π=[0.5 0.5] 试验开始我随机取了一个球因此取到红球或绿球的概率一样
A=[
0.8 0.2
0.2 0.8
] 假设我用的硬币质量不均匀，出正面的概率是0.8 背面是0.2
O={x,y} 把球落入盒子的坐标作为观察序列。

0.2 0.8 0
0 0.2 0.8
0 0 1
] (其中数值时假设的)

P（O|λ_apple）
P（O|λ_banana）
P（O|λ_kiwi）
P（O|λ_lime）
P（O|λ_orange）
P（O|λ_peach）
P（O|λ_pineapple）

## 连续大词汇量语音识别

### HMM模型

ONE w ah n
YOUNG y ah ng

A_k=[
0 1 0 0 0
0 1/2 1/2 0 0
0 0 1/2 1/2 0
0 0 0 1/2 1/2
0 0 0 0 0
]

A_ao=[
0 1 0 0 0
0 1/2 1/2 0 0
0 0 1/2 1/2 0
0 0 0 1/2 1/2
0 0 0 0 0
]

A_l=[
0 1 0 0 0
0 1/2 1/2 0 0
0 0 1/2 1/2 0
0 0 0 1/2 1/2
0 0 0 0 0
]

A_call=[
0 1 0 0 0 0 0 0 0 0 0
0 1/2 1/2 0 0 0 0 0 0 0 0
0 0 1/2 1/2 0 0 0 0 0 0 0
0 0 0 1/2 1/2 0 0 0 0 0 0
0 0 0 0 1/2 1/2 0 0 0 0 0
0 0 0 0 0 1/2 1/2 0 0 0 0
0 0 0 0 0 0 1/2 1/2 0 0 0
0 0 0 0 0 0 0 1/2 1/2 0 0
0 0 0 0 0 0 0 0 1/2 1/2 0
0 0 0 0 0 0 0 0 0 1/2 1/2
0 0 0 0 0 0 0 0 0 0 0
]

### Embedded Training

①训练用的句子的文本文件、例子： call nine one one !
②读音词典。例子：
CALL k ao l
NINE n ay n
ONE w ah n
③语音文件。sen01.wav sen02.wav
④音素列表级音素原始HMM定义
embedded training的基本思路是：读取一个句子的文本、把文本表示的一个句子转成用音素表示的句子(利用读音词典)、利用原始音素HMM定义串联起来构成句子的HMM定义(可能是非常长的HMM了)。然后把一个句子的wav文件转换成帧、提取特征，变成长的特征序列。然后把整个句子的特征序列看做是句子HMM模型的观察序列，直接却训练长的句子HMM模型(这时和普通的HMM模型一样用Baum-Welch算法)。句子训练好以后实际上各个音素也训练好了(因为句子的A、B矩阵就由音素的A、B矩阵构成)。看下图：

1.为每个训练句子建立整句HMM模型
2.初始化整句HMM模型的A矩阵，其中开始和结束状态外，每个状态只能到自己或下一个状态，概率分别是0.5
3.所有状态的B矩阵(一般用混合高斯模型)用全部训练样本的均值和方差初始化高斯模型的均值和方差
4.多次执行Baum-Welch算法。

### 语音解码

- 连续语音中词与词的边界是不知道的
- 给出一段语音信号，其中包含的额文字个数是未知的
- 搜索全部可能性是很难的，比如总共有M个词，语音长度是V个帧，全部组合是Mv$M^v$这么多可能性，无法全部遍历。

Token Passing Algorithm

### 三音素模型

HMM基础论http://www.cs.ubc.ca/~murphyk/Software/HMM/rabiner.pdf
HTK说明文档 htkbook.pdf （网上可下载）
【Speech and Language Processing】Daniel Jurafsky & James H. Martin 提取码：ceye

### 实例

PMTK是matlab中机器学习工具包。包含了很多现成的模型和算法。下面主要用他的HMM模型做一下语音识别。其中和HMM有关的代码都在pmtk3\toolbox\LatentVariableModels\hmm里面。下面是一些函数的使用说明：

创建HMM模型
function model = hmmCreate(type, pi, A, emission)
type：类型，有discrete：离散  gauss：高斯，每个状态用单个高斯分布描述  mixgausstied：混合高斯，但所有状态都用相同的分布   student：
pi:初始状态概率
A：转移矩阵
emission:发射状态描述（对应于离散时的B矩阵，但PMTK里用一个结构体描述，可以对应离散和连续多种情况，这个参数一般用其他函数生成）

function model = mkRndGaussHmm(nstates, d)
nstates：状态个数
d:一个观察符号的维度

function [observed, hidden] = hmmSample(model, len, nsamples)

HMM训练函数
function [model, loglikHist] = hmmFitEm(data, nstates, type, varargin)
data:观察序列
nstates：状态个数
type：类型，有discrete：离散  gauss：高斯  mixgausstied：混合高斯   student：
varargin：可变参数，这里可以携带很多参数，一般采用 'name','value','name','value' 的形式

model = hmmFit(data, 2, 'gauss', 'verbose', true, 'piPrior', [3 2], 'emissionPrior', prior, 'nRandomRestarts', 2, 'maxIter', 10);
【'verbose', true, 'piPrior', [3 2], 'emissionPrior', prior, 'nRandomRestarts', 2, 'maxIter', 10】就是属于varargin部分，比如【'maxIter', 10】代表最大迭代次数是10，其他都是算法中要用到的指标，可以指定也可以走默认值。

HMM训练函数
function [model, loglikHist] = hmmFit(data, nstates, type, varargin)

hmmLogprob(trueModel, observed)：计算对数概率，再模型trueModel下，观察到observed的概率是多少。出来的都是负数，绝对值越大表示概率越接近0（越小）

path=hmmMap(model, X):再模型model下，观察数列X最后可能的路径计算函数(viterbi算法)。

function [gamma, logp, alpha, beta, B] = hmmInferNodes(model, X)
% logp = log p(X | model)   对数概率
% alpha(i, t) = p(S(t)=i | X(:, 1:t)    (filtered)   α变量
% beta(i,t) propto p(X(:, t+1:T) | S(t=i))   β变量
% gamma(i,t)  = p(S(t)=i | X(:, 1:T))   (smoothed)  γ变量
% B - soft evidence  根据X推算出来的B矩阵

HMM模型用了PMTK3的库PMTK3

pmtk3 和 mfcc 自行下载并安装或把路径制定到你自己的matlab环境里，我在win7+matlab2010a 环境运行正常

clc;  clear all;  close all;

fs=8000;

Tw = 25;                % analysis frame duration (ms)  帧长
Ts = 10;                % analysis frame shift (ms)  帧移
alpha = 0.97;           % preemphasis coefficient  语音增强
M = 20;                 % number of filterbank channels
C = 12;                 % number of cepstral coefficients MFCC个数
L = 22;                 % cepstral sine lifter parameter
LF = 300;               % lower frequency limit (Hz)
HF = 3700;              % upper frequency limit (Hz)

%训练用数据
%转MFCC参数
train_feature={};
for speech=train_signals
[ MFCCs, FBEs, frames ] = mfcc( speech{1}, fs, Tw, Ts, alpha, @hamming, [LF HF], M, C+1, L );
train_feature(end+1,1)={MFCCs};
end;

%测试用数据
%转MFCC参数
test_feature={};
for speech=test_signals
[ MFCCs, FBEs, frames ] = mfcc( speech{1}, fs, Tw, Ts, alpha, @hamming, [LF HF], M, C+1, L );
test_feature(end+1,1)={MFCCs};
end;

%% HMM模型建立和训练
d = C+1; %一个观察符号的维度
nstates = 5;%状态个数
nmix    = 3; % 混合高斯分布个数

%训练apple的HMM
pi0=[1  0  0  0  0]; %初始状态概率
trans0=[1/2  1/2   0   0     0  ; %转移概率
0    1/2  1/2   0    0  ;
0    0    1/2  1/2   0  ;
0    0     0   1/2  1/2 ;
0    0     0    0    1  ];

models={};
[unique_train_labels, ~, indices] = unique(train_labels);%重复去掉，只取唯一值
for i = 1:length(unique_train_labels)
display(sprintf('modeling %s...', char(unique_train_labels(i))))

%发射概率初始值计算
stackedData = cell2mat(train_feature(indices == i)')';%训练用数据
mu = zeros(d, nmix);%均值
Sigma = zeros(d, d, nmix);%方差
for k = 1:nmix
XX             = stackedData + randn(size(stackedData));
mu(:, k)       = colvec(mean(XX));%所有训练数据的均值
Sigma(:, :, k) = cov(XX);%所有训练数据的方差
end
M = normalize(rand(nstates, nmix), 1);%混合高斯系数
emission = condMixGaussTiedCpdCreate(mu, Sigma, M);%发射状态描述结构体获得

%训练HMM
%verbose：是否打印信息
%nRandomRestarts:随机训练几次
%maxiter：最大迭代次数
%nmix：混合高斯个数
%pi0：HMM初始状态概率
%trans0：HMM初始转移概率
%emission0：初始发射状态描述
%piPrior：防止除数为0的情况出现   a/b=a+piPrior/b+piPrior  变成这种形式
%transPrior：防止除数为0的情况出现   a/b=a+piPrior/b+piPrior  变成这种形式
models(end+1,1) = {hmmFit(train_feature(indices == i), nstates, 'mixGaussTied', 'verbose', false, ...
'nRandomRestarts', 3, 'maxiter', 50, 'nmix', nmix,...
'pi0',pi0,'trans0',trans0,'emission0',emission, ...
'piPrior',pi0,'transPrior',trans0.*10)};
end;

%训练样本识别率计算
errorcount=0;
for j=1:length(train_feature)
p=zeros(length(unique_train_labels),1);
for i = 1:length(unique_train_labels)
p(i)=hmmLogprob(models{i}, train_feature(j));%计算概率值
end;
[~, i]=max(p);%取最大概率的模型作为识别
%display(sprintf('"%s" is recognized as "%s"', train_labels{j},char(unique_train_labels(i))))
if ~strcmp(train_labels{j},char(unique_train_labels(i)))
errorcount=errorcount+1;%错误累计
end;
end;
display(sprintf('train accuracy is %0.2f', (length(train_feature)-errorcount)*100/length(train_feature)));

%% 测试样本识别率计算
errorcount=0;
for j=1:length(test_feature)
p=zeros(length(unique_train_labels),1);
for i = 1:length(unique_train_labels)
p(i)=hmmLogprob(models{i}, test_feature(j));%计算概率值
end;
[~, i]=max(p);%取最大概率的模型作为识别
%display(sprintf('"%s" is recognized as "%s"', test_labels{j},char(unique_train_labels(i))))
if ~strcmp(test_labels{j},char(unique_train_labels(i)))
errorcount=errorcount+1;%错误累计
end;
end;
display(sprintf('test accuracy is %0.2f', (length(test_labels)-errorcount)*100/length(test_labels)));

function [audio_signals, word_labels] = load_audio_from_folder(audio_folder)
audio_signals = {};
word_labels = {};

for word_folder = struct2cell(dir(audio_folder))
for word_file = struct2cell(dir(sprintf('%s/%s/*.wav', audio_folder, char(word_folder(1)))))
file_path = sprintf('%s/%s/%s', audio_folder, char(word_folder(1)), char(word_file(1)));
audio_signals(end + 1) = {x(:,1)}; %#ok<AGROW>
word_labels(end + 1) = word_folder(1); %#ok<AGROW>
end
end
end

1. MFCC特征现在只用了最原始MFCC，可以把一介差分、二介差分加入构成39维度的特征
2. 连续语音识别模型建立
3. 连续语音三音素和参数共享等

## 可以自己进一步实现上述功能，不过也可以直接使用HTK工具。不过上述内容可以加深对HMM的理解，并有助于提高自己的编程能力。

2019-01-09 17:51:32 qq_14962179 阅读数 91
• ###### 自然语言处理——实战分词

742 人正在学习 去看看 魏俊
1. GMM高斯混合模型
2. HMM隐马尔科夫模型
3. EM算法
4. GMM参数估计(EM算法应用)
5. GMM-HMM模型训练
6. Baum-Welch算法(HMM前向后向算法，EM算法应用）
7. 维特比(Viterb)算法(动态规划算法的实际应用)

2019-03-07 18:55:21 weixin_43409302 阅读数 1830
• ###### 自然语言处理——实战分词

742 人正在学习 去看看 魏俊

### 一、项目简介

• 语音识别是人工智能领域的一个重要的应用场景，那么程序究竟是如何听懂语音的呢？
• 本文将用真实的音频案例，用代码呈现语音识别的基本原理和流程。
• 同时，将各种声音信号的MFCC矩阵进行可视化，“把声音的美丽画成图”。

#### 2、案例数据

• 本文采集7类声音数据（7种水果的英文读音wav文件）
• 7类分别为：apple(苹果),banana(香蕉),orange(橘子),lime(青柠檬),
kiwi(猕猴桃)，peach(桃子),pineapple(菠萝)
• 训练集各准备14个样本，测试集各一个样本

### 二、模型创建

#### 1、导入相关工具包

• 核心工具包：hmmlearn,scipy.io.wavfile,python_speech_features
import os
import warnings
import numpy as np
import scipy.io.wavfile as wf
import python_speech_features as sf
import matplotlib.pyplot as mp
import hmmlearn.hmm as hl
# 过滤代码运行中无关的警告日志
warnings.filterwarnings(
'ignore', category=DeprecationWarning)
np.seterr(all='ignore')


#### 2、定义语音识别的各类函数

# 定义声音文件的标签和路径的映射字典函数
def search_speeches(directory, speeches):
# 如果路径不是文件夹则报出异常
if not os.path.isdir(directory):
raise IOError("路径" + directory + '不是文件夹')
# 获取文件夹中的子目录
for entry in os.listdir(directory):
# 获取分类文件夹的名称作为分类标签
label = directory[directory.rfind(
os.path.sep) + 1:]
# 拼接新的文件路径
path = os.path.join(directory, entry)
# 如果路径为文件夹则继续递归向内查询
if os.path.isdir(path):
search_speeches(path, speeches)
# 如果路径为'.wav'后缀的文件名
elif os.path.isfile(path) and \
path.endswith('.wav'):
# 判断speeches中是否存在label标签
if label not in speeches:
speeches[label] = []
speeches[label].append(path)
return speeches

# 获取数据集的MFCC矩阵和标签列表
def gen_matrix(speeches):
path_x, path_y = [], []
# 获取wav文件类型标签和文件集
for label, filenames in speeches.items():
mfccs = np.array([])
# 遍历每一个wav文件
for filename in filenames:
# 提取wav文件的采样率和信号值
# 获取每个音频文件的mfcc
mfcc = sf.mfcc(sigs, sample_rate)
if len(mfccs) == 0:
mfccs = mfcc
else:
mfccs = np.append(mfccs, mfcc, axis=0)
path_x.append(mfccs)
path_y.append(label)
return path_x, path_y

# 进行模型训练,获取训练后的模型集
def model_train(path_x, path_y):
models = {}
for mfccs, label in zip(path_x, path_y):
# 利用HMM算法创建模型
model = hl.GaussianHMM(
n_components=4, covariance_type='diag',
n_iter=1000)
# 获取每个训练样本训练得到的model
models[label] = model.fit(mfccs)
return models

# 模型预测,获取样本测试的标签集
def model_pred(path_x, path_y, models):
pred_test_y = []
for mfccs in path_x:
# 初始化最优模型得分和对应的类别
best_score, best_label = None, None
# 获取模型和对应的标签
for label, model in models.items():
# 计算模型的测试得分
score = model.score(mfccs)
# 选择每个类别对应的最优模型参数
if (best_score is None) or \
best_score < score:
best_score, best_label = score, label
pred_test_y.append(best_label)
return pred_test_y

# 定义可视化函数，绘制wav文件对应的MFCC图像
def visualize(path_x, path_y):
for mfcc, label in zip(path_x, path_y):
mp.matshow(mfcc.T, cmap='jet', fignum=label)
mp.title(label, fontsize=20)
mp.xlabel('Sample', fontsize=14)
mp.ylabel('Feature', fontsize=14)
mp.tick_params(which='both', top='False',
labeltop='False', labelbottom='True',
labelsize=10)
mp.show()


#### 3、模型训练阶段

# 训练模型阶段
# 获取训练集的标签、文件字典
train_path = 'speeches/training'
train_speeches = {}
train_speeches = search_speeches(
train_path, train_speeches)
# print(train_speeches)
# 获取格式化训练样本数据集
train_x, train_y = gen_matrix(train_speeches)
# 获取训练模型集合
models = model_train(train_x, train_y)
# print(len(models))


#### 4、模型测试阶段

# 模型预测阶段
# 获取测试集的标签、文件字典
test_path = 'speeches/testing'
test_speeches = {}
test_speeches = search_speeches(
test_path, test_speeches)
# print(test_speeches)
# 获取格式化训练样本数据集
test_x, test_y = gen_matrix(test_speeches)
# 获取预测结果集
pred_test_y = model_pred(
test_x, test_y, models)
print('True Value:\n', pred_test_y)
print('Predict Value:\n', test_y)


#### 5、wav文件可视化

# 可视化各种类别的MFCC图像
visualize(test_x, test_y)


2016-03-30 19:34:54 GarfieldEr007 阅读数 2093
• ###### 自然语言处理——实战分词

742 人正在学习 去看看 魏俊

一、每个单词的读音都对应一个HMM模型，大家都知道HMM模型中有个状态集S，那么每个状态用什么来表示呢，数字？向量？矩阵？其实这个状态集中的状态没有具体的数学要求，只是一个名称而已，你可以用’1’, ’2’, ‘3’…表示，也可以用’a’, ‘b’, ’c ’表示。另外每个HMM模型中到底该用多少个状态，是通过先验知识人为设定的。

二、HMM的每一个状态都对应有一个观察值，这个观察值可以是一个实数，也可以是个向量，且每个状态对应的观察值的维度应该相同。假设现在有一个单词的音频文件，首先需要将其进行采样得到数字信息（A/D转换），然后分帧进行MFCC特征提取，假设每一帧音频对应的MFCC特征长度为39，则每个音频文件就转换成了N个MFCC向量（不同音频文件对应的N可能不同），这就成了一个序列，而在训练HMM模型的参数时（比如用Baum-Welch算法），每次输入到HMM中的数据要求就是一个观测值序列。这时，每个状态对应的观测值为39维的向量，因为向量中元素的取值是连续的，需要用多维密度函数来模拟，通常情况下用的是多维高斯函数。在GMM-HMM体系中，这个拟合函数是用K个多维高斯混合得到的。假设知道了每个状态对应的K个多维高斯的所有参数，则该GMM生成该状态上某一个观察向量（一帧音频的MFCC系数）的概率就可以求出来了。

三、对每个单词建立一个HMM模型，需要用到该单词的训练样本，这些训练样本是提前标注好的，即每个样本对应一段音频，该音频只包含这个单词的读音。当有了该单词的多个训练样本后，就用这些样本结合Baum-Welch算法和EM算法来训练出GMM-HMM的所有参数，这些参数包括初始状态的概率向量，状态之间的转移矩阵，每个状态对应的观察矩阵（这里对应的是GMM，即每个状态对应的K个高斯的权值，每个高斯的均值向量和方差矩阵）。

四、在识别阶段，输入一段音频，如果该音频含有多个单词，则可以手动先将其分割开（考虑的是最简单的方法），然后提取每个单词的音频MFCC特征序列，将该序列输入到每个HMM模型（已提前训练好的）中，采用前向算法求出每个HMM模型生成该序列的概率，最后取最大概率对应的那个模型，而那个模型所表示的单词就是我们识别的结果。

五、在建立声学模型时，可以用Deep Learning的方法来代替GMM-HMM中的GMM，因为GMM模拟任意函数的功能取决于混合高斯函数的个数，所以具有一定的局限性，属于浅层模型。而Deep Network可以模拟任意的函数，因而表达能力更强。注意，这里用来代替GMM的Deep Nets模型要求是产生式模型，比如DBN，DBM等，因为在训练HMM-DL网络时，需要用到HMM的某个状态产生一个样本的概率。

六、GMM-HMM在具体实现起来还是相当复杂的。

七、一般涉及到时间序列时才会使用HMM，比如这里音频中的语音识别，视频中的行为识别等。如果我们用GMM-HMM对静态的图片分类，因为这里没涉及到时间信息，所以HMM的状态数可设为1，那么此时的GMM-HMM算法就退化成GMM算法了。

MFCC:

MFCC的matlab实现教程可参考：张智星老师的网页教程mfcc. 最基本的12维特征。

function mfcc=frame2mfcc(frame, fs, filterNum, mfccNum, plotOpt)
% frame2mfcc: Frame to MFCC conversion.
%    Usage: mfcc=frame2mfcc(frame, fs, filterNum, mfccNum, plotOpt)
%
%    For example:
%        waveFile='what_movies_have_you_seen_recently.wav';
%        startIndex=12000;
%        frameSize=512;
%        frame=y(startIndex:startIndex+frameSize-1);
%        frame2mfcc(frame, fs, 20, 12, 1);

%    Roger Jang 20060417

if nargin<1, selfdemo; return; end
if nargin<2, fs=16000; end
if nargin<3, filterNum=20; end
if nargin<4, mfccNum=12; end
if nargin<5, plotOpt=0; end

frameSize=length(frame);
% ====== Preemphasis should be done at wave level
%a=0.95;
%frame2 = filter([1, -a], 1, frame);
frame2=frame;
% ====== Hamming windowing
frame3=frame2.*hamming(frameSize);
% ====== FFT
[fftMag, fftPhase, fftFreq, fftPowerDb]=fftOneSide(frame3, fs);
% ====== Triangular band-pass filter bank
triFilterBankPrm=getTriFilterBankPrm(fs, filterNum);    % Get parameters for triangular band-pass filter bank
% Triangular bandpass filter.
for i=1:filterNum
tbfCoef(i)=dot(fftPowerDb, trimf(fftFreq, triFilterBankPrm(:,i)));%得到filterNum个滤波系数
end
% ====== DCT
mfcc=zeros(mfccNum, 1); %DCT变换的前后个数也没有变
for i=1:mfccNum
coef = cos((pi/filterNum)*i*((1:filterNum)-0.5))'; %mfcc中的前mfccNum个系数
mfcc(i) = sum(coef.*tbfCoef');%直接按照DCT公式
end
% ====== Log energy
%logEnergy=10*log10(sum(frame.*frame));
%mfcc=[logEnergy; mfcc];

if plotOpt
subplot(2,1,1);
plot(frame, '.-');
set(gca, 'xlim', [-inf inf]);
title('Input frame');
subplot(2,1,2);
plot(mfcc, '.-');
set(gca, 'xlim', [-inf inf]);
title('MFCC vector');
end

% ====== trimf.m (from fuzzy toolbox)
function y = trimf(x, prm) %由频率的横坐标算出三角形内的纵坐标,0~1
a = prm(1); b = prm(2); c = prm(3);
y = zeros(size(x));
% Left and right shoulders (y = 0)
index = find(x <= a | c <= x);
y(index) = zeros(size(index)); %只考虑三角波内的量
% Left slope
if (a ~= b)
index = find(a < x & x < b);
y(index) = (x(index)-a)/(b-a);
end
% right slope
if (b ~= c)
index = find(b < x & x < c);
y(index) = (c-x(index))/(c-b);
end
% Center (y = 1)
index = find(x == b);
y(index) = ones(size(index));

% ====== Self demo
function selfdemo
waveFile='what_movies_have_you_seen_recently.wav';
startIndex=12000;
frameSize=512;
frame=y(startIndex:startIndex+frameSize-1);
feval(mfilename, frame, fs, 20, 12, 1);

ZCR:

过0检测，用于判断每一帧中过零点的数量情况，最简单的版本可参考：zeros cross rate.

waveFile='csNthu.wav';
frameSize=256;
overlap=0;
frameMat=enframe(y, frameSize, overlap);
frameMat(:,i)=frameMat(:,i)-mean(frameMat(:,i));    % mean justification
end
zcr=sum(frameMat(1:end-1, :).*frameMat(2:end, :)<0);
sampleTime=(1:length(y))/fs;
subplot(2,1,1); plot(sampleTime, y); ylabel('Amplitude'); title(waveFile);
subplot(2,1,2); plot(frameTime, zcr, '.-');
xlabel('Time (sec)'); ylabel('Count'); title('ZCR');

EPD:

端点检测，检测声音的起始点和终止点，可参考：EPD in Time Domain,在时域中的最简单检测方法。

waveFile='sunday.wav';
frameSize = 256;
overlap = 128;

wave=wave-mean(wave);                % zero-mean substraction
frameMat=buffer2(wave, frameSize, overlap);    % frame blocking,每一列代表一帧
frameNum=size(frameMat, 2);            % no. of frames
volume=frame2volume(frameMat);        % volume,求每一帧的能量，绝对值或者平方和,volume为行向量
volumeTh1=max(volume)*0.1;            % volume threshold 1
volumeTh2=median(volume)*0.1;            % volume threshold 2
volumeTh3=min(volume)*10;            % volume threshold 3
volumeTh4=volume(1)*5;                % volume threshold 4
index1 = find(volume>volumeTh1); %找出volume大于阈值的那些帧序号
index2 = find(volume>volumeTh2);
index3 = find(volume>volumeTh3);
index4 = find(volume>volumeTh4);
%frame2sampleIndex()为从帧序号找到样本点的序号(即每一个采样点的序号)
%endPointX长度为2,包含了起点和终点的样本点序号
endPoint1=frame2sampleIndex([index1(1), index1(end)], frameSize, overlap);
endPoint2=frame2sampleIndex([index2(1), index2(end)], frameSize, overlap);
endPoint3=frame2sampleIndex([index3(1), index3(end)], frameSize, overlap);
endPoint4=frame2sampleIndex([index4(1), index4(end)], frameSize, overlap);

subplot(2,1,1);
time=(1:length(wave))/fs;
plot(time, wave);
ylabel('Amplitude'); title('Waveform');
axis([-inf inf -1 1]);
line(time(endPoint1(  1))*[1 1], [-1, 1], 'color', 'm');%标起点终点线
line(time(endPoint2(  1))*[1 1], [-1, 1], 'color', 'g');
line(time(endPoint3(  1))*[1 1], [-1, 1], 'color', 'k');
line(time(endPoint4(  1))*[1 1], [-1, 1], 'color', 'r');
line(time(endPoint1(end))*[1 1], [-1, 1], 'color', 'm');
line(time(endPoint2(end))*[1 1], [-1, 1], 'color', 'g');
line(time(endPoint3(end))*[1 1], [-1, 1], 'color', 'k');
line(time(endPoint4(end))*[1 1], [-1, 1], 'color', 'r');
legend('Waveform', 'Boundaries by threshold 1', 'Boundaries by threshold 2', 'Boundaries by threshold 3', 'Boundaries by threshold 4');

subplot(2,1,2);
plot(frameTime, volume, '.-');
ylabel('Sum of Abs.'); title('Volume');
axis tight;
line([min(frameTime), max(frameTime)], volumeTh1*[1 1], 'color', 'm');
line([min(frameTime), max(frameTime)], volumeTh2*[1 1], 'color', 'g');
line([min(frameTime), max(frameTime)], volumeTh3*[1 1], 'color', 'k');
line([min(frameTime), max(frameTime)], volumeTh4*[1 1], 'color', 'r');
legend('Volume', 'Threshold 1', 'Threshold 2', 'Threshold 3', 'Threshold 4');

GMM:

GMM用在拟合数据分布上，本质上是先假设样本的概率分布为GMM，然后用多个样本去学习这些GMM的参数。GMM建模在语音中可用于某个单词的发音，某个人的音色等。其训练过程可参考:speaker recognition.

function [M, V, W, logProb] = gmmTrain(data, gaussianNum, dispOpt)
% gmmTrain: Parameter training for gaussian mixture model (GMM)
%    Usage: function [M, V, W, logProb] = gmm(data, gaussianNum, dispOpt)
%        data: dim x dataNum matrix where each column is a data point
%        gaussianNum: No. of Gaussians or initial centers
%        dispOpt: Option for displaying info during training
%        M: dim x meanNum matrix where each column is a mean vector
%        V: 1 x gaussianNum vector where each element is a variance for a Gaussian
%        W: 1 x gaussianNum vector where each element is a weighting factor for a Gaussian

% Roger Jang 20000610

if nargin==0, selfdemo; return; end
if nargin<3, dispOpt=0; end

maxLoopCount = 50;    % Max. iteration
minImprove = 1e-6;    % Min. improvement
minVariance = 1e-6;    % Min. variance
logProb = zeros(maxLoopCount, 1);   % Array for objective function
[dim, dataNum] = size(data);

% Set initial parameters
% Set initial M
%M = data(1+floor(rand(gaussianNum,1)*dataNum),:);    % Randomly select several data points as the centers
if length(gaussianNum)==1,
% Using vqKmeans to find initial centers
fprintf('Start KMEANS to find the initial mu...\n');
%    M = vqKmeansMex(data, gaussianNum, 0);
M = vqKmeans(data, gaussianNum, 0); %利用聚类的方法求均值,聚成gaussianNum类
%    M = vqLBG(data, gaussianNum, 0);
fprintf('Start GMM training...\n');
if any(any(~isfinite(M))); keyboard; end
else
% gaussianNum is in fact the initial centers
M = gaussianNum;
gaussianNum = size(M, 2);
end
% Set initial V as the distance to the nearest center
if gaussianNum==1
V=1;
else
distance=pairwiseSqrDist(M);%pairwiseSqrDist是dll
%distance=pairwiseSqrDist2(M);

distance(1:(gaussianNum+1):gaussianNum^2)=inf;    % Diagonal elements are inf
[V, index]=min(distance);    % Initial variance for each Gaussian
end
% Set initial W
W = ones(1, gaussianNum)/gaussianNum;    % Weight for each Gaussian,初始化时是均分权值

if dispOpt & dim==2, displayGmm(M, V, data); end
for i = 1:maxLoopCount  %开始迭代训练参数,EM算法
% Expectation step:
% P(i,j) is the probability of data(:,j) to the i-th Gaussian
% Prob为每个样本在GMM下的概率
[prob, P]=gmmEval(data, M, V, W);
logProb(i)=sum(log(prob)); %所有样本的联合概率
if dispOpt
fprintf('i = %d, log prob. = %f\n',i-1, logProb(i));
end
PW = diag(W)*P;
BETA=PW./(ones(gaussianNum,1)*sum(PW));    % BETA(i,j) is beta_i(x_j)
sumBETA=sum(BETA,2);

% Maximization step:  eqns (2.96) to (2.98) from Bishop p.67:
M = (data*BETA')./(ones(dim,1)*sumBETA');

DISTSQ = pairwiseSqrDist(M, data);                    % Distance of M to data
%DISTSQ = pairwiseSqrDist2(M, data);                    % Distance of M to data

V = max((sum(BETA.*DISTSQ, 2)./sumBETA)/dim, minVariance);    % (2.97)
W = (1/dataNum)*sumBETA;                    % (2.98)

if dispOpt & dim==2, displayGmm(M, V, data); end
if i>1, if logProb(i)-logProb(i-1)<minImprove, break; end; end
end
[prob, P]=gmmEval(data, M, V, W);
logProb(i)=sum(log(prob));
fprintf('Iteration count = %d, log prob. = %f\n',i, logProb(i));
logProb(i+1:maxLoopCount) = [];

% ====== Self Demo ======
function selfdemo
%[data, gaussianNum] = dcdata(2);
data = rand(1000,2);
gaussianNum = 8;
data=data';
plotOpt=1;
[M, V, W, lp] = feval(mfilename, data, gaussianNum, plotOpt);

pointNum = 40;
x = linspace(min(data(1,:)), max(data(1,:)), pointNum);
y = linspace(min(data(2,:)), max(data(2,:)), pointNum);
[xx, yy] = meshgrid(x, y);
data = [xx(:) yy(:)]';
z = gmmEval(data, M, V, W);
zz = reshape(z, pointNum, pointNum);
figure; mesh(xx, yy, zz); axis tight; box on; rotate3d on
figure; contour(xx, yy, zz, 30); axis image

% ====== Other subfunctions ======
function displayGmm(M, V, data)
% Display function for EM algorithm
figureH=findobj(0, 'tag', mfilename);
if isempty(figureH)
figureH=figure;
set(figureH, 'tag', mfilename);
colordef black
plot(data(1,:), data(2,:),'.r'); axis image
theta=linspace(-pi, pi, 21);
x=cos(theta); y=sin(theta);
sigma=sqrt(V);
for i=1:length(sigma)
circleH(i)=line(x*sigma(i)+M(1,i), y*sigma(i)+M(2,i), 'color', 'y');
end
set(circleH, 'tag', 'circleH', 'erasemode', 'xor');
else
circleH=findobj(figureH, 'tag', 'circleH');
theta=linspace(-pi, pi, 21);
x=cos(theta); y=sin(theta);
sigma=sqrt(V);
for i=1:length(sigma)
set(circleH(i), 'xdata', x*sigma(i)+M(1,i), 'ydata', y*sigma(i)+M(2,i));
end
drawnow
end

Speaker identification:

给N个人的语音资料，用GMM可以训练这N个人的声音模型，然后给定一段语音，判断该语音与这N个人中哪个最相似。方法是求出该语音在N个GMM模型下的概率，选出概率最大的那个。可参考:speaker recognition.

function [recogRate, confusionMatrix, speakerData]=speakerIdentify(speakerData, speakerGmm, useIntGmm)
% speakerIdentify: speaker identification using GMM parameters
%    Usage: [recogRate, confusionMatrix, speakerData]=speakerIdentify(speakerData, speakerGmm, useIntGmm)
%        speakerData: structure array generated by speakerDataRead.m
%        speakerGmm: speakerGmm(i).gmmPrm is the GMM parameters for speaker i.
%        useIntGmm: use fixed-point GMM

%    Roger Jang, 20070517, 20080726

if nargin<3, useIntGmm=0; end

% ====== Speaker identification using GMM parameters
speakerNum=length(speakerData);
for i=1:speakerNum
%    fprintf('%d/%d: Recognizing wave files by %s\n', i, speakerNum, speakerData(i).name);
for j=1:length(speakerData(i).sentence)
%        fprintf('\tSentece %d...\n', j);
%找出一个句子，看它属于哪个speaker
for k=1:speakerNum,
%            fprintf('\t\tSpeaker %d...\n', k);
%    logProb(k, :)=gmmEval(speakerData(i).sentence(j).fea, speakerGmm(k).gmmPrm);
if ~useIntGmm
%    logProb(k, :)=gmmEvalMex(speakerData(i).sentence(j).fea, gmm(k).mean, gmm(k).covariance, gmm(k).weight);
logProb(k, :)=gmmEval(speakerData(i).sentence(j).fea, speakerGmm(k).gmmPrm);
else
%    logProb(k, :)=gmmEvalIntMex(speakerData(i).sentence(j).fea, gmm(k).mean, gmm(k).covariance, gmm(k).weight);
logProb(k, :)=gmmEvalIntMex(speakerData(i).sentence(j).fea, speakerGmm(i).gmmPrm);
end
end
cumLogProb=sum(logProb, 2);
[maxProb, index]=max(cumLogProb);
speakerData(i).sentence(j).predictedSpeaker=index; %找出身份
speakerData(i).sentence(j).logProb=logProb;
end
end

% ====== Compute confusion matrix and recognition rate
confusionMatrix=zeros(speakerNum);
for i=1:speakerNum,
predictedSpeaker=[speakerData(i).sentence.predictedSpeaker];
[index, count]=elementCount(predictedSpeaker);
confusionMatrix(i, index)=count;
end
recogRate=sum(diag(confusionMatrix))/sum(sum(confusionMatrix));

GMM-HMM:

训练阶段：给出HMM的k个状态，每个状态下的观察样本的生成可以用一个概率分布来拟合，这里是采用GMM拟合的。其实，可以把GMM-HMM整体看成是一个生成模型。给定该模型的5个初始参数(结合随机和训练样本获得)，启动EM算法的E步：获得训练样本分布，即计算训练样本在各个状态下的概率。M步：用这些训练样本重新评估那5个参数。

测试阶段：(以孤立词识别为例)给定每个词发音的frame矩阵，取出某一个GMM-HMM模型，算出该发音每一帧数据在取出的GMM-HMM模型各个state下的概率，结合模型的转移概率和初始概率，获得对应的clique tree，可用图模型的方法inference出生成该语音的概率。比较多个GMM-HMM模型，取最大概率的模型对应的词。

参考资料：