完整的dtw算法

来源：互联网发布：调查报告数据分析模板编辑：程序博客网时间：2024/05/16 16:59

首先申明：这个完整的dtw算法是我自己收集并且手动敲上去的，算法不是我写的！只有读取动态语音的部分是我自己写的，我也是为了方便做语音识别的人！！

首先mfcc.m如下：

我就不用m文件了

function ccc=mfcc(x)
%归一化mel滤波器数组系数
bank =melbankm(24,256,8000,0,0.5,'m');

bank=full(bank);
%dct系数，12*24
for k=1:12
n=0:23;
dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24));

end
%归一化倒谱提升窗口
w=1+6*sin(pi*[1:12]./12);
w=w/max(w);
%预加重滤波器
xx=double(x);
xx=filter([1 -0.9375],1,xx);
%语音信号分帧
xx=enframe(xx,256,80);
%计算每帧的mfcc系数
for i=1:size(xx,1);
y=xx(1,:);
s=y'.*hamming(256);
t=abs(fft(s));
t=t.^2;
c1=dctcoef*log(bank*t(1:129));
c2=c1.*w';
m(i,:)=c2';

end
%差分参数
dtm=zeros(size(m));

for i=3:size(m,1)-2
dtm(1,:)=-2*m(i-2,:)-m(i-1,:)+m(i+1,:)+2*m(i+2,:);

end
dtm=dtm/3;
%合并mfcc参数和差分mfcc参数
ccc=[m dtm];
%去除首尾两帧，因为这两帧的一阶差分参数为0
disp('mfcc 参数‘');
ccc=ccc(1:size(m,1),:);

其次：dtw.m

function dist =dtw(t,r)
n=size(t,1);
m=size(r,1);
%帧匹配距离矩阵
d=zeros(n,m);
for i=1:n
for j=1:m
d(i,j)=sum((t(i,:)-r(j,:)).^2);

end
end
%累积距离矩阵
D=ones(n,m)*realmax;
D(1,1)=d(1,1);
%动态规划
for i=2:n
for j=1:m
D1=D(i-1,j);
if j>1
D2=D(i-1,j-1);
else
D2=realmax;
end
if j>2
D3=D(i-1,j-2);
else
D3=realmax;
end
D(i,j)=d(i,j)+min([D1,D2,D3]);

end
end
dist=D(n,m);
end

再次：testdtw.mdisp('正在计算参考模板的参数。。。');
filepath1='H:\2\990881\new\';%注册文件位置
filepath2='H:\2\990881\new1\';
for i=1:10
%fname=sprintf('H:\2\9902130\%d.wav',i-1);
str1=filepath1;
str2=num2str(i-1);
str3='.wav';
str=strcat(str1,str2,str3);
[x,fs]=audioread(str);%读取注册文件
%x=audioread(fname);
[x1,x2]=vad(x);
m=mfcc(x);

m = m(x1:x2-2,:);

ref(i).mfcc=m
end
disp('正在计算测试模板的参数');
for i=1:10
%fname=sprintf('H:\2\9902130\%d.wav',i-1);
%x=audioread(fname);
str1=filepath2;
str2=num2str(i-1);
str3='.wav';
str=strcat(str1,str2,str3);

[x,fs]=audioread(str);%读取注册文件
[x1,x2]=vad(x);
m=mfcc(x);
m=m(x1:x2-2,:);
test(i).mfcc=m;
end
disp('正在进行末班匹配');
dist=zeros(10,10);
for i=1:10
for j=1:10
dist(i,j)=dtw(test(i).mfcc,ref(j).mfcc);

end
end
disp('正在计算匹配结果。。。。');
for i=1:10
[d,j]=min(dist(i,:));
fprintf('测试模板%d的识别结果为;%d\n',i,j);

end
vad.m

function [x1,x2]=vad(x)
%幅度归一化到【-1,1】
x=double(x);
x=x/max(abs(x));
%常数设置
%参数设置
FrameLen = 256; %帧长
FrameInc = 80; %未重叠部分
amp1=10;
amp2=2;
zcr1=10;
zcr2=5;
maxsilence=3;%3*10ms=30ms
minlen=15;%15*10ms=150ms
status=0;
count=0;
silence=0;
%计算过零率
tmp1 = enframe(x(1:end-1), FrameLen, FrameInc);
tmp2 = enframe(x(2:end) , FrameLen,FrameInc);
signs = (tmp1.*tmp2)<0;
diffs = (tmp1 -tmp2)>0.01;
zcr = sum(signs.*diffs,2);
%计算短时能量
amp = sum((abs(enframe(filter([1 -0.9375], 1, x), FrameLen, FrameInc))).^2, 2);
%调整能量门限
amp1 = min(amp1, max(amp)/4);
amp2 = min(amp2, max(amp)/8);
%开始断点检测
for n=1:length(zcr)
goto=0;
switch status
case {0,1}
if amp(n)>amp1
x1=max(n-count-1,1);
status=2;
silence=0;
count=count+1;
elseif amp(n)>amp2||zcr(n)>zcr2%可能处于语音段
zcr(n)>zcr2;
status=1;
count=count+1;
else%处于静音状态
status=0;
count=0;

end
case 2
if amp(n)>amp2||zcr(n)>zcr2%保持在语音段
count=count+1;
else%语音段结束
silence=silence+1;
if silence<maxsilence%静音不够长语音还未结束
count=count+1;
elseif count<minlen%认为是噪音
status=0;
silence=0;
count=0;
else
status=3;

end
end
case 3
break;
end
end
count=count-silence/2;
%x1=round(x1);
x2=x1+count-1;
%disp('x1');
%disp(x1);

end

1 0