梅尔倒谱系数（MFCC）做为语音信号处理中的常见特征之一，在各个语音任务中都取得了不错的效果，本文将讲解一下MFCC特征的通用的提取过程。与大家一起探讨一下~
预设连续语音为 
 共含有107000个点；预加重系数 
帧长为 
 ,帧移 
 ；窗函数为 
 ，fft的个数为 
 ，采样频率为 
预加重、分帧、加窗
预加重
预加重的目的是为了去除口唇辐射的影响，增加语音的高频分辨率，其具体实现可以用差分方程实现：

预加重后，语音信号变化为
 ，数组大小还是107000。 其中预加重系数
分帧
为了保证输入信号是平稳的，我们需要将语音分成一小段(帧)，也就是分帧。帧数为 
 ,则其
numframes = 1 + int(math.ceil((1.0 * len(signal) - frame_len) / frame_step))
对连续语音 
 以
个采样点为一帧的长度，以
个采样点为相邻帧移动的距离，对不足一帧长的进行补零，最终我们可以获得 
 个帧。
加窗
加窗是为了解决由于信号的非周期截断，导致频谱在整个频带内发生了拖尾现象的泄漏问题，可以使得使全局更加连续，避免出现吉布斯效应。窗函数一般为“汉明”窗、‘’‘汉宁’窗、矩形窗。
加窗相当于把每一语音帧里面的抽样点与窗函数中对应元素的相乘。设每一帧为 
  则加窗操作为:
for i in range(frame_len):
    x[i] = x[i]*w[i]
快速傅里叶变换 FFT
离散傅里叶变化的公式如下：
 ，其中 
为满足FFT“分治”的策略，我们需要对帧长 
 进行补零，使其补零后的长度 
 为2的幂次方；比如当原始帧长为 
 ，我们对其进行补零后的长度为 
 ，其中 
 。
在采样频率 为
的情况下，
的n只是一个离散的数值，那每一点到底代表了什么频率分量？
假设采样频率是
，则
，所以事实上把
的式子里面的 
 用
替换。很明显，在特定的
点的频率分量的频率为
，从而得出
代表的是频率为
的频率“分量”
因为傅里叶变化固有的性质：实数信号变换的结果
是一组复数，里面一半数据和另一半是共轭的，因此有一半的FFT是多余的，所以经过的FFT输出的维度为 
 ,（实数和虚数不分开表示的情况下）经过FFT后，计算频谱图的功率谱。
Mel滤波器组
MFCC的分析着眼于人耳的听觉特征，人耳所听到的声音高低与声音的频率并不成线性正比关系，而用Mel频率尺度更符合人耳的听觉特性。Mel频率与实际频率的具体关系为如下：
其实现过程如下:
根据上述公式将实际频率尺度转化为Mel频率尺度
在Mel频率域上确定最低频率 
 、与最高频率 
 和Mel滤波器个数 
 。将每一个三角形滤波器的中心频率 
 在Mel频率与上等间隔分配。设 
 分别为第 
 个三角形滤波器的下限、中心频率、上限，则： 
每个三角形滤波器 
 为：
4.根据语音信号的功率谱 
 ，求每一个三角形滤波器的输出：
注意上述式子的采样点“k”与实际频率f的对应关系为： 
倒谱分析
语音信号的倒谱分析就是求信号倒谱特征参数的过程，可以通过同态处理来处理。同态处理实现了卷积关系变化为求和关系的分离操作。
上图是一个语音的频率图，其中红色细线为频谱的包络（Spectral Envelope），包络携带声音的辨识属性，特别重要，因此我们需要把它提取出来。为了将其提取出来，我可以使用同态处理的方法：
我们可以将频谱看成频谱包络（Spectral Envelope）与频谱细节（Spectral details）的乘积。因此我们对其求log，将其转化为加性关系
其中频谱包络（Spectral Envelope）可以视为“低频信号”、频谱细节（Spectral details）可以识别“高频信号”
通过DCT将其分离，我们称为“倒谱”，对DCT结果保存前若干维的低频信息即为频谱包络（Spectral Envelope）
Delta与二阶Delta
MFCC只是描述了一帧语音上的谱包络，但是语音信号似乎有一些动态上的信息，因此就引入了一阶差分(deltas)和二阶差分(deltas-deltas),其计算公式如下：
其中t表示第几帧，N通常取2，c指的就是MFCC中的某个系数。deltas-deltas就是在deltas上再计算以此deltas。
引用
furious：深入理解离散傅里叶变换(DFT)​zhuanlan.zhihu.com
Anssel：AI(I)语音(I):MFCC特征参数提取​zhuanlan.zhihu.com
倒谱分析（Cepstrum Analysis）​blog.csdn.net
《语音信号处理》赵力著

