梅尔倒谱系数MFCC由浅入深（超详细）

大家好，欢迎来到IT知识分享网。

MFCC梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients）

前言

clear all; clc; close all; y=load('su1.txt'); % 读入数据 fs=16000; nfft=1024; % 采样频率和FFT的长度 time=(0:nfft-1)/fs; % 时间刻度 figure(1), subplot 211; plot(time,y,'k'); % 画出信号波形 title('信号波形'); axis([0 max(time) -0.7 0.7]); ylabel('幅值'); xlabel(['时间/s' 10 '(a)']); grid; figure(2) nn=1:nfft/2; ff=(nn-1)*fs/nfft; % 计算频率刻度 Y=log(abs(fft(y))); % 按式(1)取实数部分 subplot 211; plot(ff,Y(nn),'k'); hold on; % 画出信号的频谱图 z=ifft(Y); % 按式(1)求取倒谱 figure(1), subplot 212; plot(time,z,'k'); % 画出倒谱图 title('信号倒谱图'); axis([0 time(512) -0.2 0.2]); grid; ylabel('幅值'); xlabel(['倒频率/s' 10 '(b)']); mcep=29; % 分离声门激励脉冲和声道冲激响应 zy=z(1:mcep+1); zy=[zy' zeros(1,1000-2*mcep-1) zy(end:-1:2)']; % 构建声道冲激响应的倒谱序列 ZY=fft(zy); % 计算声道冲激响应的频谱 figure(2), % 画出声道冲激响应的频谱，用灰线表示 line(ff,real(ZY(nn)),'color',[.6 .6 .6],'linewidth',3); grid; hold off; ylim([-4 5]); title('信号频谱（黑线）和声道冲激响频谱（灰线）') ylabel('幅值'); xlabel(['频率/Hz' 10 '(a)']); ft=[zeros(1,mcep+1) z(mcep+2:end-mcep)' zeros(1,mcep)]; % 构建声门激励脉冲的倒谱序列 FT=fft(ft); % 计算声门激励脉冲的频谱 subplot 212; plot(ff,real(FT(nn)),'k'); grid;% 画出声门激励脉冲的频谱 title('声门激励脉冲频谱') ylabel('幅值'); xlabel(['频率/Hz' 10 '(b)']); 

（2）离散余弦变换（Discrete Cosine Transform, DCT）
离散余弦变换经常用于信号处理和图像处理，用来对信号和图像进行有损数据压缩，这是由于离散余弦变换具有很强的”能量集中”特性：大多数的自然信号（包括声音和图像）的能量都集中在离散余弦变换后的低频部分，实际就是对每帧数据在进行一次降维。
对能量进行DCT变换。因为不同的Mel滤波器是有交集的，因此它们是相关的，我们可以用DCT变换去掉这些相关性，从而后续的建模时可以利用这一点(比如常见的GMM声学模型我们可以使用对角的协方差矩阵，从而简化模型)。

MFCC基本概述

MFCC计算过程

1.预处理

2.频域变换及能量求取

（1）FFT
由于信号在时域上的变换通常很难看出信号的特性，所以通常将它转换为频域上的能量分布来观察，不同的能量分布，就能代表不同语音的特性。所以在乘上汉明窗后，每帧还必须再经过快速傅里叶变换以得到在频谱上的能量分布。对分帧加窗后的各帧信号进行快速傅里叶变换得到各帧的频谱。

式中x(n)为输入的语音信号，N表示傅里叶变换的点数。

（2）谱线能量

对FFT后的频谱取模平方得到语音信号的谱线能量。

3.计算通过Mel滤波器的能量

• 三角形是低频密、高频疏的，这可以模仿人耳在低频处分辨率高的特性；
• 对频谱进行平滑化，并消除谐波的作用，突显原先语音的共振峰。（因此一段语音的音调或音高，是不会呈现在 MFCC 参数内，换句话说，以 MFCC 为特征的语音辨识系统，并不会受到输入语音的音调不同而有所影响；在每个三角形内积分，就可以消除精细结构，只保留音色的信息）
• 还可以减少数据量，降低运算量。

计算每个滤波器组输出的对数能量为：

深入内容：mel三角滤波器原理和设计

1.为什么会产生出Ｍel 这种尺度的机制呢？

Mel就是单词Melody（调，旋律）的简写，是根据人耳的听觉感知刻画出来的频率。

• 人耳朵具有特殊的功能，可以使得人耳朵在嘈杂的环境中，以及各种变异情况下仍能正常的分辨出各种语音
• 其中，耳蜗有关键作用，耳蜗实质上的作用相当于一个滤波器组，耳蜗的滤波作用是在对数频率尺度上进行的，在1000HZ以下为线性尺度，1K HZ以上为对数尺度，使得人耳对低频信号敏感，高频信号不敏感。
  例如，人们可以比较容易地发现500和1000Hz的区别，但很难发现7500和8000Hz的区别。
  

2.Mel刻度定义：

有时，还会写成：

区别也看显然，一个是log以10为底，一个是ln。

我们一般用三角滤波器来设计mel滤波器组，当然也可以用其他的滤波器，三角是mel设计时最常用的。

如下图所示，40个三角滤波器组成滤波器组，低频处滤波器密集，门限值大，高频处滤波器稀疏，门限值低。恰好对应了频率越高人耳越迟钝这一客观规律。该滤波器形式叫做等面积梅尔滤波器（Mel-filter bank with same bank area），在人声领域（语音识别，说话人辨认）等领域应用广泛。

3.mel滤波器的实现：

（1）确定最低频率（0HZ），最高频率（fs/2），Mel滤波器个数M（如23）；

（2）转换最低频率和最高频率的Mel(f)；

low_freq_mel = 0 high_freq_mel = 2595 * np.log10(1 + (fs / 2) / 700) print(low_freq_mel, high_freq_mel) 

（4）将每个中心Mel频率转化为频率f（非等间距）；

nfilt = 40 #滤波器个数 mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # 所有的mel中心点，为了方便后面计算mel滤波器组，左右两边各补一个中心点 hz_points = 700 * (10  (mel_points / 2595) - 1) 

fbank = np.zeros((nfilt, int(NFFT / 2 + 1))) # 各个mel滤波器在能量谱对应点的取值 bin = (hz_points / (fs / 2)) * (NFFT / 2) # 各个mel滤波器中心点对应FFT的区域编码，找到有值的位置 for i in range(1, nfilt + 1): left = int(bin[i-1]) center = int(bin[i]) right = int(bin[i+1]) for j in range(left, center): fbank[i-1, j+1] = (j + 1 - bin[i-1]) / (bin[i] - bin[i-1]) for j in range(center, right): fbank[i-1, j+1] = (bin[i+1] - (j + 1)) / (bin[i+1] - bin[i]) print(fbank) 

matlab版实现：

fs=8000; fl=0; fh=fs/2; bl=1125*log(1+fl/700);%将频率转换为Mel频率 bh=1125*log(1+fh/700); p=24;%滤波器个数 nfft=256;%FFT点数 B=bh-bl; y=linspace(0,B,p+2);%产生0到B之间p+2个数 Fb=700*(exp(y/1125)-1);%将Mel频率转换为频率 W2=nfft/2+1;%fs/2内对应的FFT点数 df=fs/nfft; freq=(0:W2-1)*df;%采样频率值 bank=zeros(24,W2);%生成一个24行W2列的全零数组 for k=2:p+1%why从2开始？因为k-1 f1=Fb(k-1); f2=Fb(k+1); f0=Fb(k); n1=floor(f1/df)+1;%f(m-1)在频域中的谱线索引号 n2=floor(f2/df)+1;%f(m+1)在频域中的谱线索引号 n0=floor(f0/df)+1;%f(m)在频域中的谱线索引号。f(m)是从0开始，而在MATLAB中数组的索引是从1开始，所以要加1，否则会出现index=0的错误 for i=1 : W2 if i>=n1 & i<=n0 bank(k-1,i)=(i-n1)/(n0-n1); elseif i>n0 & i<=n2 bank(k-1,i)=(n2-i)/(n2-n0); end end plot(freq,bank(k-1,:),'r','linewidth',2); hold on end grid; 

4. 对数能量取log, 计算DCT倒谱

(1)对能量取log
对于滤波器组的能量，我们对它取log。这也是受人类听觉的启发：人类对于声音大小(loudness)的感受不是线性的。为了使人感知的大小变成2倍，我们需要提高8倍的能量。这意味着如果声音原来足够响亮，那么再增加一些能量对于感知来说并没有明显区别。log这种压缩操作使得我们的特征更接近人类的听觉。为什么是log而不是立方根呢？因为log可以让我们使用倒谱均值减(cepstral mean subtraction)这种信道归一化技术（这可以归一化掉不同信道的差别）。

(2)DCT
经离散余弦变换（DCT）得到MFCC系数 :

将上述的对数能量带入离散余弦变换，求出L阶的Mel参数。L阶指MFCC系数阶数，通常取12-16。这里M是三角滤波器个数。

Deltas和Delta-Deltas特征

实现

MATLAB

注：在提取MFCC参数之前需要加载并使用VOICEBOX工具包

Df=5; fs=8000; N=fs/Df; t=0:1./fs:(N-1)./fs; x=sin(2*pi*200*t); bank=melbankm(24,256,8000,0,0.5,'t');%Mel滤波器的阶数为24，fft变换的长度为256，采样频率为8000Hz %归一化mel滤波器组系数 bank=full(bank); bank=bank/max(bank(:)); % DCT系数,12*p for k=1:12 n=0:23; dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24)); end %归一化倒谱提升窗口 w=1+6*sin(pi*[1:12]./12); %w=w/max(w); %语音信号分帧 xx=enframe(x,256,80);%对x 256点分为一帧 %计算每帧的MFCC参数 for i=1:size(xx,1) y=xx(i,:); s=y'.*hamming(256); t=abs(fft(s));%fft快速傅立叶变换 t=t.^2; c1=dctcoef*log(bank*t(1:129)); c2=c1.*w'; end plot(c2);title('MFCC'); 

