DDTW（Derivative Dynamic Time Warping）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import numpy as np

def derivative_series(x: np.ndarray) -> np.ndarray:
    """中央差分で一次微分を近似する。"""
    x = np.asarray(x, dtype=float)
    if x.size < 3:
        raise ValueError("系列の長さは 3 以上にしてください。")

    d = np.empty_like(x)
    d[0] = ((x[1] - x[0]) + (x[2] - x[0]) / 2.0) / 2.0
    d[-1] = ((x[-1] - x[-2]) + (x[-1] - x[-3]) / 2.0) / 2.0
    for i in range(1, len(x) - 1):
        d[i] = ((x[i] - x[i - 1]) + (x[i + 1] - x[i - 1]) / 2.0) / 2.0
    return d

def ddtw_distance(
    x: np.ndarray,
    y: np.ndarray,
    window: int | None = None,
) -> tuple[float, list[tuple[int, int]]]:
    dx, dy = derivative_series(x), derivative_series(y)
    n, m = len(dx), len(dy)

    if window is None:
        window = max(n, m)
    window = max(window, abs(n - m))

    cost = np.full((n + 1, m + 1), np.inf)
    cost[0, 0] = 0.0

    for i in range(1, n + 1):
        for j in range(max(1, i - window), min(m + 1, i + window + 1)):
            dist = abs(dx[i - 1] - dy[j - 1])
            cost[i, j] = dist + min(cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1])

    i, j = n, m
    path: list[tuple[int, int]] = []
    while i > 0 and j > 0:
        path.append((i - 1, j - 1))
        step = np.argmin((cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1]))
        if step == 0:
            i -= 1
        elif step == 1:
            j -= 1
        else:
            i -= 1
            j -= 1
    path.reverse()
    return cost[n, m], path

1
2

python --version  # 例: Python 3.13.0
pip install matplotlib

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import matplotlib.pyplot as plt

np.random.seed(0)
t = np.linspace(0, 4 * np.pi, 120)

series_a = np.sin(t)
series_b = np.sin(t) + 0.6
series_c = np.sin(1.2 * t + 0.3)

plt.figure(figsize=(10, 3))
plt.plot(t, series_a, label="Series A", color="black")
plt.plot(t, series_b, label="Series B (offset)", color="tab:green")
plt.plot(t, series_c, label="Series C (shape change)", color="tab:orange")
plt.legend()
plt.title("比較する 3 つの時系列")
plt.tight_layout()
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

def dtw_distance(x: np.ndarray, y: np.ndarray) -> float:
    n, m = len(x), len(y)
    cost = np.full((n + 1, m + 1), np.inf)
    cost[0, 0] = 0.0
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            dist = abs(x[i - 1] - y[j - 1])
            cost[i, j] = dist + min(cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1])
    return cost[n, m]

dtw_ab = dtw_distance(series_a, series_b)
dtw_ac = dtw_distance(series_a, series_c)
ddtw_ab, _ = ddtw_distance(series_a, series_b)
ddtw_ac, _ = ddtw_distance(series_a, series_c)

print(f"DTW (A, B):  {dtw_ab:.3f}")
print(f"DTW (A, C):  {dtw_ac:.3f}")
print(f"DDTW(A, B):  {ddtw_ab:.3f}")
print(f"DDTW(A, C):  {ddtw_ac:.3f}")

1
2
3
4

DTW (A, B):  72.248
DTW (A, C):  51.158
DDTW(A, B):   0.002
DDTW(A, C):  13.632

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

def ddtw_cost_matrix(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, list[tuple[int, int]]]:
    dx, dy = derivative_series(x), derivative_series(y)
    n, m = len(dx), len(dy)
    cost = np.full((n + 1, m + 1), np.inf)
    cost[0, 0] = 0.0

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            dist = abs(dx[i - 1] - dy[j - 1])
            cost[i, j] = dist + min(cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1])

    i, j = n, m
    path: list[tuple[int, int]] = []
    while i > 0 and j > 0:
        path.append((i - 1, j - 1))
        step = np.argmin((cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1]))
        if step == 0:
            i -= 1
        elif step == 1:
            j -= 1
        else:
            i -= 1
            j -= 1
    path.reverse()
    return cost[1:, 1:], path

cost_ac, path_ac = ddtw_cost_matrix(series_a, series_c)

plt.figure(figsize=(5, 4))
plt.imshow(cost_ac, origin="lower", cmap="viridis")
path_arr = np.array(path_ac)
plt.plot(path_arr[:, 1], path_arr[:, 0], color="white", linewidth=1.5)
plt.title("DDTW cumulative cost (Series A vs Series C)")
plt.xlabel("Series index j")
plt.ylabel("Series index i")
plt.colorbar(shrink=0.8, label="Cumulative cost")
plt.tight_layout()
plt.show()