MFCC

MFCC（Mel-frequency cepstral coefficients）：梅尔频率倒谱系数。梅尔频率是基于人耳听觉特性提出来的，它与Hz频率成非线性对应关系。
MFCC提取过程：

首先对语音进行预处理。
预处理又包括对语音进行预加重、分帧、加窗。
快速傅里叶变换
对分帧加窗后的每帧语音数据进行fft变换。
计算谱线能量
对fft变换后的每帧信号取平方。
梅尔滤波器组滤波
将求出的每帧谱线能量谱与梅尔滤波器组相乘。
DCT变换
经过梅尔滤波器滤波后的每帧信号进行DCT（离散余弦变换），得到MFCC参数。

详细过程可参加其他博客，本文主要使用matlab编程实现。

1.梅尔滤波器组实现

matlab编程实现梅尔滤波器组：

clc;
clear;
close all;
fs = 16000;  %信号采样频率
p = 16;   %梅尔滤波器个数
N = 1024
fl=0;fh=fs/2;%定义频率范围，低频和高频
bl=2595*log10(1+fl/700);%得到梅尔刻度的最小值
bh=2595*log10(1+fh/700);%得到梅尔刻度的最大值
B=bh-bl;%梅尔刻度长度
mm=linspace(0,B,p+2);%规划18个不同的梅尔刻度。但只有16个梅尔滤波器
fm=700*(10.^(mm/2595)-1);%将Mel频率转换为频率k=((N+1)*fm)/fs;%计算18个不同的k值
hm=zeros(p,N);%创建hm矩阵
df=fs/N;
freq=(0:N-1)*df;%采样频率值%绘制梅尔滤波器
for i=2:p+1%取整n0=floor(k(i-1));  %floor为向下取整，如果取整，表示三角滤波器组的各中心点被置于FFT 数据对应的的点位频率上，所以频率曲线的三角形顶点正好在FFT 数据点上，即美尔滤波器所有顶点在同一高度。n1=floor(k(i));n2=floor(k(i+1));
%     n0=k(i-1);  %不取整代表直接按精确指定的中心频率点来计算的，这时滤波器中心频率就不一定在FFT 数据点对应位置，按FFT 点位所得曲线在低频段就不是完美的三角形，即美尔滤波器所有顶点不在同一高度。
%     n1=k(i);
%     n2=k(i+1);%要知道k(i)分别代表的是每个梅尔值在新的范围内的映射，其取值范围为：0-N/2%以下实现求取三角滤波器的频率响应。for j=1:Nif n0<=j & j<=n1hm(i-1,j)=(j-n0)/(n1-n0);elseif n1<=j & j<=n2hm(i-1,j)=(n2-j)/(n2-n1);endend%此处求取hm结束。
end%绘图,且每条颜色显示不一样
c=colormap(lines(p));%定义16条不同颜色的线条
% figure()
% set(gcf,'position',[180,160,1020,550]);%设置画图的大小
for i=1:pplot(freq,hm(i,:),'-','color',c(i,:),'linewidth',2);%开始循环绘制每个梅尔滤波器
%     plot(freq,hm(i,:),'-','linewidth',2.5);%开始循环绘制每个梅尔滤波器hold on;
endgrid on;%显示方格axis([0 8000 0 inf]);%设置显示范围xlabel('频率（Hz）');ylabel('幅值');

画出的梅尔滤波器组图：

上图中，三角滤波器组的各中心点被置于FFT 数据对应的的点位频率上，所以频率曲线的三角形顶点正好在FFT 数据点上，即美尔滤波器所有顶点在同一高度。

也可以按照精确指定的中心频率点来计算，此时三角滤波器不等高。
程序：

