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:N
if n0<=j & j<=n1
hm(i-1,j)=(j-n0)/(n1-n0);
elseif n1<=j & j<=n2
hm(i-1,j)=(n2-j)/(n2-n1);
end
end
%此处求取hm结束。
end
%绘图,且每条颜色显示不一样
c=colormap(lines(p));%定义16条不同颜色的线条
% figure()
% set(gcf,'position',[180,160,1020,550]);%设置画图的大小
for i=1:p
plot(freq,hm(i,:),'-','color',c(i,:),'linewidth',2);%开始循环绘制每个梅尔滤波器
% plot(freq,hm(i,:),'-','linewidth',2.5);%开始循环绘制每个梅尔滤波器
hold on;
end
grid 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:N
if n0<=j & j<=n1
hm(i-1,j)=(j-n0)/(n1-n0);
elseif n1<=j & j<=n2
hm(i-1,j)=(n2-j)/(n2-n1);
end
end
%此处求取hm结束。
end
%绘图,且每条颜色显示不一样
c=colormap(lines(p));%定义16条不同颜色的线条
% figure()
% set(gcf,'position',[180,160,1020,550]);%设置画图的大小
for i=1:p
plot(freq,hm(i,:),'-','color',c(i,:),'linewidth',2);%开始循环绘制每个梅尔滤波器
% plot(freq,hm(i,:),'-','linewidth',2.5);%开始循环绘制每个梅尔滤波器
hold on;
end
grid 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:N2
y(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:a
C(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:N
t=abs(D(i,j));
E(i,j)=(t^2);
end
end
fl=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:N
if n0<=j & j<=n1
hm(i-1,j)=(j-n0)/(n1-n0);
elseif n1<=j & j<=n2
hm(i-1,j)=(n2-j)/(n2-n1);
end
end
%此处求取H1(k)结束。
end
%梅尔滤波器滤波
H=E*hm';
%对H作自然对数运算
%因为人耳听到的声音与信号本身的大小是幂次方关系,所以要求个对数
for i=1:a
for j=1:p
H(i,j)=log(H(i,j));
end
end
%作离散余弦变换
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
else
len = 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:N2
y(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:a
C(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:N
t=abs(D(i,j));
E(i,j)=(t^2);
end
end
fl=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:N
if n0<=j & j<=n1
hm(i-1,j)=(j-n0)/(n1-n0);
elseif n1<=j & j<=n2
hm(i-1,j)=(n2-j)/(n2-n1);
end
end
%此处求取H1(k)结束。
end
%梅尔滤波器滤波
H=E*hm';
%对H作自然对数运算
%因为人耳听到的声音与信号本身的大小是幂次方关系,所以要求个对数
for i=1:a
for j=1:p
H(i,j)=log(H(i,j));
end
end
%作离散余弦变换
Fbank = H';
Fbank1 = dct(Fbank);
mfcc = Fbank1';
J=mfcc(:,(1:p));
feat=mfcc;
%求一阶差分系数
dtfeat=0;
dtfeat=zeros(size(feat));%默认初始化
for i=3:a-2
dtfeat(i,:)=-2*feat(i-2,:)-feat(i-1,:)+feat(i+1,:)+2*feat(i+2,:);
end
dtfeat=dtfeat/3;
%求取二阶差分系数
%二阶差分系数就是对前面产生的一阶差分系数dtfeat再次进行操作。
dttfeat=0;
dttfeat=zeros(size(dtfeat));%默认初始化
for i=3:a-2
dttfeat(i,:)=-2*dtfeat(i-2,:)-dtfeat(i-1,:)+dtfeat(i+1,:)+2*dtfeat(i+2,:);
end
dttfeat=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