%%

%input: data (ch*time*trials)
%trialsIndex: the trials to be used
%frequency: tragets frequency (in vector)
%harmonics: number of harmonics


function [out] = myCCA2_1(input, trialsIndex, resolution, frequency, harmonics)
%column major

template = zeros(size(input,2),harmonics*2);
time = size(input,2);
T = 1/time * resolution;                     % Sample time
L = size(input,2);                     % Length of signal
t = (0:L-1)*T;                % Time vector
index=1;
corrVector = zeros(length(trialsIndex),length(frequency));

for trials = 1:length(trialsIndex)
    data = zeros(size(input,2),size(input,1));
    data = squeeze(input(:,:,trialsIndex(trials))');
    for f = 1:length(frequency)
        for harm=1:(harmonics)
            template(:,index) = sin(2*harm*pi*frequency(f)*t);
            template(:,index+1) = cos(2*harm*pi*frequency(f)*t);
            index= index+2;
        end
        [A,B,r,U,V] = canoncorr(data,template);
        corrVector(trials,f) = r(1);
        index=1;
        template=[];
    end
end
out = mean(corrVector,1);

end