function [ MS, D ] = ITD( X, N, d )

% function ITD
%
% outputs:  MS - Modal shape matrix
%           D - f_r matrix
% inputs:   X - trajectory matrix from n-channel sensor
%           N - amount of loaded samples
%           d - delta (shift in ITD)
%%
% size of input matrix
sizeX = size(X);

% cycle while: shifting input displacement data
i = 1;
while i <= sizeX(1,1)
    X_1(i, :) = X(i, 1 : N);
    X_2(i, :) = X(i, 1+d : N+d);
    X_3(i, :) = X(i, 1+2*d : N+2*d);
    i = i + 1;
end

% X_bar, X_hat
X_bar = [X_1; X_2];
X_hat = [X_2; X_3];

% ITD eigenvalue decomposition for time periods longer then 2n
A_bar = (X_hat * X_bar') * (X_bar * X_bar')^-1;

[V, D] = eig(A_bar);    % V - mode shapes
                        % D - natural frequencies, damping

% size of matrix V
sizeV = size(V);

% matrix of mode shapes
MS = V(1 : sizeV(1, 1)/2, :);

end

