## Distance Time Warp (DTW) function implementation in mathematica

7

3

The DTW (Distance time warp) is given here http://en.wikipedia.org/wiki/Dynamic_time_warping

There is an implementation of DTW in R as well as MATLAB. However i am unable to find something for DTW in Mathematica.

Can someone help me how to implement dtw in Mathematica ??

The DTW function in R is given here:

function (x, y = NULL, dist.method = "Euclidean", step.pattern = symmetric2,
window.type = "none", keep.internals = FALSE, distance.only = FALSE,
open.end = FALSE, open.begin = FALSE, ...)
lm <- NULL
if (is.null(y)) {
if (!is.matrix(x))
stop("Single argument requires a global cost matrix")
lm <- x
}
else if (is.character(dist.method)) {
x <- as.matrix(x)
y <- as.matrix(y)
lm <- proxy::dist(x, y, method = dist.method)
}
else if (is.function(dist.method)) {
stop("Unimplemented")
}
else {
stop("dist.method should be a character method supported by proxy::dist()")
}
wfun <- .canonicalizeWindowFunction(window.type)
dir <- step.pattern
norm <- attr(dir, "norm")
if (!is.null(list(...)$partial)) { warning("Argument partial' is obsolete. Use open.end' instead") open.end <- TRUE } n <- nrow(lm) m <- ncol(lm) if (open.begin) { if (is.na(norm) || norm != "N") { stop("Open-begin requires step patterns with 'N' normalization (e.g. asymmetric, or R-J types (c)). See papers in citation().") } lm <- rbind(0, lm) np <- n + 1 precm <- matrix(NA, nrow = np, ncol = m) precm[1, ] <- 0 } else { precm <- NULL np <- n } gcm <- globalCostMatrix(lm, step.matrix = dir, window.function = wfun, seed = precm, ...) gcm$N <- n
gcm$M <- m gcm$call <- match.call()
gcm$openEnd <- open.end gcm$openBegin <- open.begin
gcm$windowFunction <- wfun lastcol <- gcm$costMatrix[np, ]
if (is.na(norm)) {
}
else if (norm == "N+M") {
lastcol <- lastcol/(n + (1:m))
}
else if (norm == "N") {
lastcol <- lastcol/n
}
else if (norm == "M") {
lastcol <- lastcol/(1:m)
}
gcmjmin <- m if (open.end) { if (is.na(norm)) { stop("Open-end alignments require normalizable step patterns") } gcmjmin <- which.min(lastcol)
}
gcm$distance <- gcm$costMatrix[np, gcm$jmin] if (is.na(gcm$distance)) {
stop("No warping path exists that is allowed by costraints")
}
if (!is.na(norm)) {
gcm$normalizedDistance <- lastcol[gcm$jmin]
}
else {
gcm$normalizedDistance <- NA } if (!distance.only) { mapping <- backtrack(gcm) gcm <- c(gcm, mapping) } if (open.begin) { gcm$index1 <- gcm$index1[-1] - 1 gcm$index2 <- gcm$index2[-1] lm <- lm[-1, ] gcm$costMatrix <- gcm$costMatrix[-1, ] gcm$directionMatrix <- gcm$directionMatrix[-1, ] } if (!keep.internals) { gcm$costMatrix <- NULL
gcm$directionMatrix <- NULL } else { gcm$localCostMatrix <- lm
if (!is.null(y)) {
gcm$query <- x gcm$reference <- y
}
}
class(gcm) <- "dtw"
return(gcm)


}

There is some here http://forums.wolfram.com/mathgroup/archive/2003/Jul/msg00544.html

– Nasser – 2014-10-25T05:39:44.837

Already saw it. It is not helping me that much. I am dealing with numbers itself and not the string – curiousman – 2014-10-25T05:57:35.720

3I believe WarpingDistance is the function you're looking for which is new in v11. – rfrasier – 2017-04-02T11:53:39.007

6

Dynamic time warping (DTW) is a built-in function by the name of WarpingDistance as of Mathematica 11.

6

Here is a direct (recursive) implementation of the DTW. You need three pieces: a distance function, construction of the DTW matrix, and a function to find the best path through the DTW matrix.

(*distance function*)
dist[s_, t_] := Abs[s - t];
(*boundary conditions*)
dtw[1, 1] = dist[x[], y[]];
dtw[1, j_] := dtw[1, j] = dist[x[], y[[j]]] + dtw[1, j - 1];
dtw[i_, 1] := dtw[i, 1] = dist[x[[i]], y[]] + dtw[i - 1, 1];
(*main recursion*)
dtw[i_, j_] := dtw[i, j] = dist[x[[i]], y[[j]]] +
Min[dtw[i - 1, j - 1], dtw[i - 1, j], dtw[i, j - 1]];
(*finding best path through dtwMatrix*)
pathFind[{i_, j_}] := Module[{nbhd},
nbhd = {{i, Min[j + 1, nY]}, {Min[i + 1, nX], j}, {Min[i + 1, nX],
Min[j + 1, nY]}};
nbhd[[First[Ordering[Map[dtwMat[[#[], #[]]] &, nbhd]]]]]];


To use this, set up inputs x and y (here we demonstrate with sinusoids)

x = Table[Sin[2 Pi 6 t], {t, 0, 3, 0.01}];
y1 = Table[Sin[2 Pi 7 t], {t, 0, 1.5, 0.005}];
y2 = Table[Sin[2 Pi 3 t], {t, 1.5, 3, 0.005}];
y = Flatten[{y1, y2}]; nX = Length[x]; nY = Length[y];


Distance matrix distMat & cumulative distance matrix dtwMat:

distMat = Outer[dist, x, y];
dtwMat = dtwPath = Outer[dtw, Range[nX], Range[nY]];
bestPath = NestWhileList[pathFind, {1, 1}, (#[] < nX) && (#[] < nY) &];
Map[(dtwPath[[#[], #[]]] = None) &, bestPath];


Display distance matrix and highlight best path:

GraphicsRow[{ArrayPlot[distMat], ArrayPlot[dtwPath, Background -> Blue]}] Find the length of the shortest path:

dtw[nX, nY]
66.1387


Of course, it's probably easier to use the built in version, as suggested by C. E.!