地质时间序列的动态时间扭曲，Matlab答案

【问题标题】：Dynamic Time Warping for geology time series, Matlab地质时间序列的动态时间扭曲，Matlab
【发布时间】：2012-02-02 15:55:39
【问题描述】：

背景：基本上，我使用像语音识别中使用的动态时间扭曲算法来尝试扭曲地质数据（从环境条件中过滤掉噪声）这两个问题之间的主要区别是 dtw 打印一个扭曲函数，允许输入的两个向量被扭曲，而对于我试图解决的问题，我需要保持一个参考向量不变，同时拉伸和收缩测试变量向量以适应。

这是matlab中的dtw：

function [Dist,D,k,w]=dtw()
%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
[t,r,x1,x2]=randomtestdata();
[rows,N]=size(t);
[rows,M]=size(r);
%for n=1:N
%    for m=1:M
%        d(n,m)=(t(n)-r(m))^2;
%    end
%end
d=(repmat(t(:),1,M)-repmat(r(:)',N,1)).^2; %this replaces the nested for loops from         above Thanks Georg Schmitz 

D=zeros(size(d));
D(1,1)=d(1,1);

for n=2:N
    D(n,1)=d(n,1)+D(n-1,1);
end
for m=2:M
    D(1,m)=d(1,m)+D(1,m-1);
end
for n=2:N
    for m=2:M
        D(n,m)=d(n,m)+min([D(n-1,m),D(n-1,m-1),D(n,m-1)]);
    end
end

Dist=D(N,M);
n=N;
m=M;
k=1;
w=[];
w(1,:)=[N,M];
while ((n+m)~=2)
    if (n-1)==0
        m=m-1;
    elseif (m-1)==0
        n=n-1;
    else 
      [values,number]=min([D(n-1,m),D(n,m-1),D(n-1,m-1)]);
      switch number
      case 1
        n=n-1;
      case 2
        m=m-1;
      case 3
        n=n-1;
        m=m-1;
      end
  end
    k=k+1;
    w=cat(1,w,[n,m]);
end
w=flipud(w)

%w is a matrix that looks like this:

%    1 1
%    1 2
%    2 2
%    3 3
%    3 4
%    3 5
%    4 5
%    5 6
%    6 6

所以这就是说第二个向量的第一个和第二个点都应该映射到第一个向量的第一个点。即 1 1 1 2 并且第一个向量上的第五个和第六个点应该在第六个点映射到第二个向量。等等，所以 w 包含扭曲数据的 x 坐标。

通常我会说

X1=w(:,1);
X2=w(:,2);
for i=1:numel(reference vector)
 Y1(i)=reference vector(X1(i));
 Y2(i)=test vector(X2(i));
end

但我不需要拉伸参考向量，因此我需要使用 X1 中的重复来了解如何缩小 Y2，并使用 X2 中的重复来了解如何拉伸 Y2，而不是使用 X1 中的重复来拉伸 Y1 并重复X2 拉伸 Y2。

我尝试使用 find 方法查找 X1 和 X2 中的重复，然后根据需要平均（收缩）或线性插值（拉伸），但代码变得非常复杂且难以调试。

这真的不清楚吗？我很难解释这个问题，但我只需要知道如何获取 w 并创建一个相应地拉伸和收缩的 Y2。

【问题讨论】：

我不确定这个广泛的解释是否相关。假设这是您唯一的问题，您也许可以只展示 W 和 Y2 现在的一两个示例以及 Y2 应该变成的样子。
你能展示一些关于你拥有什么和你想要什么的图表吗？
“这真的不清楚吗？”如果你问，你自己可能也有这种感觉。我当然很难追随你最终的追求。背景解释很好，但您需要简洁地陈述您的问题，理想情况下，背景信息几乎是多余的。
注意：我在这里回答了一个关于时间扭曲的问题：stackoverflow.com/questions/15283816/…，但我认为这与您的问题无关。

标签： matlab

【解决方案1】：

首先，这是从 wikipedia 上的伪代码翻译而来的 Matlab 中的 DTW：

t = 0:.1:2*pi;
x0 = sin(t) + rand(size(t)) * .1;
x1 = sin(.9*t) + rand(size(t)) * .1;

figure
plot(t, x0, t, x1);
hold on

DTW = zeros(length(x0), length(x1));
DTW(1,:) = inf;
DTW(:,1) = inf;
DTW(1,1) = 0;

for i0 = 2:length(x0)
    for i1 = 2:length(x1)
        cost = abs(x0(i0) - x1(i1));
        DTW(i0, i1) = cost + min( [DTW(i0-1, i1) DTW(i0, i1-1) DTW(i0-1, i1-1)] );
    end
end

无论您是将x_0 变形为x_1、x_1 变形为x_0，还是将它们彼此变形，您都可以从矩阵DTW 中得到答案。在你的情况下：

[cost, path] = min(DTW, [], 2);
plot(t, x1(path));
legend({'x_0', 'x_1', 'x_1 warped to x_0'});

【讨论】：

【解决方案2】：

我没有答案，但我一直在使用从维基百科文章中的伪代码实现的 @tokkot 代码。它有效，但我认为它缺少 DTW 的三个要求：

两个序列的第一个和最后一个点必须匹配，使用 min()，其中一个序列的一些（或许多）第一个和结束点丢失。
输出序列不是单调递增的。我用 x1(sort(path)) 代替，但我不相信这是真正的最小距离。
此外，由于我尚未找到的原因，扭曲序列的一些中间点丢失了，我认为这与 DTW 不兼容。

我仍在寻找像 DTW 这样的算法，其中一个序列是固定的（不是扭曲的）。我需要将等距温度测量的时间序列与另一个序列进行比较。第一个不能时移，没有意义。

【讨论】：