《语音信号处理试验教程》（梁瑞宇等）的代码主要是Matlab实现的，现在Python比较热门，所以把这个项目大部分内容写成了Python实现，大部分是手动写的。使用CSDN博客查看帮助文件：
Python语音基础操作–2.1语音录制，播放，读取
Python语音基础操作–2.2语音编辑
Python语音基础操作–2.3声强与响度
Python语音基础操作–2.4语音信号生成
Python语音基础操作–3.1语音分帧与加窗
Python语音基础操作–3.2短时时域分析
Python语音基础操作–3.3短时频域分析
Python语音基础操作–3.4倒谱分析与MFCC系数
Python语音基础操作–4.1语音端点检测
Python语音基础操作–4.2基音周期检测
Python语音基础操作–4.3共振峰估计
Python语音基础操作–5.1自适应滤波
Python语音基础操作–5.2谱减法
Python语音基础操作–5.4小波分解
Python语音基础操作–6.1PCM编码
Python语音基础操作–6.2LPC编码
Python语音基础操作–6.3ADPCM编码
Python语音基础操作–7.1帧合并
Python语音基础操作–7.2LPC的语音合成
Python语音基础操作–10.1基于动态时间规整(DTW)的孤立字语音识别试验
Python语音基础操作–10.2隐马尔科夫模型的孤立字识别
Python语音基础操作–11.1矢量量化(VQ)的说话人情感识别
Python语音基础操作–11.2基于GMM的说话人识别模型
Python语音基础操作–12.1基于KNN的情感识别
Python语音基础操作–12.2基于神经网络的情感识别
Python语音基础操作–12.3基于支持向量机SVM的语音情感识别
Python语音基础操作–12.4基于LDA，PCA的语音情感识别
代码可在Github上下载：busyyang/python_sound_open
隐马尔科夫模型（Hidden Markov Models, HMM）作为语音信号的一种统计模型，在语音处理中得到广泛应用。
一个用于语音识别的HMM通常用三组模型参数M={A,B,π}\bold{M=\{A,B,\pi\}}M={A,B,π}来定义，假设某HMM一共有N个状态{Si}i−1N\{S_i\}_{i-1}^N{Si}i−1N,那么参数的定义为：
A\bold{A}A:状态转移概率矩阵；
A=[a11...a1N.........aN1...aNN]A=\begin{bmatrix} a_{11}& ...& a_{1N}\\...&...&...\\a_{N1}&...&a_{NN} \end{bmatrix}A=⎣⎡a11...aN1.........a1N...aNN⎦⎤
其中aija_{ij}aij是从状态SiS_iSi到状态SjS_jSj转移时的转移概率，并且0⩽aij⩽1,∑i=1Naij=10\leqslant a_{ij}\leqslant 1,\sum_{i=1}^Na_{ij}=10⩽aij⩽1,∑i=1Naij=1
π\piπ:系统初始状态概率的集合，{π}i=1N\{\pi\}_{i=1}^N{π}i=1N表示初始状态是sis_isi的概率，即：πi=P[S1=si](1⩽i⩽N),∑i=1Nπij=1\pi_i=P[S_1=s_i](1\leqslant i\leqslant N),\sum_{i=1}^N\pi_{ij}=1πi=P[S1=si](1⩽i⩽N),∑i=1Nπij=1
B\bold{B}B:处处观测值概率集合，B={bij(k)}\bold{B}=\{b_{ij}(k)\}B={bij(k)},其中bij(k)b_{ij}(k)bij(k)是从状态SiS_iSi到状态SjS_jSj转移时，观测值为k的输出概率。根据观察集合XXX的取值可将HMM分为连续型和离散型。
前向-后向算法
用来计算给定观测值序列O=o1o2...oT\bold{O}=o_1o_2...o_TO=o1o2...oT以及一个模型M={A,B,π}\bold{M=\{A,B,\pi\}}M={A,B,π}时，由模型MMM产生出OOO的概率P(O∣M)\bold{P(O|M)}P(O∣M),设S1S_1S1是初始状态，SNS_NSN是终了状态，则前向-后向算法描述如下：
（1）前向算法
根据输出观察值序列重前向后递推计算输出序列；
αt(j)\bold{\alpha_t(j)}αt(j)可以由递推计算：初始化为α0(1)=1,α0(j)=0(j≠1)\bold{\alpha_0(1)}=1,\bold{\alpha_0(j)}=0(j\neq 1)α0(1)=1,α0(j)=0(j=1)
递推公式为：
αt(j)=∑iαt−1(i)aijbij(ot),(t=1,2,...,T;i,j=1,2,...,N)\bold{\alpha_t}(j)=\sum_i\bold{\alpha_{t-1}}(i)a_{ij}b_{ij}(o_t),(t=1,2,...,T;i,j=1,2,...,N)αt(j)=i∑αt−1(i)aijbij(ot),(t=1,2,...,T;i,j=1,2,...,N)
最后结果：
P(O∣M)=αT(N)\bold{P(O|M)}=\bold{\alpha_T}(N)P(O∣M)=αT(N)
t时刻的αt(j)\bold{\alpha_t}(j)αt(j)等于t-1时刻的所有状态的αt−1(i)aijbij(ot)\bold{\alpha_{t-1}}(i)a_{ij}b_{ij}(o_t)αt−1(i)aijbij(ot)的和，如果状态SiS_iSi到状态SjS_jSj没有转移时，aij=0a_{ij}=0aij=0.前向算法计算量大大减小，为N(N+1)(T−1)N(N+1)(T-1)N(N+1)(T−1)次乘法和N(N−1)(T+1)N(N-1)(T+1)N(N−1)(T+1)次加法。
（2）后向算法
定义βt(i)\bold{\beta}_t(i)βt(i)为后向概率，即从状态SiS_iSi开始到状态SNS_NSN结束输出部分符号序列为ot+1,ot+2,...,oTo_{t+1},o_{t+2},...,o_{T}ot+1,ot+2,...,oT的概率，初始化：βT(N)=1,βT(j)=0,(j≠N)\bold{\beta_T(N)}=1,\bold{\beta_T(j)}=0,(j\neq N)βT(N)=1,βT(j)=0,(j=N)
递推公式：
βt(i)=∑jβt+1(j)aijbij(ot+1),(t=T,T−1,...,1;i,j=1,2,...,N)\bold{\beta}_t(i)=\sum_j\beta_{t+1}(j)a_{ij}b_{ij}(o_{t+1}),(t=T,T-1,...,1;i,j=1,2,...,N)βt(i)=j∑βt+1(j)aijbij(ot+1),(t=T,T−1,...,1;i,j=1,2,...,N)
所以：P(O∣M)=∑i=1Nβ1(i)πi=β0(1)\bold{P(O|M)}=\sum\limits_{i=1}^N\beta_1(i)\pi_i=\beta_0(1)P(O∣M)=i=1∑Nβ1(i)πi=β0(1)
后向计算的计算量为N2TN^2TN2T，根据定义可以知道：P(O∣M)=∑i=1N∑j=1Nαt(i)aijbij(ot+1)βt+1(j)\bold{P(O|M)}=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\alpha_t(i)a_{ij}b_{ij}(o_{t+1})\beta_{t+1}(j)P(O∣M)=i=1∑Nj=1∑Nαt(i)aijbij(ot+1)βt+1(j)
维特比（Viterbi）算法
Viterbi解决的是给定观察符号序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT和模型M={A,B,π}\bold{M=\{A,B,\pi\}}M={A,B,π},求出状态序列S=s1s2...sTS=s_1s_2...s_TS=s1s2...sT的问题。最佳意义上的状态序列是使P(S,O∣M)P(S,O|M)P(S,O∣M)最大时确定的状态序列，即HMM输出一个观察值序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT时，可能通过的状态序列有多种，这里面使输出概率最大的序列状态。
初始化：α0′(1)=1,α0′(j)=0(j≠1)\alpha_0'(1)=1,\alpha_0'(j)=0(j\neq 1)α0′(1)=1,α0′(j)=0(j=1)
递推公式：αt′(j)=maxiαt−1′(j−1)aijbij(ot),(t=1,2,...,T;i,j=1,2,...,N)\alpha_t'(j)=\underset{i}{\max}\alpha_{t-1}'(j-1)a_{ij}b_{ij}(o_t),(t=1,2,...,T;i,j=1,2,...,N)αt′(j)=imaxαt−1′(j−1)aijbij(ot),(t=1,2,...,T;i,j=1,2,...,N)
最后：Pmax(S,O∣M)=αT′(N)P_{\max}(S,O|M)=\alpha_T'(N)Pmax(S,O∣M)=αT′(N)
每一次使αt′(j)\alpha_t'(j)αt′(j)最大的状态i组成的状态序列就是所求的最佳状态序列。求最佳状态序列的方式为：
1）给定每个状态准备一个数组变量αt′(j)\alpha_t'(j)αt′(j),初始化时，令初始状态S1S_1S1的数组变量α0′(1)=1\alpha_0'(1)=1α0′(1)=1,其他状态α0′(j)=0\alpha_0'(j)=0α0′(j)=0
2）根据t时刻输出的观察符号oto_tot计算αt′(j)=maxiαt−1′aijbij(ot)=maxi{αt−1′a1jb1j(ot),αt−1′a2jb2j(ot),...,αt−1′aNjbNj(ot)}\alpha_t'(j)=\underset{i}{\max}\alpha_{t-1}'a_{ij}b_{ij}(o_t)=\underset{i}{\max}\{\alpha_{t-1}'a_{1j}b_{1j}(o_t),\alpha_{t-1}'a_{2j}b_{2j}(o_t),...,\alpha_{t-1}'a_{Nj}b_{Nj}(o_t)\}αt′(j)=imaxαt−1′aijbij(ot)=imax{αt−1′a1jb1j(ot),αt−1′a2jb2j(ot),...,αt−1′aNjbNj(ot)}，当状态SiS_iSi到状态SjS_jSj没有转移时，aij=0a_{ij}=0aij=0。设计一个符号数组变量，称为追加状态序列寄存器，利用这个最佳状态序列寄存器吧每次使得αt′(j)\alpha_t'(j)αt′(j)最大的状态保存下来。
3）当t≠Tt\neq Tt=T时转移到2），否则转移到4）。
4）把最终的状态寄存器αT′(N)\alpha_T'(N)αT′(N)内的值取出，则Pmax(S,O∣M)=αT′(N)P_{\max}(S,O|M)=\alpha_T'(N)Pmax(S,O∣M)=αT′(N)为最佳状态序列寄存器的值，就是所求的最佳状态序列。
Baum-Welch算法
Baum-Welch算法解决的是HMM训练，即HMM参数估计问题，给定一个观察序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT，该算法能确定一个M={A,B,π}M=\{A,B,\pi\}M={A,B,π}，使P(O∣M)P(O|M)P(O∣M)最大，这是一个泛函极值问题。由于给定的训练序列有限，因而不存在一个最佳的方法来估计模型，Baum-Welch算法也是利用递归思想，使P(O∣M)P(O|M)P(O∣M)局部放大后，最后得到优化的模型参数M={A,B,π}M=\{A,B,\pi\}M={A,B,π}。利用Baum-Welch算法重估公式得到的新模型M^\hat MM^,一定有P(O∣M^)>P(O∣M)P(O|\hat M)>P(O|M)P(O∣M^)>P(O∣M),重复估计过程，直到P(O∣M^)P(O|\hat M)P(O∣M^)收敛，不再明显增大，此时的M^\hat MM^即为所求模型。
给定一个观察值符号序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT，以及一个需要通过训练进行重估计参数的HMM模型M={A,B,π}M=\{A,B,\pi\}M={A,B,π},按前向-后向算法，设对于符号序列，在时刻t从状态SiS_iSi到状态SjS_jSj的转移概率为γt(i,j)\gamma_t(i,j)γt(i,j):
γt(i,j)=αt−1(i)aijbij(ot)βt(j)αT(N)=αt−1(i)aijbij(ot)βt(j)∑iαt(i)βt(i)\gamma_t(i,j)=\frac{\alpha_{t-1}(i)a_{ij}b_{ij}(o_t)\beta_t(j)}{\alpha_T(N)}=\frac{\alpha_{t-1}(i)a_{ij}b_{ij}(o_t)\beta_t(j)}{\sum_i\alpha_t(i)\beta_t(i)}γt(i,j)=αT(N)αt−1(i)aijbij(ot)βt(j)=∑iαt(i)βt(i)αt−1(i)aijbij(ot)βt(j)
同时，对于符号序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT，在t时刻的Markov链处于状态SiS_iSi的概率为：
∑j=1Nγt(i,j)=αt(i)βt(i)∑iαt(i)βt(i)\sum\limits_{j=1}^N\gamma_t(i,j)=\frac{\alpha_t(i)\beta_t(i)}{\sum_i\alpha_t(i)\beta_t(i)}j=1∑Nγt(i,j)=∑iαt(i)βt(i)αt(i)βt(i)
这时，状态SiS_iSi到状态SjS_jSj转移次数的期望为∑tγt(i,j)\sum_t\gamma_t(i,j)∑tγt(i,j)，而从状态SiS_iSi转移出去的次数期望为∑j∑tγt(i,j)\sum_j\sum_t\gamma_t(i,j)∑j∑tγt(i,j)，所以重估公式为：
a^ij=∑tγt(i,j)∑j∑tγt(i,j)=∑tαt−1(i)aijbij(ot)βt(j)∑tαt(i)βt(i)\hat a_{ij}=\frac{\sum_t\gamma_t(i,j)}{\sum_j\sum_t\gamma_t(i,j)}=\frac{\sum_t\alpha_{t-1}(i)a_{ij}b_{ij}(o_t)\beta_t(j)}{\sum_t\alpha_t(i)\beta_t(i)}a^ij=∑j∑tγt(i,j)∑tγt(i,j)=∑tαt(i)βt(i)∑tαt−1(i)aijbij(ot)βt(j)
b^ij=∑t:ot=kγt(i,j)∑tγt(i,j)=∑t:ot=kαt−1(i)aijbij(ot)βt(j)∑tαt−1(i)aijbij(ot)βt(j)\hat b_{ij}=\frac{\sum\limits_{t:o_t=k}\gamma_t(i,j)}{\sum_t\gamma_t(i,j)}=\frac{\sum\limits_{t:o_t=k}\alpha_{t-1}(i)a_{ij}b_{ij}(o_t)\beta_t(j)}{\sum_t\alpha_{t-1}(i)a_{ij}b_{ij}(o_t)\beta_t(j)}b^ij=∑tγt(i,j)t:ot=k∑γt(i,j)=∑tαt−1(i)aijbij(ot)βt(j)t:ot=k∑αt−1(i)aijbij(ot)βt(j)
得到的新模型就是M^={A^,B^,π^}\hat M=\{\hat A,\hat B,\hat\pi\}M^={A^,B^,π^}。
具体的实现步骤为：
1)适当地选择aija_{ij}aij和bij(k)b_{ij}(k)bij(k)的初始值，常用的设定方式为：给予从状态i转移出去的每条弧相等的转移概率，即
aij=1从状态i转移出去的弧的条数a_{ij}=\frac{1}{从状态i转移出去的弧的条数}aij=从状态i转移出去的弧的条数1
给予每个输出观察符号相等的输出概率初始值，即：
bij(k)=1码本中码字的个数b_{ij}(k)=\frac{1}{码本中码字的个数}bij(k)=码本中码字的个数1
并且每条弧上给予相同的输出概率矩阵。
2）给定一个（训练）观察值符号序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT，由初始模型计算γt(i,j)\gamma_t(i,j)γt(i,j),并且由重估公式，计算a^ij\hat a_{ij}a^ij和b^ij(k)\hat b_{ij}(k)b^ij(k).
3)再给定一个（训练）观察值序列O=o1o2...oTO=o_1o_2...o_TO=o1o2...oT，吧前一次的a^ij\hat a_{ij}a^ij和b^ij(k)\hat b_{ij}(k)b^ij(k)作为初始模型计算γt(i,j)\gamma_t(i,j)γt(i,j)，重新计算a^ij\hat a_{ij}a^ij和b^ij(k)\hat b_{ij}(k)b^ij(k).
4）直到a^ij\hat a_{ij}a^ij和b^ij(k)\hat b_{ij}(k)b^ij(k)收敛为止。
语音识别一般采用从左到右的HMM，所以初始状态概率πi\pi_iπi不需要顾及，总设定为:
π1=1,πi=0,(i=2,...,N)\pi_1=1,\pi_i=0,(i=2,...,N)π1=1,πi=0,(i=2,...,N)
from chapter3_分析实验.mel import Nmfccfrom scipy.io import wavfile, loadmatfrom hmmlearn import hmmfrom sklearn.externals import joblibimport numpy as npimport os"""代码来自：/chinatelecom08/article/details/82901480并进行了部分更改"""def gen_wavlist(wavpath):"""得到数据文件序列:param wavpath::return:"""wavdict = {}labeldict = {}for (dirpath, dirnames, filenames) in os.walk(wavpath):for filename in filenames:if filename.endswith('.wav'):filepath = os.sep.join([dirpath, filename])fileid = filename.strip('.wav')wavdict[fileid] = filepathlabel = fileid.split('_')[1]labeldict[fileid] = labelreturn wavdict, labeldictdef compute_mfcc(file):"""读取数据并计算mfcc:param file: 文件名:return: mfcc系数""""""有手动修改wavfile.read()函数的返回值，添加了bit_depth的返回，如果报错，修改调用方式为：fs, audio = wavfile.read(file)-3-20 Jie Y."""fs, audio, bits = wavfile.read(file)"""由于部分信号太短而报错，所以fs//2了"""mfcc = Nmfcc(audio, fs // 2, 12, frameSize=int(fs // 2 * 0.025), inc=int(fs // 2 * 0.01))return mfcc'''&usage:搭建HMM-GMM的孤立词识别模型参数意义：CATEGORY:所有标签的列表n_comp:每个孤立词中的状态数n_mix:每个状态包含的混合高斯数量cov_type:协方差矩阵的类型n_iter:训练迭代次数'''class Model:def __init__(self, CATEGORY=None, n_comp=3, n_mix=3, cov_type='diag', n_iter=1000):super(Model, self).__init__()self.CATEGORY = CATEGORYself.category = len(CATEGORY)self.n_comp = n_compself.n_mix = n_mixself.cov_type = cov_typeself.n_iter = n_iter# 关键步骤，初始化models，返回特定参数的模型的列表self.models = []for k in range(self.category):model = hmm.GMMHMM(n_components=self.n_comp, n_mix=self.n_mix, covariance_type=self.cov_type,n_iter=self.n_iter)self.models.append(model)def train(self, tdata):for i in range(tdata.shape[1]):model = self.models[i]for x in range(tdata[0, i].shape[1]):data = tdata[0, i][0, x].squeeze()mfcc = Nmfcc(data, 8000, 24, 256, 80)model.fit(mfcc)def test(self, pdata):label = []result = []for k in range(pdata.shape[1]):for i in range(pdata[0, k].shape[1]):label.append(str(k + 1))data = pdata[0, k][0, i].squeeze()mfcc = Nmfcc(data, 8000, 24, 256, 80)result_one = []for m in range(self.category):model = self.models[m]re = model.score(mfcc)result_one.append(re)result.append(self.CATEGORY[np.argmax(np.array(result_one))])print('识别得到结果：\n', result)print('原始标签类别：\n', label)# 检查识别率，为：正确识别的个数/总数totalnum = len(label)correctnum = 0for i in range(totalnum):if result[i] == label[i]:correctnum += 1print('识别率:', correctnum / totalnum)def save(self, path="models.pkl"):joblib.dump(self.models, path)def load(self, path="models.pkl"):self.models = joblib.load(path)tdata = loadmat('tra_data.mat')['tdata']pdata = loadmat('rec_data.mat')['rdata']CATEGORY = [str(i + 1) for i in range(tdata.shape[1])]# 进行训练models = Model(CATEGORY=CATEGORY)models.train(tdata)models.test(tdata)models.test(pdata)
