【重大修改】动态时间规整（Dynamic Time Warping）

本文只是简单的介绍DTW算法的目的和实现。具体的DTW可以参考一下文献：

离散序列的一致性度量方法：动态时间规整（DTW）  http://blog.csdn.net/liyuefeilong/article/details/45748399

动态时间归整/规整/弯曲(Dynamic time warping,DTW)   http://www.cnblogs.com/flypiggy/p/3603192.html

【更新日志】2017-7-17

评论区一位同学指出文章有误，笔者能力有限，暂未推导出错误之处，这里提供一个一般不会出错的理论介绍网址：维基百科，毕竟能上维基百科的都是大家公认的。DTW相关维基百科戳这里，有兴趣同学可以在评论区一起讨论文章错误之处，谢谢

DTW是干什么的？

    动态时间规整算法，故名思议，就是把两个代表同一个类型的事物的不同长度序列进行时间上的“对齐”。比如DTW最常用的地方，语音识别中，同一个字母，由不同人发音，长短肯定不一样，把声音记录下来以后，它的信号肯定是很相似的，只是在时间上不太对整齐而已。所以我们需要用一个函数拉长或者缩短其中一个信号，使得它们之间的误差达到最小。

     再来看看运动捕捉，比如当前有一个很快的拳击动作，另外有一个未加标签的动作，我想判断它是不是拳击动作，那么就要计算这个未加标签的动作和已知的拳击动作的相似度。但是呢，他们两个的动作长度不一样，比如一个是100帧，一个是200帧，那么这样直接对每一帧进行对比，计算到后面，误差肯定很大，那么我们把已知拳击动作的每一帧找到无标签的动作的对应帧中，使得它们的距离最短。这样便可以计算出两个运动的相似度，然后设定一个阈值，满足阈值范围就把未知动作加上“拳击”标签。

DTW怎么计算？

下面我们来总结一下DTW动态时间规整算法的简单的步骤：

1. 首先肯定是已知两个或者多个序列，但是都是两个两个的比较，所以我们假设有两个序列A={a1,a2,a3,...,am}  B={b1,b2,b3,....,bn}，维度m>n

2. 然后用欧式距离计算出每序列的每两点之间的距离，D(ai,bj) 其中1≤i≤m，1≤j≤n

   画出下表：

                                                          

3.  接下来就是根据上图将最短路径找出来。从D(a1,a2)沿着某条路径到达D(am,bn)。找路径满足：假如当前节点是D(ai,bj)，那么下一个节点必须是在D(i+1,j)，D(i,j+1)，D(i+1,j+1)之间选择，并且路径必须是最短的。计算的时候是按照动态规划的思想计算，也就是说在计算到达第(i,j)个节点的最短路径时候，考虑的是左上角也即第(i-1,j)、(i-1,j-1)、(i,j-1)这三个点到(i,j)的最短距离。

4. 接下来从最终的最短距离往回找到那条最佳的输出路径， 从D(a1,b1)到D(am,bn)。他们的总和就是就是所需要的DTW距离

【注】如果不回溯路径，直接在第3步的时候将左上角三个节点到下一个节点最短的点作为最优路径节点，就是贪婪算法了。DTW是先计算起点到终点的最小值，然后从这个最小值回溯回去看看这个最小值都经过了哪些节点。

举个栗子：

已知：两个列向量a=[8 9 1]'，b=[ 2 5 4 6]'，其中'代表转置，也就是把行向量转换为列向量了

求：两个向量利用动态时间规整以后的最短距离

第一步：计算对应点的欧式距离矩阵d【注意开根号】

     6     3     4     27     4     5     31     4     3     5

第二步：计算累加距离D【从6出发到达5的累加距离】

     6     9    13    1513    10    14    1614    14    13    18

计算方法如下：

D(1,1)=d(1,1)=6

D(1,2)=D(1,1)+d(1,2)=9

...

D(2,2)=min(D(1,2),D(1,1),D(2,1))+d(2,2)=6+4=10

...

