matlab实现MFCC

article/2025/10/20 22:16:50

MFCC

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

  1. 首先对语音进行预处理。
    预处理又包括对语音进行预加重、分帧、加窗。
  2. 快速傅里叶变换
    对分帧加窗后的每帧语音数据进行fft变换。
  3. 计算谱线能量
    对fft变换后的每帧信号取平方。
  4. 梅尔滤波器组滤波
    将求出的每帧谱线能量谱与梅尔滤波器组相乘。
  5. 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


http://chatgpt.dhexx.cn/article/vwlNAbIO.shtml

相关文章

理解MFCC

文章目录 提取音频的整体步骤预加重分帧加窗FFT(快速傅里叶变换)声谱图&#xff08;Spectrogram&#xff09;梅尔频谱和梅尔倒谱 倒谱&#xff08;cepstrum&#xff09;就是一种信号的傅里叶变换经对数运算后再进行傅里叶反变换得到的谱记住一句话&#xff0c;在梅尔频谱上做倒…

MFCC详细步骤及解析

MFCC(Mel-frequency cepstral coefficients):梅尔频率倒谱系数。梅尔频率是基于人耳听觉特性提出来的&#xff0c; 它与Hz频率成非线性对应关系。梅尔频率倒谱系数(MFCC)则是利用它们之间的这种关系&#xff0c;计算得到的Hz频谱特征。主要有 以下几个步骤&#xff1a;预加重&a…

MFCC理解

MFCC 在语音识别&#xff08;SpeechRecognition&#xff09;和话者识别&#xff08;SpeakerRecognition&#xff09;方面&#xff0c;最常用到的语音特征就是梅尔倒谱系数&#xff08;Mel-scaleFrequency Cepstral Coefficients&#xff0c;简称MFCC&#xff09;。根据人耳听觉…

MFCC特征介绍

MFCC特征介绍 在语音识别技术中&#xff0c;需要提取音频的特征&#xff0c;然后就可以使用该音频进行模型的训练或者是进行识别&#xff0c;目前很常用的一种特征叫做MFCC特征&#xff0c;又叫做梅尔倒谱系数特征。MFCC特征保留了语义相关的一些内容&#xff0c;过滤掉了诸如…

深入理解MFCC(梅尔频率倒谱系数)

从倒谱图出发 MFCC是Mel Frequency Cepstral Coefficient的简称&#xff0c;要理解MFCC特征&#xff0c;就需要先明白这里引入的一个新的概念——Cepstral&#xff0c;这个形容词的名词形式为Cepstrum&#xff0c;即倒谱图&#xff08;频谱图Spectrum前四个字母倒着拼&#xf…

MFCC特征提取

在语音识别方面&#xff0c;最常用到的语音特征就是梅尔倒谱系数&#xff08;Mel-scaleFrequency Cepstral Coefficients&#xff0c;简称MFCC&#xff09;。 MFCC的提取过程包括预处理、快速傅里叶变换、Mei滤波器组、对数运算、离散余弦变换、动态特征提取等步骤。 1.预处理 …

MFCC算法讲解及实现(matlab)

史上最详细的MFCC算法实现&#xff08;附测试数据&#xff09; 1.matlab安装voicebox语音包2.MFCC原理讲解3.MFCC算法设计实现&#xff08;matlab&#xff09;3.1 .wav格式语音文件提取【x(200000*1)】3.2 预加重【x(200000*1)】3.3 分帧{S(301*1103)}3.4 加窗{C(301*1103)}3.5…

Parquet encoding

Dictionary encoding

Parquet原理剖析

行存VS列存 广义的数据分析系统大致分为可以分为计算层、数据格式层和存储层。 计算层主要负责数据查询的介入和各种逻辑计算&#xff0c;如&#xff1a;MR、Spark、Flink。 存储层承载数据持久化存储&#xff0c;以文件语义或类似文件语义(对象存储)对接计算层。 数据格式层&…

Spark 实战 - 3.一文搞懂 parquet

