mfcc svm



大家好，又见面了，我是你们的朋友全栈君。

史上最详细的MFCC算法实现（附测试数据）

这里该包的安装我直接附上我们师姐写过的一篇文章，里边的介绍很详细：

整个MFCC过程大致可以分为以下几步：

1.音频文件读取（最好是.wav文件） 2.预先加重 3.分帧 4.加窗 5.傅里叶变换（当是2的N次方时，可以使用FFT快速傅里叶变换） 6.梅尔滤波器组 7.离散余弦变换DCT

uad 在matlab中,使用函数audioread函数来读取本地wav文件，这里要注意的是，采样频率一般为8000Hz和16000Hz，采样频率需要大于真实信号最大频率的2倍，才不会导致频谱混叠。

这里我们用于测试的数据的采样频率 f s f_s fs​ 44100,这个是由 audioread函数默认决定的。

uad 为了避免在后边的FFT操作中出现数值问题，我们需要加强一下高频信息，因为一般高频能量比低频小。其预加重函数如下所示： y ( n ) = x ( n ) − α ⋅ x ( n − 1 ) y(n)=x(n)-alpha cdot x(n-1) y(n)=x(n)−α⋅x(n−1) 其中 α alpha α一般取值为 0.97 、 0.95 0.97、0.95 0.97、0.95

然后我们再来对比一下原始文件和预加重后的数据差异

我们可以看到整个的数据其幅度范围是有所减小的，但是可以看得出来，高频部分的缩小倍数是小于低频部分的缩小倍数的。

uad 我们要对语音数据做傅里叶变换，将信息从时域转化为频域。但是如果对整段语音做FFT，就会损失时序信息。因此，我们假设在很短的一段时间t内的频率信息不变，对长度为t的帧做傅里叶变换，就能得到对语音数据的频域和时域信息的适当表达。例如我们这里的采样点数为个点，如果真的这样做的话，就很麻烦了，于是我们在语音分析中引入分帧的概念，将原始语音信号分成大小固定的N段语音信号，这里每一段语音信号都被称为一帧。 uad 但是，如果我们这样分帧的话，帧与帧之间的连贯性就会变差，于是我们每一帧的前N个采样点数据与前一帧的后N个采样点数据一样。其原理图大致如下所示：

uad 对于整个采样点数据可以分为多少帧以及帧与帧之间交叉的采样点个数N，不是随便分的，一般来说帧长设置为 25 m s 25ms 25ms,帧移设置为 10 m s 10ms 10ms,对于我这次的仿真，其帧数和帧长数值如下： 帧 数 = f s ⋅ 0.025 = 44100 ⋅ 0.025 = 1103 ( 个 采 样 点 ) 帧 移 = f s ⋅ 0.001 = 44100 ⋅ 0.01 = 441 ( 个 采 样 点 ) quad 帧数=f_s cdot 0.025=44100 cdot0.025=1103(个采样点)\ quad \ quad 帧移=f_s cdot 0.001=44100 cdot 0.01=441(个采样点) 帧数=fs​⋅0.025=44100⋅0.025=1103(个采样点)帧移=fs​⋅0.001=44100⋅0.01=441(个采样点) uad 在这里我们要调用matlab的enframe函数来进行分帧操作，要知道这个函数是包含在voicebox工具箱里边，首先确保其已经安装成功。

uad 将信号分帧后,我们将每一帧代入窗函数，窗外的值设定为0，其目的是消除各个帧两端可能会造成的信号不连续性（即谱泄露 spectral leakage）。常用的窗函数有方窗、汉明窗和汉宁窗等，根据窗函数的频域特性，常采用汉明窗（hamming window）。接下来我来讲解一下怎么加窗：我们需要做的就是为每一帧数据，也就是301帧数据都加入大小为1103的汉明窗。其汉明窗的表达公式如下所示： W ( n ) = ( 1 − a ) − a ⋅ c o s ( 2 ⋅ π ⋅ n / N ) 1 ≤ n ≤ N W(n)=(1-a)-a cdot cos(2cdot pi cdot n/N) uad 1 leq n leq N W(n)=(1−a)−a⋅cos(2⋅π⋅n/N)1≤n≤N 对于a的取值不同，将会产生不同的汉明窗，一般情况下， a = 0.46 a=0.46 a=0.46

