基于动态时间规整（DTW）的孤立字语音识别
注意！这个博客里给出的文件是 在matlab2012版本里使用的，而我学习时使用的时matlab2018b，因此发现了很多新版本不兼容的问题，但我没有给出修改后能直接在新版本中用的代码。
另外，matlab要安装voicebox才可以正常进行实验！
还有一件事，录音文件要修改好文件名和文件目录，不然会检测不到。
这个实验在没有同学和老师的指导和帮助下完成，一共花了大概两周多（也没有每天做啦），请珍惜劳动成果，不要直接拿去装逼。
有个关键函数提前说一下：
DTW函数：myDTW
功能：利用DTW完成模式匹配
调用格式：
[cost] = myDTW(F,R)
参数说明：输入参数F为模板的MFCC特征；R为当前待识别语音的MFCC。输出参数cost为F和R之间的匹配距离。
看看目录
下面是具体的程序清单：
function [cost] = myDTW(F,R) [r1,c1] = size(F); %模板的维度 [r2,c2] = size(R); %当前待识别语音的维度 distance = zeros(r1,r2); for n = 1:r1 for m=1:r2 FR = F(n,:)-R(m,:)； FR = FR.^2; distance(n,m) = sqrt(sum(FR))/c1; %采用欧式距离度量F和R的相似程度 end end D=zeros(r1+1,r2+1); D(1,:)=inf; D(:,1)=inf; D(1,1)=0; D(2:(r1+1),2:(r2+1)) =distance; %寻找整个过程的最短匹配距离 for i=1:r1; for j= 1:r2; [dmin] = min([D(i,j),D(i,j+1),D(i+1,j)]); D(i+1,j+1) =D(i+1,j+1) +dmin; end end cost =D(r1+1,r2+1); %返回最终距离,没有反向回溯
这个基于动态时间规整（DTW）的孤立字语音识别的实验，具体步骤有很多，都写在下面了：
（1）训练过程：运行setTemplates.m文件
1）实验的语音来自三个人，每人说一遍0~9这十个数字，记录下相应的语音，保存在 目录“SpeechData”中。以第一个人的语音“1”为例，首先，调用[speechIn1,FS1] =audioread(q)读入语音“1”，保存在speechIn1中，并获得采样率FS1。
2）对speechIn1进行预处理，首先调用my_vad函数进行端点检测，提取出有效语音段，然后调用mfccf函数对有效语音段提取mfcc特征。
3）以同样的方式处理其他训练语音，把第一/二/三个说话人的十个数字0~9语音分别保存s1/s2/s3结构体中，并最终保存在Vectors1.mat/Vectors2.mat/Vectors3.mat文件中。
（2）识别过程：运行matchTemplates.m文件
1） 首先程序提示“Press any key to start 2 seconds of speech recording…”,在MATLAB主窗口按任意键开始录入待识别的语音，程序将录入的语音保存在speechIn变量中。
2） 对speechIn进行预处理，包括端点检测和mfcc特征提取，这一过程和训练中完全相同，得到（NP）矩阵。其中N为帧数，P为维度，接着调用CMN函数，对每一维度分量进行标准化，结果保存在rMatrix中。
3） 调用DTWScores函数，将当前待识别语音的rMatrix与训练过程中保存的s1、s2、s3中 的0~9各数字匹配，此函数中调用myDTW的匹配过程，返回rMarix和每个说话人说的每一个数字的匹配分值，保存在Sco中（（3*10）的矢量）。
4） 将Sco从大到小进行排序，首先得到每一个模板匹配的两个最低分值对应的数字序号（有2*3=6个序号），而后返回出现频率最高的序号作为判决结果，从而完成过程。
但是讲这么多，还是很难把实验完成对不对，现在我把全部用到的都放出来
识别过程matchTemplates.m程序清单：
clear all;
close all;
ncoeff = 12; %MFCC参数阶数
N = 10; %10个数字
fs=16000; %采样频率
duration2 = 2; %录音时长
k = 3; %训练样本的人数
speech = audiorecorder(fs,16,1);
disp('Press any key to start 2 seconds of speech recording...');
pause
disp('Recording speech...');
recordblocking(speech,duration2) % duration*fs 为采样点数
speechIn=getaudiodata(speech);
disp('Finished recording.');
disp('System is trying to recognize what you have spoken...');
speechIn = my_vad(speechIn); %端点检测
rMatrix1 = mfccf(ncoeff,speechIn,fs); %采用MFCC系数作为特征矢量
rMatrix = CMN(rMatrix1); %归一化处理
Sco = DTWScores(rMatrix,N); %计算DTW值
[SortedScores,EIndex] = sort(Sco,2); %按行递增排序，并返回对应的原始次序
Nbr = EIndex(:,1:2) %得到每个模板匹配的2个最低值对应的次序
[Modal,Freq] = mode(Nbr(:)); %返回出现频率最高的数Modal及其出现频率Freq
Word = char('zero','One','Two','Three','Four','Five','Six','Seven','Eight','Nine');
if mean(abs(speechIn)) < 0.01
fprintf('No microphone connected or you have not said anything.\n');
else if (Freq <2) %频率太低不确定
fprintf('The word you have said could not be properly recognised.\n');
else
fprintf('You have just said %s.\n',Word(Modal,:));
end
训练过程setTemplates----程序清单：
%设置模板
ncoeff=12; %mfcc系数的个数
fMatrix1 = cell(1,10);
fMatrix2 = cell(1,10);
fMatrix3 = cell(1,10);
for i = 1:10
q = ['F:\ SpeechData\p1\' num2str(i-1) '.wav'];
[speechIn1,FS1] = audioread(q);
speechIn1 = my_vad(speechIn1);
fMatrix1(1,i) = {mfccf(ncoeff,speechIn1,FS1)};
end
for j = 1:10
q = ['F: \SpeechData\p2\' num2str(j-1) '.wav'];
[speechIn2,FS2] = audioread(q);
speechIn2 = my_vad(speechIn2);
fMatrix2(1,j) = {mfccf(ncoeff,speechIn2,FS2)};
end
for k = 1:10
q = ['F: \SpeechData\p2\' num2str(k-1) '.wav'];
[speechIn3,FS3] = audioread(q);
speechIn3 = my_vad(speechIn3);
fMatrix3(1,k) = {mfccf(ncoeff,speechIn3,FS3)};
end
%将数据保存为mat文件
fields = {'zero','One','Two','Three','Four','Five','Six','Seven','Eight','nine'};
s1 = cell2struct(fMatrix1, fields, 2); %fields项作为行
save Vectors1.mat -struct s1;
s2 = cell2struct(fMatrix2, fields, 2);
save Vectors2.mat -struct s2;
s3 = cell2struct(fMatrix3, fields, 2);
save Vectors3.mat -struct s3;
端点检测my_vad（）函数----程序清单：
function trimmed_X = my_vad(x) %端点检测；输入为录入语音，输出为有用信号 Ini = 0.1; %初始静默时间 Ts = 0.01; %窗的时长 Tsh = 0.005; %帧移时长 Fs = 16000; %采样频率 counter1 = 0; %以下四个参数用来寻找起始点和结束点 counter2 = 0; counter3 = 0; counter4 = 0; ZCRCountf = 0; %用于存储过零率检测结果 ZCRCountb = 0; ZTh = 40; %过零阈值 w_sam = fix(Ts*Fs); %窗口长度 o_sam = fix(Tsh*Fs); %帧移长度 lengthX = length(x); segs = fix((lengthX-w_sam)/o_sam)+1; %分帧数 sil = fix((Ini-Ts)/Tsh)+1; %静默时间帧数 win = hamming(w_sam); Limit = o_sam*(segs-1)+1; %最后一帧的起始位置 FrmIndex = 1:o_sam:Limit; %每一帧的起始位置 ZCR_Vector = zeros(1,segs); %记录每一帧的过零点数 %短时过零点 for t = 1:segs ZCRCounter = 0; nextIndex = (t-1)*o_sam+1; for r = nextIndex+1:(nextIndex+w_sam-1) if (x(r) >= 0) && (x(r-1) >= 0) elseif (x(r) > 0) && (x(r-1) < 0) ZCRCounter = ZCRCounter + 1; elseif (x(r) < 0) && (x(r-1) < 0) elseif (x(r) < 0) && (x(r-1) > 0) ZCRCounter = ZCRCounter + 1; end end ZCR_Vector(t) = ZCRCounter; end %短时平均幅度 Erg_Vector = zeros(1,segs); for u = 1:segs nextIndex = (u-1)*o_sam+1; Energy = x(nextIndex:nextIndex+w_sam-1).*win; Erg_Vector(u) = sum(abs(Energy)); end IMN = mean(Erg_Vector(1:sil)); %静默能量均值（噪声均值） IMX = max(Erg_Vector); %短时平均幅度的最大值 I1 = 0.03 * (IMX-IMN) + IMN; %I1，I2为初始能量阈值 I2 = 4 * IMN; ITL = 100*min(I1,I2); %能量阈值下限,前面系数根据实际情况更改得到合适结果 ITU = 10* ITL; %能量阈值上限 IZC = mean(ZCR_Vector(1:sil)); stdev = std(ZCR_Vector(1:sil)); %静默阶段过零率标准差 IZCT = min(ZTh,IZC+2*stdev); %过零率阈值 indexi = zeros(1,lengthX); indexj = indexi; indexk = indexi; indexl = indexi; %搜寻超过能量阈值上限的部分 for i = 1:length(Erg_Vector) if (Erg_Vector(i) > ITU) counter1 = counter1 + 1; indexi(counter1) = i; end end ITUs = indexi(1); %第一个能量超过阈值上限的帧 %搜寻能量超过能量下限的部分 for j = ITUs:-1:1 if (Erg_Vector(j) < ITL) counter2 = counter2 + 1; indexj(counter2) = j; end end start = indexj(1)+1; %第一级判决起始帧 Erg_Vectorf = fliplr(Erg_Vector);%将能量矩阵关于中心左右对称，若是一行向量相当于逆序 %重复上面过程相当于找结束帧 for k = 1:length(Erg_Vectorf) if (Erg_Vectorf(k) > ITU) counter3 = counter3 + 1; indexk(counter3) = k; end end ITUf = indexk(1); for l = ITUf:-1:1 if (Erg_Vectorf(l) < ITL) counter4 = counter4 + 1; indexl(counter4) = l; end end finish = length(Erg_Vector)-indexl(1)+1;%第一级判决结束帧 %从第一级判决起始帧开始进行第二判决（过零率）端点检测 BackSearch = min(start,25); for m = start:-1:start-BackSearch+1 rate = ZCR_Vector(m); if rate > IZCT ZCRCountb = ZCRCountb + 1; realstart = m; end end if ZCRCountb > 3 start = realstart; end FwdSearch = min(length(Erg_Vector)-finish,25); for n = finish+1:finish+FwdSearch rate = ZCR_Vector(n); if rate > IZCT ZCRCountf = ZCRCountf + 1; realfinish = n; end end if ZCRCountf > 3 finish = realfinish; end x_start = FrmIndex(start); %最终的起始位置 x_finish = FrmIndex(finish-1); %最终的结束位置 trimmed_X = x(x_start:x_finish);
特征提取mfccf（）函数----程序清单
function FMatrix=mfccf(num,s,Fs) %计算并返回信号s的mfcc参数及其一阶和二阶差分参数 %参数说明： num 是mfcc系数阶数，Fs为采样频率 N=512; % FFT 数 Tf=0.02; %窗口的时长 n=Fs*Tf; %每个窗口的长度 M=24; %M为滤波器组数 l=length(s); Ts=0.01; %帧移时长 FrameStep=Fs*Ts; %帧移 a=1; b=[1, -0.97]; %预加重处理的一阶FIR高通滤波器系数 noFrames=fix((l-n)/FrameStep)+1; %帧数 FMatrix=zeros(noFrames-2, num); %倒谱系数初始矩阵，一阶差分的首尾2帧为0, lifter=1:num; %倒谱加权初始矩阵 lifter=1+floor(num/2)*(sin(lifter*pi/num));%倒谱加权函数 if mean(abs(s)) > 0.01 s=s/max(s); %初始化 end %计算 MFCC 系数 for i=1:noFrames frame=s((i-1)*FrameStep+1:(i-1)*FrameStep+n); %frame用于存储每一帧数据 framef=filter(b,a,frame); %预加重滤波器 F=framef.*hamming(n); %加汉明窗 FFTo=fft(F,n); %计算FFT melf=melbankm(M,N,Fs); %计算mel滤波器组系数 halfn=1+floor(N/2); spectr=log(melf*(abs(FFTo(1:halfn)).^2)+1e-22);%三角滤波器进行滤波加上1e-22防止进行对数运算后变成无穷 %进行DCT变换 c=zeros(1,num); for p=1:num for m=1:M c(p)=c(p)+spectr(m)*cos(p*(m-0.5)*pi/M); end end ncoeffs=c.*lifter; %对MFCC参数进行倒谱加权 FMatrix(i, :)=ncoeffs; end %调用deltacoeff函数计算MFCC差分系数 d=deltacoeff(FMatrix); %计算一阶差分系数 d1=deltacoeff(d); %计算二阶差分系数 FMatrix=[FMatrix,d,d1]; %将三组数据作为特征向量
melbankm（）函数---程序清单:
function H=melbankm(M,N,fs,fl,fh) %Mel滤波器组 % 输入参数: M 滤波器组数量 % N FFT长度 % fs 采样频率 % fl -fh 为线性功率谱的有用频带(默认为0-0.5*fs) % %输出参数: H返回滤波器组，每一行为一个滤波器 if nargin<4 fl=0*fs; fh=0.5*fs; end %计算每个滤波器的中心频率 f=zeros(1,M+1); for m=1:M+2 f(m)=floor((N/fs)*mel2freq(freq2mel(fl)... +(m-1)*(freq2mel(fh)-freq2mel(fl))/(M+1))); end %求滤波器组H c=floor(N/2)+1; y=zeros(1,c); H=zeros(M,c); for m=2:M+1 for k=1:c %由于fh最高为fs/2，那么最多需c位就能存储H if f(m-1)<=k&&k<=f(m) y(k)=(k-f(m-1))/(f(m)-f(m-1)); elseif f(m)<=k&&k<=f(m+1) y(k)=(f(m+1)-k)/(f(m+1)-f(m)); else y(k)=0; end end H(m,:)=y; end
freq2mel（）函数----程序清单：
function b=freq2mel(f) %%实现线性频域到mel频域的转换 %输入参数：f为线性频域频率 %输出参数：b为mel频域频率 b=2595*log10(1+f/700);
mel2freq（）函数----程序清单：
function f=mel2freq(b) %%实现mel频域转换到线性频域 %输入参数：b为mel频域内频率 %输出参数：f为线性频域内频率 f=700*(10^(b/2595)-1); deltacoeff（）函数----程序清单： function diff = deltacoeff(x) %计算MFCC差分系数 [nr,nc]=size(x); N=2; diff=zeros(nr,nc); for t=3:nr-2 for n=1:N diff(t,:)=diff(t,:)+n*(x(t+n,:)-x(t-n,:)); end diff(t,:)=diff(t,:)/10; %10=2*(1^2+2^2) end
CMN（）函数----程序清单
function NormMatrix = CMN(Matrix) %归一化处理 [r,c]=size(Matrix); NormMatrix=zeros(r,c); for i=1:c MatMean=mean(Matrix(:,i)); NormMatrix(:,i)=Matrix(:,i)-MatMean; end
DTWScores（）函数----程序清单
function AllScores = DTWScores(rMatrix,N)
%动态时间规整（DTW）寻找最小失真
%输入参数：rMatrix为当前读入语音的MFCC参数矩阵,N为每个模板数量词汇数
%输出参数：AllScores
%初始化DTW判别矩阵
Scores1 = zeros(1,N);
Scores2 = zeros(1,N);
Scores3 = zeros(1,N);
%加载模板数据
s1 = load('Vectors1.mat');
fMatrixall1 = struct2cell(s1);
s2 = load('Vectors2.mat');
fMatrixall2 = struct2cell(s2);
s3 = load('Vectors3.mat');
fMatrixall3 = struct2cell(s3);
%计算DTW
for i = 1:N
fMatrix1 = fMatrixall1{i,1};
fMatrix1 = CMN(fMatrix1);
Scores1(i) = myDTW(fMatrix1,rMatrix);
end
for j = 1:N
fMatrix2 = fMatrixall2{j,1};
fMatrix2 = CMN(fMatrix2);
Scores2(j) = myDTW(fMatrix2,rMatrix);
end
for k= 1:N
fMatrix3 = fMatrixall3{k,1};
fMatrix3 = CMN(fMatrix3);
Scores3(k) = myDTW(fMatrix3,rMatrix);
end
AllScores = [Scores1;Scores2;Scores3];