DTW（Dynamic Time Warping / 动态时间归整） python实现

[python] view plain copy

1. from math import *  
2. import matplotlib.pyplot as plt  
3. import numpy  
4.   
5. def print_matrix(mat) :  
6.     print '[matrix] width : %d height : %d' % (len(mat[0]), len(mat))  
7.     print '-----------------------------------'  
8.     for i in range(len(mat)) :  
9.         print mat[i]#[v[:2] for v in mat[i]]  
10.   
11. def dist_for_float(p1, p2) :  
12.     dist = 0.0  
13.     elem_type = type(p1)  
14.     if  elem_type == float or elem_type == int :  
15.         dist = float(abs(p1 - p2))  
16.     else :  
17.         sumval = 0.0  
18.         for i in range(len(p1)) :  
19.             sumval += pow(p1[i] - p2[i], 2)  
20.         dist = pow(sumval, 0.5)  
21.     return dist  
22.   
23. def dtw(s1, s2, dist_func) :  
24.     w = len(s1)  
25.     h = len(s2)  
26.       
27.     mat = [([[0, 0, 0, 0] for j in range(w)]) for i in range(h)]  
28.       
29.     #print_matrix(mat)  
30.       
31.     for x in range(w) :  
32.         for y in range(h) :              
33.             dist = dist_func(s1[x], s2[y])  
34.             mat[y][x] = [dist, 0, 0, 0]  
35.               
36.     #print_matrix(mat)  
37.   
38.     elem_0_0 = mat[0][0]  
39.     elem_0_0[1] = elem_0_0[0] * 2  
40.   
41.     for x in range(1, w) :  
42.         mat[0][x][1] = mat[0][x][0] + mat[0][x - 1][1]  
43.         mat[0][x][2] = x - 1  
44.         mat[0][x][3] = 0  
45.   
46.     for y in range(1, h) :  
47.         mat[y][0][1] = mat[y][0][0] + mat[y - 1][0][1]              
48.         mat[y][0][2] = 0  
49.         mat[y][0][3] = y - 1  
50.   
51.     for y in range(1, h) :  
52.         for x in range(1, w) :  
53.             distlist = [mat[y][x - 1][1], mat[y - 1][x][1], 2 * mat[y - 1][x - 1][1]]  
54.             mindist = min(distlist)  
55.             idx = distlist.index(mindist)  
56.             mat[y][x][1] = mat[y][x][0] + mindist  
57.             if idx == 0 :  
58.                 mat[y][x][2] = x - 1  
59.                 mat[y][x][3] = y  
60.             elif idx == 1 :  
61.                 mat[y][x][2] = x  
62.                 mat[y][x][3] = y - 1  
63.             else :  
64.                 mat[y][x][2] = x - 1  
65.                 mat[y][x][3] = y - 1  
66.   
67.     result = mat[h - 1][w - 1]  
68.     retval = result[1]  
69.     path = [(w - 1, h - 1)]  
70.     while True :  
71.         x = result[2]  
72.         y = result[3]  
73.         path.append((x, y))  
74.   
75.         result = mat[y][x]  
76.         if x == 0 and y == 0 :  
77.             break  
78.           
79.     #print_matrix(mat)  
80.   
81.     return retval, sorted(path)  
82.   
83. def display(s1, s2) :  
84.   
85.     val, path = dtw(s1, s2, dist_for_float)  
86.       
87.     w = len(s1)  
88.     h = len(s2)  
89.       
90.     mat = [[1] * w for i in range(h)]  
91.     for node in path :  
92.         x, y = node  
93.         mat[y][x] = 0  
94.   
95.     mat = numpy.array(mat)  
96.       
97.     plt.subplot(2, 2, 2)  
98.     c = plt.pcolor(mat, edgecolors='k', linewidths=4)  
99.     plt.title('Dynamic Time Warping (%f)' % val)  
100.   
101.     plt.subplot(2, 2, 1)  
102.     plt.plot(s2, range(len(s2)), 'g')  
103.       
104.     plt.subplot(2, 2, 4)  
105.     plt.plot(range(len(s1)), s1, 'r')  
106.           
107.   
108.     plt.show()  
109.       
110.   
111. s1 = [1, 2, 3, 4, 5, 5, 5, 4]  
112. s2 = [3, 4, 5, 5, 5, 4]  
113. s2 = [1, 2, 3, 4, 5, 5]  
114. s2 = [2, 3, 4, 5, 5, 5]  
115. #val, path = dtw(s1, s2, dist_for_float)  
116. display(s1, s2)  

运行结果如下：