HTK以及My_htk数据链接： https://pan.baidu.com/s/1Ajo7d-odrRiAwmCB_CQTzQ 提取码：hqnv
一：文件准备

HTK 和 HTK–samples

下载 HTK 和 HTK–samples 两个压缩文件，保存至 F 盘根目录下。 下载地址：http://htk.eng.cam.ac.uk/download.shtml（注：HTK 官方 网站需要先注册，才能下载）
Visual Studio

安装 Visual Studio 软件, 我安装的是 VS2013

二： HTK 编译
1.创建 bin.win32 和运行 VS2013
进入到解压后的 HTK-3.4.1 文件夹下的 htk 文件夹中，在本目录下创 建一个名为 bin.win32 的文件夹，这个文件夹是用来存放 htk 生成的各个 exe 程序的。找到自己 VS 安装目录下 VC 的子目录下的 bin 文件夹，将该路径添加到 path 变量中。

在cmd中输入 ：path=%path%;F:\VS2013\VC\BIN即可 （ 注意不要加入不必要的空格 ） ， 添加完成后继续输入 path， 看看有没有添加成功 （ 末尾是否有刚添加的路径 ） ， 如果添加成功 ， 输入VCVARS32


2.进行编译

编译 HTKLib

进入到HTKLib目录 ：cd HTKLib， 在命令行里输入nmake/f htk_htklib_nt.mkf all 进行编译 ， 输入后 ，DOS窗口会出现如下情况 （ 截取其中部分 ）：


编译 HTKTools

在命令行中输入 cd…， 表示返回上一层目录 ， 再输入 cd HTKTools， 进入HTKTools目录 ， 输入nmake/f htk_htktools_nt.mkf all 编译该目录下的文件 。




编译 HLMLib

输入cd…， 再输入cd HLMLib， 进入HLMLib目录 ， 输入 ：nmake/f htk_hlmlib_nt.mkf all 编译该目录下的文件 。


编译 HLMTools

输入cd…， 再输入cd HLMTools， 进入HLMTools目录 ， 输入 ：nmake/f htk_hlmtools_nt.mkf all 编译该目录下的文件 。


3.编译验证
这时候我们已经编译完所有 exe 程序，我们打开 bin. win32 文件夹看 一下，如果如下图所示，则证明已经编译成功。



4.工具路径加入到 path

为了能在DOS命令行中使用我们编译生成的工具 ， 要将生成的工具路径加入到 path中。 在cmd中输入path=%path%;F:\htk\bin.win32回车 ， 再输入path 查看添加后的结果如下 ：

【 特别提醒 】： 最好在同一个DOS窗口实现所有步骤 ， 否则可能不成功 。 如果下次需要特征提取 ， 需要将path=%path%;F:\htk\bin.win32 回车添加进去 。


三： HTK 工具箱的使用（将 WAV 生成 MFCC 文件）
1. htk 文件准备

需要准备如下两个文件（再同一文件夹下新建两个 txt 文件，把文件内 容写进去后，保存，再修改文件名即可）：



hcopy.conf

其中：hcopy.conf 文件是配置文件, 将下列内容粘贴进去即可，有可以 根据要求修改参数：

SOURCEFORMAT   = WAV  
TARGETKIND     = MFCC_0_D_A  
TARGETRATE     = 100000.0     ##10000 = 10000*100ns = 1ms
WINDOWSIZE     = 250000.0    
NUMCEPS        = 12
PREEMCOEF      = 0.97
NUMCHANS       = 27                  #定义美尔频谱的频道数量 
CEPLIFTER      = 22      #定义倒谱所用到的滤波器组内滤波器个数。



hcopy.scp

第二个文件是输入输出文件：hcopy.scp （标明语音文件的地址和对应 输出 mfcc 文件的地址, 提前建好一个 mfcc 文件夹)。hcopy.scp 文件 内容如下：


ubm

ubm 文件里面就是需要提取特征的.wav 文件

2.mfcc 特征提取
在 DOS 窗口利用 HCopy 文件进行 mfcc 特征提取
命令：HCopy -A -D -C hcopy.conf -S hcopy.scp

我图片演示的是第二次 HTK 的用法，如果是第一次就可以在同一个 DOS 窗口下输入该命令

