7
$\begingroup$

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)
}
gcm$jmin <- m
    if (open.end) {
        if (is.na(norm)) {
            stop("Open-end alignments require normalizable step patterns")
        }
        gcm$jmin <- 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)

}

$\endgroup$
3
  • $\begingroup$ There is some here forums.wolfram.com/mathgroup/archive/2003/Jul/msg00544.html $\endgroup$
    – Nasser
    Oct 25, 2014 at 5:39
  • $\begingroup$ Already saw it. It is not helping me that much. I am dealing with numbers itself and not the string $\endgroup$
    – curiousman
    Oct 25, 2014 at 5:57
  • 3
    $\begingroup$ I believe WarpingDistance is the function you're looking for which is new in v11. $\endgroup$
    – rfrasier
    Apr 2, 2017 at 11:53

2 Answers 2

6
$\begingroup$

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

$\endgroup$
6
$\begingroup$

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[[1]], y[[1]]];
dtw[1, j_] := dtw[1, j] = dist[x[[1]], y[[j]]] + dtw[1, j - 1];
dtw[i_, 1] := dtw[i, 1] = dist[x[[i]], y[[1]]] + 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[[#[[1]], #[[2]]]] &, 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}, (#[[1]] < nX) && (#[[2]] < nY) &];
Map[(dtwPath[[#[[1]], #[[2]]]] = None) &, bestPath];

Display distance matrix and highlight best path:

GraphicsRow[{ArrayPlot[distMat], ArrayPlot[dtwPath, Background -> Blue]}]

enter image description here

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.!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.