动态时间规整DTW
好吧，最近耳机坏了，日剧计划暂停一下，好好更新一下学术方面的东西吧。。。
今天介绍的是动态时间规整算法（Dynamic Time Warping），我对这个算法有着非常深的感情，嗯，你知道为什么吗？今年上半年还是大四的我，显然是到了要做毕设的时候了，然后指导老师在2月20才交代给我要做什么，然后我就盲目的乱弄了一个多月，因为是算法设计的类型问题，想不出算法你就别指望可以写得出什么东西来，比起那些只有工作量没有太多思维量的毕设题目，我对我的前景堪忧啊，多少个晚上躺在床上，突然突发奇想想到一个东西，然后就下床开电脑验证一下，再然后，就是失望地回去睡觉了。。直到那一天，我遇到了你，DTW，是你拯救了我！！让我产生了一个神奇的灵感，然后花没几天就把毕设给做完了。。。不仅出了不错的结果，还拿了个校级优秀论文。。
我。。。好像又TM扯远了。。。
好吧，其实最开始DTW适用于语音识别的，假设有这样一个问题，我们有一个数据库，里面存了几个简单的词汇的语音数据，语音数据嘛，就是一个一维向量，也不一定是几个，有很多个也行，然后我们现在采集到一个新的语音数据了，然后问我们可不可以通过和数据库里面的数据进行匹配，然后分别出这个新的语音说的是什么？？
所谓的匹配就是做差，差越小说明越相似。如下图所示，首先我们有一个从数据库里面拿出来的参考数据R，一个是刚得到的测试数据T，第二排第一个图就是所谓的直接做差！灰色区域就是差的绝对值。
发现问题没？两个模板居然不一样长！！！但是你也必须承认，这是肯定的嘛！！我们说话假设采样率是3k，你两次说同样一个词，得到的数据长度你能保证一样吗？？不可能的嘛。更何况不同人说的？？对吧。
那要怎么办呢？？好，不一样长我们就把它们用采样插值之类的方法变到一样长，说白了就是缩放。上图的D2所示。
好吧，又发现一个问题，假如我们说的单词是“你好”两个字，A可能“你”这个字说的慢一点，“好”这个字说的快一点，但是B却正好相反，这要怎么办呢？？嗯，没错，这样就引进了第三种做差方式，非线性匹配，不同段采用不用的“缩放率”。
你觉得很难想象，很难实现？？且听我慢慢废话。。。
首先我们把问题重述一遍，我们现在有两个以为数据，R和T，一个长度为M，一个长度为N，现在我们要找到一个匹配方式，非线性的，使得他们的差最小。
既然有差这个概念，首先就要有距离这个概念，什么叫做两个东西的差，用空间的角度来看，就是两个的距离。所以我们先定义R中的第i个数据与T中的第j个数据之间的距离用d(Ri, Tj)表示。这语音这个实例中，其实就是R的第i个点的值减去T的第j个点的值，再做绝对值。
先看下图：
我们现在把两个模板分别定义为x轴和y轴，显然一个长度是M，一个长度是N，然后画出网格子来,像上面那样。然后我们定义网格上的每一个点就是R和T对应点的距离，也就是说点(i,j)的值是d(Ri, Tj)。非线性匹配是什么意思呢？就是以左下角点(1,1)为起始点【非要说(0,0)的这辈子都是码农的命！！】，右上角(M,N)为终点，寻找一条路径，是的路径上经过的点的值的和最小！！
如果要用数学来表达的话，假设我们的匹配方法是 Im=Φ(In)，也就是说路径上经过(n,m)这个点，那么我们的目标就是：
D = min∑d(Ri,TΦ(i))
我觉得嘛，有时候看公式会对东西理解比较清晰，但是有时候非公式的言语说的更清楚。
再补充一下说明吧，非线性映射就是照上面一条路径，而我们一开始讨论的直接匹配找到的路劲就是y=x找条线，而线性匹配呢，找到的是y= Nx/M这条线。我们要找的，是一条弯曲的线！！
在找这条线之前，我们要先约定一下下面几个条件，不然出来的结果就不符合常理了：
- Φ(1) = 1；言下之意就是起始点必须是左下角那个点！！
- Φ(N) = M；言下之意就是重点必须是右上角那个点！！
- Φ(i+1) ≥ Φ(i)，就是说路径必须是一直往上走的，不可以往下走，往下走什么意思呢，就是i和j点匹配，i+1点却和j-1点匹配，这不符合语音的常理嘛。
- Φ(i+1)-Φ(i) ≤ 1；这个说明我们必须考虑所有的数据，不可以跳过任何一个点。
有了上述约定，我们找到的路径就是我们最想要的路径了。假设我们找到了这条路径，使得一路上的累加的距离和min∑d(Ri,TΦ(i))最小，我们称这个距离为DTW距离。对于我们一开始提出来的语音匹配问题，我们只要计算测试模板与所有数据库中的模板的DTW距离，最小那个，就是和它相匹配的那个结果。这就是最简单的语音识别算法。
问题来了，这条路径怎么走呢？？搞过ACM的人一眼就可以看出来要怎么弄，稍微了解一点算法，加上知道DTW的中文名是动态时间规整，也应该想到了，没错，就是动态规划！！锵锵锵锵！！！
好吧，我错了，我发现上面4个约定不够，再加一个约定吧。。。囧。
约定5：我们假设路径上可以到达(i,j)这个点的前一个点，只能是(i-1,j),(i-1,j-1),(i,j-1)三个点中的一个。其实为什么把这个约定和上面的约定专门分开出来说呢？
因为上面4个约定是不管什么DTW算法都必须要遵守的，但是现在我们的约定5呢，是可以改变的，也就是说根据实际应用，不一定要按照我说的做，比如说，可以改成：
约定5版本二：我们假设路径上可以到达(i,j)这个点的前一个点，只能是(i-1,j),(i-1,j-1),(i-1,j-2)三个点中的一个。懂了么？
现在我们可以用动态规划来找路径了。
为了方便说明，现在我们定义一个大写的字母D的距离，D(i,j)表示点(i,j)开始到点(M,N)这个点的最小累积距离，懂什么意思么？就是说假设我们只匹配右上角的某一小块，这一小块左下角那个点是原图的(i,j)那个点，右上角还是(N,M)这个点，这一小块的DTW距离成为D(i,j)。
假设我们要计算D(1,1)，所谓D(1,1)是什么，就是我们要求的DTW距离，要计算这个，我们只要计算下面三个值：
D(1,2),D(2,2),D(2,1),如果可以计算出这三个来，那么:
D(1,1) = d(1,1) + min{D(1,2),D(2,2),D(2,1)}
懂什么意思了么？点(1,1)的下一个点只能是(1,2),(2,2),(2,1)三个中的一个，那么这三个点对应的D中最小的那个，加上d(1,1)本身的值，那么就是(1,1)这个点到终点所要累加的距离了。
你说D(1,2),D(2,2),D(2,1)这三个是多少你不知道？？
这篇破文章你都看到这里了，难道你就提炼不出下面这个公式？
D(i,j) = d(i,j) + min{D(i+1,j),D(i+1,j+1),D(i,j+1)}
懂了吗。。。了吗。。。吗。。。
就是个递归吗，先算出D(N,M) = d(N,M)
然后从右上往左下慢慢推，最后推到D(1,1)来。。
唉，好吧，写的我好累。。
下面给个matlab代码吧。。
code:
function [dist path]= DTW(t,r,test_model)
if(nargin == 2)
test_model = 0;
end
if(size(t,1) == 1)
t = t';
end
if(size(r,1) == 1)
r = r';
end
n = size(t,1);
m = size(r,1);
% 帧匹配距离矩阵
d = zeros(n,m);
for i = 1:n
for j = 1:m
d(i,j) = sqrt(sum((t(i,:)-r(j,:)).^2));
end
end
% 累积距离矩阵
D = ones(n,m) * realmax;
D(1,1) = d(1,1);
map = zeros(n,m,2);
map(1,1,:) = [1 1];
% 动态规划
for i = 1:n
for j = 1:m
if(i == 1 && j == 1)
continue;
end
if(i > 1)
D1 = D(i-1,j);
else
D1 = realmax;
end
if(j > 1 && i > 1)
D2 = D(i-1,j-1);
else
D2 = realmax;
end
if(j > 1)
D3 = D(i,j-1);
else
D3 = realmax;
end
if(test_model==1)
index = find([D1,D2,D3] == min([D1,D2,D3]));
index = index(1);
if(index == 1)
map(i,j,1) = i-1;
map(i,j,2) = j;
end
if(index == 2)
map(i,j,1) = i-1;
map(i,j,2) = j-1;
end
if(index == 3)
map(i,j,1) = i;
map(i,j,2) = j-1;
end
end
D(i,j) = d(i,j) + min([D1,D2,D3]);
end
end
dist = D(n,m);
% map
path = [n m];
if(test_model == 1)
while(path(1,1) >1 || path(1,2) > 1)
path = [map(path(1,1),path(1,2),1) map(path(1,1),path(1,2),2); path];
end
end
【完】
本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com
求详解！！！！