【MFCC梅尔倒频谱参数】浅谈语音特征参数MFCC提取过程

一、概述

• 在语音识别（Speech Recognition）和话者识别（Speaker
  Recognition）方面，最常用到的语音特征就是梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients，简称MFCC）。根据人耳听觉机理的研究发现，人耳对不同频率的声波有不同的听觉敏感度。从200Hz到5000Hz的语音信号对语音的清晰度影响对大。两个响度不等的声音作用于人耳时，则响度较高的频率成分的存在会影响到对响度较低的频率成分的感受，使其变得不易察觉，这种现象称为掩蔽效应。由于频率较低的声音在内耳蜗基底膜上行波传递的距离大于频率较高的声音，故一般来说，低音容易掩蔽高音，而高音掩蔽低音较困难。在低频处的声音掩蔽的临界带宽较高频要小。所以，人们从低频到高频这一段频带内按临界带宽的大小由密到疏安排一组带通滤波器，对输入信号进行滤波。将每个带通滤波器输出的信号能量作为信号的基本特征，对此特征经过进一步处理后就可以作为语音的输入特征。由于这种特征不依赖于信号的性质，对输入信号不做任何的假设和限制，又利用了听觉模型的研究成果。因此，这种参数比基于声道模型的LPCC相比具有更好的鲁邦性，更符合人耳的听觉特性，而且当信噪比降低时仍然具有较好的识别性能。

• 梅尔倒谱系数（Mel-scale Frequency Cepstral
  Coefficients，简称MFCC）是在Mel标度频率域提取出来的倒谱参数，Mel标度描述了人耳频率的非线性特性，它与频率的关系可用下式近似表示：
  
  下图是梅尔频率关系图
  

• 倒谱(Cepstrum): 信号的傅里叶变化经过对数运算后在进行逆傅里叶变换得到的谱
  

• 对滤波器组的输出使用离散余弦变换(DCT)去除相关性的倒谱系数c(n):
  

二、提取过程

上图是MFCC的流程图。
以下的代码叙述都是以Python语言作为例子

提取语音信号

1. 提取语音信号
   一般python提取音频信号可以使用标准库scipy拿到语音信号值。

def read_wav(audio_path):
    """
    使用scipy库下的方法读取音频文件wav
    时间序列y / 采样率sr = 音频时长
    
    :param audio_path: 音频路径
    :return: 返回时间坐标，时间序列y（数据类型和声道数由文件本身决定）和 采样率sr(赫兹)
    """
    from scipy.io.wavfile import read
    from numpy import arange
    sr, y = read(filename=audio_path)  # 读取音频文件，返回音频采样率和时间序列
    x = arange(0, len(y)/sr, 1/sr)  # 总秒数 = 总采样点/采样频率
    return x, y, sr

预加重

2. 预加重
   预加重处理其实是将语音信号通过一个高通滤波器：

   Y(n) = X(n) - μX(n-1)
   y(n) 指输入的是离散信号
   H(z) = 1 - μz^-1
   H(z) 指输入的是连续信号

   上式中μ的值介于0.9-1.0之间，我们通常取0.97。预加重的目的是提升高频部分，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱。同时，也是为了消除发生过程中声带和嘴唇的效应，来补偿语音信号受到发音系统所抑制的高频部分，也为了突出高频的共振峰。

def pre_emphasis(sinal, coefficient=0.97):
	import numpy
	return numpy.append(signal[0], signal[1:] - coefficient * signal[:-1])

分帧

3. 分帧
   先将N个采样点集合成一个观测单位，称为帧。通常情况下N的值为256或512，涵盖的时间约为20~30ms左右。为了避免相邻两帧的变化过大，因此会让两相邻帧之间有一段重叠区域，此重叠区域包含了M个取样点，通常M的值约为N的1/2或1/3。通常语音识别所采用语音信号的采样频率为8KHz或16KHz，以8KHz来说，若帧长度为256个采样点，则对应的时间长度是256/8000×1000=32ms。

