%
%  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 = M_and_R_QRS(name,alpha,N,M,fraction);
%
% load the data
%
  fs = 1000;
  data = load(name);
%
%  remove the baseline wander
%
  B = [1 -1]*(1+alpha)/2;
  A = [1 -alpha];
%
%  Now filter the signal
%
  x = filter(B,A,data);
%
%  normalize the signal
%
  n = floor(log(0.02)/log(alpha)+1);
  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 = arg+abs(x(n-i+1)-x(n-i)).^2*(N-i+1);
      end;
      g1(n) = arg;
  end;
%  
%  apply a moving average filter to smooth out the signal
%
  A = zeros(1,M); 
  A(1) = M;
  B = ones(1,M);
%
%  Now filter the signal
%
  g = filter(B,A,g1);
%
%  set the threshold
%
  th = fraction*max(g);
%
  blank = 20;  % a blanking period
%
  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
%
%
%  now plot the beats and the trigger locations
%
  yy = data(beat);
  figure;
  orient tall
  subplot(2,1,1); 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)]);
     title(['M&R Method:  N = ' num2str(N) ', M = ' num2str(M), ', fraction = ', num2str(fraction) ', \alpha = ', num2str(alpha)]);
  subplot(2,1,2)
  plot(time_x,data,'-',(beat-1)/fs,yy,'o'); grid; axis('tight')
    xlabel('Time (sec)');