本文将讲解一下Kaldi的提取MFCC的源码，MFCC特征作为语音信号处理技术的常用特征之一，主要包含以下几个部分：
其中kaldi的提取的模块架构图如下
接口函数
featbin/http://compute-mfcc-feats.cc
输入：waveform---音频信号，wave_data.SampFreq() ----音频采样率，vtln_warp_local---vtln参数  
输出：features--- MFCC特征
分帧、加窗、预加重
分帧、加窗、预加重的函数接口在feat/feature-common-inl.h ,
其具体实现在feat/http://feature-window.cc 
分帧 
ExtractWindow（）
输入： sample_offset ----偏移量（一般为零） wave---音频信号  f---第几帧（帧数的下标）opts---分帧的配置参数（设置帧长、帧移动等参数） window_function---窗函数类型
输出：window---加窗后的帧信号  log_energy_pre_window---每个窗函数的log能量
作用：通过opts中设置的frame_length，frame_shift，计算每帧信号的起始点与终止点，并送入加窗函数ProcessWindow中。
预加重、加窗
输入：opts---分帧的配置参数（设置帧长、帧移动等参数） window_function---窗函数类型  window---加窗前的帧信号
输出：window---加窗后的帧信号  log_energy_pre_window---每个窗函数的log能量
作用：在ProcessWindow函数，kaldi先通过Preemphasize函数对加窗前的信号进行预加重，并通过MulElements函数，计算加窗后的帧信号window 并输出。
预加重
输入：waveform---加窗前的帧信号，preemph_coeff---预加重系数
输出：waveform---预加重后的帧信号 
作用：waveform[i]=waveform[i]-a*waveform[n-1] , a = 预加重系数;
          对于kaldi 其每帧的初始信号加重后为waveform[0]=waveform[0]-a*waveform[0]。 
 FFT  功率谱
feat/http://feature-mfcc.cc
FFT与计算功率谱的方法，各个函数介绍的很清楚，这边就在不在细说了~
其中RealFft函数为计算FFT函数，ComputePowerSpectrum函数为计算功率谱 
生成Mel banks,并计算Mel频谱
其中红框内GetMelBanks函数为生成Mel滤波器组的函数，红框内Compute函数为功率谱函数经过滤波器后的输出
生成Mel banks
feat/http://mel-computations.cc
输入：opts---mel滤波器相关参数（滤波器组的个数等）  frame_opts----分帧的配置参数（设置帧长、帧移动等参数）                                 vtln_warp_local---vtln参数  
输出：返回一个滤波器组
作用：通过对opts.low_freq与opts.high_freq进行梅尔坐标转换后，计算梅尔值域。根据滤波器组的个数num_bins，在梅尔至于建立若干个三角滤波器
其中每个滤波器中的每个点值为（如下）每个滤波器共有fft个点
 if (mel <= center_mel)
          weight = (mel - left_mel) / (center_mel - left_mel);
 else
         weight = (right_mel-mel) / (right_mel-center_mel);
计算Mel频谱
输入：power_spectrum---功率谱 
输出;mel_energies_out---mel频谱输出
作用:mel滤波器组举证与功率频谱相乘
求倒谱
Log
Dct

	

点击上方“AI大道理”，选择“置顶”公众号
重磅干货，细致入微AI大道理
 ——————
1 特征提取流程
    在语音识别和说话者识别方面，最常用到的语音特征就是梅尔倒谱系数(Mel-scaleFrequency Cepstral Coefficients，简称MFCC)。
     MFCC提取过程包括预处理、快速傅里叶变换、Mei滤波器组、对数运算、离散余弦变换、动态特征提取等步骤。
 2 快速傅里叶变换
    快速傅里叶变换即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称，简称FFT。
    FFT不是Fast FT，而是Fast DFT 
    FT的种类很多，以最简单的基于2的FFT为例。
    FFT实际上一种分治算法。FFT将长度为的信号分解成两个长度为信号进行处理，这样分解一直到最后，每一次的分解都会减少计算的次数。理解FFT分以下三个步骤进行：
