%
%  this routine computes the (normalized) cross correlation of a template
%  with a signal. 
%
%  template is the template we are to use
%  signal is the signal we are comparing the template to
%  namex is the file the template was constructed from (x)
%  namey is the file we use to get (y)
%
function Rxy = normalized_cross_correlation(template,namex,namey,start)
%
%  load the signal file
%
  signalx = load(namex);
  signalx = signalx(start:end);
  signaly = load(namey);
  signaly = signaly(start:end);
%
%  get the length of the signal and the length of the template
%
  Nt = length(template)-1;
  Nlag = 100;
%
  for k = 1:Nlag
%
%       estimate the numerator for the cross correlation
%
%         num = sum(template.*signaly());
%
%       estimate the numerators for the autocorrelations
%
%         num_xx = sum(template.*signalx();
%         num_yy = sum(signaly().*signaly());
%
%       use the following denominator for computing the cross/autocorrelation 
%    
%         den = ;
%         Rxy(k) = num/den;
%         Rxx(k) = num_xx/den;
%         Ryy(k) = num_yy/den;
  end;
%
%  Normalize to a maximum of 1 for cross correlation plots
%
%   scale = ;
%
%  set the sampling rate
%
  fs = 100;
%
  figure; orient tall;
  
  subplot(5,1,1); plot([0:Nt]/fs,template); grid; xlabel('Time (sec)'); axis('tight'); ylabel('Template (x)'); 
  title([' Crosscorrelation using files ' namex ' and ' namey ', starting at ' , num2str(start)]);
  
  subplot(5,1,2); plot([0:Nt+Nlag-1]/fs, signalx(1:Nt+Nlag)); grid; xlabel('Time (sec)'); axis('tight');
  ylabel('Signal for R_{xx}');
  
  subplot(5,1,3); plot([0:Nt+Nlag-1]/fs, signaly(1:Nt+Nlag)); grid; xlabel('Time (sec)'); axis('tight');
  ylabel('Signal for R_{yy}');
  
%   subplot(5,1,4); plot([0:Nlag-1]/fs,Rxy*scale); grid;
%      ylabel('\rho_{xy} ');
%      xlabel('Delay in Seconds');
%      axis([0 0.5 -1.0 1.0]);
% %
%  Now try and find the power spectrum
%
%   [Xmag,Xphase,df,frq] = baf_fft();
%   Xmag = Xmag/max(Xmag);
%   XdB = 20*log10(Xmag);
%   subplot(5,1,5); plot(frq,XdB); grid; axis([0 max(frq) min(XdB) max(XdB)]);
%      ylabel('S_{xy}(in dB)');
%      xlabel('Frequency in Hz ');