音频入门: 最全面详细的Mel频谱和MFCC讲解

前言

在音频领域，mel频谱和mfcc是非常重要的特征数据，在深度学习领域通常用此特征数据作为网络的输入训练模型，来解决音频领域的各种分类、分离等业务，如端点侦测、节奏识别、和弦识别、音高追踪、乐器分类、音源分离、回声消除等相关业务。

当然，针对深度学习音频领域的业务，不是用下这两个特征、选几个网络、打个标签，放数据训练就完事了， 仅仅基于mel频谱和mfcc这两个特征，解决好上述业务某些情况下还是远远不够的，熟悉这些特征的内在逻辑性、衍生细节和延展，才能更好的结合深度学习解决业务问题。

下面讲解mel频谱和mfcc特征的算法流程和一些细节、延展，这些细节从局部角度来看，都会影响到最终特征呈现的细节差异，这些差异放大到模型训练结果的准确性、鲁棒性上怎么样是非常值得研究的，某些情况下可能会有质的变化，质的变化无论正向还是负向都是值得关注的，最怕的是没变化；同时，一些问题的延展从广义角度来看，带来不同的特征组合、网络结构设计思考等也是解决业务问题非常重要的思想源泉。

注： 文中一些图例来自开源项目audioFlux: https://github.com/libAudioFlux/audioFlux

算法流程

设 sr 为采样率，fftLength 为帧长度，slideLength 为滑动长度

下面是一张mel频谱和mfcc的大概算法流程图。

1. 预加重（Pre-emphasis）

如流程图所示的第1步，属于信号的预处理，补偿高频分量损失，提升高频分量，一般情况下可以忽略此步骤，属于信号的简单增强，对特征有一定的提升效果。公式如下

$$x[n]=x[n]-\alpha x[n-1]$$

$\alpha$一般取0.97
公式属于差分一阶即高通滤波器。

2. 分帧（Frame）

现实中大多数信号都是非平稳的，但大多数短时间内可以近似看做是平稳的，可以用短时傅里叶变换表现非平稳信号频域特征。 一般语音中采用10ms~30ms左右，乐音中可以更长一些64ms~256ms。

分帧涉及到前后重叠（overlap），一般情况下以滑动帧长的1/4或1/2（前后重叠3/4或1/2）进行，即$slideLength=\frac{fftLength}{4}或\frac{fftLength}{2}$。

3. 加窗（Window）

加窗目的是减少频谱泄露，降低泄漏频率干扰，提升频谱效果，默认不处理即加矩形窗（Rect），干扰泄漏较严重，一般情况下加Hann窗，针对大多数信号都有不错的效果。公式如下
$$w(n)=0.5\left(1-\cos \left(2\pi \frac{n}{N}\right)\right),0\le n\le N$$

Hann为余弦窗，N表示阶数。

4. 短时傅里叶变换（STFT）

分帧加窗傅里叶变换即短时傅里叶变换。公示如下
$$X(\tau ,f)={\int }_{-\infty }^{\infty }x(t)w(t-\tau )e^{-j2\pi ft}dt$$
$$X(m,k)=\sum _{n=0}^{N-1}x[n]W[n-m]e^{\frac{-j2\pi kn}{N}}$$