def audio2frame(signal, frame_length, frame_step, winfunc=lambda x: numpy.ones((x,))):
    """将音频信号转化为帧。
	参数含义：
	signal:原始音频型号
	frame_length:每一帧的长度(这里指采样点的长度，即采样频率乘以时间间隔)
	frame_step:相邻帧的间隔（同上定义）
	winfunc:lambda函数，用于生成一个向量
    """
    signal_length = len(signal)  # 信号总长度
    frame_length = int(round(frame_length))  # 以帧帧时间长度
    frame_step = int(round(frame_step))  # 相邻帧之间的步长
    if signal_length <= frame_length:  # 若信号长度小于一个帧的长度，则帧数定义为1
        frames_num = 1
    else:  # 否则，计算帧的总长度
        frames_num = 1+int(math.ceil((1.0*signal_length-frame_length)/frame_step))
    pad_length = int((frames_num-1)*frame_step+frame_length)  # 所有帧加起来总的铺平后的长度, 最后一帧长度当满长度处理,帧与帧之间的相同部分不算
    zeros = numpy.zeros((pad_length-signal_length,))  # 不够的长度使用0填补，类似于FFT中的扩充数组操作
    pad_signal = numpy.concatenate((signal, zeros))  # numpy的矩阵拼接, 填补后的信号记为pad_signal
    # numpy.tile(x, y) 生成新数组，新数组为复制y次x数组->一行中有y个x的数组； y=[x0, y0]可为数组时->有x0行的x数组，y0列的x数组
    indices = numpy.tile(numpy.arange(0, frame_length), (frames_num, 1))+numpy.tile(numpy.arange(0, frames_num * frame_step, frame_step), (frame_length, 1)).T  # 相当于对所有帧的时间点进行抽取，得到frames_num*frame_length长度的矩阵
    indices = numpy.array(indices, dtype=numpy.int32)  # 将indices转化为矩阵
    frames = pad_signal[indices]  # 得到帧信号
    win = numpy.tile(winfunc(frame_length), (frames_num, 1))  # window窗函数，这里默认取1
    return frames*win   # 返回帧信号矩阵

• 下面来详解以下这个算法, 本人会以表格形式来展现信号，设一个信号长度为11，时间序列值为1-11的信号，并且令帧长度为4，帧间隔为2.
  ∵帧间隔一般是帧长度的1/2或1/3,在这里取1/2
  (1) 依照算法，123行代码所得到的值依次是11、4、2
  (2) 接下来通过判断信号总长度与帧长度大小：
  如果信号长度小于帧长度，说明一个帧就足够了，帧总数赋值为1。否则计算帧。
  

  前面的1的作用是将最后一帧的长度若不够帧长度的则按满长度处理，因为除法得到的结果是小数，则进行取整，这会导致缺少一帧，固先把最后一帧当作满处理，再通过减一个帧长度，就得到的总帧数。在这里得到5
  (3) 计算处理后的总长度，就是每帧按满帧处理，其中帧与帧的重复部分只取一次,所以在这里得到值12
  (4) 通过计算处理后的长度与原信号长度的差，拿到需要补0的数目。在这里长度差为12 - 11 = 1，需要补一个0，其中在这里先将信号转为矩阵，再向信号矩阵后面添加0矩阵。
  
  (5) indices的生成代码:
  第一part: 数组[0, …, 帧长度] 的数组扩充至 [帧总数，1]的矩阵。 在这里第一部分的矩阵为:
  
  第二部分：数组生成规则是，从0到帧总数*帧步长的范围内，每隔一个帧步长就取一个数。并将该数组复制至[帧长度，1]的矩阵最后取逆矩阵。在这里第二部分的矩阵为:
  
  取逆矩阵后为：
  
  最后第一部分与第二部分进行矩阵加法:
  
  在这里就显而易见了，每一行代表一个帧序列，一共有5个帧，且每一个帧为4个帧长度，帧的值为信号矩阵的索引值，即signal[0]， signal[1]， …
  (6) 将帧生成二维数组转化为矩阵
  (7) 将拿到的帧矩阵通过signal的值转化为帧信号
  
  (8) 窗矩阵生成 (一般选用的窗函数是汉明窗，下面有详解介绍)
  默认不使用窗函数，在这里即窗函数为[1, 1, 1, 1]
  再将窗函数复制扩充至[帧总数行，帧长度列]
  在这里即为
  
  若窗函数为[2, 2, 2, 2]，则意义是将每一帧的信号扩大1倍
  (9) 使用窗矩阵点乘帧信号矩阵，对每一个帧都通过窗函数进行处理。
  (10) 最后返回处理完的帧信号矩阵
  在这里，最终返回的信号是
  