Python实现方法

import numpy as np from scipy import signal from scipy.fftpack import dct import pylab as plt def enframe(wave_data, nw, inc, winfunc): '''将音频信号转化为帧。 参数含义： wave_data:原始音频型号 nw:每一帧的长度(这里指采样点的长度，即采样频率乘以时间间隔) inc:相邻帧的间隔（同上定义） ''' wlen=len(wave_data) #信号总长度 if wlen<=nw: #若信号长度小于一个帧的长度，则帧数定义为1 nf=1 else: #否则，计算帧的总长度 nf=int(np.ceil((1.0*wlen-nw+inc)/inc)) pad_length=int((nf-1)*inc+nw) #所有帧加起来总的铺平后的长度 zeros=np.zeros((pad_length-wlen,)) #不够的长度使用0填补，类似于FFT中的扩充数组操作 pad_signal=np.concatenate((wave_data,zeros)) #填补后的信号记为pad_signal indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T #相当于对所有帧的时间点进行抽取，得到nf*nw长度的矩阵 indices=np.array(indices,dtype=np.int32) #将indices转化为矩阵 frames=pad_signal[indices] #得到帧信号 win=np.tile(winfunc,(nf,1)) #window窗函数，这里默认取1 return frames*win #返回帧信号矩阵 Df=5 fs=8000 N=fs/Df t = np.arange(0,(N-1)/fs,1/fs) wave_data=np.sin(2*np.pi*200*t) #预加重 #b,a = signal.butter(1,1-0.97,'high') #emphasized_signal = signal.filtfilt(b,a,wave_data) #归一化倒谱提升窗口 lifts=[] for n in range(1,13): lift =1 + 6 * np.sin(np.pi * n / 12) lifts.append(lift) #print(lifts)  #分帧、加窗  winfunc = signal.hamming(256) X=enframe(wave_data, 256, 80, winfunc) #转置的原因是分帧函数enframe的输出矩阵是帧数*帧长 frameNum =X.shape[0] #返回矩阵行数18，获取帧数 #print(frameNum) for i in range(frameNum): y=X[i,:] #fft yf = np.abs(np.fft.fft(y)) #print(yf.shape) #谱线能量 yf = yf2 #梅尔滤波器系数 nfilt = 24 low_freq_mel = 0 NFFT=256 high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700)) # 把 Hz 变成 Mel mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # 将梅尔刻度等间隔 hz_points = (700 * (10(mel_points / 2595) - 1)) # 把 Mel 变成 Hz bin = np.floor((NFFT + 1) * hz_points / fs) fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1)))) for m in range(1, nfilt + 1): f_m_minus = int(bin[m - 1]) # left f_m = int(bin[m]) # center f_m_plus = int(bin[m + 1]) # right for k in range(f_m_minus, f_m): fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1]) for k in range(f_m, f_m_plus): fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m]) filter_banks = np.dot(yf[0:129], fbank.T) filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks) # 数值稳定性 filter_banks = 10 * np.log10(filter_banks) # dB  filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8) #print(filter_banks) #DCT系数 num_ceps = 12 c2 = dct(filter_banks, type=2, axis=-1, norm='ortho')[ 1 : (num_ceps + 1)] # Keep 2-13 c2 *= lifts print(c2) plt.plot(c2) plt.show() 

总结

因此，MFCC的全部组成其实是由： N维MFCC参数（N/3 MFCC系数+ N/3 一阶差分参数+ N/3 二阶差分参数）+帧能量（此项可根据需求替换）。

这里的帧能量是指一帧的音量（即能量），也是语音的重要特征，而且非常容易计算。因此，通常再加上一帧的对数能量（定义：一帧内信号的平方和，再取以10为底的对数值，再乘以10）使得每一帧基本的语音特征就多了一维，包括一个对数能量和剩下的倒频谱参数。另外，解释下最开始说的40维是怎么回事，假设离散余弦变换的阶数取13，那么经过一阶二阶差分后就是39维了再加上帧能量总共就是40维，当然这个可以根据实际需要动态调整。

相关资料：
http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/