【语音信号处理四】DTW算法
【参考资料】
【1】《基于DTW距离的时序相似性方法提取水稻遥感信息》
【2】《动态时间规整—DTW算法》 https://blog.csdn.net/qq_39516859/article/details/81705010
【3】《DTW算法Python实现》 https://www.cnblogs.com/ningjing213/p/10502519.html
【4】《动态规划讲解+例子》https://wenku.baidu.com/view/29ffed3e974bcf84b9d528ea81c758f5f61f2989.html#
1. 算法描述
DTW主要用户评价不同长度的时序序列相似度。
1.1 假设存在
长度为n的序列
S
=
{
s
1
,
s
2
,
…
s
n
}
S = \{s_1, s_2, \dots s_n \}
S={s1,s2,…sn}
长度为m的序列
T
=
{
t
1
,
t
2
,
…
t
m
}
T = \{t_1, t_2, \dots t_m \}
T={t1,t2,…tm}
其中
n
≠
m
n \ne m
n=m
1.2 构造距离矩阵
A n × m A_{n \times m} An×m，矩阵A的每个元素 a i j = d ( s i , t j ) a_{ij} = d(s_i, t_j) aij=d(si,tj)代表两个序列对应点的距离。这个距离计算的方法通常就是欧式距离。
构造这个矩阵的目的实际上是未来寻求两个不长度曲线的最佳匹配。如下图所示：
如果时间长度序列一样这一次匹配各点欧式距离即可得到相似度。但如果时序长度不同
这需要存在一种匹配关系，比如在上图的例子中红色序列的第一个点同时匹配绿色序列的二三两个点。
备注：这张图由于绘制的时候由于时间间隔画比较大所以不那么直观。
那么这种匹配关系的查找在算法中可以体现为在矩阵A中寻找一条路径从 a 1 1 a_11 a11移动到 a n m a_nm anm，然后根据这个匹配原则计算距离，使得距离和最小【最佳匹配】。
下图中的红、黑、蓝轨迹实际就代表三种匹配规则。
1.3 算法逻辑
DTW算法的本质是用动态规划的思路去搜索这个最佳匹配。从这个角度讲应该用什么算法是无所谓的，只要能搜索到这个最佳匹配 ????
动态规划算法备注：
动态规划是每一轮都基于上一轮的最优解，dijkstra的最短路径算法应该是最经典的，需要注意的是，如上图所示。
在计算B轮的时候，并不因为 { C 1 , D 1 , E } \{C_1,D_1,E\} {C1,D1,E}是最有解，就只搜索 { B 1 , C 1 } \{B_1,C_1\} {B1,C1}、 { B 2 , C 1 } \{B_2,C_1\} {B2,C1}、 { B 3 , C 1 } \{B_3,C_1\} {B3,C1} 而是搜索全部从B到C的9条可选路径，只是对于上一轮,如果从 C 1 C_1 C1出发这取上一轮 C 1 C_1 C1的最优解，同理 C 2 C_2 C2、 C 3 C_3 C3
DTW算法备注：
2. python代码实现
s = [1,2,3,4,5]
t = [1,1,2,3,4,5]
#构筑矩阵
# [[0, 1, 4, 9, 16], [0, 1, 4, 9, 16], [1, 0, 1, 4, 9],
# [4, 1, 0, 1, 4], [9, 4, 1, 0, 1], [16, 9, 4, 1, 0]]
def _get_min_in_three(x,y,z):
min = x if x < y else y
min = min if min < z else z
return min
def _do_dtw(s, t):
#1. 构筑距离矩阵A
A = [[0 for i in range(len(s))] for j in range(len(t))]
for i in range(len(s)):
for j in range(len(t)):
A[j][i] = abs(s[i] - t[j])
#2. 采用动态规划所搜最短路径
#移动的路径是从矩阵的左上角移动到右下角
D = [[0 for i in range(len(s))] for j in range(len(t))]
for i in range(len(s)):
for j in range(len(t)):
if 0 == i and 0 == j:
D[j][i] = 0
elif j == 0: # 贴列
D[j][i] = A[j][i] + D[j][i - 1]
elif i == 0: #贴行
D[j][i] = A[j][i] + D[j - 1][i]
else:
# D[j-1,i-1] D[j-1, i]
# D[j,i-1] D[j, i ]
# 这里 D[j, i ] 为 D[j-1][i-1] D[j-1[i] D[j[i-1] 中最小值 + A[j][i]
D[j][i] = _get_min_in_three(D[j-1][i-1],D[j-1][i],D[j][i-1]) + A[j][i]
return A, D
A, D = _do_dtw(s, t)