D(m,n)=min(D(m-1,n),D(m-1,n-1),D(m,n-1))+d(m,n)

网上有很多代码，但是大部分代码有误，比如网上流传比较多的错误版本：http://www.cnblogs.com/luxiaoxun/archive/2013/05/09/3069036.html

正确的代码可以自己写，挺简单，但是我找了一个可视化的代码，大家可以参考一下：

dtw.m

 

function [Dist,D,k,w,rw,tw]=dtw(r,t,pflag)
%
% [Dist,D,k,w,rw,tw]=dtw(r,t,pflag)
%
% Dynamic Time Warping Algorithm
% Dist is unnormalized distance between t and r
% D is the accumulated distance matrix
% k is the normalizing factor
% w is the optimal path
% t is the vector you are testing against
% r is the vector you are testing
% rw is the warped r vector
% tw is the warped t vector
% pflag  plot flag: 1 (yes), 0(no)
%
% Version comments:
% rw, tw and pflag added by Pau Mic[row,M]=size(r); if (row > M) M=row; r=r'; end;
[row,N]=size(t); if (row > N) N=row; t=t'; end;
d=sqrt((repmat(r',1,N)-repmat(t,M,1)).^2); %this makes clear the above instruction Thanks Pau Mic
d
D=zeros(size(d));
D(1,1)=d(1,1);for m=2:MD(m,1)=d(m,1)+D(m-1,1);
end
for n=2:ND(1,n)=d(1,n)+D(1,n-1);
end
for m=2:Mfor n=2:ND(m,n)=d(m,n)+min(D(m-1,n),min(D(m-1,n-1),D(m,n-1))); % this double MIn construction improves in 10-fold the Speed-up. Thanks Sven Mensingend
endDist=D(M,N);
n=N;
m=M;
k=1;
w=[M N];
while ((n+m)~=2)if (n-1)==0m=m-1;elseif (m-1)==0n=n-1;else [values,number]=min([D(m-1,n),D(m,n-1),D(m-1,n-1)]);switch numbercase 1m=m-1;case 2n=n-1;case 3m=m-1;n=n-1;endendk=k+1;w=[m n; w]; % this replace the above sentence. Thanks Pau Mic
end% warped waves
rw=r(w(:,1));
tw=t(w(:,2));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if pflag% --- Accumulated distance matrix and optimal pathfigure('Name','DTW - Accumulated distance matrix and optimal path', 'NumberTitle','off');main1=subplot('position',[0.19 0.19 0.67 0.79]);image(D);cmap = contrast(D);colormap(cmap); % 'copper' 'bone', 'gray' imagesc(D);hold on;x=w(:,1); y=w(:,2);ind=find(x==1); x(ind)=1+0.2;ind=find(x==M); x(ind)=M-0.2;ind=find(y==1); y(ind)=1+0.2;ind=find(y==N); y(ind)=N-0.2;plot(y,x,'-w', 'LineWidth',1);hold off;axis([1 N 1 M]);set(main1, 'FontSize',7, 'XTickLabel','', 'YTickLabel','');colorb1=subplot('position',[0.88 0.19 0.05 0.79]);nticks=8;ticks=floor(1:(size(cmap,1)-1)/(nticks-1):size(cmap,1));mx=max(max(D));mn=min(min(D));ticklabels=floor(mn:(mx-mn)/(nticks-1):mx);colorbar(colorb1);set(colorb1, 'FontSize',7, 'YTick',ticks, 'YTickLabel',ticklabels);set(get(colorb1,'YLabel'), 'String','Distance', 'Rotation',-90, 'FontSize',7, 'VerticalAlignment','bottom');left1=subplot('position',[0.07 0.19 0.10 0.79]);plot(r,M:-1:1,'-b');set(left1, 'YTick',mod(M,10):10:M, 'YTickLabel',10*rem(M,10):-10:0)axis([min(r) 1.1*max(r) 1 M]);set(left1, 'FontSize',7);set(get(left1,'YLabel'), 'String','Samples', 'FontSize',7, 'Rotation',-90, 'VerticalAlignment','cap');set(get(left1,'XLabel'), 'String','Amp', 'FontSize',6, 'VerticalAlignment','cap');bottom1=subplot('position',[0.19 0.07 0.67 0.10]);plot(t,'-r');axis([1 N min(t) 1.1*max(t)]);set(bottom1, 'FontSize',7, 'YAxisLocation','right');set(get(bottom1,'XLabel'), 'String','Samples', 'FontSize',7, 'VerticalAlignment','middle');set(get(bottom1,'YLabel'), 'String','Amp', 'Rotation',-90, 'FontSize',6, 'VerticalAlignment','bottom');% --- Warped signalsfigure('Name','DTW - warped signals', 'NumberTitle','off');subplot(1,2,1);set(gca, 'FontSize',7);hold on;plot(r,'-bx');plot(t,':r.');hold off;axis([1 max(M,N) min(min(r),min(t)) 1.1*max(max(r),max(t))]);grid;legend('signal 1','signal 2');title('Original signals');xlabel('Samples');ylabel('Amplitude');subplot(1,2,2);set(gca, 'FontSize',7);hold on;plot(rw,'-bx');plot(tw,':r.');hold off;axis([1 k min(min([rw; tw])) 1.1*max(max([rw; tw]))]);grid;legend('signal 1','signal 2');title('Warped signals');xlabel('Samples');ylabel('Amplitude');end


测试代码

 

test.m

 

clear
clc
a=[8 9 1 9 6 1 3 5]';
b=[2 5 4 6 7 8 3 7 7 2]';
[Dist,D,k,w,rw,tw] = DTW(a,b,1);
fprintf('最短距离为%d\n',Dist)
fprintf('最优路径为')
w


测试结果：

d =6     3     4     2     1     0     5     1     1     67     4     5     3     2     1     6     2     2     71     4     3     5     6     7     2     6     6     17     4     5     3     2     1     6     2     2     74     1     2     0     1     2     3     1     1     41     4     3     5     6     7     2     6     6     11     2     1     3     4     5     0     4     4     13     0     1     1     2     3     2     2     2     3最短距离为27
最优路径为
w =1     11     21     31     41     51     62     63     74     85     96    107    108    10


规整以后可视化如下

 

 

【更新日志：2019-9-11】

为了验证结果正确性，在python中也找到了一个工具库`dtw`

安装方法

pip install dtw

接收参数与返回参数

def dtw(x, y, dist, warp=1, w=inf, s=1.0):"""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:param int warp: how many shifts are computed.:param int w: window size limiting the maximal distance between indices of matched entries |i,j|.:param float s: weight applied on off-diagonal moves of the path. As s gets larger, the warping path is increasingly biased towards the diagonalReturns the minimum distance, the cost matrix, the accumulated cost matrix, and the wrap path."""

输入参数：前面x和y是输入的两个序列点，dist是计算距离的方法，因为有很多距离计算的方法，只不过通常使用欧式距离，直接使用numpy里面的linalg.norm即可。

返回值：这个d有点蒙，但是如果想得到上面的27那个结果，可以看accumulated cost matrix这个返回值，wrap path就是两条路径的对应点。

调用方法(官方案例):

import numpy as np# We define two sequences x, y as numpy array
# where y is actually a sub-sequence from x
x = np.array([2, 0, 1, 1, 2, 4, 2, 1, 2, 0]).reshape(-1, 1)
y = np.array([1, 1, 2, 4, 2, 1, 2, 0]).reshape(-1, 1)from dtw import dtweuclidean_norm = lambda x, y: np.abs(x - y)d, cost_matrix, acc_cost_matrix, path = dtw(x, y, dist=euclidean_norm)print(d)

本文代码和文章已经同步到微信公众号中，公众号与本博客将持续同步更新运动捕捉、机器学习、深度学习、计算机视觉算法，敬请关注

 

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/news/246750.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

从机器学习谈起

从机器学习谈起

很好的一篇文章，转载自博客园：http://www.cnblogs.com/subconscious/p/4107357.html 在本篇文章中，我将对机器学习做个概要的介绍。本文的目的是能让即便完全不了解机器学习的人也能了解机器学习，并且上手相关的实践。这篇文档也算…
阅读更多...
核函数

核函数

由于下一篇要学机器学习的另外一种模型——核模型，里面涉及到核函数，所以先找了一下核函数的相关知识。 在知乎上看到了一些比较好的解答，详细参考：http://www.zhihu.com/question/24627666 首先举一个核函数把低维空间映射到高…
阅读更多...
关于Matlab编程的思考（待续）

关于Matlab编程的思考（待续）

Matlab编程的规范化思考 1.并行化 2.释放内存 3.需要调参的变量太多，可考虑将变量都放到一个结构体里面。 4.find(y),就是要找到y中那些非零项的指引 5.代码运行出现问题的时候，在命令行输入why就可以得到答案 6.输入bench可以给电脑跑分。 7.hom…
阅读更多...
拉普拉斯锐化图像

拉普拉斯锐化图像

在图像增强中，平滑是为了消除图像中噪声的干扰，或者降低对比度，与之相反，有时为了强调图像的边缘和细节，需要对图像进行锐化，提高对比度。 图的边缘是指在局部不连续的特征。 简要介绍一下原理&#xff1…
阅读更多...
运动捕捉数据的描述ASF/AMC

运动捕捉数据的描述ASF/AMC

运动捕捉数据有多种格式：ASF/AMC，BVH，C3D等，这三个是比较常用的，一般的matlab实验用的是ASF/AMC，其次就是BVH。 ASF/AMC文件格式是Acclaim Games公司设计开发的，全称是Acclaim Skeleton File/A…
阅读更多...
应用深度学习(台大陈蕴侬李宏毅) Part1

应用深度学习(台大陈蕴侬李宏毅) Part1

History of Deep Learning Big Data & GPU 端到端 Universality Theorem Core Factors for Applied Deep Learning 参考文献 http://v.qq.com/vplus/578e2d6f5e1fadc1/foldervideos/8n1000201qzzkx5 Deep Learning ◦Goodfellow, Bengio, and Courville, “Deep Learning…
阅读更多...
世界坐标

世界坐标

世界坐标是最直观反映人体在世界坐标系下运动位置的变化信息，对分析运动行为有重要的作用。下面介绍如何根据ASF/AMC文件计算人体各个关节的世界坐标。 根据前面讲的ASF/AMC文件的格式，可以知道人体运动可以看做是通过根节点root的平移以及其他关节绕其父…
阅读更多...
人工神经网络——径向基函数(RBF)神经网络

人工神经网络——径向基函数(RBF)神经网络

此博客排版不好，重新用Markdown写了一篇，同时附上了代码，戳这里 本文摘自：《模式识别与智能计算——matlab技术实现第三版》与《matlab神经网络43个案例分析》 【注】蓝色字体为自己的理解部分 径向基函数神经网络的优点&#xf…
阅读更多...
李宏毅机器学习课程－Transfer Learning

李宏毅机器学习课程－Transfer Learning

深度学习 -> 强化学习 ->迁移学习（杨强教授报告) 李宏毅机器学习课程－Transfer Learning 迁移学习－吴恩达 freeze 待处理的 理解深层神经网络中的迁移学习及TensorFlow实现 Transfer Learning模式 Similar domain, different task…
阅读更多...
matlab实现RBF的相关函数

matlab实现RBF的相关函数

摘自《matlab神经网络43个案例分析》 (1)newrb() 该函数可以用来设计一个近似径向基网络(approximate RBF)。调用格式为： [net,tr]newrb(P,T,GOAL,SPREAD,MN,DF) 其中P为Q组输入向量组成的R*Q位矩阵，T为Q组目标分类向量组成的S*Q维矩阵。GOAL为均方误…
阅读更多...
李宏毅机器学习课程－Structured Learning

李宏毅机器学习课程－Structured Learning

Simple structured learning framework for python pystruct-github Slides for explaining structured prediction and PyStruct -github 一、Structured Learning-Unifed Framework 之前的input and output 都是vectors Training阶段就是找一个Ｆ来评估&#xff…
阅读更多...
matlab之norm函数

matlab之norm函数

简单点说就是用来计算范数的一个函数。 假设A是一个矩阵，那么norm(A)或者norm(A,2)计算的就是A的2范数；同理norm(A,1)计算的就是1范数了. 2范数：计算步骤是先计算A*A‘（这里A’代表转置，也就是原矩阵*（原…
阅读更多...
matlab之unwrap函数

matlab之unwrap函数

网上的说法： 要计算一个系统相频特性，就要用到反正切函数，计算机中反正切函数规定，在一、二象限中的角度为0～pi，三四象限的角度为0～-pi。 若一个角度从0变到2pi，但实际得到的结果是…
阅读更多...
Python xrange与range的区别

Python xrange与range的区别

xrange 与 range xrange 用法与 range 完全相同，所不同的是生成的不是一个list对象，而是一个生成器。 要生成很大的数字序列的时候，用xrange会比range性能优很多，因为不需要一上来就开辟一块很大的内存空间。 xrange 和 range 这…
阅读更多...
受限玻尔兹曼机准备知识——蒙特卡洛方法

受限玻尔兹曼机准备知识——蒙特卡洛方法

先了解几个基本概率知识，不急着看蒙特卡洛方法的定义，具体的MC方法参考网上各种资料。 两个比较好的学习MC方法的文章：蒙特卡洛方法入门 (结合了实例)和 蒙特卡洛方法 (推荐，非常详细) 更新日志：2016-11-19&#xff…
阅读更多...
受限玻尔兹曼机准备知识——MCMC方法和Gibbs采样

受限玻尔兹曼机准备知识——MCMC方法和Gibbs采样

先点明几个名词 MCMC方法：马尔可夫链-蒙特卡洛方法 (千万别叫成梅特罗波利斯蒙特卡罗方法了) Metropolis-Hastings采样：梅特罗波利斯-哈斯廷斯采样 Gibbs采样：吉布斯采样 还是介绍一下学习MCMC和Gibbs采样比较好的一个资料：随机采…
阅读更多...
受限玻尔兹曼机——简单证明

受限玻尔兹曼机——简单证明

花了很久看了一下玻尔兹曼机，感觉水有点深，总之一步一步来嘛~~~~ 先说一下一个非常好的参考资料： 受限玻尔兹曼机（RBM）学习笔记 ，有兴趣的可以再看看这篇文章的参考文献或者博客，写的也非常好&…
阅读更多...
受限玻尔兹曼机RBM实现及能量值思考——matlab实现

受限玻尔兹曼机RBM实现及能量值思考——matlab实现

网址：http://www.cs.toronto.edu/~hinton/MatlabForSciencePaper.html 这个代码主要是在mnist上做手写数字识别的代码，贴出来的目的主要是想研究一下在迭代过程中能量的变化情况。 1. 标准能量函数 标准的能量函数的表达式为： 那么就将这个…
阅读更多...
Cheat_Sheet ---Keras、Matlab、Matplotlib、Numpy、Pandas、Scikit-Learn、SciPy

Cheat_Sheet ---Keras、Matlab、Matplotlib、Numpy、Pandas、Scikit-Learn、SciPy

Cheat_Sheet ---KerasCheat_Sheet ---MatlabCheat_Sheet ---MatplotlibCheat_Sheet ---NumpyCheat_Sheet ---PandasCheat_Sheet ---Scikit-LearnCheat_Sheet ---SciPy参考文献 http://ddl.escience.cn/f/IDkq#path%2F8215264
阅读更多...
WPF中DataContext的绑定技巧-粉丝专栏

WPF中DataContext的绑定技巧-粉丝专栏

（关注博主后，在“粉丝专栏”，可免费阅读此文） 先看效果： 上面的绑定值都是我们自定义的属性，有了以上的提示，那么我们可以轻松绑定字段，再也不用担心错误了。附带源码。 …
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023