%% leave-one-trial-out

%data(normal): (ch*time*trial)
%excludeTrials: template should not include these trials
%testData(sick): (ch*time*trial)
%excludeTrials2: template should not include these trials

% output: power of target frequency

function [ out ,out2] = myTRCA( data ,excludeTrials, data2, excludeTrials2)

[V,D]  = trca_matrix(data(:,:,excludeTrials));
numTrials = (1:size(data,3));
template = squeeze(mean(data(:,:,excludeTrials),3));
newData = squeeze(mean(data(:,:,numTrials(~ismember(numTrials,numTrials(excludeTrials)))),3));
newData2 = squeeze(mean(data2(:,:,excludeTrials2),3));
template2 = squeeze(mean(data2(:,:,numTrials(~ismember(numTrials,numTrials(excludeTrials2)))),3));
freq = (8:0.2:11.8);

for i=1:20
    pNewData{i} = (V'*newData)';
    pTemplateData{i} = (V'*template)';
    powerPnewdata = 10*log(abs(fft(pNewData{i}(:,1:size(pNewData{i},2)))));
    newR = powerPnewdata(round(freq(i)*5+1),:);
    meanNoiseOfCh = mean(powerPnewdata(round(freq(~ismember(freq,freq(i)))*5+1),:),1);
    out(i) = mean(newR(1));%snr = newR./meanNoiseOfCh;snr(1)
    
    pNewData2{i} = (V'*newData2)';
    pTemplateData2{i} = (V'*template2)';
    
    powerPnewdata2 = 10*log(abs(fft(pNewData2{i}(:,1:size(pNewData2{i},2)))));
    newR2 = powerPnewdata2(round(freq(i)*5+1),:);
    meanNoiseOfCh2 = mean(powerPnewdata2(round(freq(~ismember(freq,freq(i)))*5+1),:),1);
    out2(i) = mean(newR2(1));%snr2 = newR4./meanNoiseOfCh2;snr2(1)
end



% 
% numTrials = (1:size(data,3));
% template = squeeze(mean(data(:,:,excludeTrials),3))';
% newData = squeeze(mean(data(:,:,numTrials(~ismember(numTrials,numTrials(excludeTrials)))),3))';
% 
% newData2 = squeeze(mean(data2(:,:,excludeTrials2),3))';
% template2 = squeeze(mean(data2(:,:,numTrials(~ismember(numTrials,numTrials(excludeTrials2)))),3))';
% 
% 
% if(isempty(numTrials(~ismember(numTrials,numTrials(excludeTrials)))))
%     template = newData;
%     template2 = newData2;
% end
% [c, wx, wy] = myCCA2__(template,5,8:0.2:11.8,harmonics);
% freq = (8:0.2:11.8);
% debug = zeros(20,harmonics*2,2);
% for i=1:20
%     pNewData{i} = newData*wx{i};
%     pTemplateData{i} = template*wx{i};
%     %for rr = 1:size(pNewData{i},2)
%     %[A,B,r,U,V]  = canoncorr(pNewData{i}(:,1:size(pNewData{i},2)),pTemplateData{i}(:,1:size(pNewData{i},2)));
%     %[A,B,r,U,V]  = canoncorr(pNewData{i},pTemplateData{i});
%     powerPnewdata = 10*log(abs(fft(pNewData{i}(:,1:size(pNewData{i},2)))));
%     newR = powerPnewdata(round(freq(i)*5+1),:);
%     %newR(rr) = r(1);
%     %end
%     out(i) = mean(newR(1));
% %     out(i) = r(1);
%     %debug(i,:,1) = newR;
%     
%     pNewData2{i} = newData2*wx{i};
%     pTemplateData2{i} = template2*wx{i};
%     
%     powerPnewdata2 = 10*log(abs(fft(pNewData2{i}(:,1:size(pNewData2{i},2)))));
%     newR2 = powerPnewdata2(round(freq(i)*5+1),:);
%     %for rr = 1:size(pNewData2{i},2)
% %     [A,B,r2,U,V]  = canoncorr(pNewData2{i}(:,1:size(pNewData2{i},2)),pTemplateData{i}(:,1:size(pNewData2{i},2)));
%     %[A,B,r,U,V]  = canoncorr(pNewData2{i},pTemplateData{i});
%     %    newR2(rr) = r(1);
%     %end
%     out2(i) = mean(newR2(1));
% %     out2(i) = r2(1);
%     %debug(i,:,2) = newR2;
%     if(length(excludeTrials)>=51)
% %         figure;subplot(2,1,1);plot(freq,powerPnewdata(round(freq(:)*5+1),:));
% %         hold on;plot(freq,mean(powerPnewdata(round(freq(:)*5+1),:),2),'r','LineWidth',2);hold off;
% %         subplot(2,1,2);plot(freq,powerPnewdata2(round(freq(:)*5+1),:))
% %         hold on;plot(freq,mean(powerPnewdata2(round(freq(:)*5+1),:),2),'r','LineWidth',2);hold off;
% %         title([num2str(freq(i)) ' hz']);
%     end
% end
 %figure;plot(freq,out,'b',freq,out2,'r');
end

function [corrVector, Wx, Wy] = myCCA2__(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(1,length(frequency));

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(input,template);
    corrVector(f) = r(1);
    Wx{f} = A;
    Wy{f} = B;
    index=1;
    template=[];
end
end

