clf;
clear all;

load e:\mepps mepp

time = linspace(0,.1,1000)';

%build a filter for estimating the rising phase
[b,a]=butter(4,.1);

plot(time,mepp(:,100),'r','linewidth',2) %reference mepp
hold on
plot(time,mepp(:,1),'g','linewidth',1)

%plot(time,filter(b,a,mepp(:,100)))
%align the mepps
%--smooth the peak to get the amp
%--get 1/2 point on rise
%--shift waveform
goodmepplist = [];

for i=1:100
   %get the rough maximum
   [Vmax,indexMax] = max(mepp(:,i));
   %95%range
   gt95indexs = find(mepp(:,i)>0.95*Vmax);
   gt95first = gt95indexs(1);
   gt95last = gt95indexs(end);
   %clip out the 95% range
   peakTimes = time(gt95first:gt95last);
   peakVolts = mepp(gt95first:gt95last, i);
   %fit a quadratic to the range
   quadfit = polyfit(peakTimes, peakVolts, 2);
   %and get the maximum of the FITTED quadratic
   [Vmax,index] = ...
      max(polyval(quadfit,peakTimes));
   %find the 50% rise point
   halfptlist = ...
      find(filter(b,a,mepp(100:gt95first,i))>=0.5*Vmax);
   if (length(halfptlist)>=1)
      halfpt = halfptlist(1);
      %shift mepp to align 50% points
      mepp(:,i) = ...
         [mepp(halfpt:end,i);mepp(1:halfpt-1,i)];
      goodmepplist = [goodmepplist i];
   else
      i %display incorrect mepp ids
   end
   
end

avgmepp = sum(mepp(:,goodmepplist),2)/100;
plot(time,avgmepp,'b','linewidth',2)
%plot(time,filter(b,a,avgmepp),'b','linewidth',2)
title('Average of 100 mepps, 15% std.dev. noise')
legend('no-noise reference',...
   'example noisy mepp',...
   'signal averaged'...
   );