首页 > 模式识别 > 动态时间规整DTW

动态时间规整DTW

2012年11月21日 发表评论 阅读评论

好吧,最近耳机坏了,日剧计划暂停一下,好好更新一下学术方面的东西吧。。。

今天介绍的是动态时间规整算法(Dynamic Time Warping),我对这个算法有着非常深的感情,嗯,你知道为什么吗?今年上半年还是大四的我,显然是到了要做毕设的时候了,然后指导老师在2月20才交代给我要做什么,然后我就盲目的乱弄了一个多月,因为是算法设计的类型问题,想不出算法你就别指望可以写得出什么东西来,比起那些只有工作量没有太多思维量的毕设题目,我对我的前景堪忧啊,多少个晚上躺在床上,突然突发奇想想到一个东西,然后就下床开电脑验证一下,再然后,就是失望地回去睡觉了。。直到那一天,我遇到了你,DTW,是你拯救了我!!让我产生了一个神奇的灵感,然后花没几天就把毕设给做完了。。。不仅出了不错的结果,还拿了个校级优秀论文。。

我。。。好像又TM扯远了。。。

好吧,其实最开始DTW适用于语音识别的,假设有这样一个问题,我们有一个数据库,里面存了几个简单的词汇的语音数据,语音数据嘛,就是一个一维向量,也不一定是几个,有很多个也行,然后我们现在采集到一个新的语音数据了,然后问我们可不可以通过和数据库里面的数据进行匹配,然后分别出这个新的语音说的是什么??

所谓的匹配就是做差,差越小说明越相似。如下图所示,首先我们有一个从数据库里面拿出来的参考数据R,一个是刚得到的测试数据T,第二排第一个图就是所谓的直接做差!灰色区域就是差的绝对值。
DTW1发现问题没?两个模板居然不一样长!!!但是你也必须承认,这是肯定的嘛!!我们说话假设采样率是3k,你两次说同样一个词,得到的数据长度你能保证一样吗??不可能的嘛。更何况不同人说的??对吧。

那要怎么办呢??好,不一样长我们就把它们用采样插值之类的方法变到一样长,说白了就是缩放。上图的D2所示。

好吧,又发现一个问题,假如我们说的单词是“你好”两个字,A可能“你”这个字说的慢一点,“好”这个字说的快一点,但是B却正好相反,这要怎么办呢??嗯,没错,这样就引进了第三种做差方式,非线性匹配,不同段采用不用的“缩放率”。

你觉得很难想象,很难实现??且听我慢慢废话。。。

首先我们把问题重述一遍,我们现在有两个以为数据,R和T,一个长度为M,一个长度为N,现在我们要找到一个匹配方式,非线性的,使得他们的差最小。

既然有差这个概念,首先就要有距离这个概念,什么叫做两个东西的差,用空间的角度来看,就是两个的距离。所以我们先定义R中的第i个数据与T中的第j个数据之间的距离用d(Ri, Tj)表示。这语音这个实例中,其实就是R的第i个点的值减去T的第j个点的值,再做绝对值。

先看下图:
DTW2
我们现在把两个模板分别定义为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) = 1;言下之意就是起始点必须是左下角那个点!!
  2. Φ(N) = M;言下之意就是重点必须是右上角那个点!!
  3. Φ(i+1) ≥ Φ(i),就是说路径必须是一直往上走的,不可以往下走,往下走什么意思呢,就是i和j点匹配,i+1点却和j-1点匹配,这不符合语音的常理嘛。
  4. Φ(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

  • 澈子

    求详解!!!!