计算两个时间序列的soft dtw的matlab代码

以下是两个时间序列的soft dtw的matlab代码：

function [D, gamma] = softdtw(x, y, gamma)
% SOFTDTW Dynamic Time Warping with Soft-DTW loss
%
% [D, gamma] = softdtw(x, y, gamma)
%
% Input:
% x: NxD array, a time series (D dimensions)
% y: MxD array, another time series (D dimensions)
% gamma: regularization parameter for Soft-DTW (gamma >= 0)
% Output:
% D: Soft-DTW distance
% gamma: optimal gamma
%
% Soft-DTW was originally proposed in the following paper:
% M. Cuturi, "Fast Global Alignment Kernels," ICML'11.
% https://arxiv.org/abs/1105.5716

if nargin < 3
gamma = 1.0;
end

% initialization
N = size(x, 1);
M = size(y, 1);
D = zeros(N+2, M+2) + inf;
D(1,1) = 0;

% compute distances
for i=1:N
for j=1:M
cost = sum((x(i,:)-y(j,:)).^2);
D(i+1,j+1) = cost + min([D(i+1,j), D(i,j+1), D(i,j)]);
end
end

% compute optimal gamma
if gamma > 0

% initialization
log_r = (-log(N+M-2)-psi(gamma))/gamma;
log_u = (-log(N)-psi(gamma))/gamma;
log_v = (-log(M)-psi(gamma))/gamma;
log_K = -inf(N+1,M+1);
log_K(1,1) = log(1);

% forward recursion
for i=1:N
for j=1:M
log_K(i+1,j+1) = logsumexp([
log_K(i,j+1)+log_r+(D(i+1,j+1)-D(i+1,j)-D(i,j+1))/gamma,
log_K(i,j)+log_u+(D(i+1,j+1)-D(i,j)-D(i+1,j))/gamma,
log_K(i+1,j)+log_v+(D(i+1,j+1)-D(i,j+1)-D(i+1,j))/gamma]);
end
end

% backward recursion
path = [N,M];
i = N;
j = M;
while i > 1 || j > 1
c = [D(i,j),D(i+1,j),D(i,j+1)];
if i == 1
m = 2;
elseif j == 1
m = 3;
else
[~,m] = min(c);
end
if m == 1 % diagonal
path = [i-1,j-1; path];
i = i - 1;
j = j - 1;
elseif m == 2 % up
path = [i-1,j; path];
i = i - 1;
else % left
path = [i,j-1; path];
j = j - 1;
end
end

gamma = (D(N+1,M+1)-sum(exp(log_K(N+1,M+1))))/(N+M-2);

end

D = sqrt(D(N+1,M+1));

function s = logsumexp(a, dim)
% Computes log(sum(exp(a),dim)) while avoiding numerical underflow.
% Default is dim = 1 (columns).

% compute along the first dimension by default
if nargin < 2
dim = 1;
end

% subtract the largest in each column
[y, i] = max(a,[],dim);
dims = ones(1,ndims(a));
dims(dim) = size(a,dim);
a = a - repmat(y, dims);

% compute log(sum(exp(a),dim)) while avoiding numerical overflow
s = y + log(sum(exp(a),dim));
i = find(~isfinite(y));
if ~isempty(i)
s(i) = y(i);
end

% replace NaN with -Inf
s(isnan(s)) = -Inf;

end

function y = psi(x)
% Digamma function of x (psi(x) = d/dx log(gamma(x)))
% https://en.wikipedia.org/wiki/Digamma_function
%
% This implementation has been adapted from the GPy package.
% https://github.com/SheffieldML/GPy/tree/devel/GPy/util

y = zeros(size(x));
for i=1:numel(x)
if x(i) <= 0
y(i) = NaN;
elseif x(i) <= 6
y(i) = psi_0_6(x(i));
else
y(i) = psi(x(i) - 1) + 1/(x(i) - 1);
end
end

function y = psi_0_6(x)
% Digamma function of x for 0 < x <= 6
% Coefficients were obtained using the Remez algorithm

p = [0.0083297439961316088,...
-0.040665510375325418,...
0.042909613820387576,...
-0.074550332888769239,...
0.12164618014696546,...
-0.50000000002530351];
y = polyval(p, x - 3)/((x - 1)*(x - 2)) - 1/(2*x - 1) - 1/(2*x);

end

end