%
%  this routine signal averages to determine a template
%
function averaged_template = signal_average(trigger,lefti,righti,template,threshold,name);
%
  s = load(name);
%
%  for now, format the data
%
  temp = s';
  data = temp(:);
  N = length(data);
%
%  loop over all beats, and compute the correlation coefficient
%  don't use the beat used to make the template
%
  value = [];
  for i = 1:length(trigger);
      if (((trigger(i)-lefti)>=0) && (trigger(i)+righti < N))
%           temp = corrcoef(data(trigger(i)-lefti:trigger(i)+righti),template)
%           cc = temp(1,2)
          x = data(trigger(i)-lefti:trigger(i)+righti);
          y = template;
          cc = sum((x-mean(x)).*(y-mean(y)))/sqrt(sum((y-mean(y)).^2)*sum((x-mean(x)).^2));
          if ((cc > threshold) && (cc < 0.999))  
              value = [value trigger(i)]; 
          end;
      end;
  end;
%
%  determine how many beats to average over
%
  good = length(value);
  n = good+1;
  while (n > good)
    fprintf(1,'There are %d beats with correlation coefficients that exceed the threshold \n', length(value));
    nn = input('Average over how many beats? ','s');
    n = str2num(nn);
  end;
%
%  Now average
%
  temp = template;
  for i = 1:n-1
      temp = temp+data(value(i)-lefti:value(i)+righti);
  end;
  averaged_template = temp/n;
%
%  Now plot the the original and final templates
%
  fs = 125; % sample rate = 125 Hz
 % fs = 1000; % sample rate = 1000 Hz
  time = [0:length(template)-1]/fs;
  orient tall
  plot(time*1000,template,'--',time*1000,averaged_template,'-'); grid;
  axis('tight'); legend('Original','Averaged');
  title(['Averaged over ' num2str(n) ' beats']);
  xlabel('Time (msec)');