步骤1：将信号分解成两个子信号
偶数样本点信号：；
奇数样本点信号：；
步骤2：将两个求和项理解成两个长度为的DFT
步骤3：FFT的具体计算过程
对于任意都要进行次加法操作，所以DFT共有次乘法操作。
对于任意都要进行次加法操作，DFT共有次加法操作。
FFT共有次乘法操作和次加法操作。
语音信号是有限长的离散信号。
预处理后的语音信号：
FFT后效果：
3 Mel滤波器组将能量谱通过一组Mel尺度的三角形滤波器组，定义一个有M个滤波器的滤波器组(滤波器的个数和临界带的个数相近)，采用的滤波器为三角滤波器，中心频率为 。M通常取22-26。各f(m)之间的间隔随着m值的减小而缩小，随着m值的增大而增宽，如图所示：
三角带通滤波器有两个主要目的：
(1)三角形是低频密、高频疏的，这可以模仿人耳在低频处分辨率高的特性；
(2)对频谱进行平滑化，并消除谐波的作用，突显原先语音的共振峰。频谱有包络和精细结构，分别对应音色与音高。对于语音识别来讲，音色是主要的有用信息，音高一般没有用。在每个三角形内积分，就可以消除精细结构，只保留音色的信息。
(3)傅里叶变换得到的序列很长(一般为几百到几千个点)，把它变换成每个三角形下的能量，可以减少数据量
Mel频率和频率f的对应关系：
或者
Mel滤波器实现过程：
(1)确定最低频率(0HZ)，最高频率(fs/2)，Mel滤波器个数M(23)；
(2)转换最低频率和最高频率的Mel(f)；
(3)计算相连两个Mel滤波器中心Mel频率的距离，在Mel频率上，两两之间的中心频率是等间距的；
(4)将各种中心Mel频率转化为频率f(非等间距)；
(5)计算频率所对应的FFT中点的下标；
灵魂的拷问：为什么有些Mel滤波器组不等高，我设计的是等高的？这样有影响吗？有优势吗？
AI大语音：不等高的原因是乘了一个递减的系数，就是实现上一些细节的差别，保证了每个滤波器的能量和一样。横轴指的频率，低频的系数高，就是对低频更加的关注 。没有太大的影响，一般主要用等高的。
经过梅尔滤波器组后的Fbank特征：
4 对数运算
将原语音信号经过傅里叶变换得到频谱：                                                        X[k]=H[k]E[k]
只考虑幅度就是：
                            |X[k] |=|H[k]||E[k] |
两边取对数：
              log||X[k] ||= log ||H[k] ||+ log ||E[k] ||
再在两边取逆傅里叶变换得到：
                                    x[k]=h[k]+e[k]
灵魂的拷问：为什么要进行对数运算？它在干嘛？
对数运算包括取绝对值和log运算。取绝对值是仅使用幅度值，忽略相位的影响，因为相位信息在语音识别中作用不大。
log运算是为了分别包络和细节，包络代表音色，细节带包音高，显然语音识别就是为了识别音色。另外，人的感知与频率的对数成正比，正好使用log模拟。
FFT变换后，卷积变成了乘法，取对数后，乘法变成了加法，把卷积信号转换成加性信号。
5 离散余弦变换(DCT)
再在两边取逆傅里叶变换得到：
                                    x[k]=h[k]+e[k]
在上一步中，我们成功地把基音信息与声道信息变成了加性的。那么如何分离呢？它们有如下性质：
频谱图中(注意是一帧FFT变换内)
(1)基音信息在频域是快速变化的。
(2)声道信息在频域是缓慢变化的。
因此再做一次DCT可以将其分离。我们称之为"倒谱域"。因此倒谱域的低频部分刻画了声道信息，高频部分刻画了基音信息。   
由此得到12维的MFCC特征：
 由于许多要处理的信号都是实信号，在使用DFT时由于傅里叶变换时由于实信号傅立叶变换的共轭对称性导致DFT后在频域中有一半的数据冗余。
    将DFT式子拆开，抽出实数部分：
则实数部分：
虚数部分：
又有：
而当x[n]是实偶信号时：
把DFT写成：
    但是实际中并没有那么多实偶信号，我们就认为造出来。将信号长度扩大成原来的两倍,并变成2N，又为了让造出来的信号关于0对称，把整个延拓的信号向右平移 0.5 个单位，最终DCT变换公式：
6 动态特征提取
标准的倒谱参数MFCC只反映了语音参数的静态特性，语音的动态特性可以用这些静态特征的差分谱来描述。实验证明：把动、静态特征结合起来才能有效提高系统的识别性能。差分参数的计算可以采用下面的公式：
 式中,dt表示第t个一阶差分，Ct表示第t个倒谱系数，Q表示倒谱系数的阶数，K表示一阶导数的时间差，可取1或2。将上式的结果再代入就可以得到二阶差分的参数。
因此，MFCC的全部组成其实是由：N维MFCC参数(N/3 MFCC系数+ N/3 一阶差分参数+ N/3 二阶差分参数)+帧能量(此项可根据需求替换)。
这里的帧能量是指一帧的音量(即能量)，也是语音的重要特征。
d_mfcc_feat = delta(wav_feature, 1)                     d_mfcc_feat2 = delta(wav_feature, 2) 
feature = np.hstack((wav_feature, d_mfcc_feat, d_mfcc_feat2))
最终39维MFCC图：
附录(魔鬼写手)
  
——————
浅谈则止，细致入微AI大道理
扫描下方“AI大道理”，选择“关注”公众号
           欢迎加入！
    
▼下期预告▼AI大语音(五)——声学模型▼往期精彩回顾▼AI大语音(一)——语音识别基础AI大语音(二)——语音预处理AI大语音(三)——傅里叶变换家族
                                                                          留你心，言你想