%% 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 ] = mySNR( 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

