function EMGdat = loadEMG(infile,path1,gait,chans,EMGleg,FPbounds)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  loadEMG.m loads the EMG data stored in a .mat file produced by
%  the Harvard/Bertec setup.  If requested it also determines the
%  heel strike and toe off times for both feet.  Its output is a
%  data structure.  
%
%  Inputs:
%  infile   - name of .mat file containing data
%  path1    - subject initials and date
%  gait     - 1 if gait events required, 0 else
%  chans    - vector of EMG channels to be processed
%  EMGleg   - L or R; refers to leg with EMG electrodes
%  Fpbounds - if exists, automatically determines FP noise cutoff
%
%  Outputs:
%  EMGdat - data structure containing EMG info, data and gait events
%
%  Jared Markowitz
%  Summer 2011
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

path_full = ['~/Documents/MATLAB/data/' path1 '/raw/'];

cmd1 = ['load ''' path_full 'mat/' infile '.mat'';'];
eval(cmd1);
cmd2 = ['analog = ' infile '.Analog;'];
eval(cmd2);

EMGdat.subject  = path1;
EMGdat.trial    = infile;
EMGdat.EMGleg   = EMGleg;
EMGdat.Fs       = analog.Frequency;

if gait == 1
  cmd3 = ['FP1 = ' infile '.Force(1);'];
  eval(cmd3);
  cmd4 = ['FP2 = ' infile '.Force(2);'];
  eval(cmd4);  
  close all;

  if exist('FPbounds') == 1
    thresh1 = FPbounds;
  else  
    figure(1)
    plot(FP1.Force(3,:));
    grid on;
    axis([0 length(FP1.Force(3,:)) 0 100]);
    figure(2)
    plot(FP2.Force(3,:),'r');
    grid on;
    axis([0 length(FP2.Force(3,:)) 0 100]);
    thresh1 = input('Please input vector of force plate thresholds. Blue is FP1, Red is Fp2.\n:');
    close all;
  end;
    
  on1 = find(FP1.Force(3,:) > thresh1(1));
  on2 = find(FP2.Force(3,:) > thresh1(2));
  
  EMGdat.FPthresh = thresh1;
  
  EMGdat.Lto = on1(find(diff(on1)>100));
  EMGdat.Lhs = on1(find(diff(on1)>100)+1);
  EMGdat.Rto = on2(find(diff(on2)>100));
  EMGdat.Rhs = on2(find(diff(on2)>100)+1);
end;

% This may miss one GC at the beginning, but worth throwing one out
% if it ensures not having a partial one.

A1 = cellstr(analog.Labels);

for i = 1:length(chans)
  if (strcmp(path1,'AD052611')+strcmp(path1,'AM062411')+strcmp(path1,'EM072611')+strcmp(path1,'DS100611'))>0   
    cmd = ['ind = (find(strcmp(analog.Labels,''EMG' num2str(chans(i)) ''')==1));'];
    eval(cmd);
  else
    ind = chans(i)+18;  % EMG channels went from 19 -> 34
  end;  
  cmd = ['EMGdat.chan' num2str(chans(i)) ' = analog.Data(ind,:);' ];
  eval(cmd);
end;
