function [corrVector, Wx, Wy] = myFBCCA(input, resolution, frequency, harmonics)
%column major
template = zeros(size(input,1),harmonics*2);
time = size(input,1);
T = 1/time * resolution;                     % Sample time
L = size(input,1);                     % Length of signal
t = (0:L-1)*T;                % Time vector
index=1;
corrVector = zeros(length(frequency),harmonics);
for subBand = 1:harmonics
    filterInput = bp40_FB1(double(input),512,subBand);
    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(filterInput,template);
        corrVector(f,subBand) = r(1);
        Wx{f,subBand} = A;
        Wy{f,subBand} = B;
        index=1;
        template=[];
    end
end
end