uad 由上边的公式我们可以得到汉明窗矩阵 C ( 301 , 1103 ) C_{(301,1103)} C(301,1103)​,其矩阵大小为(301,1103)，由于汉明窗矩阵和分帧后的矩阵S具有相同大小，所以在matlab中使这两个矩阵的对应位置相乘，即可得到加窗后的矩阵 S C ( 301 , 1103 ) SC_{(301,1103)} SC(301,1103)​，其大小为(301,1103)。接下来我将随便选取一帧数据来展示一下汉明窗、原始数据、加窗后的数据。其matlab代码如下所示：

uad 在上边的图示中我们就可以看到，在每一帧的低频部分和高频部分都被汉明窗相乘后起了较大抑制作用，使其结果接近于0。

uad 对于加窗后的矩阵 S C ( 301 , 1103 ) SC_{(301,1103)} SC(301,1103)​,它是一个301*1103的矩阵，也就是说，它有301帧数据，且每一帧数据都有1103个采样点，那么我们接下来就要对这301帧的每一帧都要进行N=4096的FFT快速傅里叶变换,得到一个大小为(301,4096)大小的矩阵 D ( 301 , 4096 ) D_{(301,4096)} D(301,4096)​,其帧数还是301帧，对每一帧的4096个数据点分别取模再取平方，然后除以4096；便得到能量谱密度 E ( 301 , 4096 ) E_{(301,4096)} E(301,4096)​,其大小为301×4096,然后再对每一帧得到的能量进行相加，即得到一个301×1的矩阵 F ( 301 , 1 ) F_{(301,1)} F(301,1)​,其中的每个元素代表每一帧能量的总和。

uad 首先我要讲一下什么是梅尔值，这是一个新的量度，相比于正常的频率机制，梅尔值更加接近于人耳的听觉机制，其在低频范围内增长速度很快，但在高频范围内，梅尔值的增长速度很慢。每一个频率值都对应着一个梅尔值，其对应关系如下 m = 2595 ⋅ l o g 10 ( 1 + f 700 ) m=2595 cdot log_{10}(1+frac{f}{700}) m=2595⋅log10​(1+700f​) 在matlab上画出其对应图像如下所示：

对于该函数图像的显示的matlab代码如下所示:

但是如果我们要将梅尔频率m转换为频率f呢，我们对上式整理即可得到： f = 700 ⋅ ( 1 0 m / 2595 − 1 ) f=700 cdot (10^{m/2595}-1) f=700⋅(10m/2595−1) uad 好了介绍到这里，对于如何得到梅尔滤波器组我们还是无从下手，于是，我在这里描述了一下获得梅尔滤波器的几个简单步骤。然后接下来的操作我们也就将会按照这个步骤来实现。

其中过程1、2、3、4的实现代码如下所示：

上边几步都比较好理解，但是对于接下来谱线索引号k的定义，或许就需要一些理解了，其定义公式如下所示： k = ( 1 + N ) ⋅ f m f s k=frac{(1+N)cdot f_m}{f_s} k=fs​(1+N)⋅fm​​ uad 其中 N N N为FFT点数， f s f_s fs​为抽样频率， f m f_m fm​为之前那28个梅尔刻度转化为频率后的值，最后我们得到的 k k k值为一个1*28的矩阵。且k值范围为 0 − N / 2 0-N/2 0−N/2。这个式子是把频率对应到频谱中2048个频率分量的某个。 以下则是k值得求解过程：

