function [V,D] = trca_matrix(X)
% Task-related component analysis (TRCA)
%   X : eeg data (Num of channels * num of sample points * number of trials)

nChans  = size(X,1);
dur = size(X,2);
nTrials = size(X,3);

S = zeros(nChans, nChans);
% Computation of correlation matrices
for ch_i = 1:1:nChans
    for ch_j = 1:1:nChans
        for trial_i = 1:1:nTrials
            for trial_j = 1:1:nTrials
                if trial_i ~= trial_j
                    x_i = X(ch_i, :, trial_i);
                    x_j = X(ch_j, :,trial_j);
                    S(ch_i,ch_j) = S(ch_i,ch_j) + (x_i-mean(x_i,2))*(x_j-mean(x_j,2))';
                end % if
            end % trial_j
        end % trial_i
    end % ch_j
end % ch_i

X1 = X(1:nChans,:);
X1 = X1 - repmat(mean(X1,2),1,size(X1,2));
Q = X1*X1';

% TRCA eigenvalue algorithm
[V,D] = eig(Q\S);
%Y = V'*X;