动态时间规整-DTW算法

大家好，欢迎来到IT知识分享网。

作者：桂。

时间：2017-05-31  16:17:29

链接：动态时间规整-DTW算法 

前言

动态时间规整（Dynamic Time Warping，DTW）是孤立词识别的早期技术，梳理一下，主要包括：

　　1）孤立词识别操作步骤;

　　2）DTW原理;

内容基本就是两个博文的整合，最后一并给出链接。

一、孤立词识别操作步骤

基本原理：

基本操作是预加重、分帧，端点检测技术又叫有话帧检测（Voice activity detection, VAD）技术。特征提取参考之前的博文。例如:

检测语音→特征提取。

多说一句，倒谱就是将乘性关系变为加性关系：xy→logx+logy，一般的谱分析我们都是采用频谱，或者小波这样与频谱的区别只是不同量度，这些都是解决加性噪声的滤波问题，倒谱是一种为了滤除乘性噪声的谱方法，简单的说就是对功率谱求log，再反傅里叶变换，公式如 。

特征提取之后就是特征的模板匹配，也就是DTW算法。

二、DTW算法思路

　　A-DTW必要性

语音识别的匹配需要解决的一个关键问题是说话人对同一个词的两次发音不可能完全相同，这些差异不仅包括音强的大小、频谱的偏移，更重要的是发音时音节的长短不可能完全相同，而且两次发音的音节往往不存在线性对应关系。

设参考模板有M帧矢量{R(1)，R(2)，…R(m)，…，R(M)}，R(m)为第m帧的语音特征矢量，测试模板有N帧矢量{T(1)，T(2)，…T(n)，…，T(N)}，T(n)是第n帧的语音特征矢量。d(T(i_n),R(i_m))表示T中第i_n帧特征与R中i_m帧特征之间的欧几里得距离。直接匹配是假设测试模板和参考模板长度相等，即i_n=i_m；线性时间规整技术假设说话速度是按不同说话单元的发音长度等比例分布的，即。这两种假设其实都不符合实际语音的发音情况，我们需要一种更加符合实际情况的非线性时间规整技术，也就是DTW算法。三种匹配模式的对比：

 　　B-DTW思路

　　首先还是介绍下DTW的思想：假设现在有一个标准的参考模板R，是一个M维的向量，即R={R（1），R（2），……，R（m），……，R（M）}，每个分量可以是一个数或者是一个更小的向量。现在有一个才测试的模板T，是一个N维向量，即T={T（1），T（2），……，T（n），……，T（N）}同样每个分量可以是一个数或者是一个更小的向量，注意M不一定等于N，但是每个分量的维数应该相同。

     由于M不一定等于N，现在要计算R和T的相似度，就不能用以前的欧式距离等类似的度量方法了。那用什么方法呢？DTW就是为了解决这个问题而产生的。

首先我们应该知道R中的一个分量R（m）和T中的一个分量T（n）的维数是相同的，它们之间可以计算相似度（即距离）。在运用DTW前，我们要首先计算R的每一个分量和T中的每一个分量之间的距离，形成一个M*N的矩阵。（为了方便，行数用将标准模板的维数M，列数为待测模板的维数N）。

然后下面的步骤该怎么计算呢？用个例子来看看。

