%
%  This program implements a simple turns count algorithm
%
%  N_lp is the number of points to use in the lowpass filter
%  cutoff is the cutoff frequency of the lowpass filter (in Hz)
%  delay is the estimted delay of the lowpass filter and derivative
%
function turns_count(N_lp,cutoff,delay,delta)
%
  x = importdata('emg_dog2.dat');
  fs = 10000;
  xrange = [14000:15000];
 % xrange = [11000:13000];
  x = x(xrange);
%
%  first compute the derivative and smooth it
%
  dx = zeros(1,length(x));
  dx(2:end) = (x(2:end)-x(1:end-1))*fs;
  [B,A] = butter(N_lp,2*cutoff/fs);
  dx = filter(B,A,dx);
%
%  Now look for sign changes in dx
%
%  signs = sign(dx(2:end).*dx(1:end-1));
  dsign = sign(dx(2:end).*dx(1:end-1))<0;
%
%  Find the indices that go with the sign changes
%
  index = find(dsign)-delay;
  while( index(1) < 1 )
      index = index(2:end);
  end;
  Ni = length(index);
%
  k = 1;
  change = [];
  while( k < Ni )
      j = k+1;
      while( (j <= Ni) && (abs(x(index(k))-x(index(j)))< delta)) 
          j = j+1; 
      end;
      if ( j <= Ni )  % must have found a good point
          change = [change index(k)];
          k = j;
      else
          k = k+1;
      end
  end;
%  
%  plot stuff
%
  orient tall
  time = xrange/fs;
  subplot(3,1,1); plot(time,x,'-',time(index),x(index),'o'); grid; xlabel('Time (sec)'); ylabel('x(t) and sign changes'); axis('tight');
  title([' N_{lp} = ' num2str(N_lp) ', f_{lp} = ' num2str(cutoff) ' Hz, delay = ' num2str(delay) ' \Delta = ' num2str(delta) ', Turning points = ' num2str(length(change))]); 
  subplot(3,1,2); plot(time,dx); grid; xlabel('Time (sec)'); ylabel('dx/dt'); axis('tight');
  subplot(3,1,3); plot(time,x,'-',time(change),x(change),'o'); grid; xlabel('Time (sec)'); ylabel('Turning Points'); axis('tight');
