%% leave-one-trial-out, apply training based CCA then calculate the SNR

%data: (ch*time*trial)
%excludeTrials: the trials index of input data i.e. template should not include these
%        trials

%output: one dimension snr vector over frequency
%ex: out(8:0.2:11.8), out(i) =  (power of frequency(i)) / mean(power of freqnecy(!i))

function [ snr ] = myCCA5( data , excludeTrials, harmonics, resolution,freq, varargin)
numTrials = (1:size(data,3));
newData = squeeze(mean(data(:,:,excludeTrials),3))';
template = squeeze(mean(data(:,:,numTrials(~ismember(numTrials,numTrials(excludeTrials)))),3))';
if(isempty(numTrials(~ismember(numTrials,numTrials(excludeTrials)))))
    template = newData;
end
[c, wx, wy] = myCCA2__(template,resolution,freq,harmonics);
frequency = freq;
for i=1:length(freq)
    pNewData{i} = newData*wx{i};
    tempData(i,:) = pNewData{i}(:,1);
end

for ch =1:size(tempData,1)
    snr(ch) = mySNR1(tempData(i,:),freq,ch,resolution);
end

end

function [ out ] = mySNR1( data,freq,signalFreqIndex,resolution )
out = zeros(size(data,1),1);

noiseIndex = freq(~ismember(freq,freq(signalFreqIndex)))*resolution+1;
signalIndex = freq(signalFreqIndex)*resolution+1;
for i=1:size(data,1)%for each channel
    tempPower = abs(fft(data(i,:)));
    noisepower = 0;
    for j=1:length(noiseIndex)
        noisepower = noisepower+tempPower(1,floor(noiseIndex(j)));
    end
    noisepower = noisepower/length(noiseIndex);
    out(i) = tempPower(1,floor(signalIndex))/ noisepower;
end

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



%% snr power ratio of target/nonTarget
% for f = 1:length(frequency)
%     nonTargetIndex = int16(freq2Index_2(frequency(~ismember(frequency,frequency(f)))));
%     
%     out(f) = tempPower(f,int16(freq2Index(frequency(f))))/...
%         mean(tempPower(f,int16(freq2Index(frequency(nonTargetIndex)))));
%     
% end


function [index] = freq2Index_2(hz)
index = (hz-9)*4+1;
end

% function [hz] = index2Freq(index)
% hz = (index-1)/5;
% end

function [index] = freq2Index(hz)
index = (hz)*4+1;
end



%% test function
function [template] = getSin(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
template(:,1) = sin(2*pi*frequency*t);

end
