耳蜗实质上相当于一个滤波器组，耳蜗的滤波作用是在对数频率尺度上进行的，在1000HZ下，人耳的感知能力与频率成线性关系；而在1000HZ以上，人耳的感知能力与频率不构成线性关系，而更偏向于对数关系，这就使得人耳对低频信号比高频信号更敏感。Mel频率的提出是为了方便人耳对不同频率语音的感知特性的研究。频率与Mel频率的转换公式为：
MFCC在一定程度上模拟了人耳对语音的处理特点，应用了人耳听觉感知方面的研究成果，采用这种技术语音识别系统的性能有一定提高。
MFCC参数的提取
1、 预加重处理
预加重处理其实是一个高通滤波器，该高通滤波顺的传递函数为：
其中的取值为0.97，该高通滤波器作用是滤去低频，使语音信号的高频特性更加突现。
2、 分帧及加窗处理
由于语音信号只在较短的时间内呈现平稳性（一般认为10-30ms），因此将语音信号划分为一个一个的短时段即一帧。同时为避免丢失语音信号的动态信息，相邻帧之间要有一段重叠区域，重叠区域一段为帧长的1/2或1/3。然后再将每帧乘上窗函数，以增加每帧左端和右端的连续性。
3、 各帧信号的FFT变换
对分帧加窗后的各帧信号进行FFT变换得到各帧的频谱。并对语音信号的频谱取模平方得到语音信号的功率谱。
4、 三角滤波器系数的求取
定义若干个带通三角滤波器(k)，0<=m<=M，M为滤波器个数，其中心频率为f(m),每个带通三角滤波器的频率响应为
且满足Mel(f(m))-Mel(f(m-1))=Mel(f(m+1))-Mel(f(m))
求得滤波系数为m(i),i=1,…,p，p为滤波器阶数
5、 三角滤波并进行离散余弦变换DCT
C(i)即为所要求提取的特征参数。
特征参数的识别
特征参数的识别主要采用BP神经网络算法进行预测，而在预测前需要用一定数量的样本对网络进行训练，使网络具有联想记忆和预测能力。
网络训练步骤如下：
（1） 网络初始化。确定网络输入层、隐层、输出层数目，输出层到隐层的连接权值及隐层到输出层的连接权值
，同时初始化隐层阈值a和输出层阈值b;
（2） 隐层的输出计算。隐层输出式中，l为隐含层节点数；f为隐含层激励函数，本实验选取函数为:
（3） 输出层输出计算。根据隐含层输出H，连接权值和阈值b，计算BP神经网络预测输出O.
（4） 误差计算。根据网络预测输出O和期望输出Y，计算网络预测误差e.
（5） 权值更新。根据网络预测误差e更新网络连接权值,
其中，为学习速率。
（6） 阈值更新。根据网络预测误差e更新网络节点阈值a,b.
（7） 判断算法迭代是否结束，若没有结束，返回步骤（2）。
BP神经网络分类
用训练好的BP神经网络分类语音特征信号，根据分类结果分析BP神经网络分类能力。
将四种音乐风格的数据各1500组共6000组MEL特征数据随机抽取4000组进行训练网络，剩下的2000组特征向量进行辨识，得到各类风格的正确率如下：
rightridio =
1.0000 0.9261 0.9129 0.8399
下图为BP网络误差：
MATLAB程序如下：
clear;
clc;
%%%%%%%%采集的四种音乐各500000个数据%%%%%%%%%%%%%%%%%%%%%
load f:课程voice_recoc1 x1
load f:课程voice_recoc2 x2
load f:课程voice_recoc3 x3
load f:课程voice_recoc4 x4
%%%%%%%%%%%%%%%对语音信号进行预加重处理%%%%%%%%%%%%%%%%%%
len=length(x1);
heigt=0.98;
for i=2:len
x1(i)=x1(i)-heigt*x1(i-1);
end
for i=2:len
x2(i)=x2(i)-heigt*x2(i-1);
end
for i=2:len
x3(i)=x3(i)-heigt*x3(i-1);
end
for i=2:len
x4(i)=x4(i)-heigt*x4(i-1);
end
%%%%%%%%%%%%%%MEL三角滤波参数%%%%%%%%%%%%%%%%%%%%%%%%%%%
fh=20000;
melf=2595*log(1+fh/700);
M=24;
i=0:25;
f=700*(exp(melf/2595*i/(M+1))-1);
N=256;
for m=1:24
for k=1:256
x=fh*k/N;
if (f(m)<=x)&&(x<=f(m+1))
F(m,k)=(x-f(m))/(f(m+1)-f(m));
else if (f(m+1)<=x)&&(x<=f(m+2))
F(m,k)=(f(m+2)-x)/(f(m+2)-f(m+1));
else
F(m,k)=0;
end
end
end
end
m=N/2;
for k=1:12
n=0:23;
dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24));
end
count=floor(length(x1)/m);
%%%%%%%%%%%%%%%四种语音的特征参数的求取%%%%%%%%%%%
c1=zeros(count,12);
for i=1:count-2
x_frame=x1(m*(i-1)+1:m*(i-1)+N);
Fx=abs(fft(x_frame));
s=log(Fx.^2*F');
c1(i,:)=s*dctcoef';
end
c1=zeros(count,12);
for i=1:count-2
x_frame=x2(m*(i-1)+1:m*(i-1)+N);
Fx=abs(fft(x_frame));
s=log(Fx.^2*F');
c2(i,:)=s*dctcoef';
end
c3=zeros(count,12);
for i=1:count-2
x_frame=x3(m*(i-1)+1:m*(i-1)+N);
Fx=abs(fft(x_frame));
s=log(Fx.^2*F');
c3(i,:)=s*dctcoef';
end
c4=zeros(count,12);
for i=1:count-2
x_frame=x4(m*(i-1)+1:m*(i-1)+N);
Fx=abs(fft(x_frame));
s=log(Fx.^2*F');
c4(i,:)=s*dctcoef';
end
%save c1 c1
%save c2 c2
%save c3 c3
%save c4 c5
%四个特征信号矩阵合成一个矩阵
data(1:1500,:)=c1(1:1500,:);
data(1501:3000,:)=c2(1:1500,:);
data(3001:4500,:)=c3(1:1500,:);
data(4501:6000,:)=c4(1:1500,:);
%%%%%%%%%%%特征信号第一列为所属类别%%%%%%%%%%%%%%
for i=1:6000
if (i>=1)&&(i<=1500)
data(i,1)=1;
else if(i>=501)&&(i<=3000)
data(i,1)=2;
else if (i>=1001)&&(i<=4500)
data(i,1)=3;
else
data(i,1)=4;
end
end
end
end
%从1到2000间随机排序
k=rand(1,6000);
[m,n]=sort(k);
%输入输出数据
input=data(:,2:12);
output1 =data(:,1);
%把输出从1维变成4维
for i=1:6000
switch output1(i)
case 1
output(i,:)=[1 0 0 0];
case 2
output(i,:)=[0 1 0 0];
case 3
output(i,:)=[0 0 1 0];
case 4
output(i,:)=[0 0 0 1];
end
end
%随机提取1500个样本为训练样本，500个样本为预测样本
input_train=input(n(1:4000),:)';
output_train=output(n(1:4000),:)';
input_test=input(n(4001:6000),:)';
output_test=output(n(4001:6000),:)';
%输入数据归一化
[inputn,inputps]=mapminmax(input_train);
%% 网络结构初始化
innum=11;
midnum=12;
outnum=4;
%权值初始化
w1=rands(midnum,innum);
b1=rands(midnum,1);
w2=rands(midnum,outnum);
b2=rands(outnum,1);
w2_1=w2;w2_2=w2_1;
w1_1=w1;w1_2=w1_1;
b1_1=b1;b1_2=b1_1;
b2_1=b2;b2_2=b2_1;
%学习率
xite=0.1;
alfa=0.01;
%% 网络训练
for ii=1:10
E(ii)=0;
for i=1:1:4000
%% 网络预测输出
x=inputn(:,i);
% 隐含层输出
for j=1:1:midnum
I(j)=inputn(:,i)'*w1(j,:)'+b1(j);
Iout(j)=1/(1+exp(-I(j)));
end
% 输出层输出
yn=w2'*Iout'+b2;
%% 权值阀值修正
%计算误差
e=output_train(:,i)-yn;
E(ii)=E(ii)+sum(abs(e));
%计算权值变化率
dw2=e*Iout;
db2=e';
for j=1:1:midnum
S=1/(1+exp(-I(j)));
FI(j)=S*(1-S);
end
for k=1:1:innum
for j=1:1:midnum
dw1(k,j)=FI(j)*x(k)*(e(1)*w2(j,1)+e(2)*w2(j,2)+e(3)*w2(j,3)+e(4)*w2(j,4));
db1(j)=FI(j)*(e(1)*w2(j,1)+e(2)*w2(j,2)+e(3)*w2(j,3)+e(4)*w2(j,4));
end
end
w1=w1_1+xite*dw1';
b1=b1_1+xite*db1';
w2=w2_1+xite*dw2';
b2=b2_1+xite*db2';
w1_2=w1_1;w1_1=w1;
w2_2=w2_1;w2_1=w2;
b1_2=b1_1;b1_1=b1;
b2_2=b2_1;b2_1=b2;
end
end
%% 语音特征信号分类
inputn_test=mapminmax('apply',input_test,inputps);
for ii=1:1
for i=1:200000
%隐含层输出
for j=1:1:midnum
I(j)=inputn_test(:,i)'*w1(j,:)'+b1(j);
Iout(j)=1/(1+exp(-I(j)));
end
fore(:,i)=w2'*Iout'+b2;
end
end
%% 结果分析
%根据网络输出找出数据属于哪类
for i=1:2000
output_fore(i)=find(fore(:,i)==max(fore(:,i)));
end
%BP网络预测误差
error=output_fore-output1(n(4001:6000))';
%画出预测语音种类和实际语音种类的分类图
figure(1)
plot(output_fore,'r')
hold on
plot(output1(n(4001:6000))','b')
legend('预测语音类别','实际语音类别')
%画出误差图
figure(2)
plot(error)
title('BP网络分类误差','fontsize',12)
xlabel('语音信号','fontsize',12)
ylabel('分类误差','fontsize',12)
%print -dtiff -r600 1-4
k=zeros(1,4);
%找出判断错误的分类属于哪一类
for i=1:2000
if error(i)~=0
[b,c]=max(output_test(:,i));
switch c
case 1
k(1)=k(1)+1;
case 2
k(2)=k(2)+1;
case 3
k(3)=k(3)+1;
case 4
k(4)=k(4)+1;
end
end
end
%找出每类的个体和
kk=zeros(1,4);
for i=1:2000
[b,c]=max(output_test(:,i));
switch c
case 1
kk(1)=kk(1)+1;
case 2
kk(2)=kk(2)+1;
case 3
kk(3)=kk(3)+1;
case 4
kk(4)=kk(4)+1;
end
end
%正确率
rightridio=(kk-k)./kk