clc;
clear;
close all;
fs = 16000;  %信号采样频率
p = 16;   %梅尔滤波器个数
N = 1024
fl=0;fh=fs/2;%定义频率范围，低频和高频
bl=2595*log10(1+fl/700);%得到梅尔刻度的最小值
bh=2595*log10(1+fh/700);%得到梅尔刻度的最大值
B=bh-bl;%梅尔刻度长度
mm=linspace(0,B,p+2);%规划18个不同的梅尔刻度。但只有16个梅尔滤波器
fm=700*(10.^(mm/2595)-1);%将Mel频率转换为频率k=((N+1)*fm)/fs;%计算18个不同的k值
hm=zeros(p,N);%创建hm矩阵
df=fs/N;
freq=(0:N-1)*df;%采样频率值%绘制梅尔滤波器
for i=2:p+1%取整%n0=floor(k(i-1));  %floor为向下取整，如果取整，表示三角滤波器组的各中心点被置于FFT 数据对应的的点位频率上，所以频率曲线的三角形顶点正好在FFT 数据点上，即美尔滤波器所有顶点在同一高度。%n1=floor(k(i));%n2=floor(k(i+1));n0=k(i-1);  %不取整代表直接按精确指定的中心频率点来计算的，这时滤波器中心频率就不一定在FFT 数据点对应位置，按FFT 点位所得曲线在低频段就不是完美的三角形，即美尔滤波器所有顶点不在同一高度。n1=k(i);n2=k(i+1);%要知道k(i)分别代表的是每个梅尔值在新的范围内的映射，其取值范围为：0-N/2%以下实现求取三角滤波器的频率响应。for j=1:Nif n0<=j & j<=n1hm(i-1,j)=(j-n0)/(n1-n0);elseif n1<=j & j<=n2hm(i-1,j)=(n2-j)/(n2-n1);endend%此处求取hm结束。
end%绘图,且每条颜色显示不一样
c=colormap(lines(p));%定义16条不同颜色的线条
% figure()
% set(gcf,'position',[180,160,1020,550]);%设置画图的大小
for i=1:pplot(freq,hm(i,:),'-','color',c(i,:),'linewidth',2);%开始循环绘制每个梅尔滤波器
%     plot(freq,hm(i,:),'-','linewidth',2.5);%开始循环绘制每个梅尔滤波器hold on;
endgrid on;%显示方格axis([0 8000 0 inf]);%设置显示范围xlabel('频率（Hz）');ylabel('幅值');

画出的图：

上图是直接按精确指定的中心频率点来计算的，这时滤波器中心频率就不一定在FFT 数据点对应位置，按FFT 点位所得曲线就不是完美的三角形，即美尔滤波器所有顶点不在同一高度。

2. MFCC主程序

matlab实现MFCC的整个过程：

frameSize=1024;  %帧长
N = frameSize;
inc=512;  %帧移
p = 32;  %梅尔滤波器个数
[x,fs]=audioread('jian_2_1000.wav');%读取语音wav文件，本文语音长度3秒，单声道，采样频率16000
% 本文语音信号的fs为16000Hz
N2=length(x);  %语音信号长度%预加重
for i=2:N2y(i)=x(i)-0.97*x(i-1);
end
y=y';%对y取转置
S=enframe(y,frameSize,inc);%分帧,对语音信号x进行分帧，
[a,b]=size(S);  %a为矩阵行数，b为矩阵列数
%创建汉明窗矩阵C
C=zeros(a,b);
ham=hamming(b);
for i=1:aC(i,:)=ham;
end
%将汉明窗C和S相乘得SC
SC=S.*C;
F=0;N=frameSize;
for i=1:a%对SC作N=1024的FFT变换D(i,:)=fft(SC(i,:),N);%以下循环实现求取谱线能量for j=1:Nt=abs(D(i,j));E(i,j)=(t^2);end
endfl=0;fh=fs/2;%定义频率范围，低频和高频
bl=2595*log10(1+fl/700);%得到梅尔刻度的最小值
bh=2595*log10(1+fh/700);%得到梅尔刻度的最大值
%梅尔坐标范围
% p=26;%滤波器个数
B=bh-bl;%梅尔刻度长度
mm=linspace(0,B,p+2);%规划34个不同的梅尔刻度。但只有32个梅尔滤波器
fm=700*(10.^(mm/2595)-1);%将Mel频率转换为频率
W2=N/2+1;%fs/2内对应的FFT点数,1024个频率分量k=((N+1)*fm)/fs;%计算34个不同的k值
hm=zeros(p,N);%创建hm矩阵
df=fs/N;
freq=(0:N-1)*df;%采样频率值%绘制梅尔滤波器
for i=2:p+1%取整n0=floor(k(i-1));  %floor为向下取整，如果取整，表示三角滤波器组的各中心点被置于FFT 数据对应的的点位频率上，所以频率曲线的三角形顶点正好在FFT 数据点上，即美尔滤波器所有顶点在同一高度。n1=floor(k(i));n2=floor(k(i+1));
%     n0=k(i-1);  %不取整代表直接按精确指定的中心频率点来计算的，这时滤波器中心频率就不一定在FFT 数据点对应位置，按FFT 点位所得曲线在低频段就不是完美的三角形，即美尔滤波器所有顶点不在同一高度。
%     n1=k(i);
%     n2=k(i+1);%要知道k(i)分别代表的是每个梅尔值在新的范围内的映射，其取值范围为：0-N/2%以下实现求取三角滤波器的频率响应。for j=1:Nif n0<=j & j<=n1hm(i-1,j)=(j-n0)/(n1-n0);elseif n1<=j & j<=n2hm(i-1,j)=(n2-j)/(n2-n1);endend%此处求取H1(k)结束。
end%梅尔滤波器滤波H=E*hm';%对H作自然对数运算%因为人耳听到的声音与信号本身的大小是幂次方关系，所以要求个对数for i=1:afor j=1:pH(i,j)=log(H(i,j));endend%作离散余弦变换
Fbank = H';
Fbank1 = dct(Fbank);
mfcc = Fbank1';
plot(mfcc);

