function D = dtw(x,y,dst)
% D = dtw(x,y,dst)
%
% x first sequence
% y second sequence
% dst distance indicator
%     1:    edit distance
%     2:    Euclidean DTW
%     else  Derivative DTW (default)
%
% Returns DTW distance between x and y.

if nargin < 2
    error('Too few input parameters.');
end
if nargin < 3
    dst = 0;
end

switch dst
    case 1
        mydist = @(x,y,i,j) (x(i) ~= y(j))*1.0; % function handle syntax: 
                                                % handle = @(arglist) funktion
                                                % handle = @funktion
    case 2
        mydist = @(x,y,i,j) (x(i) -  y(j))^2;
    otherwise
        mydist = @(x,y,i,j) (dx(x,i) - dx(y,j))^2;
end

m = length(x);
n = length(y);

gamma = zeros(m,n);
gamma(1,1) = mydist(x,y,1,1);
for i = 2:m
    gamma(i,1) = mydist(x,y,i,1) + gamma(i-1,1);
end
for j = 2:n
    gamma(1,j) = mydist(x,y,1,j) + gamma(1,j-1);
end
for i = 2:m
    for j = 2:n
        h = [gamma(i-1,j-1) gamma(i,j-1) gamma(i-1,j)];
        gamma(i,j) = mydist(x,y,i,j) + min(h);
    end
end
D = gamma(m,n);


function delta_x = dx(x,i)

m = length(x);
if i < 1 || i > m
    delta_x = 0;
    return
end


if i == 1
    if i == m
        delta_x = 0;
    else
        delta_x = x(i + 1) - x(i);
    end
else
    if i == m
        delta_x = x(i) - x(i - 1);
    else
        delta_x = 1/2*(x(i) - x(i - 1) + (x(i + 1) - x(i - 1))/2);
    end
end