加窗

4. 加窗(默认是汉明窗)
   将每一帧乘以汉明窗，以增加帧左端和右端的连续性。假设分帧后的信号为S(n), n=0,1…,N-1, N为帧数目，那么乘上汉明窗后,W(n)形式如下：
   
   不同的a值会产生不同的汉明窗，一般情况下a取0.46.

def hanming(n, a, N):
    import numpy as np
    return (1 - a) - a * np.cos((2 * np.pi * n) / (N - 1))

# 汉明窗展示
from matplotlib.pyplot import plot, show


def hanming(n, a, N):
    import numpy as np
    return (1 - a) - a * np.cos((2 * np.pi * n) / (N - 1))


N = 100
a = 0.46
x = linspace(0, 1, N)	# 频率为N
y = [hanming(i, a, N) for i in range(N)]
plot(x, y)
show()

• N越大，曲线越平滑。随着a的变化，趋势如下图.
  

  汉明窗的特性是窗内中间的数据极大程度被保留，而边缘的两边被削弱。所以有上面的帧步长为1/2的帧长度，将窗内边缘的部分通过下个帧表现出来。
  而为什么要加窗呢？
  在信号处理中，可以说加窗处理是一个必经的过程，因为我们的计算机只能处理有限长度的信号，因此原始信号X(t)要以T（采样时间）截断，即有限化，成为XT(t)后再进一步处理，这个过程序就是加窗处理，但什么时候用什么窗呢？这时我们就要对所需用到的函数窗做一定的了解。在平时，我们用得最多的是矩形窗。实际的信号处理过程中，我们用的矩形窗在边缘处将信号突然截断时，窗外时域信息全部消失，导致在频域增加了频率分量的现象，即频谱泄漏。避免泄漏的最佳方法是满足整周期采样条件，但实际中是不可能做到的。对于非整周期采样的情况，必须考虑如何减少加窗时造成的泄漏误差，主要的措施是使用合理的加窗函数，使信号截断的锐角钝化，从而使频谱的扩散减到最少。

  首先介绍一下为什么要用函数窗：函数窗的主要用于对截断处的不连续变化进行平滑，减少泄漏。此外，加窗处理还有很多其它的原因，如减少噪声干扰、限定测试的持续时间、从频率接近的信号中分离出幅值不同的信号……

  常见的几种窗的基本指标：
  
  
  下面是通过生成一个不同长度的余弦函数，对余弦函数进行快速傅里叶变换以及对加了汉明窗的余弦函数进行快速傅里叶变换，进行对比。在这里，长度取值N我分别取值了10、100、500、1000，步长为0.1.(在这里我并没有对x坐标进行频率变换)
  每一行代表的是同一个长度的余弦函数，第一列是原函数，第二列是原函数的快速傅里叶变换图，第三列是加汉明窗后的函数，第四列是加汉明窗后函数的快速傅里叶变换。
  
  在这里x坐标轴我没有做频率的变换，这是演示一下图像在加窗前后的傅里叶变换的差别。
  通过上图可以看到，在n等于10的时候，图像已经完全失去了cos的频率特性，所以在原函数的FFT图像中，频谱已经失真了，但是经过加窗后，频谱竟然有些修复。(cos函数的频谱图不清楚的建议去看一下信号与系统这本书或者百度一下)。n=100时，其频谱图也有些失真，但是经过加窗后，失真现象也缓和了。到n=500, n=1000时，因为时间序列已经比较长了，所以原函数的FFT的图像与原FFT图像尽管有些误差，但是误差已经可以算是可以允许了。这说明，汉明窗的确是能缓和频谱泄露现象。