一.引用 parquet 文件常见于 Spark、Hive、Streamin、MapReduce 等大数据场景&#xff0c;通过列式存储和元数据存储的方式实现了高效的数据存储与检索&#xff0c;下面主要讲 parquet 文件在 spark 场景下的存储&#xff0c;读取与使用中可能遇到的坑。 二.Parquet 加载方式 …

Spark Parquet使用

Spark SQL下的Parquet使用最佳实践和代码实战 分类&#xff1a; spark-sql&#xff08;1&#xff09; 一、Spark SQL下的Parquet使用最佳实践 1&#xff09;过去整个业界对大数据的分析的技术栈的Pipeline一般分为以下两种方式&#xff1a; a&#xff09;Data Source -> HD…

Arrow 之 Parquet

Parquet-format 左边是文件开头及具体的数据&#xff0c; 右边是文件结尾的 Footer Metadata There are three types of metadata: file metadata, column (chunk) metadata and page header metadata. All thrift structures are serialized using the TCompactProtocol. Co…

parquet存入mysql_解密列存 parquet

在做数据分析的时候,相对于传统关系型数据库,我们更倾向于计算列之间的关系。在使用传统关系型数据库时,基于此的设计,我们会扫描很多我们并不关心的列,这导致了查询效率的低下,大部分数据库 io 比较低效。因此目前出现了列式存储。Apache Parquet 是一个列式存储的文件格…

Parquet原理

在互联网大数据应用场景下&#xff0c;通常数据量很大且字段很多&#xff0c; 但每次查询数据只针对其中的少数几个字段&#xff0c;这时候列式存储是极佳的选择。 列式存储要解决的问题&#xff1a; 把IO只给查询需要用到的数据 只加载需要被计算的列空间节省 列式的压缩效…

parquet--golang使用

github 其实如果不适用一些可视化工具解析parquet文件&#xff0c;不太好看parquet文件内部正常应该是什么样的。但是使用一些可视化工具的话&#xff0c;可以发现&#xff0c;parquet文件会像表格&#xff0c;如excel文件&#xff0c;csv文件那样&#xff0c;排列数据。通过结…

Parquet

动机 创建Parquet是利用压缩性,高效的列式存储来在Haddop生态圈任何项目中应用. 记住Parquet是构建在复杂嵌套的数据结构, 并且使用记录分解和集成的算法在Dremely论文中描述.我们相信这种方法是更强大的的可以非常简单的使嵌套命令空间的扁平化. Parquet构建可以非常高效的…

Parquet 存储格式

1.介绍 Apache Parquet 是 Hadoop 生态圈中一种新型列式存储格式&#xff0c;它可以兼容 Hadoop 生态圈中大多数计算框架(Mapreduce、Spark 等)&#xff0c;被多种查询引擎支持&#xff08;Hive、Impala、Drill 等&#xff09;&#xff0c;并且它是语言和平台无关的。 2.特点…

parquet 简介

参考文章&#xff1a;parquet 简介 Parquet原理 【2019-05-29】Parquet 简介 Apache Parquet是一种能够有效存储嵌套数据的列式存储格式。 面向分析型业务的列式存储格式 由 Twitter 和 Cloudera 合作开发&#xff0c;2015 年 5 月从 Apache 的孵化器里毕业成为 Apache 顶…

Parquet文件详解

1、parquet文件简介 Apache Parquet是Apache Hadoop生态系统的一种免费的开源面向列的数据存储格式。 它类似于Hadoop中可用的其他列存储文件格式&#xff0c;如RCFile格式和ORC格式。 Apache Parquet 是由 Twitter 和 Cloudera 最先发起并合作开发的列存项目&#xff0c;也是…

Gson解析json数据

gson是谷歌推出的&#xff0c;除此之外还有阿里的FastJson&#xff0c;官方json和jackjson。下面通过一个实例来讲解使用gson来解析json数据&#xff1a; 1.先做好准备工作&#xff0c;在网上下载Gson的jar包&#xff0c;放到工程的libs(没有此目录的话自己建一个)目录下: ht…