设数据长度为 dataLength，
$t=\{\begin{array}{ll}\frac{(dataLength-fftLength)}{slideLength}+1,& 无填充\\ \frac{dataLength}{slideLength}+1,& 填充fftLength\end{array}$

STFT属于标准的数学变换，为复数域，尺寸为 $t\ast fftLength$，为表示区分，一般的如$|X(m,k)|$取模，尺寸为$t\ast (fftLength/2+1)$ 表示为STFT频谱，有以下类型频谱。

$|X(m,k)|$，STFT幅值频谱
$|X(m,k){|}^2$，STFT功率频谱
$\log (|X(m,k)|)$，STFT dB（分贝）频谱
$20\log \left(\frac{|X(m,k)|}{fftLength}\right)$，STFT标准dB频谱

注：

1. dB频谱属于相对谱，加减乘除对频率相对dB值没影响，标准dB频谱相当于建立一个基准参考线，方便分析频谱和量化，大多数频谱示例都是此种类型。
2. 深度学习中使用dB谱训练大多数要优于其它数值类型的频谱。

5. 滤波器组过程（Filter bank）

此过程是计算mel频谱关键部分和mfcc的重要一步。流程图如下

mel刻度（scale）是一种基于人耳听觉设计的log压缩刻度，人耳针对低频比较敏感，高频不太敏感，比如110hz和116hz一般人都能区分出来，但4000hz和4100hz大部分分不出来。mel刻度和hz的转换公式如下
$$\{\begin{array}{l}mel=2595{\log }_{10}(1+\frac{hz}{700})\\ hz=700(10^{\frac{mel}{2595}}-1)\end{array}$$

图中前三个步骤即根据业务fre边界和num频带个数计算mel刻度下所映射的freBandArr。
接下来就是STFT频带如何映射到mel刻度的频带，用频带加三角窗进行计算mel刻度的filterBank matrix，三角窗公式如下
$$w(n)=\{\begin{array}{ll}\frac{2n}{N},& 0\le n\le \frac{N}{2}\\ 2-\frac{2n}{N},& \frac{N}{2}\le n\le N\end{array}$$

如下图

当然这个是对称的，线性频带用对称的没问题，mel刻度产生的频带log比例，左右不会对称，这时候使用三角窗就是非对称的，下面是mel刻度的filterBank matrix公式
$$w(k,h)=\{\begin{array}{ll}\frac{h-f_{k-1}}{f_k-f_{k-1}},& f_{k-1}<h\le f_k\\ \frac{f_{k+1}-h}{f_{k+1}-f_k},& f_k<h\le f_{k+1}\\ 0,& other\end{array}$$

下面给出一个简单的示意图

最后，一般使用STFT功率频谱和filterBank matrix做矩阵乘法运算即得出mel功率频谱。

6. 非线性校正（Rectification）

对上一步mel功率频谱取log运算，即mel dB频谱就是通常所使用的“mel频谱”，一般情况下5、6步作为一个整体看待，合到一块可以理解为mel频谱的计算。

下面是STFT频谱和mel频谱（num=128）的效果对比图。

注： 图例来自开源项目audioFlux: https://github.com/libAudioFlux/audioFlux

7. 离散余弦变换（DCT）

离散余弦变换，即数据为偶函数的实数傅里叶的变换，有去相关和能量集中特点。公式如下
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{\frac{-j2\pi kn}{N}}=\sum _{n=0}^{N-1}x[n]\left[\cos (\frac{2\pi kn}{N})-j\sin (\frac{2\pi kn}{N})\right]$$

当$\sin (\theta )=0$，即DCT变换， $\cos (\theta )=0$，即DST变换。

DCT变换，属于数据高度去相关后的特征反映，由于能量集中的特点，在数据压缩领域广泛常用，一般的，DCT考虑到数据解析延拓时中间点选取、边界、方式等情况有8种对应变换，常用的为DCT-II变换，公式如下
$$X(k)=\sum _{n=0}^{N-1}x[n]\cos \left[\frac{\pi }{N}(n+\frac{1}{2})k\right]$$

mel频谱经过DCT变换后得到倒谱系数（cepstral coeffs）即MFCCs。

8. 能量和delta

能量和delta的计算属于mfcc特征体系下的可选操作。
能量特征相当于给mfcc加上bias偏置，具有一定抗噪作用，公式如下
$$energy=\sum _{n=1}^N{x}^2[n]=\frac{1}{N}\sum _{m=1}^N|X[m]{|}^2$$

根据公式，可以从时域计算，也可以基于频域计算，很明显，这样计算出的结果动态范围太小，一般针对结果取log运算即$energy=\log ({\sum }_{n=1}^N{x}^2[n])$，又称log能量，可以替换首个mfcc直流分量值或首位置追加。

delta是计算数据的变化，基于当前点区域的局部斜率最小二乘近似值，公式如下
$$delta=\frac{{\sum }_{k=-M}^{M}kx[k]}{{\sum }_{k=-M}^M{k}^2}$$

M为阶数，为奇数，一般取9

针对mfcc计算其delta，然后再计算delta的delta，可以侦测mfcc状态的变化，变化的变化，可以作为mfcc的两组辅助特征参与网络模型的训练，某些情况下起到更好的准确性和泛化能力。

mfcc相关效果图如下

注： 图例来自开源项目audioFlux: https://github.com/libAudioFlux/audioFlux

综上所有，详细描述解释了算法流程中每一步计算流程，下面将对一些步骤流程涉及到的细节思考点和延伸点做深入的展开。

相关细节

1. weight-A计权

针对上面预加重方式而言，本质上是一个高通滤波器，是衰减低频的，本身并不会提升高频，“补偿高频分量损失，提升高频分量”并不太严谨但比较形象，衰减低频相对来说就是变相提升高频，是缩减高低频动态范围的。

为什么要缩减高低频动态范围，有没有更好的方式，答案就是weight-A计权。

上面有提到“人耳对低频比较敏感，高频不太敏感”，这个敏感是针对频率分辨而言的，针对频率的强度而言，这句话是反过来的即人耳对高频比较敏感，低频不太敏感，以吉他乐器为例，用同样力度拨最粗和最细的弦发出的音，即物理上两者强度（振幅）是一样的，但听觉上高频比低频要响的多，人们用weight-A计权来量化这种“响度”的心理指标，数学上表示不同频段的log函数加减，效果如下图。

如图所示，整个频段不仅有衰减、还有真正的提升区域。

针对一些业务的深度学习模型训练，加weight-A计权能提升4%~5%左右精度。

2. overlap重叠

数据分帧涉及到帧长和重叠两个问题，帧长决定频域的频率分辨率和时域的时间分辨率，帧长越长，频域分辨率越精确，时域分辨率越模糊，但受限大多数信号本身非平稳特点不可能无限长，帧长越短，时域分辨率越精确，频域分辨率越模糊。

重叠问题，就是相对当前帧滑动的问题，如上面分帧流程所述一般情况下滑动帧长的1/4或1/2，当然，滑动长度也可以等同帧长（前后重叠为0），甚至超过帧长（没有重叠，前后跳跃）。

重叠多少还是不重叠还是跳跃，本身并不会提升时域分辨率，可以理解为频谱t维度的不同时间间隔采样，滑动小相当于频谱图的插值升采样，滑动大相当于频谱图的抽取降采样，针对端点侦测业务频域的相关算法，滑动太小或跳动过大都不会有好的效果，如下面效果对比。

如图所见，紫色和绿色是滑动极小和跳动过大效果。

针对端点侦测相关业务，选择合适帧长后，前后滑动一般没有重叠或少许跳动效果相对好一些。

3. window窗函数

时域加窗目的是减少频谱泄露，上述算法流程描述中，一般情况下加Hann窗，但窗函数有很多，如Triang，Hann，Hamm，Guass，Kaiser，Flattop，Blackman等等。

不同窗如何选择，在深度学习一些业务中，不同窗的选择对模型的训练和结果影响是怎样的，或者哪些窗在业务中是值得做尝试的。

Guass，Kaiser非常值得尝试，公式如下
$$\{\begin{array}{l}guass(n)=e^{-n^2/2{\sigma }^2}=e^{-\frac{1}{2}{\left(\alpha \frac{n}{(N-1)/2}\right)}^2},-(N-1)/2\le n\le (N-1)/2\\ kaiser(n)=\frac{I_0\left(\beta \sqrt{1-{\left(\frac{n-N/2}{N/2}\right)}^2}\right)}{I_0(\beta )},0\le n\le N\end{array}$$

针对Guass：

$\sigma =(N-1)/(2\alpha )$, 两者之间互为反比

默认$\alpha =2.5$

针对Kaiser：

$I_0(\beta )$为零阶第一类修正贝塞尔函数，可有下面公式级数计算，
$I_0(\beta )=1+{\sum }_{k=1}^{\infty }{\left[\frac{1}{k!}{\left(\frac{\beta }{2}\right)}^k\right]}^2$，一般取$15$项左右
默认$\beta =5$

Gauss，Kaiser相比其它大多数窗而言，属于窗口可变窗，一般窗口宽度越大，分辨率越细，过渡带越小，但阻带衰减偏大；一般窗口宽度越窄，分辨率越大，过渡带越大，但阻带衰减偏小。

针对Guass和Kaiser这两种可以调整宽度的窗，在一些业务中可以选择最大化旁瓣衰减，用相应的频谱特征训练深度学习业务模型时，可能会有不错的表现。

4. rectification非线性校正

上述算法流程描述中，一般而言5、6步合到一起是标准的频谱呈现结果，但从功率频谱非线性校正角度而言，
不止log函数一种，当然log是最重要最常用的一种，可以用cubic root（三次开方），类Relu等其它方式。

从深度学习角度来看，可以把类mel频谱当做一种网络层计算，log、cubic root等非线性操作当做激活函数。

从这个角度来看，针对不同的业务，激活函数的不同选择设计，对模型的训练和结果的影响有可能是天翻地覆的，在研发业务中是很值得关注测试的一个点。

各种刻度

1. 其它刻度

mel频谱如此知名，难道做音频领域深度学习就只有mel频谱？当然不是，有些情况mel频谱不一定是最优的。

如上面算法流程描述中mel刻度的解释，是一种基于人耳听觉设计的log压缩刻度，后面发展出更准确人耳听觉模型刻度，Bark和ERB刻度，公式如下
$$bark=\frac{26.81hz}{1960+hz}-0.53,hz=1960(\frac{bark+0.53}{26.81-bark})$$
$$erb=A{\log }_{10}(1+0.00437hz),hz=\frac{10^{\frac{erb}{A}}-1}{0.00437}$$

$A=\frac{1000\ln (10)}{(24.7)(4.37)}$

2. Gammatone filter

$$g(t)=a{t}^{n-1}{e}^{-2\pi bt}\cos (2\pi f_{c}t+\varphi )$$

公式中 $a{t}^{n-1}{e}^{-2\pi bt}$ 部分为 $\Gamma$ 形式函数，$\cos$ 理解为tone，称之为gammatone。
一般情况下，人耳听觉模型中的ERB刻度和gammatone filter的bandwidth相关联，即b为ERB刻度的bandwidth。

基于这个滤波器求频响非常复杂，只能给出近似公式，论文公式推导足足有20多页，已超过大部分人的研究上限，即使拿论文的结果公式直接编程实现也不轻松。幸运的是audioFlux开源项目有比较标准的实现，感兴趣的朋友可以研究一下。

3. 刻度频谱

除了上述刻度以外，还有基于乐音八度的octave刻度，更通用的log刻度等等，基于上面算法流程，mel刻度产生mel频谱和mfcc特征，同样的流程，bark/erb等刻度产生对应的bark/erb频谱和相应倒谱系数，一些不同刻度频谱的对比图如下

注： 图例来自开源项目audioFlux: https://github.com/libAudioFlux/audioFlux

更多的不同刻度频谱类型可尝试使用audioFlux库测试。

4. 可能的问题

在一些业务中，比如乐器相关业务中，上面所有的刻度中，可能都会有一些小问题，针对低频的频带，可能有些已经小于STFT的最小分辨率，这时候显然没意义，额外的造成干涉发散效果；针对中频的频带，有些可能过大，分辨不够细；针对高频的频带，可能还不够大，分辨还是相对过小。

上述情况下，在一些业务中，特定的一些网络结构中，中低频的分辨率差可能会导致训练模型的欠拟合，很难达到一定准确度，或者高频的分辨率过小可能会导致训练模型的欠拟合，泛化能力差，出现这种情况后可以尝试以下方式：

1，适当的增加数据集、调整网络结构和参数等通用方式。
2，增加不同维度的特征输入，让网络找到最优解或跳出局部最优解。
3，自定义刻度，既然明确原因，从特征本身优化上入手。

综上，不同scale下产生的频谱数据高低频细节、能量聚集、对比度都会有不同的差异，这种差异放大到一些业务实践中，模型的准确度和鲁棒性怎么样，就非常有尝试研究价值。

在深度学习一些业务中，如果mel刻度下的相关特征可以出结果，使用bark相关特征替换，往往有一定的效果提升。

滤波器组

1. 概念

数字信号中，滤波器有基础的高通滤波器和低通滤波器， 带通滤波器可以有低通和高通串联而成，Filter bank可以理解为多个带通滤波器。

滤波器组是一个映射矩阵，表示STFT线性频带和不同刻度下的频带映射关系，至于怎么映射，就是算法流程描述第5步中的三角窗函数法，三角窗的Filter bank如下图所示

2. 归一化

如上图所示，一般刻度产生的频带，低频带宽小，高频越来越大，意味着高频区域很长的一段频带都会参与当前频率分量的映射计算，相对低频而言显然不太合理，这时候需要对Filter bank进行归一化处理。

归一化方式有以带宽和面积区域两种方式，以带宽方式归一化效果如下图

3. 不同的窗

针对Filter bank的计算，难道就只有三角窗函数法，当然不是，可以使用Rect，Hann，Hamm等等这些窗函数，Gammatone也是一种特殊的窗函数，某种角度上讲甚至可以不计算，下面是几种不同窗下Filter bank的对比图

下面是不同窗下的频谱对比图

注： 图例来自开源项目audioFlux: https://github.com/libAudioFlux/audioFlux

图中Slaney和ETSI是两种不同类型的三角窗，其它更多窗类型测试可尝试使用audioFlux库。

4. 加窗方式

窗函数重叠处处理有很多细化方式，频带前后窗都是有重叠的，一般的处理方式是重叠点后面不再参与之前计算，重叠点之前不再参与之后计算，如下面图所示

如果各个频带窗函数可以等比例缩放，这样的话对不同刻度下的频带都变成可调整的，窗函数不同的宽度选择，对应时频分辨率不同的变化，结合PWT（伪小波变换），就是另外一种CWT效果的同等变体，不同于CWT对波函数的限定要求，可以用更广泛的窗函数研究CWT。

综上，针对Filter bank的计算，**加什么窗？怎么加窗？何种归一化方式？**使用这些组合产生出更多细粒度的不同特征，训练出模型准确性和鲁棒性如何，是非常值得尝试研究的。

倒谱系数

1. 概念

什么是倒谱系数，为什么DCT计算后是倒谱系数？

倒谱的标准公式定义如下
C ( r ) = F { log ⁡ ( F { f ( t ) } ) } C(r)=F\{\log(F\{f(t)\})\} C(r)=F{log(F{f(t)})}

数据经过FFT变换后取log后，再次FFT变换一般称为倒谱，更广义的讲即时域数据经过FFT变为频域数据，频域数据再次
FFT变换即倒谱变换。

DCT是DFT的特例，所以DCT针对之前FFT结果再次变换符合上述公式概念定义，也称倒谱系数；相比DFT，DCT能量较聚集，相当于频谱数据的再压缩，用小部分系数就能反映出数据的高度去相关的特征。

2. 倒谱应用

倒谱可以干什么？

可以估算音高，可以作为解卷积（deconv）的一种方式，分离信号，计算包络（envelope）/共振峰（Formant）等等，
解卷积推导公式如下
{ x ( t ) = g ( t ) ∗ h ( t ) F { x ( t ) } = X ( ω ) = G ( ω ) H ( ω ) log ⁡ X ( ω ) = X ^ ( ω ) = G ^ ( ω ) + H ^ ( ω ) F − 1 { X ^ ( ω ) } = x ^ ( k ) = g ^ ( k ) + h ^ ( k ) \begin{cases} x(t)=g(t)*h(t) \\ F\{x(t)\}=X(\omega)=G(\omega)H(\omega) \\ \log X(\omega)=\hat X(\omega)=\hat G(\omega)+\hat H(\omega) \\ F^{-1}\{\hat X(\omega)\}=\hat x(k)=\hat g(k)+\hat h(k) \end{cases} ⎩ ⎨ ⎧​x(t)=g(t)∗h(t)F{x(t)}=X(ω)=G(ω)H(ω)logX(ω)=X^(ω)=G^(ω)+H^(ω)F−1{X^(ω)}=x^(k)=g^​(k)+h^(k)​

3. 延展

分数傅里叶变换（FRFT）

根据倒谱的定义，延展一下，能否多次FFT，多次FFT的意义是什么，多次变换引申如下
$$F^2=F(F(x))$$
$$F^n(x)=F^{(n-1)}(F(x))$$
设$n=\frac{2\varphi }{\pi }$，换元后定义
$$F_{\varphi }(x)=F^{2\varphi /\pi }(x)$$

引用$\varphi$后，FFT变换的次数可以非整数即分数傅里叶变换，物理意义的表现就是旋转频谱一定角度，可以多一个维度分析频谱，同时也带来其它概念的延伸，如分数卷积等等。

总结

本文大体上分四个部分，算法流程详细描述解释了每一步计算流程，相关细节算是简单步骤流程的细节点思考和延展，各种刻度、滤波器组、倒谱系数算是重要概念和重要步骤的细节点思考和延展。