快速傅里叶变换

5. 快速傅里叶变换
   由于信号在时域上的变换通常很难看出信号的特性，所以通常将它转换为频域上的能量分布来观察，不同的能量分布，就能代表不同语音的特性，这一转换可以通过傅里叶变换进行。所以在乘上汉明窗后，每帧还必须再经过快速傅里叶变换以得到在频谱上的能量分布，快速傅里叶变换与傅里叶变换的区别可以看作是运算量的降低。对分帧加窗后的各帧信号进行快速傅里叶变换得到各帧的频谱后，由于傅里叶变换后的值是含有虚数部分的，所以对语音信号的频谱取模平方得到语音信号的功率谱。设语音信号的DFT为：
   式中x(n)为输入的语音信号，N表示傅里叶变换的点数。
   代码为：

from scipy import fft
result = fft(signal)	# signal是指信号的时间序列

梅尔滤波器组

6. 三角带通滤波器
   将能量谱通过一组Mel尺度的三角形滤波器组，定义一个有M个滤波器的滤波器组（滤波器的个数和临界 带的个数相近），采用的滤波器为三角滤波器，中心频率为f(m),m=1,2,…,M。M通常取22-26。各f(m)之间的间隔随着m值的减小而缩小，随着m值的增大而增宽，如图所示：
   
   可以从三角波函数可以看出，随着m的增大，幅值越来越小，这从频域上看，就是显低频，越到高频的部分越衰减的厉害。
   三角滤波器的频率响应定义为：
   
   三角带通滤波器有两个主要目的：
   对频谱进行平滑化，并消除高次谐波的作用(高次谐波也就是高频分量)，突显原先语音的共振峰。（因此一段语音的音调或音高，是不会呈现在 MFCC 参数内，换句话说，以 MFCC 为特征的语音辨识系统，并不会受到输入语音的音调不同而有所影响） 此外，还可以降低运算量。

对数能量

7. 计算每个滤波器组输出的对数能量
   
   X(t)是信号序列， Hm是梅尔滤波器组
   对数能量：一帧内信号的平方和，再取以10为底的对数值，再乘以10
   此外，一帧的音量（即能量），也是语音的重要特征。

   注：若要加入其它语音特征以测试识别率，也可以在此阶段加入，这些常用的其它语音特征包含音高、过零率以及共振峰等。

离散余弦变换（DCT）

8. 离散余弦变换（DCT）得到MFCC系数
   
   将上述的对数能量带入离散余弦变换，求出L阶的Mel-scale Cepstrum参数。L阶指MFCC系数阶数，通常取12-16。这里m是三角滤波器的阶数。

动态差分参数的提取（包括一阶差分和二阶差分）

9. 动态差分参数的提取（包括一阶差分和二阶差分）
   标准的倒谱参数MFCC只反映了语音参数的静态特性，语音的动态特性可以用这些静态特征的差分谱来描述。实验证明：把动、静态特征结合起来才能有效提高系统的识别性能。差分参数的计算可以采用下面的公式：
   
   式中,dt表示第t个一阶差分，Ct表示第t个倒谱系数，Q表示倒谱系数的阶数，K表示一阶导数的时间差，可取1或2。将上式的结果再代入就可以得到二阶差分的参数。

参考文献或博客：
[1]说话人识别中改进特征提取算法的研究[J].宋乐.太原理工大学
[2]基于时频分布与MFCC的说话人识别[J]. 金银燕,于凤芹,何艳. 计算机系统应用. 2012(04)
[3]语音识别之——mfcc什么是汉明窗，为什么加汉明窗[CSDN博客].gxiaoyaya
[4]语音识别的第一步MFCC特征提取代码（Python）[ITKeyword社区].断桥残雪断桥残雪
[5]语音识别之MFCC特征提取[百度文库].汤旭国
[6]MFCC [百度文库]