%
%  this script implements the Murthy and Rangaraj method for QRS detection
%
%  alpha is the pole of the filter to remove the baseline wander
%  N is the length of the window within which first order differences are
%  computed
%  a MA filter of length M is used to smooth g1 to get g
%  threshold is a fraction of max(g) to set the beat threshold at
%
function beat = laptop4(alpha,N,M,threshold,start,samples);
%
% load the data
%
  data = load('ecgpvc.dat');
  data = data(start:start+samples);
%
%  set the sample rate
%
  fs = 200;
%
%  remove the baseline wander
%
%   B = 
%   A = 
%
%  Now filter the signal
%
%   x = filter(B,A,data);
%
%  normalize the signal
%
%   n = 
%   x = x/max(x(n:end));
%
%  Now filter the signal
%
%   Nx = length(x);
%   g1 = zeros(1,Nx);
%   for n=N+1:Nx
%       arg = 0;
%       for i=1:N
%           arg = 
%       end;
%       g1(n) = arg;
%   end;
%  
%  apply a Mth order moving average filter to smooth out the signal
%
%   A = 
%   B = 
%
%  Now filter the signal
%
%   g = filter(B,A,g1);
%
%  set the threshold
%
%   th = threshold*max(g);
%
  blank = 20;  % a blanking period
%
%  find the peaks
%
%   beat = [];
%   n = M+1;
%   while (n < Nx-M)
%       if ( (g(n) > th) & (g(n)>g(n+1:n+M)) & (g(n)>(g(n-M:n-1))))
%           beat = [beat n];
%           n = n + blank;
%       else
%           n = n+1;
%       end;
%   end;
%
%  get the time array
%
  time_x = [0:length(x)-1]/fs;
%
%  plot the original and filtered signal
%
  figure;
  orient landscape
  subplot(3,1,1); plot(time_x,x); grid;   
    xlabel('Time (sec)'); ylabel('x'); axis([0 max(time_x) -1 1]);
    title(['M&R Method:  N = ' num2str(N) ', M = ' num2str(M), ', fraction = ', num2str(threshold) ', \alpha = ', num2str(alpha)]);
%   subplot(3,1,2); plot(time_x,g1); grid; 
%     xlabel('Time (sec)'); ylabel('g_1');  axis([0 max(time_x) min(g1) max(g1)]);
%   subplot(3,1,3); plot(time_x,g,'-',time_x,th*ones(1,length(g)),':'); grid; 
%     xlabel('Time (sec)'); ylabel('g'); axis([0 max(time_x) min(g) max(g)]);
% %
% %  now plot the beats and the trigger locations
% %
%   yy = data(beat);
%   figure;
%   orient landscape
%   plot(time_x,data,'-',(beat-1)/fs,yy,'o'); grid;
%     xlabel('Time (sec)');
%     title(['M&R Method:  N = ' num2str(N) ', M = ' num2str(M), ', fraction = ', num2str(threshold) ', \alpha = ', num2str(alpha)]);