好了，现在我们只剩下最后一步了，创建Hm矩阵，这个矩阵得定义如下所示： H m ( k ) = { 0 k < f ( m − 1 ) 2 ( k − f ( m − 1 ) ) ( f ( m + 1 ) − f ( m − 1 ) ) ( f ( m ) − f ( m − 1 ) ) f ( m − 1 ) ≤ k ≤ f ( m ) 2 ( f ( m + 1 ) − k ) ( f ( m + 1 ) − f ( m − 1 ) ) ( f ( m ) − f ( m − 1 ) ) f ( m ) ≤ k ≤ f ( m + 1 ) 0 k ≥ f ( m + 1 ) H_m(k)=begin{cases} 0& ext{k}<f(m-1)\ frac{2(k-f(m-1))}{(f(m+1)-f(m-1))(f(m)-f(m-1))}& f(m-1)le k le f(m)\ frac{2(f(m+1)-k)}{(f(m+1)-f(m-1))(f(m)-f(m-1))}& f(m)le k le f(m+1)\ 0& kge f(m+1) end{cases} Hm​(k)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​0(f(m+1)−f(m−1))(f(m)−f(m−1))2(k−f(m−1))​(f(m+1)−f(m−1))(f(m)−f(m−1))2(f(m+1)−k)​0​k<f(m−1)f(m−1)≤k≤f(m)f(m)≤k≤f(m+1)k≥f(m+1)​ uad 上式中的m代表第几个滤波器，k为横坐标0-2048。26个滤波器就是算2049个频率分量分别属于26个频率带的概率.根据上式计算26个滤波器的二维数组，也就是26*4096二维数组Hm. 以下贴出步骤6 Hm矩阵的求解

接下来将要进行最后一步，输出Hm矩阵，并且将梅尔滤波器组画出来。

uad 画出来之后，我们就会发现该梅尔滤波器，在频率很小的时候，滤波器宽度很窄，随着其频率的增大，滤波器的宽度也会逐渐增大，但其幅值也会逐渐减小。

uad 在进行离散余弦变换之前，我们还需要做的就是把第3.5节得到的二维矩阵能量谱 E ( 301 , 4096 ) E_{(301,4096)} E(301,4096)​，乘以第3.6节得到的二维数组梅尔滤波器 H m ( 26 ∗ 4096 ) Hm_{(26*4096)} Hm(26∗4096)​的转置，矩阵的转置可得到301*26的矩阵，然后满足矩阵乘法定律，得到参数H,其定义如下： H = E ⋅ H m T H=E cdot H_m^T H=E⋅HmT​ 此处的H其实是301×26的二维矩阵。 uad 由于滤波器组得到的系数是相关性很高的，因此我们用离散余弦变换（Discrete Cosine Transform）来去相关并且降维。一般来说，在自动语音识别（Automatic Speech Recognition）领域，因为大部分信号数据一般集中在变换后的低频区，所以对每一帧只取前13个数据就好了。 好了接下来我们就要进行离散余弦变换了，但是在开始之前我感觉还是先讲解一下其具体的步骤流程吧。

根据mfcc的定义，我们需要对能量的对数作离散余弦变换，即可得到MFCC参数： m f c c ( i , n ) = ∑ m = 1 M l o g [ H ( i , m ) ] ⋅ c o s [ π ⋅ n ⋅ ( 2 m − 1 ) 2 M ] mfcc(i,n)=sum_{m=1}^{M}log[H(i,m)] cdot cos[frac{pi cdot n cdot(2m-1)}{2M}] mfcc(i,n)=m=1∑M​log[H(i,m)]⋅cos[2Mπ⋅n⋅(2m−1)​] uad 其中H为我们上边得到的矩阵H,M代表梅尔滤波器个数，i代表第几帧数据（取值为1-301），n代表第i帧的第n列(n取值范围为1-26)。那么根据上述公式我们可以写出求取mfcc的代码如下：

uad 接下来我们就要根据公式进行升倒谱的计算了，前边我们已经讲到了，因为大部分的信号数据一般集中在变换后的低频区，所以对每一帧只取前13个数据就好了。其公式表达如下： K ( i ) = 1 + ( L 2 ) ⋅ s i n ( π ⋅ i L ) i = 1 , 2 , 3… , 13 K(i)=1+(frac{L}{2})cdot sin(frac{pi cdot i}{L}) uad i=1,2,3…,13 K(i)=1+(2L​)⋅sin(Lπ⋅i​)i=1,2,3...,13 uad 其中L为升倒谱系数，一般取值为22，我们将其带入即可求得矩阵K,这是一个1×13大小的矩阵，这一部分的升倒谱的其实现代码如下：

接下来我们就要求取MFCC的三个参数了，首先我们先获取mfcc的第一组数据，根据公式： f e a t ( i , j ) = J ( i , j ) ⋅ K ( j ) i = 1 , 2 , 3 , . . . , 301 ; j = 1 , 2 , . . . , 13 feat(i,j)=J(i,j) cdot K(j)\ quad \ quad i=1,2,3,…,301; quad j=1,2,…,13 feat(i,j)=J(i,j)⋅K(j)i=1,2,3,...,301;j=1,2,...,13 根据公式我们可以实现如下代码：

uad 接下来就是求取MFCC的第二组，第三组参数了。第二组参数其实就是在已有的基础参数下作一阶微分操作，第三组参数在第二组参数下作一阶微分操作。那么表达公式事什么样的呢，别急，等我慢慢道来： d t f e a t [ i ] [ j ] = f e a t [ i + 1 ] [ j ] − f e a t [ i − 1 ] [ j ] + 2 ⋅ f e a t [ i + 2 ] [ j ] − 2 ⋅ f e a t [ i − 2 ] [ j ] i = 1 ， 2 ， 3 ， . . . , 301 ; j = 1 , 2 , 3… , 13 dtfeat[i][j]=feat[i+1][j]-feat[i-1][j]+2cdot feat[i+2][j]-2 cdot feat[i-2][j] \ quad \ i=1，2，3，…,301 quad; quad j=1,2,3…,13 dtfeat[i][j]=feat[i+1][j]−feat[i−1][j]+2⋅feat[i+2][j]−2⋅feat[i−2][j]i=1，2，3，...,301;j=1,2,3...,13 按照上边的公式，我们可以使用代码实现一阶求导和二阶求导的计算

uad 好了到这里我们就完成了，MFCC三组参数的求解，我们现在就只需要将这三组数据拼接到一起形成一个301×39的矩阵即可。其实现代码如下：

uad 前面导图中，我们也有讲到过，由于一阶求导和二阶求导的前两帧和后两帧数据都为0，于是我们就要对在mfcc_final中去除这四帧数据。而后我们再选取每一帧的mfcc系数的第一个数得到 M F C C 0 MFCC_0 MFCC0​,这是一个297×1的数据，对 M F C C 0 MFCC_0 MFCC0​来进行绘图，并与原始信号进行比对。

比对结果如下所示：

uad 好了，到了这里我们就可以看到了，原始信号之前是20000个采样点的数据，而现在的 M F C C 0 MFCC_0 MFCC0​参数图形大致与原始信号一致，并且其点数只有297个点，这也就说明通过此方法 M F C C 0 MFCC_0 MFCC0​，我们可以提取出语音信号的特点以及走向趋势，也就是说在某个程度上我们可以用这297个点来代替 2 ⋅ 1 0 5 2 cdot 10^5 2⋅105点的数据。

上边图示 M F C C 0 MFCC_0 MFCC0​仅为语音MFCC参数的第一维参数，因其（ M F C C 0 MFCC_0 MFCC0​）包含了语音信号的时域能量信息，所以一般被用来用作语音信号的 端点检测。另外还有很多同学对变量的理解不是那么明白，这里特此说明一下：

最后还是要说明一下，此篇博客只是为了让大家深入理解MFCC特征的求取过程，所以其在 运行速度 和 代码规范性 方面未作过多优化，如果大家是用于做实验研究，建议大家可以直接调用matlab自带的mfcc函数（matlab2019版以后已经自带mfcc的函数），

那如果你的matlab版本较低，则可以参考4.1节代码。

本次训练是在参考了很多资料的前提下完成的，为了防止自己忘记，所以特此写了本篇文章。

很多同学想要这个mfcc版本的完整版，再加之上边的讲解是只针对自己语音，一些语音长度，包括语音帧数使得部分同学难以理解，所以在这里我又将程序进行了一下修改，使同学们可以自己输入自己的音频文件，并返回mfccs参数矩阵：

上边封装的函数是我根据上边分步讲解的内容，封装为Mymfcc函数，调用方法也很简单：

这里返回的参数就是一个197*39的mfcc参数矩阵，当然这个矩阵的大小还是由输入的语音长度决定的。

实验过程中，有些同学如果想要实验过程中的 diguashao.wav,我已经将该文件上传至百度网盘，提取码为：fjq6 希望可以帮助到大家。

1. 2. 3. 4.

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/151165.html原文链接：https://javaforall.cn

版权声明：

本文来源网络，所有图片文章版权属于原作者，如有侵权，联系删除。

本文网址：https://www.mushiming.com/mjsbk/7284.html