本案例基于HMM-GMM模型对单音源语音进行识别,主要数据是7种英文单词的单音源语音,每种单词包括15个文件。我们采用的是峰值特征,后续将会尝试使用MFCCs特征提取方法提取更优的特征进行识别,注意:本案例中的特征提取手段对于多音源及连续语音识别不适用。
本案例主要包括以下内容:
1. 语音识别框架
1.1 特征提取
1.2 解码
2. 数据探索
2.1 语音文件包含的内容
2.2 音频波形的可视化
3. 特征提取
3.1 短时傅里叶变换STFT
3.2 提取峰值
4. GMM-HMM模型解码
5. 语音识别模型训练及预测
6. 总结
语音识别是自然语言处理中最具挑战性的工作之一,也是近几年快速发展的一个领域。自动化的语音识别(Automatic Speech Recognition,ASR)可以使得机器能够理解人类所发出的声音(如各国的语言),它有着广泛的应用,比如翻译、人机交互设计和语音信息检索(Speech Information Retrieval,SIR)等。现代ASR技术将信号处理、模式识别、网络和通信融合到统一的框架中。目前,ASR中采用的技术可以分为三类:
在过去的一段时间里,HMM已经成为了ASR中最成功的语音识别模型。HMM对于时间序列的分析非常有效,几乎现在所有的超大词汇连续语音识别(Large Vocabulary Continuous Speech Recognition,LVCSR)系统都是基于HMM的。
在本案例中,我们首先简述基于HMM的LVCSR系统的核心框架,然后利用HMM模型(具体为Gaussian Mixture Model HMM, GMM-HMM)具体实现单一说话者的词汇识别。这个案例并不会完整的实现一个语音识别系统,我们只是希望通过这个案例,使得读者能够深入理解HMM在语音识别上的一些核心思想。
除了本案例实现HMM的方法,HMM的实现还可以使用Python工具包 hmmlearn。
整个基于HMM的语音识别系统分为特征提取和基于HMM的解码过程。
这一阶段从音频波形中提取便于进行语音识别的特征,相当于Encoding部分。最简单且最广泛应用的方法之一是采用梅尔倒频谱系数(mel-frequency cepstral coefficients,MFCCs)分析的方法。
Mel频率(mel-frequency)是基于人耳听觉特性提出来的(类似分贝概念),它与赫兹频率成非线性对应关系。倒谱(cepstrum)就是一种信号的傅里叶变换经对数运算后再进行傅里叶反变换得到的谱。经过MFCCs特征提取后,波形图就编程了mel倒频谱图,如下所示:
整个特征提取过程即把音频波形转换为长度固定的声矢量(acoustic vectors)$\boldsymbol{Y}_{1:T}=\boldsymbol{y}_1,\ldots,\boldsymbol{y}_T$的过程。
简而言之,解码过程就是要找到最有可能产生声矢量$\boldsymbol{Y}$的单词序列$\boldsymbol{w}_{1:L}=w_1,\ldots,w_L$,即找到
$$\hat{\boldsymbol{w}}=\arg\max_{\boldsymbol{w}}\{p(\boldsymbol{w}|\boldsymbol{Y})\}$$但是判别式的$P(\boldsymbol{w}|\boldsymbol{Y})$不容易通过直接建模得到。故需要通过贝叶斯定理将上述问题转换为另一种等价形式:
$$\hat{\boldsymbol{w}}=\arg\max_{\boldsymbol{w}}\{p(\boldsymbol{Y}|\boldsymbol{w})p(\boldsymbol{w})\}$$似然$p(\boldsymbol{Y}|\boldsymbol{w})$由声学模型(acoustic model)确定,先验概率$p(\boldsymbol{w})$由语言模型(language model)确定。在声学模型中,声音(sound)的基本单位为音素(phone),如$/b/$、$/ae/$、$/t/$,声学模型的核心即是HMM模型。这些音素由发音字典(pronunciation dictionary)确定音素所组成的单词,见基于HMM的语音识别框架图。
本案例所使用的数据集(来源为 Google Code project by Hakon Sandsmark)是单音源(single speaker)的数据,并不是多个说话者(multispeaker)的数据,所以在模型上对鲁棒性要求没有特别高。本案例所使用的特征提取方法对于多说话者的语音识别不适用。
先来看看文件夹中到底包含了哪些单词的语音:
fpaths = [] # 文件路径
labels = [] # 每一个语音文件所包含的单词内容
spoken = [] # 所有语音文件所包含的单词内容列表
import os
for f in os.listdir('./audio')[1:]:
for w in os.listdir('./audio/' + f):
fpaths.append('./audio/' + f + '/' + w)
labels.append(f)
if f not in spoken:
spoken.append(f)
print('Words spoken:', spoken)
为了便于后续的分析,我们用一个标签向量all_labels
来表示每一个语音文件的单词内容。
from scipy.io import wavfile
import numpy as np
data = np.zeros((len(fpaths), 32000)) # 每一行为一个语音样本,每一列为样本的一个特征
maxsize = -1 # 最大的列数,每一个音频文件的时长不同,所以在循环中需要不断更改
for n,file in enumerate(fpaths):
_, d = wavfile.read(file)
data[n, :d.shape[0]] = d
if d.shape[0] > maxsize:
maxsize = d.shape[0]
data = data[:, :maxsize]
#Each sample file is one row in data, and has one entry in labels
print('Number of files total:', data.shape[0])
all_labels = np.zeros(data.shape[0])
for n, l in enumerate(set(labels)):
all_labels[np.array([i for i, _ in enumerate(labels) if _ == l])] = n
print('data',data.shape)
print('Labels and label indices', all_labels)
共有105个文件,有7种不同的单词,每一个单词一共读了15遍,每一遍用单独的wav格式的语音文件存储。数据集data
总共为105行6966列。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns # 导入seaborn包
sns.set_context("notebook",font_scale=1.3) # notebook格式,放大横纵坐标标记
sns.set_palette('Set2') # 配色使用Set2
# 导入音频处理及DCT(离散余弦变换)转换所需要使用的包
import scipy.io.wavfile
from scipy.fftpack import dct
# 内嵌模式
%matplotlib inline
%config InlineBackend.figure_format = 'retina' # retina高清显示
为了方便后续画出频谱图,这里我们定义一个画图的函数spectrum_plot
:
"""
函数 spectrum_plot():画出频谱图
输入参数 sequence:时间序列,以毫秒(ms)为单位
输出参数 amplitude:幅值
"""
def spectrum_plot(sequence,amplitude):
plt.figure(figsize=(8,4))
plt.plot(time_sequence,signal)
plt.xlabel('miliseconds')
plt.ylabel('Amplitude(signed 16 bit)')
选取文件./audio/apple/apple01.wav
文件,观察其音频波形。
sample_rate, signal = scipy.io.wavfile.read('./audio/apple/apple01.wav') # 导入数据
print 'sample_rate:',sample_rate,'Hz'
print 'signal time_interval:',1000*signal.shape[0]/float(sample_rate),'ms'
整个音频采样频率为8000 Hz,这个音频持续时长为336.75 ms。定义time_sequence
设置画出波形的时长,这里需要根据采样频率计算转换为时间毫秒量纲。
time_sequence = [item*1000/float(sample_rate) for item in range(signal.shape[0])]
spectrum_plot(time_sequence,signal)
传统音频特征的提取方法的有两种,一种方法为简单的共振峰提取,另一种为应用广泛并在实际使用时证明效果显著的梅尔倒频谱系数MFCCs提取的方法。目前已经有相当一部分人开始尝试使用深度神经网络进行特征提取,我们之后在其他的案例中再考虑阐述这个方法。特征提取的方法直接决定了语音识别系统的性能。在本案例中我们的重点在于HMM模型应用部分,所以我们使用简单的共振峰提取方法提取特征。
现在我们已经将语音数据导入到了数组signal
中,接下来就需要从原始的数据当中提取出特征。我们使用简单的频率峰值检测的方法提取特征。
为了找到频谱中的峰值,我们需要用到一种叫做短时傅里叶变换(Short Time Fourier Transform, STFT)的技术。STFT是许多数字信号处理(Digital Signal Processor, DSP)过程中的关键技术,并且由许多高效的计算方式(参考基于numpy的实现,在matplotlib中有直接实现STFT的specgram()函数,由于我们的数据比较简单,所以我们选择自己实现一个实现STFT的函数。
该方法的原理其实很简单:首先语音将会被分成几个小块,然后对每一小块进行傅里叶变换FFT,得到一个二维的FFT“图片”,这个“图片”即我们通常所说的频谱。选择不同的FFT的大小会得到不同频率分辨率下的结果,而重叠这些窗口可以让我们以增加数据大小为代价来控制时间分辨率。
简单来说,假设原始数据X
是一个长度为20的向量,如果FFT的大小为10,重叠窗口为5,即FFT的一半,那么我们将得到二维的STFT_X
结果如下(用伪代码的形式展示):
STFT_X[0, :] = FFT(X[0:9])
STFT_X[1, :] = FFT(X[5:14])
STFT_X[2, :] = FFT(X[10:19])
这样我们就从原始数据X
得到了3个FFT帧。接下来,我们就需要从每一行的STFT_X
中找到峰值。
import scipy
"""
函数 stft():实现短时傅里叶变换STFT
输入参数 x:原始数据,ndarray;
fftsize:FFT大小,默认为64
overlap_pct:重叠窗口大小(FFT大小的倍数),默认为0.5
输出参数 raw[:,:(fftsize//2)]:二维STFT谱
"""
def stft(x, fftsize=64, overlap_pct=.5):
hop = int(fftsize * (1 - overlap_pct))
w = scipy.hanning(fftsize + 1)[:-1] # 转换为汉明窗口
raw = np.array([np.fft.rfft(w * x[i:i + fftsize]) for i in range(0, len(x) - fftsize, hop)])
return raw[:, :(fftsize // 2)]
定义STFT函数后,我们就来看看其中一个语音经过STFT转换后的频谱图。
plt.figure(figsize=(8,4))
plt.plot(data[0, :])
plt.title('Timeseries example for %s'%labels[0])
plt.xlim(0, 3000)
plt.xlabel('Time (samples)')
plt.ylabel('Amplitude (signed 16 bit)')
plt.figure()
# 画出STFT功率谱PSD(power spectral density)
# 为了避免计算log0,在后面加1
log_freq = 20 * np.log(np.abs(stft(data[0, :])) + 1) # stft结果取对数
print(log_freq.shape)
plt.imshow(log_freq, cmap='winter', interpolation=None)
plt.colorbar()
plt.xlabel('Freq (bin)')
plt.ylabel('Time (overlapped frames)')
plt.ylim(log_freq.shape[1])
plt.title('PSD of %s example'%labels[0])
有了STFT结果之后,我们就可以使用滑窗提取峰值,算法的主要步骤如下:
x
的窗口,在这个例子中,我们使用的窗口大小为9;data[0:9]
到data[1:10]
;我们定义一个peakfind
函数提取峰值。
from numpy.lib.stride_tricks import as_strided
"""
函数 peakfind():STFT的峰值检测
输入参数 x:STFT谱;
n_peaks:最终选择峰值数
l_size:窗口左半部分大小,默认为3
r_size:窗口中间部分大小,默认为3
c_size:窗口右半部分大小,默认为3
f:对每部分采用的统计函数,默认为均值np.mean
输出参数 heights:波峰对应的高度(幅值)
top[:n_peaks]:前n_peaks个波峰的频率值
"""
def peakfind(x, n_peaks, l_size=3, r_size=3, c_size=3, f=np.mean):
win_size = l_size + r_size + c_size
shape = x.shape[:-1] + (x.shape[-1] - win_size + 1, win_size)
strides = x.strides + (x.strides[-1],)
xs = as_strided(x, shape=shape, strides=strides)
def is_peak(x):
centered = (np.argmax(x) == l_size + int(c_size/2))
l = x[:l_size]
c = x[l_size:l_size + c_size]
r = x[-r_size:]
passes = np.max(c) > np.max([f(l), f(r)])
if centered and passes:
return np.max(c)
else:
return -1
r = np.apply_along_axis(is_peak, 1, xs)
top = np.argsort(r, None)[::-1]
heights = r[top[:n_peaks]]
top[top > -1] = top[top > -1] + l_size + int(c_size / 2.)
return heights, top[:n_peaks]
下面我们分别同样data[0,:]
来看一下STFT谱在同一帧下,不同窗口大小下的峰值检测结果。
plt.figure(figsize=(8,4))
plot_data = np.abs(stft(data[0, :]))[20,:]
values, locs = peakfind(plot_data, l_size=1,c_size=1,r_size=1,n_peaks=6)
fp_3 = locs[values > -1]
fv_3 = values[values > -1]
values, locs = peakfind(plot_data, l_size=3,c_size=3,r_size=3,n_peaks=6)
fp_9 = locs[values > -1]
fv_9 = values[values > -1]
f1 = plt.plot(plot_data)
f2 = plt.plot(fp_3, fv_3, '<')
f3 = plt.plot(fp_9, fv_9, '>')
plt.legend(['stft data','window size 3','window size 9'])
plt.xlabel('Frequency (bins)')
plt.ylabel('Amplitude')
plt.title('Peak location example')
我们观察到,窗口为9时,满足条件的峰值只有3个,并不能取到完整的6个峰值。所以我们的峰值搜索并不是完美的。因为对于大小为64的FFT,大小为9的窗口相对来说比较大,容易过滤掉一些峰值。这对后续的模型有一些负面影响。如果为了找到更多的峰值,过大的FFT可以做到,但是在同等时间分辨率的情况下,后续模型的效果也会受到负面影响。所以就这就需要一些高级的特征提取方法,如MFCCs,这里我们不会涉及到。
在特征提取的最后一部分,我们提取所有语音样本的峰值数据作为后续GMM-HMM模型的输入。
all_obs = []
for i in range(data.shape[0]):
d = np.abs(stft(data[i, :]))
n_dim = 6
obs = np.zeros((n_dim, d.shape[0]))
for r in range(d.shape[0]):
_, t = peakfind(d[r, :], n_peaks=n_dim)
obs[:, r] = t.copy()
if i % 10 == 0:
print("Processed obs %s" % i)
all_obs.append(obs)
all_obs = np.atleast_3d(all_obs)
print all_obs.shape
最终得到的是一个$105\times6\times216$的三维数组,最终每一个语音波形都得到216个重叠帧,每一帧用最大的6个峰值对应的频率来描述。
HMM模型主要用于三种情形中:
我们使用Baum-Welch算法训练HMM模型。同时,这一部分使用的是scipy 0.14用于多变量正态密度计算。接下来,我们直接定义GMM-HMM模型的类GmmHmm
,然后将前向后向算法,EM过程,模型训练、预测等过程都通过定义函数的方法封装在类GmmHmm
中。
import scipy.stats as st
import numpy as np
class GmmHmm:
def __init__(self, n_states):
self.n_states = n_states
self.random_state = np.random.RandomState(0)
# 随机初始化
self.prior = self._normalize(self.random_state.rand(self.n_states, 1))
self.A = self._stochasticize(self.random_state.rand(self.n_states, self.n_states))
self.mu = None
self.covs = None
self.n_dims = None
def _forward(self, B):
log_likelihood = 0.
T = B.shape[1]
alpha = np.zeros(B.shape)
for t in range(T):
if t == 0:
alpha[:, t] = B[:, t] * self.prior.ravel()
else:
alpha[:, t] = B[:, t] * np.dot(self.A.T, alpha[:, t - 1])
alpha_sum = np.sum(alpha[:, t])
alpha[:, t] /= alpha_sum
log_likelihood = log_likelihood + np.log(alpha_sum)
return log_likelihood, alpha
def _backward(self, B):
T = B.shape[1]
beta = np.zeros(B.shape);
beta[:, -1] = np.ones(B.shape[0])
for t in range(T - 1)[::-1]:
beta[:, t] = np.dot(self.A, (B[:, t + 1] * beta[:, t + 1]))
beta[:, t] /= np.sum(beta[:, t])
return beta
def _state_likelihood(self, obs):
obs = np.atleast_2d(obs)
B = np.zeros((self.n_states, obs.shape[1]))
for s in range(self.n_states):
#N 这一部分需要 0.14版本以上的Scipy
np.random.seed(self.random_state.randint(1))
B[s, :] = st.multivariate_normal.pdf(
obs.T, mean=self.mu[:, s].T, cov=self.covs[:, :, s].T)
return B
def _normalize(self, x):
return (x + (x == 0)) / np.sum(x)
def _stochasticize(self, x):
return (x + (x == 0)) / np.sum(x, axis=1)
def _em_init(self, obs):
if self.n_dims is None:
self.n_dims = obs.shape[0]
if self.mu is None:
subset = self.random_state.choice(np.arange(self.n_dims), size=self.n_states, replace=False)
self.mu = obs[:, subset]
if self.covs is None:
self.covs = np.zeros((self.n_dims, self.n_dims, self.n_states))
self.covs += np.diag(np.diag(np.cov(obs)))[:, :, None]
return self
def _em_step(self, obs):
obs = np.atleast_2d(obs)
B = self._state_likelihood(obs)
T = obs.shape[1]
log_likelihood, alpha = self._forward(B)
beta = self._backward(B)
xi_sum = np.zeros((self.n_states, self.n_states))
gamma = np.zeros((self.n_states, T))
for t in range(T - 1):
partial_sum = self.A * np.dot(alpha[:, t], (beta[:, t] * B[:, t + 1]).T)
xi_sum += self._normalize(partial_sum)
partial_g = alpha[:, t] * beta[:, t]
gamma[:, t] = self._normalize(partial_g)
partial_g = alpha[:, -1] * beta[:, -1]
gamma[:, -1] = self._normalize(partial_g)
expected_prior = gamma[:, 0]
expected_A = self._stochasticize(xi_sum)
expected_mu = np.zeros((self.n_dims, self.n_states))
expected_covs = np.zeros((self.n_dims, self.n_dims, self.n_states))
gamma_state_sum = np.sum(gamma, axis=1)
#在除法运算之前,将0改为1
gamma_state_sum = gamma_state_sum + (gamma_state_sum == 0)
for s in range(self.n_states):
gamma_obs = obs * gamma[s, :]
expected_mu[:, s] = np.sum(gamma_obs, axis=1) / gamma_state_sum[s]
partial_covs = np.dot(gamma_obs, obs.T) / gamma_state_sum[s] - np.dot(expected_mu[:, s], expected_mu[:, s].T)
# 对称
partial_covs = np.triu(partial_covs) + np.triu(partial_covs).T - np.diag(partial_covs)
#加上对角元素以保证半正定
expected_covs += .01 * np.eye(self.n_dims)[:, :, None]
self.prior = expected_prior
self.mu = expected_mu
self.covs = expected_covs
self.A = expected_A
return log_likelihood
def fit(self, obs, n_iter=15):
# obj支持2D或者3D数组
# 2D数组的形式为 n_features, n_dims
# 3D数组的形式为 n_examples, n_features, n_dims
# 例如一个语音片段若包含6个特征,并且由105个不同的样本
# 那么这个obj数组为三维,大小为(105,6,X),其中X为frames的个数
# 对于只包含一个样本的数据obj大小为(6,X)
if len(obs.shape) == 2:
for i in range(n_iter):
self._em_init(obs)
log_likelihood = self._em_step(obs)
elif len(obs.shape) == 3:
count = obs.shape[0]
for n in range(count):
for i in range(n_iter):
self._em_init(obs[n, :, :])
log_likelihood = self._em_step(obs[n, :, :])
return self
def transform(self, obs):
# obj支持2D或者3D数组
# 2D数组的形式为 n_features, n_dims
# 3D数组的形式为 n_examples, n_features, n_dims
# 例如一个语音片段若包含6个特征,并且由105个不同的样本
# 那么这个obj数组为三维,大小为(105,6,X),其中X为frames的个数
# 对于只包含一个样本的数据obj大小为(6,X)
if len(obs.shape) == 2:
B = self._state_likelihood(obs)
log_likelihood, _ = self._forward(B)
return log_likelihood
elif len(obs.shape) == 3:
count = obs.shape[0]
out = np.zeros((count,))
for n in range(count):
B = self._state_likelihood(obs[n, :, :])
log_likelihood, _ = self._forward(B)
out[n] = log_likelihood
return out
下面我们任意选取两个(4,40)的数组放入类中进行训练。
if __name__ == "__main__":
rstate = np.random.RandomState(0)
t1 = np.ones((4, 40)) + .001 * rstate.rand(4, 40)
t1 /= t1.sum(axis=0) # 归一化
t2 = rstate.rand(*t1.shape)
t2 /= t2.sum(axis=0) # 归一化
m1 = GmmHmm(2)
m1.fit(t1)
m2 = GmmHmm(2)
m2.fit(t2)
m1t1 = m1.transform(t1)
m2t1 = m2.transform(t1)
print("Likelihoods for test set 1")
print("M1:", m1t1)
print("M2:", m2t1)
print("Prediction for test set 1")
print("Model", np.argmax([m1t1, m2t1]) + 1)
print()
m1t2 = m1.transform(t2)
m2t2 = m2.transform(t2)
print("Likelihoods for test set 2")
print("M1:", m1t2)
print("M2:", m2t2)
print("Prediction for test set 2")
print("Model", np.argmax([m1t2, m2t2]) + 1)
对于样本t1
而言,它更符合模型M1
;对于样本t2
而言,它更符合M2
。
现在所有的准备工作都已经就绪,就只差利用实际语音样本训练模型并进行预测了。
接下来,我们使用原始105个样本中的80%作为训练集,20%作为测试集,进行模型训练和预测。为了保证测试集和训练集中的样本单词类别均衡,我们需要在划分数据集时,向参数stratify
传入每个样本单词的标记all_labels
。
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(all_obs,
all_labels,
train_size=0.8,
test_size=0.2,
random_state =666,
stratify=all_labels)
# 由于HMM得到的结果是概率,所以我们需要对特征进行归一化
for n,i in enumerate(all_obs):
all_obs[n] /= all_obs[n].sum(axis=0)
print('Size of training matrix:', X_train.shape)
print('Size of testing matrix:', X_test.shape)
为了预测每个语音中的单词,我们需要训练7个独立的GMM-HMM模型,每一个模型对应一个单词。然后我们将测试语音样本中的特征提取出来,“喂”给7个模型,选择输出相似度最高的GMM-HMM对应的单词作为预测结果。
ys = set(all_labels)
ms = [GmmHmm(6) for y in ys]
_ = [m.fit(X_train[y_train == y, :, :]) for m, y in zip(ms, ys)]
ps = [m.transform(X_test) for m in ms]
res = np.vstack(ps)
predicted_labels = np.argmax(res, axis=0)
missed = (predicted_labels != y_test)
print('Test accuracy: %.2f percent' % (100 * (1 - np.mean(missed))))
看起来准确率较高,不算太好,毕竟这是一个简单的单音源任务。下面我们来看看混淆矩阵的情况:
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_test, predicted_labels)
f,ax = plt.subplots(figsize=(8,6))
sns.heatmap(cm,
annot=True,
cmap='summer',
linewidths=1.3)
ax.set_xticklabels(spoken)
ax.set_yticklabels(spoken,rotation='horizontal')
plt.title('Confusion matrix, single speaker')
plt.ylabel('True label')
plt.xlabel('Predicted label')
经过了艰辛的探索,我们终于可以完整的提取语音波形的峰值特征,并用于GMM-HMM模型的训练,并最终实现单音源语音的识别,你可以尝试用自己的声音录几个词语的语音文件,然后使用本案例的步骤生成一个识别你自己语音的语音识别模型。在未来,我们可以尝试用MFCCs或者使用深度神经网络来处理语音识别的问题!