梅尔倒谱系数
MFCC
梅尔频率倒谱系数的分析是基于人的听觉特性机理，即根据人的听觉实验结果来分析语音的频谱。因为人耳所能听到的声音高低与声音的频率并不成线性正比关系，所以用mel频率尺度更符合人耳的听觉特性。
预加重部分用一个一阶高通滤波器，目的是为了补偿高频分量的损失，提升高频分量。滤波器常设为：
分帧：把每一帧信号当作稳定信号来处理。
加窗：减少频谱泄漏。
FFT:把时域转化为频域，计算谱线能量。
mel滤波：
MEL滤波器滤波：
梅尔频率尺度和实际频率的对应关系
美尔滤波器的传递函数
在梅尔频率下，每个滤波器直接的间隔梅尔滤波是一样的，实际频率不一样。
构造梅尔滤波器组
**设计梅尔滤波器思路：**自己的理解,不一定对。
（1）确定要设计滤波器的个数p，最高频率fh，最低频率fl，fft采样点数。
（2）将最高频率，最低频率转化为梅尔频率尺度bl,bh.
（3）在梅尔尺度中，从0到(bh-bl)分为p+1段，一共有p+2个点。
（4）将得到的p+2个点转化为实际频率。
（5）计算实际频率对应的fft点数，2-p+1个点为每个滤波器的中心频率。
（6）根据fft点数，和梅尔滤波器传递函数H(k)画出三角窗。
clc;
clear;
%Mel滤波
fs = 8000; %采样频率
fl = 0; %fl是设计滤波器的最低频率
fh = fs / 2; %fh是设计滤波器的最高频率
bl = 1125 * log(1 + fl/700); %最低频率对应的Mel频率
bh = 1125 * log(1 + fh/700); %最高频率对应的Mel频率
p = 24; %在fl和fh之间设计Mel滤波器的个数
N = 256; %FFT点数,N越大,频率分辨率越高,误差越小
MelF = linspace(0, bh-bl, p+2); %在0至bh-bl的Mel频率范围内产生p+2个Mel频率值
F = 700 * (exp(MelF/1125) - 1); %将上一步产生的p+2个Mel频率值转化为p+2个实际频率值
df = fs / N; %计算频率分辨率
n = N/2 + 1; %计算fs/2内对应的FFT点数
f = (0:n-1) * df; %计算频率序列,共有n个频率点
bank = zeros(24, n); %生成24行n列的全零矩阵
for m1 = 2:p+1 %循环起始数为2，因为下方语句中用到m-1的值,其中m-1必须为正整数,所以m为大于等于2的正整数
%F为实际频率,分别找到三角滤波器左中右对应的频率
F_left = F(m1-1);
F_mid = F(m1);
F_right = F(m1+1);
%使用ceil（）函数向上取整，计算实际频率对应的FFT点数
n_left = ceil(F_left / df); %频率F_left对应的FFT点数
n_mid = ceil(F_mid / df); %频率F_mid对应的FFT点数
n_right = ceil(F_right / df); %频率F_right对应的FFT点数
%用三角窗的频率响应函数作图
for k = 1:n %根据FFT点数画三角窗
if k>=n_left & k<=n_mid
bank(m1-1,k) = (k-n_left)/(n_mid-n_left); %画第m-1个滤波器的左半部分
elseif k>n_mid & k<=n_right
bank(m1-1,k) = (n_right-k)/(n_right-n_mid); %画第m-1个滤波器的右半部分
end
end
plot(f,bank(m-1,:))
hold on;
end
经过滤波器滤波之后，输出和滤波器一样的维数，这里我用了24个滤波器，所有输出的维数也是24维
每一帧信号经过滤波器，输入的是一个值。一个数。
clc;
clear;
[x,fs]=audioread('C3_4_y_1.wav');
%Mel滤波
fs = 8000; %采样频率
fl = 0; %fl是设计滤波器的最低频率
fh = fs / 2; %fh是设计滤波器的最高频率
bl = 1125 * log(1 + fl/700); %最低频率对应的Mel频率
bh = 1125 * log(1 + fh/700); %最高频率对应的Mel频率
p = 24; %在fl和fh之间设计Mel滤波器的个数
N = 256; %FFT点数,N越大,频率分辨率越高,误差越小
MelF = linspace(0, bh-bl, p+2); %在0至bh-bl的Mel频率范围内产生p+2个Mel频率值
F = 700 * (exp(MelF/1125) - 1); %将上一步产生的p+2个Mel频率值转化为p+2个实际频率值
df = fs / N; %计算频率分辨率
n = N/2 + 1; %计算fs/2内对应的FFT点数
f = (0:n-1) * df; %计算频率序列,共有n个频率点
bank = zeros(24, n); %生成24行n列的全零矩阵
for m1 = 2:p+1 %循环起始数为2，因为下方语句中用到m-1的值,其中m-1必须为正整数,所以m为大于等于2的正整数
%F为实际频率,分别找到三角滤波器左中右对应的频率
F_left = F(m1-1);
F_mid = F(m1);
F_right = F(m1+1);
%使用ceil（）函数向上取整，计算实际频率对应的FFT点数
n_left = ceil(F_left / df); %频率F_left对应的FFT点数
n_mid = ceil(F_mid / df); %频率F_mid对应的FFT点数
n_right = ceil(F_right / df); %频率F_right对应的FFT点数
%用三角窗的频率响应函数作图
for k = 1:n %根据FFT点数画三角窗
if k>=n_left & k<=n_mid
bank(m1-1,k) = (k-n_left)/(n_mid-n_left); %画第m-1个滤波器的左半部分
elseif k>n_mid & k<=n_right
bank(m1-1,k) = (n_right-k)/(n_right-n_mid); %画第m-1个滤波器的右半部分
end
end
end
% 归一化mel滤波器组系数
bank=full(bank);
bank=bank/max(bank(:));
% plot(f,bank(m-1,:))
% hold on;
x=filter([1 -0.9375],1,x);
%分帧
S=enframe(x,256,160);
n2=fix(256/2)+1;
for k=1:24
n=0:23;
dctxs(k,:)=cos((2*n+1)*k*pi/(2*24));
end
for i=1:size(x,1)
y = x(i,:);
s = y' .* hamming(256);
t = abs(fft(s));
t = t.^2;
bank=melbankm(24,256,fs,0,0.5,'m');
% 归一化Mel滤波器组系数
bank=full(bank);
bank=bank/max(bank(:));
c=log(bank*t(1:n2));
c1=dctxs*c;
w=1+12*sin(pi*[1:24]./24);
w=w/max(w);
c2=c1.*w';%大部分的信号数据一般集中在变换后的低频区，所以对每一帧只取前13个数据就好了。
m(i,:)=c2';
end
%差分系数 %%%N维MFCC参数（N/3MFCC系数+ N/3一阶差分参数+ N/3二阶差分参数）+帧能量（此项可根据需求替换)
dtm = zeros(size(m));
for i=3:size(m,1)-2
dtm(i,:) = -2*m(i-2,:) - m(i-1,:) + m(i+1,:) + 2*m(i+2,:);%标准的倒谱参数MFCC只反映了语音参数的静态特性，语音的动态特性可以用这些静态特征的差分谱来描述。把动、静态特征结合起来才能有效提高系统的识别性能
end
dtm = dtm / 3;
%合并MFCC参数和一阶差分MFCC参数
ccc = [m dtm];%N维MFCC参数（N/3MFCC系数+ N/3一阶差分参数+ N/3二阶差分参数）+帧能量（此项可根据需求替换）
%去除首尾两帧，因为这四帧的一阶差分参数为0
ccc = ccc(3:size(m,1)-2,:);
取对数是因为人耳听到的声音与信号本身的大小是幂次方关系，所以要求个对数.
bank输出的矩阵形式为：滤波器个数*fft点数/2-1。t为fft点数，只取fft点数/2-1。
预加重前的信号
预加重后的信号