1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要 在语音识别（Speech Recognition）和话者识别（Speaker Recognition）方面，最常用到的语音特征就是梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients，简称MFCC）。根据人耳听觉机理的研究发现，人耳对不同频率的声波有不同的听觉敏感度。从200Hz到5000Hz的语音信号对语音的清晰度影响对大。两个响度不等的声音作用于人耳时，则响度较高的频率成分的存在会影响到对响度较低的频率成分的感受，使其变得不易察觉，这种现象称为掩蔽效应。由于频率较低的声音在内耳蜗基底膜上行波传递的距离大于频率较高的声音，故一般来说，低音容易掩蔽高音，而高音掩蔽低音较困难。在低频处的声音掩蔽的临界带宽较高频要小。所以，人们从低频到高频这一段频带内按临界带宽的大小由密到疏安排一组带通滤波器，对输入信号进行滤波。将每个带通滤波器输出的信号能量作为信号的基本特征，对此特征经过进一步处理后就可以作为语音的输入特征。由于这种特征不依赖于信号的性质，对输入信号不做任何的假设和限制，又利用了听觉模型的研究成果。因此，这种参数比基于声道模型的LPCC相比具有更好的鲁邦性，更符合人耳的听觉特性，而且当信噪比降低时仍然具有较好的识别性能。 梅尔倒谱系数（Mel-scale Frequency Cepstral Coefficients，简称MFCC）是在Mel标度频率域提取出来的倒谱参数，Mel标度描述了人耳频率的非线性特性，它与频率的关系可用下式近似表示：
美尔尺度是建立从人类的听觉感知的频率——Pitch到声音实际频率直接的映射。人耳对于低频声音的分辨率要高于高频的声音。通过把频率转换成美尔尺度，我们的特征能够更好的匹配人类的听觉感知效果。从频率到美尔频率的转换公式如下：
M(f)=1125ln(1+f/700)M(f)=1125ln(1+f/700)M(f)=1125ln(1+f/700)M(f)=1125ln(1+f/700)
而从美尔频率到频率的转换公式为：
M−1(m)=700(em/1125−1)M−1(m)=700(em/1125−1)M−1(m)=700(em/1125−1)M−1(m)=700(em/1125−1)
流程图：
预处理包括预加重、分帧、加窗。 （1）预加重 预加重的目的是提升高频部分，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求频谱。同时，也是为了消除发生过程中声带和嘴唇的效应，来补偿语音信号受到发音系统所抑制的高频部分，也为了突出高频的共振峰。预加重处理其实是将语音信号通过一个高通滤波器。
（2）分帧 由于语音信号的非平稳特性 和 短时平稳特性，将语音信号分分帧。一帧有N个采样点，如N的值为256或512，涵盖的时间约为20~30ms左右。为了避免相邻两帧的变化过大，平缓过度，因此会让两相邻帧之间有一段重叠区域，此重叠区域包含了M个取样点，通常M的值约为N的1/2或1/3。通常语音识别所采用语音信号的采样频率为8KHz或16KHz，以8KHz来说，若帧长度为256个采样点，则对应的时间长度是256/8000×1000=32ms。 （3）加窗 为了缓解频谱泄漏。将每一帧乘以一个窗函数，如汉明窗，海宁窗。假设分帧后的信号为S(n), n=0,1…,N-1, N为帧的大小。
语音信号的倒谱分析就是求信号倒谱特征参数的过程，可以通过同态处理来处理。同态处理实现了卷积关系变化为求和关系的分离操作。
由于即使同一个人不同时间发出同一个声音，也不可能具有相同的长度，因此就需要用到动态时间归正（DTW）算法。把时间归正和距离测度计算结合起来的一种非线性归正技术。DTW 本质上是一个简单的动态规划算法，是用来计算两个维数不同的向量之间的相似度的问题，即计算向量 M1 和 M2 的最短距离。是一种非常常用的语音匹配算法。
对两个不同维数的语音向量 m1 和 m2进行匹配（m1 和 m2 的每一维也是一个向量，是语音每一帧的特征值，这里利用的是 MFCC 特征）。设两个向量的长度为 M1 和 M2，则距离可以表示为：
那么，就可以这样进行匹配：
每一条从（1，1）到（M1,M2）路径都有一个累计距离称为路径的代价；
每一条路径都代表一种对齐情况；
代价最小的路径就是所求的对准路径。
- 定义一个代价函数 ， 表示从起始点(1,1)出发，到达（i，j）点最小代价路径的累计距离。有：
3.MATLAB核心程序
melNo=25; %number of mel filters
syms mel(freq)
mel(freq)=1127*log(1+freq/700); % mel scale
center=zeros(1,(melNo+2));
freq1=7000;
freq0=50;
center(1)=mel(freq0);
hbandwidth=(mel(freq1)-mel(freq0))/(melNo+1);% half of bandwidth of each mel filter in mel scale
for i=2:melNo+2
center(i)=center(i-1)+hbandwidth; % center of each filter in mel scale
end
for i=1:melNo+2
center(i)=(exp(center(i)/1127)-1)*700; %moving centers to normal frequency scale
end
mj=zeros(frameNo,melNo);
for i=1:melNo
for j=1:frameNo
for k=1:(points/2)+1
if frq(k)>center(i) && frq(k)<center(i+1)
melresult=STFf(j,k)*((1/(center(i+1)-center(i)))*(frq(k)-center(i)));
elseif frq(k)>center(i+1) && frq(k)<center(i+2)
melresult=STFf(j,k)*((-1/(center(i+2)-center(i+1)))*(frq(k)-center(i+2)));
else
melresult=0;
end
mj(j,i)=melresult^2+mj(j,i);
end
end
end
mj=log10(mj);