%% (Internal) Calculate an estimate of the SNR of ECG signal, considering signal the ensemble average, and noise each realization minus the ensemble average.
%   
%   corr_gain = calc_correlation_gain(signal, heasig, references, limits, bRobust)
% 
% Arguments:
% 
%      + signal: signal to use
% 
%      + heasig: header information about the signal.
% 
%      + references: temporal samples where the synch occurs (QRS complex locations)
% 
%      + limits: the amount of samples w.r.t. references(i) to consider a a
%                time window around references(i), from references(i) -
%                limits(1) to references(i) + limits(2)
% 
%      + bRobust: Calculate the esemble average using mean or median.
%                 Default: use median.
% 
% Output:
% 
%      + corr_gain: Estimate of the SNR considering signal the ensemble
%                   average, and noise each realization minus the ensemble
%                   average.
% 
% Example:
% 
% See also CalcRRserieQuality
% 
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
% 
function corr_gain = calc_correlation_gain(signal, heasig, references, limits, bRobust)

corr_gain = nan;
[nsamp, nsig] = size(signal);
pattern_size = sum(limits);
lreferences = length(references);

sig_pack = pack_signal(signal, references(references > limits(1) & references < (nsamp-limits(2) )), limits, true);    

if( isempty(sig_pack) )
    noise_power = nan;
    return
end

if( nargin < 5 || isempty(bRobust) )
    pattern_avg = squeeze(median(sig_pack,3));
else
    pattern_avg = squeeze(mean(sig_pack,3));
end

sig_pack_sub_avg = bsxfun(@minus, sig_pack, pattern_avg);

corr_gain = 10*log10( mean(pattern_avg.^2) ./ 10.^squeeze(mean(log10(mean(sig_pack_sub_avg.^2,1)),3)) );
 
% % Obsolete: estimate noise from a random sampling simulating guessing the QRS locations. 
% 
% if( nargin < 4 || isempty(noise_power) )
%     noise_power = nan(30,nsig);
%     noise_disp = nan(30,nsig);
%     for ii = 1:30
%         noise_pack = pack_signal(signal, sort(randsample(limits(1):(nsamp-limits(2)), lreferences)), limits);    
% %         noise_avg = flipud(reshape(squeeze(median(noise_pack,3)), pattern_size, nsig));
%         noise_limits = prctile(noise_pack, [5 50 95], 3);
%         
%         noise_down = squeeze(noise_limits(:,:,1));
%         noise_avg = squeeze(noise_limits(:,:,2));
%         noise_up = squeeze(noise_limits(:,:,3));
%         
%         noise_power(ii,:) = mean(noise_avg.^2);
%         noise_disp(ii,:) = mean(noise_up - noise_down);
%     end
%     noise_power = nanmedian(noise_power);
%     noise_disp = nanmedian(noise_disp);
% end
% 
% corr_gain = 10*log10(mean(pattern_avg.^2) ./ noise_power );


% pattern_limits = prctile(sig_pack, [5 50 95], 3);
% 
% pattern_down = squeeze(pattern_limits(:,:,1));
% pattern_up = squeeze(pattern_limits(:,:,3));
% mean_pattern_disp = mean(pattern_up - pattern_down);

% corr_gain_2 = 20*log10( noise_disp ./ mean_pattern_disp );


% debug
% aux_val = [];
% for ii  = 1:100
%     noise_pack = pack_signal(signal, sort(randsample(limits(1):(nsamp-limits(2)), length(references))), limits);    
%     noise_avg = flipud(reshape(squeeze(median(noise_pack)), pattern_size, nsig));
%     aux_val = [aux_val; mean(noise_avg.^2)];
% end
% 
% figure
% hist( aux_val, 100)
% legend()