cross_corr_error(cross_corr,lags,len_sig,bias_flag)

function [error_corr x_axis] = cross_corr_error(cross_corr,lags,len_sig,bias_flag)
%This function estimates the error inculcated while calculating the
%correlation function of two signals. Since most of the times only a single
%realization of the signal is avalaible there is a variance assosciated
%with the calculated correlation function. Please refer to Digital Signal
%Processing Principles, Algorithms and Applications, Third Edition. By John
%G. Proakis and Dimitris G. Manolakis. Lags is the shift indices
%corresponding to the autocorrelation function. The inbuilt function
%xcorr/xcov gives this variable as an output. It is to be noted that the
%biased case has a lesser variance in the answer.
%Author - Ambarish Jash
%Date - 29th May, 2009

const_power = sum(cross_corr.^2);
%Now the second part of the summation formula given in the reference will
%be evaluated.
num_cores = feature('numCores');
%This has to be done for single core machines
if num_cores>1
    matlabpool (num_cores-1);
end
error_corr = zeros(1,max(size(cross_corr)));
parfor m_count = 1:max(size(cross_corr)) %These are the m values
    m_value = lags(m_count);
    init_count = 0; %Summation variable initialized.
   for n_count = 1:1:max(size(cross_corr)) %These are the n values
       n_value = lags(n_count);
       corr_pos1 = find(lags == n_value - m_value);
       corr_pos2 = find(lags == n_value + m_value);
       %The condition statement is necessary since there may be situations
       %where there is no data for the lag. In that case the init_count
       %value does not change
       if (size(corr_pos1,2) == 0 || size(corr_pos2,2) == 0)
           init_count = init_count + 0;
       else
           init_count = cross_corr(corr_pos1)*cross_corr(corr_pos2) + init_count;
       end
   end
   if (strcmp(bias_flag,'unbiased'))
       error_corr(m_count) = len_sig.*(init_count + const_power)./((len_sig - (abs(m_value)))^2);
   elseif (strcmp(bias_flag,'biased'))
       error_corr(m_count) = (init_count + const_power)./(len_sig);
   end
end
x_axis = lags;
matlabpool close