dtw.m

Dynamic Time Warping Function in MATLAB

function [Y,dis,B,D,d]=DTW(w1,w2)
% calculate the distance matrix
N=size(w1,2);
M=size(w2,2);
d=zeros(N,M);
for k=1:N
    for l=1:M
        d(k,l)=sqrt(dot(w1(:,k)-w2(:,l),(w1(:,k)-w2(:,l))'));
    end
end
%calculate the optimum path
D=zeros(N,M);
B=zeros(N,M);
D(1,1)=d(1,1);
B(1,1)=1+i*1;
D(1,2)=D(1,1)+d(1,2);
D(1,2)=inf;
D(1,3:M)=inf;
% for l=2:M
% D(1,l)=D(1,l-1)+d(1,l);
% end
B(1,2:M)=1+i*(1:M-1);
D(2,1)=D(1,1)+d(2,1);
D(2,1)=inf;
D(3:N,1)=inf;
% for k=2:N
% D(k,1)=D(k-1,1)+d(k,1);
% end
 B(2:N,1)=(1:N-1)+1*i;
for k=2:N
    for l=2:M
         if (imag(B(k-1,l))==l && real(B(k,l-1))==k)
             [D(k,l),B(k,l)]=min([inf D(k-1,l-1)+d(k,l) inf]);
         elseif imag(B(k-1,l))==l
             [D(k,l),B(k,l)]=min([inf D(k-1,l-1)+d(k,l) D(k,l-1)+d(k,l)]);
         elseif real(B(k,l-1))==k
             [D(k,l),B(k,l)]=min([D(k-1,l)+d(k,l) D(k-1,l-1)+d(k,l) inf]);
         elseif (imag(B(k-1,l))==l && real(B(k,l-1))==k)
             [D(k,l),B(k,l)]=min([inf D(k-1,l-1)+d(k,l) inf]);
         elseif (imag(B(k-1,l))~=l && real(B(k,l-1))~=k)
             [D(k,l),B(k,l)]=min([D(k-1,l)+d(k,l) D(k-1,l-1)+d(k,l) D(k,l-1)+d(k,l)]);
         end
%             [D(k,l),B(k,l)]=min([D(k-1,l)+d(k,l) D(k-1,l-1)+d(k,l) D(k,l-1)+d(k,l)]);
        if B(k,l)==1
            B(k,l)=(k-1)+i*l;
        elseif B(k,l)==2
            B(k,l)=(k-1)+i*(l-1);
        else
            B(k,l)=k+i*(l-1);
        end
    end
end


%back tracking and termination
k=N;
l=M;
X=[];
while((k~=1)||(l~=1))
    t=k;
    k=real(B(k,l));
    l=imag(B(t,l));
      if(l~=imag(B(k,l)))
        X=[w1(:,k) X];
      end
end
if (imag(B(N,M))==M)
    Y=[w1(:,1) X(:,1:end-1) w1(:,N)];
else
    Y=[w1(:,1) X w1(:,N)];
end
dis=D(N,M);