最后得到的mfcc就是MFCC参数。

MFCC程序求取过程中的分帧函数：
分帧函数为enframe函数，matlab程序为：

function frameout=enframe(x,win,inc)nx=length(x(:));            % 取数据长度
nwin=length(win);           % 取窗长
if (nwin == 1)              % 判断窗长是否为1，若为1，即表示没有设窗函数len = win;               % 是，帧长=win
elselen = nwin;              % 否，帧长=窗长
end
if (nargin < 3)             % 如果只有两个参数，设帧inc=帧长inc = len;
end
nf = fix((nx-len+inc)/inc); % 计算帧数
frameout=zeros(nf,len);            % 初始化
indf= inc*(0:(nf-1)).';     % 设置每帧在x中的位移量位置
inds = (1:len);             % 每帧数据对应1:len
frameout(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));   % 对数据分帧
if (nwin > 1)               % 若参数中包括窗函数，把每帧乘以窗函数w = win(:)';            % 把win转成行数据frameout = frameout .* w(ones(nf,1),:);  % 乘窗函数
end

画出的语音的mfcc参数数据图：

图中横轴代表语音分帧的帧数，纵轴代表mfcc参数的幅值。

mfcc参数图：
在mfcc程序后面加上：

imagesc(mfcc');
axis xy;   %y轴从下往上
colorbar;
xlabel('语音分帧帧数');
ylabel('mfcc参数维数');

画出的图：

3. 加上一阶差分和二阶差分的MFCC

加上一阶差分和二阶差分的MFCC-D-A参数的matlab程序：

frameSize=1024;  %帧长
N = frameSize;
inc=512;  %帧移
p = 32;  %梅尔滤波器个数
[x,fs]=audioread('jian_2_1000.wav');%读取语音wav文件，本文语音长度3秒，单声道，采样频率16000
% 本文语音信号的fs为16000Hz
N2=length(x);  %语音信号长度%预加重
for i=2:N2y(i)=x(i)-0.97*x(i-1);
end
y=y';%对y取转置
S=enframe(y,frameSize,inc);%分帧,对语音信号x进行分帧，
[a,b]=size(S);  %a为矩阵行数，b为矩阵列数
%创建汉明窗矩阵C
C=zeros(a,b);
ham=hamming(b);
for i=1:aC(i,:)=ham;
end
%将汉明窗C和S相乘得SC
SC=S.*C;
F=0;N=frameSize;
for i=1:a%对SC作N=1024的FFT变换D(i,:)=fft(SC(i,:),N);%以下循环实现求取谱线能量for j=1:Nt=abs(D(i,j));E(i,j)=(t^2);end
endfl=0;fh=fs/2;%定义频率范围，低频和高频
bl=2595*log10(1+fl/700);%得到梅尔刻度的最小值
bh=2595*log10(1+fh/700);%得到梅尔刻度的最大值
%梅尔坐标范围
% p=26;%滤波器个数
B=bh-bl;%梅尔刻度长度
mm=linspace(0,B,p+2);%规划34个不同的梅尔刻度。但只有32个梅尔滤波器
fm=700*(10.^(mm/2595)-1);%将Mel频率转换为频率
W2=N/2+1;%fs/2内对应的FFT点数,2049个频率分量k=((N+1)*fm)/fs;%计算28个不同的k值
hm=zeros(p,N);%创建hm矩阵
df=fs/N;
freq=(0:N-1)*df;%采样频率值%绘制梅尔滤波器
for i=2:p+1%取整n0=floor(k(i-1));  %floor为向下取整，如果取整，表示三角滤波器组的各中心点被置于FFT 数据对应的的点位频率上，所以频率曲线的三角形顶点正好在FFT 数据点上，即美尔滤波器所有顶点在同一高度。n1=floor(k(i));n2=floor(k(i+1));
%     n0=k(i-1);  %不取整代表直接按精确指定的中心频率点来计算的，这时滤波器中心频率就不一定在FFT 数据点对应位置，按FFT 点位所得曲线在低频段就不是完美的三角形，即美尔滤波器所有顶点不在同一高度。
%     n1=k(i);
%     n2=k(i+1);%要知道k(i)分别代表的是每个梅尔值在新的范围内的映射，其取值范围为：0-N/2%以下实现求取三角滤波器的频率响应。for j=1:Nif n0<=j & j<=n1hm(i-1,j)=(j-n0)/(n1-n0);elseif n1<=j & j<=n2hm(i-1,j)=(n2-j)/(n2-n1);endend%此处求取H1(k)结束。
end%梅尔滤波器滤波H=E*hm';%对H作自然对数运算%因为人耳听到的声音与信号本身的大小是幂次方关系，所以要求个对数for i=1:afor j=1:pH(i,j)=log(H(i,j));endend%作离散余弦变换
Fbank = H';
Fbank1 = dct(Fbank);
mfcc = Fbank1';
J=mfcc(:,(1:p));
feat=mfcc;%求一阶差分系数dtfeat=0;dtfeat=zeros(size(feat));%默认初始化for i=3:a-2dtfeat(i,:)=-2*feat(i-2,:)-feat(i-1,:)+feat(i+1,:)+2*feat(i+2,:); enddtfeat=dtfeat/3;
%求取二阶差分系数
%二阶差分系数就是对前面产生的一阶差分系数dtfeat再次进行操作。dttfeat=0;dttfeat=zeros(size(dtfeat));%默认初始化for i=3:a-2dttfeat(i,:)=-2*dtfeat(i-2,:)-dtfeat(i-1,:)+dtfeat(i+1,:)+2*dtfeat(i+2,:); enddttfeat=dttfeat/10;%这里的10是根据数据确定的，默认为2%将得到的mfcc的三个参数feat、dtfeat、dttfeat拼接到一起mfcc_final=[feat,dtfeat,dttfeat];%拼接完成imagesc(mfcc_final');
axis xy;   %y轴从下往上
colorbar;
xlabel('语音分帧帧数');
ylabel('MFCC-D-A参数维数');

最后得到的mfcc就是MFCC参数。
MFCC程序求取过程中的分帧函数和上面的enframe函数一样。
画出的加上一阶差分和二阶差分的MFCC-D-A参数图：

以上就是所有MFCC及MFCC-D-A参数的matlab求取过程。

4.MFCC调包

matlab中还自带了mfcc的程序包（本文matlab版本为2019a，较低版本中可能不含有如下的mfcc包），其调用方法为：

[mfcc,delta_mfcc,delta_delta_mfcc]= mfcc(sph,Fs,'WindowLength',512,'OverlapLength',256 ,'NumCoeffs',32,'LogEnergy','Ignore');  % 对语音提取mfcc特征参数

其中等号左边mfcc为mfcc参数，delta_mfcc为mfcc的一阶差分参数，delta_delta_mfcc为mfcc的二阶差分参数。
等号右边括号中
sph为输入的语音信号
Fs为采样频率
‘WindowLength’,512代表分帧的帧长为512
‘OverlapLength’,256代表重叠部分为256
‘NumCoeffs’,32表示m梅尔滤波器个数取32个
‘LogEnergy’,‘Ignore’表示不要能量，如果为’LogEnergy’,‘Append’表示将能量放在mfcc参数前面，如果为’LogEnergy’,'Replace’表示用能量代替mfcc的第一维系数
使用方法：

[sph,Fs]=audioread('jian_2_1000.wav');%读取语音wav文件，本文语音长度3秒，单声道，采样频率16000
[mfcc,delta_mfcc,delta_delta_mfcc]= mfcc(sph,Fs,'WindowLength',512,'OverlapLength',256 ,'NumCoeffs',32,'LogEnergy','Ignore');  % 对语音提取mfcc特征参数

有关于MFCC的问题可联系qq：1479115257

matlab实现MFCC相关推荐

【语音识别】基于matlab GUI MFCC+VAD端点检测智能语音门禁系统【含Matlab源码 451期】
⛄一.MFCC简介 1 引言语音识别是一种模式识别, 就是让机器通过识别和理解过程把语音信号转变为相应的文本或命令的技术.语音识别技术主要包括特征提取技术.模式匹配准则及模型训练技术3个方面.目前一 ...
梯度倒谱matlab程序,MFCC梅尔倒谱参数及matlab代码
右端就会是连续的,那就可以不需要乘上汉明窗了.但是在实作上,由于基本周期的计算会需要额外的时间,而且也容易算错,因此我们都用汉明窗来达到类似的效果. 5.三角带通滤波器(Triangular Band ...
基于节拍谱的语音音乐分类模型
节拍谱的取得方式为按照论文<THE BEAT SPECTRUM:A NEW APPROACH TO RHYTHM ANALYSIS>获得,但是在分类中,有点阴差阳错的计算出不同的结果. 按 ...
MFCC特征参数提取（一）（基于MATLAB和Python实现）
1.MFCC概述在语音识别(Speech Recognition)和话者识别(Speaker Recognition)方面,最常用到的语音特征就是梅尔倒谱系数(Mel-scale Frequency ...
语音识别中的MFCC的提取原理和MATLAB实现
一.首先让我们借用并澄清几个语音学中的概念 1.临界频带与听觉掩蔽听觉临界频带:设纯音频率为,用噪声(设频率为)掩蔽纯音时,在噪声湮没的纯音的过程中,起作用的是频率在以内的噪声,称为临界频带.即当噪 ...
TensorFlow1.14或TensorFlow2内部获取mfcc原理探索（matlab复现或python复现）
MFCC概述声音是因为物体振动而产生的声波,是可以被人或动物的听觉器官所感知的波动现象.声音有多种特性,比如音色.音调.响度.频率.人耳是能够通过这些特征区分出声音的种类,那么如何让机器去区别不同种 ...
【语音识别】基于MFCC的小波变换DTW实现说话人识别matlab代码
1 简介小波变换的发展为语音信号提供了新的处理方法与技术,从而使语音处理技术取得了较快的发展.说话人识别提取说话人的语音特征对说话人的身份进行确认或辨认.语音识别研究领域的一个重要研究方向,就是从语 ...
【语音识别】基于MFCC和MEL倒频系数实现声纹识别附matlab代码
1 内容介绍提出了以Mel频率倒谱系数(Mel Frequency Cepstrum Coefficients,MFCC)和MEL倒频系数作为特征提取技术,以KNN作为分类器的语音识别方法,实验结果 ...
Matlab mfcc函数参数详解（英文附例）
Matlab mfcc函数参数详解其实可以直接打开源代码看哈. %MFCC Extract the mfcc, log-energy, delta, and delta-delta of audio ...

matlab实现MFCC