Динамическое временное деформирование для временных рядов геологии, Matlab

Question

Фон: В основном я использую динамический алгоритм деформирования времени, используемый в распознавании речи, чтобы попытаться деформировать геологические данные (отфильтровать шум из условий окружающей среды). Основное различие между этими двумя проблемы в том, что DTW печатает функцию коробления, что позволяет оба вектора, которые вводятся быть искаженно, в то время как для проблемы я пытаюсь решить я должен держать один опорный вектор постоянная при растяжении и усадки переменного вектора теста, чтобы соответствовать.

здесь Dtw в MATLAB:

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 и что пятая и шестая точки на первом векторе должны быть отображены во второй вектор в точке шесть. и т. д., так что 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, чтобы узнать, как уменьшить У2 и повторы в X2, чтобы узнать, как чтобы растянуть Y2, а не использовать повторы в X1, чтобы растянуть Y1 и повторить в X2, чтобы растянуть Y2.

Я попытался использовать метод поиска, чтобы найти повторы как в X1, так и в X2, а затем усреднить (сжать) или интерполировать линейно (растянуть) по мере необходимости, но код стал очень сложным и трудным для отладки.

Было ли это действительно непонятно? Мне было трудно объяснить эту проблему, но мне просто нужно знать, как взять w и создать Y2, растянутый и сжатый соответственно.

Answer 1

Во-первых, здесь DTW в Matlab переводе с псевдокода на wikipedia:

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'});