%
%  Example 3: Evaluation of the autocorrelation function 
%  using ensemble averages  .............
%

clear all;

N = 32 ;     % N samples in each realization
ff = 0.15 ;   %  Normalized frequency
ite = 1000;   % number of realizations (ensembles)
ite = 1 ;

snr = 20.0 ;  % SNR in dB

for n=1:5
   for k=1:5
     r(n, k) = 0.0 ;
   end;
end;

for jj=1:ite
   
v = (randn(1,N) + j * randn(1,N))/sqrt(2.0) ;    % this is complex noise
sf = exp(-log(10.0)*snr/20) ; 

s =  exp(j*2*pi*ff*[0:N-1]) + sf*v ;    % signal + noise 

%
% .... Evaluate autocorrelation using ensembel avarages
%

for n=1:5
   for k=1:5
     r(n, k) = r(n, k) + s(n)*s(n+k-1)' / ite ;
   end;
end;

end ;  % end of iteration loop jj

r     % Print the autocorrelation function r(n, n-k) ...a 5X5 matrix...

%
%  ... from results, note that r(n, n-k) is independent of n.
%      This due to stationarity. Therefore, the estimates can be impoved
%      by averaging in the n direction .... as shown below
%      NOTE: this is time avaeraging

auto = sum(r)/5

% Note : r matrix can be easily evaluated using a data correlation matrix.
%  as shown i.e. evaluate E(s*s') ....

%
%  Example 4: Evaluation of the data correlation matrix 
%  
%

dat = zeros(N,N) ;

for jj=1:ite
   
v = (randn(1,N) + j * randn(1,N))/sqrt(2.0) ;    % this is complex noise
sf = exp(-log(10.0)*snr/20) ; 

s =  exp(j*2*pi*ff*[0:N-1]) + sf*v ;    % signal + noise 

%
% .... Evaluate  ensembel expectations ...
%

dat = dat + s'*s/ite ;

end ;

dat(1:5,1:5)  % print only the the 5X5 sub-matrix

% Compare the results with previous example ....
% Due to stationarity the dat matrix is Toeplitz.
% The estimates can be improved by averaging along the diagonals
% to obtain the autocorrelation function, auto ...
%      NOTE: this is time avaeraging

%
%  Example 5: 
%      The ensemble  averaging and time averaging can be interchanged....
%  *** Evaluation of the autocorrelation function using time averages....
%

ite = 100 ;
ac = zeros(1,2*N-1) ;
snr = -5.0 ;

for jj=1:ite
   
v = (randn(1,N) + j * randn(1,N))/sqrt(2.0) ;    % this is complex noise
sf = exp(-log(10.0)*snr/20) ; 

s =  exp(j*2*pi*ff*[0:N-1]) + sf*v ;    % signal + noise 

%
% .... Evaluate  autocorrelation function ...
%

ac = ac + xcorr(s) ;
  if(jj == 1) 
    ac_1 = ac ;
  end;
  if(jj == 10) 
    ac_10 = ac/10 ;
  end;
  if(jj == 100) 
    ac_100 = ac/100 ;
  end;

end ;
tt = [-N+1:N-1];
figure(1); plot(tt,ac_1,tt,ac_10,tt,ac_100) ; grid ;

xlabel('Lag (samples)');
ylabel('Arbitory Units');
legend('avg = 1', 'avg = 10', 'avg = 100');

%
% Fourier transform of the autocorrelation gives the PSD ...
%

tt2 = tt/(2*N) ;
pw_1 = abs(fft(ac_1))/(N*N);
pw_10 = abs(fft(ac_10))/(N*N) ;
pw_100 = abs(fft(ac_100))/(N*N) ;
figure(2); plot(tt2,fftshift(pw_1),tt2,fftshift(pw_10),tt2,fftshift(pw_100)) ; grid ;

xlabel('Normalized Frequency');
ylabel('Power');
legend('avg = 1', 'avg = 10', 'avg = 100');