%% data = ch*time
%  snr = power(signalFreq)/ avg(power(freq w/o signalfreqIndex))
% example mySNR(data,[10 11 12 13], 2, 5)
% return SNR of power(11)/avg.power(10 12 13)

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

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

end
function [out] = myCCA(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;
for harm=1:(harmonics)
    template(:,index) = sin(2*harm*pi*frequency*t);
    template(:,index+1) = cos(2*harm*pi*frequency*t);
    index= index+2;
end
[A,B,r,U,V] = canoncorr(input,template);
out = U;
end