这个例子中假设标准模板R为字母ABCDEF(6个)，测试模板T为1234(4个)。R和T中各元素之间的距离已经给出。如下：

     既然是模板匹配，所以各分量的先后匹配顺序已经确定了，虽然不是一一对应的。现在题目的目的是要计算出测试模板T和标准模板R之间的距离。因为2个模板的长度不同，所以其对应匹配的关系有很多种，我们需要找出其中距离最短的那条匹配路径。现假设题目满足如下的约束：当从一个方格((i-1,j-1)或者(i-1,j)或者(i,j-1))中到下一个方格(i,j)，如果是横着或者竖着的话其距离为d(i,j)，如果是斜着对角线过来的则是2d(i,j).其约束条件如下图像所示：

     其中g(i,j)表示2个模板都从起始分量逐次匹配，已经到了M中的i分量和T中的j分量，并且匹配到此步是2个模板之间的距离。并且都是在前一次匹配的结果上加d(i,j)或者2d(i,j),然后取最小值。

     所以我们将所有的匹配步骤标注后如下：

     怎么得来的呢？比如说g(1,1)=4, 当然前提都假设是g(0,0)=0,就是说g(1,1)=g(0,0)+2d(1,1)=0+2*2=4.

     g(2,2)=9是一样的道理。首先如果从g(1,2)来算的话是g(2,2)=g(1,2)+d(2,2)=5+4=9,因为是竖着上去的。

     如果从g(2,1)来算的话是g(2,2)=g(2,1)+d(2,2)=7+4=11,因为是横着往右走的。

     如果从g(1,1)来算的话，g(2,2)=g(1,1)+2*d(2,2)=4+2*4=12.因为是斜着过去的。

     综上所述，取最小值为9. 所有g(2,2)=9.

     当然在这之前要计算出g(1,1),g(2,1),g(1,2).因此计算g(I,j)也是有一定顺序的。

其基本顺序可以体现在如下：

     计算了第一排，其中每一个红色的箭头表示最小值来源的那个方向。当计算了第二排后的结果如下：

     最后都算完了的结果如下：

     到此为止，我们已经得到了答案，即2个模板直接的距离为26. 我们还可以通过回溯找到最短距离的路径，通过箭头方向反推回去。如下所示：

     到这里，估计大家动手算一下就会明白了，其实这个就是动态规划的思路：

主要code:

from numpy import array, zeros, argmin, inf, equal, ndim
from scipy.spatial.distance import cdist

def dtw(x, y, dist):
    """
    Computes Dynamic Time Warping (DTW) of two sequences.
    :param array x: N1*M array
    :param array y: N2*M array
    :param func dist: distance used as cost measure
    Returns the minimum distance, the cost matrix, the accumulated cost matrix, and the wrap path.
    """
    assert len(x)
    assert len(y)
    r, c = len(x), len(y)
    D0 = zeros((r + 1, c + 1))
    D0[0, 1:] = inf
    D0[1:, 0] = inf
    D1 = D0[1:, 1:] # view
    for i in range(r):
        for j in range(c):
            D1[i, j] = dist(x[i], y[j])
    C = D1.copy()
    for i in range(r):
        for j in range(c):
            D1[i, j] += min(D0[i, j], D0[i, j+1], D0[i+1, j])
    if len(x)==1:
        path = zeros(len(y)), range(len(y))
    elif len(y) == 1:
        path = range(len(x)), zeros(len(x))
    else:
        path = _traceback(D0)
    return D1[-1, -1] / sum(D1.shape), C, D1, path


def _traceback(D):
    i, j = array(D.shape) - 2
    p, q = [i], [j]
    while ((i > 0) or (j > 0)):
        tb = argmin((D[i, j], D[i, j+1], D[i+1, j]))
        if (tb == 0):
            i -= 1
            j -= 1
        elif (tb == 1):
            i -= 1
        else: # (tb == 2):
            j -= 1
        p.insert(0, i)
        q.insert(0, j)
    return array(p), array(q)

放在音频里就是：模板/测试样本的STFT取幅度后，每一帧存在一个sum(dist)的距离d_ij,i为模板的第i帧，j表示测试样本的第j帧，这样就得出了路径，从而计算最短路径，完成识别：

三、其他

其实DTW是模板匹配的一般思路，这是理论层面。孤立词识别是它的一个应用，对于尺寸不同的匹配问题，DTW都应该有一定效果，比如数字识别.首先建立模板：

然后是测试样本：

测试样本有大有小，直接模板匹配是不合适的。实际应用中数字可能形状存在差异、且有倾斜，需要把字符切割出来，而不考虑白框。这里只是梳理思路，仅考虑最理想的情况。

例如利用DTW：将图片二值化im2bw→拉成向量img(:)→将测试数据利用DTW，分别与模板匹配，找出最短距离的最小值，即是对应的数字。

参考

• http://www.cnblogs.com/ChengQH/p/2dc8272d6b045b9cee3a02d221662251.html
• http://www.cnblogs.com/tornadomeet/archive/2012/03/23/2413363.html