function [mean_act,std_act,mean_stim,std_stim] = comb_avg_click(listfile,path1,proc,numpts,plotit,enter_gc,enter_th,HAMin1,HABin1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% comb_avg_click.m produces means and standard deviations of the
% stimulation and activation profiles of the trials contained in
% listfile. These are normalized to numpts per gait cycle.
% The script provides three different options for discarding
% anomalous gait cycles, including the default "clicking" method.
% It provides two different options for thresholding stimulation
% signals, including the clicking method and reading thresholds
% from a file.
%
% Inputs:
% listfile  - file containing list of activation files to be read.
% path1     - data directory (subject name and date).
% proc      - processing method (e.g. 'SangZaj' or 'lowpass').
% numpts    - number of points in normalized gait cycle.
% plotit    - 1 if just act plots desired, 2 if act and stim plots desired.
% enter_gc  - means to enter discarded GC.  Default is to select by clicking;
%             1 for reading from file, 2 for entering manually.
% enter_th  - means to enter thresholding levels.  Default is to select by clicking;
%             1 for reading from file.
% HAMin1    - 0 to use SMEM for HAM est, 1 to use BFLH.
% HABin1    - 0 to use GLUU for HAB est, 1 to use GMED.
%
% Outputs: 
% Saves all normalized activation profiles as well as mean and
% standard deviations to an output file with the same name as listfile
% but in a different directory. Also saves files produced by
% thresholding stimulation estimations in a different directory
% (but with same variable names).
%
% Example Usage:
% comb_avg_click('walk125','JB011613','low100_mvc_sh',100,[1 2],1,0,0,0)
%
% Jared Markowitz
% Spring 2013
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath('~/Dropbox/harvard/processMOT');

path_in   = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/act_est/combined_new'];
path_list = ['~/Documents/MATLAB/data/' path1 '/processed/mot/simm_out/'];

fid1  = fopen([path_list '/' listfile '.txt']);
tline = fgetl(fid1);
list1 = {};

while ischar(tline)
  list1 = [list1; {tline}];  
  tline = fgetl(fid1);
end; 

fclose('all');

% Ensure that enter_th, enter_gc, plotit, save options are consistent:
if exist('enter_th')>0 && length(find(enter_th ~= 1)) > 0
  save_thr = 1;
  plotit   = [1 2];
end;
if exist('enter_gc')>0 && enter_gc ~= 1
  save_rem = 1;
  enter_th = 0;
end;

if plotit > 0
  close all;
end;

% Load Perry muscles for reference:
load ~/Documents/MATLAB/data/Perry_muscles/stim_act.mat
p_sVAS  = p_sVLAT;   clear p_sVLAT;
p_sADDL = p_sAL;     clear p_sAL;
p_sADDM = p_sAM;     clear p_sAM;

% Initialize cells:
load([path_in '/' list1{1}]);  
musc = fieldnames(muscles);

for l=1:size(musc,1)
  eval(['a' musc{l} ' = {};']);
  eval(['s' musc{l} ' = {};']);
end;

% Load all gait cycles into cells:
for i = 1:length(list1)
  clear rsang* rlow* act* muscles;
  load([path_in '/' list1{i}]);

  if strcmp(proc,'low100_mvc_lo') == 1
    act1   = act100_mvc_lo;
    stim1  = rlow100_mvc;    
  elseif strcmp(proc,'low100_mvc_sh') == 1
    act1   = act100_mvc_sh;
    stim1  = rlow100_mvc;    
  elseif strcmp(proc,'low100_fw_lo') == 1
    act1   = act100_fw_lo;
    stim1  = rlow100_fw;    
  elseif strcmp(proc,'low100_fw_sh') == 1
    act1   = act100_fw_sh;
    stim1  = rlow100_fw;    
  elseif strcmp(proc,'sang_mvc_lo') == 1
    act1   = act_mvc_lo;
    stim1  = rsang_mvc;    
  elseif strcmp(proc,'sang_mvc_sh') == 1
    act1   = act_mvc_sh;
    stim1  = rsang_mvc;    
  elseif strcmp(proc,'sang_fw_lo') == 1
    act1   = act_fw_lo;
    stim1  = rsang_fw;    
  elseif strcmp(proc,'sang_fw_sh') == 1
    act1   = act_fw_sh;
    stim1  = rsang_fw;    
  else
    fprintf('Faulty proc.\n');
  end;  
  
  for nn = 1:size(musc,1)
    eval(['chan_no = muscles.' musc{nn} ';']);
    eval(['stim.' musc{nn} ' = stim1.chan' num2str(chan_no) ';']);
  end;    

  eval(['heels = act.' act.EMGleg 'hs;']);

  for m = 1:size(musc,1)
    for j = 1:(length(heels)-1)
      start1 = heels(j);
      stop1  = heels(j+1);
      eval(['a' musc{m} '=[a' musc{m} '; {act1.' musc{m} '(start1:stop1)}];']);
      eval(['s' musc{m} '=[s' musc{m} '; {stim.' musc{m} '(start1:stop1)}];']);    
    end; 
  end;  
end;
  
% Convert to percent gait cycle:

for n = 1:size(musc,1)
  eval(['n_a' musc{n} ' = [];']);
  eval(['n_s' musc{n} ' = [];']);
end;

Fs = act.Fs;
dt = 1/Fs;
numpts = numpts+1;

eval(['len_k=size(a' musc{1} ',1);']);

for k = 1:len_k
  eval(['x1 = 0:dt:(length(a' musc{1} '{k})-1)*dt;']);
  x2 = linspace(0,x1(end),numpts);
  for p = 1:size(musc,1)
    eval(['i_a' musc{p} '=interp1(x1,a' musc{p} '{k},x2,''spline'');']);
    %eval(['i_a' musc{p} '=i_a' musc{p} '-min(i_a' musc{p} ');']);
    eval(['n_a' musc{p} '=[n_a' musc{p} '; i_a' musc{p} '];']);
    eval(['i_s' musc{p} '=interp1(x1,s' musc{p} '{k},x2,''spline'');']);
    %eval(['i_s' musc{p} '=i_s' musc{p} '-min(i_s' musc{p} ');']);
    eval(['n_s' musc{p} '=[n_s' musc{p} '; i_s' musc{p} '];']);
  end;
end;

% Remove outlier gait cycles.  Note that these may be different for
% each muscle:

if strcmp(proc(1:3),'low') == 1; 
  procr = 'low100';
elseif strcmp(proc(1:4),'sang') == 1; 
  procr = 'sang';
else
  fprintf('Faulty proc input.\n');
  return;
end;

if exist('enter_gc') == 1 && enter_gc == 1
  path_rem  = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/removed/' procr ];
  eval(['load ' path_rem '/' listfile '.mat;']);
  clear HAMin; clear HABin;            
elseif exist('enter_gc') == 1 && enter_gc == 2
  for q = 1:size(musc,1)
    eval(['r' musc{q} '= [];']);
    eval(['r' musc{q} '=input(''Enter vector of ' musc{q} ' gait cycles to be removed:\n'');']);
  end;
else
  eval(['input_x = ones(size(n_a' musc{1} ',1),1)*(0:(numpts-1));']);
  for kk = 1:size(musc,1)
    eval(['r' musc{kk} '=select_curves(input_x,n_a' musc{kk} ',10,''' musc{kk} ' Act. Est.'',''Pct. GC'',''Act. Est.'');']);
  end;
  clear input_x gco;
end;

%if exist('rGLUU') == 0
%  input_x = ones(size(n_aGLUU,1),1)*(0:(numpts-1));
%  rGLUU =select_curves(input_x,n_aGLUU,10,'GLUU Act. Est.','Pct. GC','Act. Est.');
%  save_rem = 1;
%end;

% Find arrays, means, sds with anomalous GCs removed:
for ll = 1:size(musc,1)
  eval(['a' musc{ll} '(r' musc{ll} ') = [];']);
  eval(['s' musc{ll} '(r' musc{ll} ') = [];']);  
  eval(['n_a' musc{ll} '(r' musc{ll} ',:) = [];']);
  eval(['n_s' musc{ll} '(r' musc{ll} ',:) = [];']);
  eval(['m_a' musc{ll} ' = mean(n_a' musc{ll} ',1);']);  
  eval(['m_s' musc{ll} ' = mean(n_s' musc{ll} ',1);']);  
  eval(['sd_a' musc{ll} ' = std(n_a' musc{ll} ',0,1);']);  
  eval(['sd_s' musc{ll} ' = std(n_s' musc{ll} ',0,1);']);  
end;  

if length(find(plotit==1)) > 0
  addpath '~/Dropbox/harvard/processMOT/';  
  for ii = 1:size(musc,1)
    eval(['figure(' num2str(ii+50) ');']);
    eval(['niceBars(0:(numpts-1),m_a' musc{ii} ',sd_a' musc{ii} ',''b'',0.5);']);
    hold on;
    eval(['plot(0:(numpts-1),m_a' musc{ii} ',''k'',''linewidth'',4);']);
    eval(['title(''' musc{ii} ' Activation Profile'',''fontsize'',12);']);
    xlabel('Percent GC','fontsize',12);
    ylabel('Estimated Activation','fontsize',12);
    grid on; axis tight;
    %eval('print -dpng ''' musc{ii} '_a.png'';']);
  end;
end;
if length(find(plotit==2)) > 0
  for jj = 1:size(musc,1)
    eval(['figure(' num2str(jj+30) ');']);    
    eval(['niceBars(0:(numpts-1),m_s' musc{jj} ',sd_s' musc{jj} ',''b'',0.5);']);
    hold on;
    eval(['plot(0:(numpts-1),m_s' musc{jj} ',''k'',''linewidth'',4);']);
    eval(['title(''' musc{jj} ' Stimulation Profile'',''fontsize'',12);']);
    xlabel('Percent GC','fontsize',12);
    ylabel('Estimated Stimulation','fontsize',12);
    grid on; axis tight;   
    %eval('print -dpng ''' musc{jj} '_s.png'';']);
  end;
end; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Thresholding

% Can likely do one per stimulation estimation method (Sang/LP)
% and use this across speed and with scaling for MVC->FW.

% Determine threshold levels:
if exist('enter_th') == 1 && enter_th == 1
  eval(['load ~/Documents/MATLAB/data/' path1 '/processed/EMG/thresh_values/' procr '/' listfile '.mat']);
  % Would add support for FW normalization here....
elseif exist('enter_th') == 1 && enter_th ~= 1
  for pp = 1:size(musc,1)  
    eval(['figure(' num2str(pp+30) ');']); 
    hold on;
    if strcmp(musc{pp},'BFLH')
      plot(p_sSMEM,'k','linewidth',4);
    elseif strcmp(musc{pp},'GLU')
      plot(p_sGLU_low,'k','linewidth',4);
      plot(p_sGLU_upp,'y','linewidth',4);        
    elseif strcmp(musc{pp},'GLUU')
      plot(p_sGLU_upp,'k','linewidth',4);        
    else  
      eval(['plot(p_s' musc{pp} ',''k'',''linewidth'',4);']);
    end;
    axis tight;
    resp1 = 'n';
    h1 = plot([0 100],[0 0],'r','linewidth',3);
    h2 = plot([0 100],[0 0],'m','linewidth',3);
    while sum(strcmp(resp1,{'Y';'y'})) == 0
      delete(h1);
      delete(h2);
      fprintf('Click vertical level for threshold:\n');
      [xthresh,ythresh] = ginput(1);
      hold on;
      h1 = plot([0 100],[ythresh ythresh],'r','linewidth',4);
      eval(['q1 = m_s' musc{pp} ';']);
      ind = find(q1<=ythresh);
      q1(ind) = 0;
      h2 = plot(0:(numpts-1),q1,'m','linewidth',4);
      resp1 = input('Good enough? (y/n):\n','s');
    end;  
    eval(['thresh_values.' musc{pp} '=ythresh;']);
  end;  
else
  fprintf('Invalid enter_th input.\n');
end;

% Perform thresholding:
for ppp = 1:size(musc,1)
  eval(['num_gc = size(n_a' musc{ppp} ',1);']);
  eval(['full_stim_' musc{ppp} ' = [];']);
  for qq = 1:num_gc
    eval(['full_stim_' musc{ppp} ' = cat(1,full_stim_' musc{ppp} ',s' musc{ppp} '{qq}) ;']);
  end;  
  eval(['ind=find(full_stim_' musc{ppp} '<=thresh_values.' musc{ppp} ');']);
  eval(['thresh_stim_' musc{ppp} ' = full_stim_' musc{ppp} ';']);
  eval(['thresh_stim_' musc{ppp} '(ind) = 0;']);  
end;

% Pass remaining stim through Raasch bilinear dynamics:
for rr = 1:size(musc,1)
  eval(['tt = (0:1/Fs:(length(thresh_stim_' musc{rr} ')-1)/Fs)'';']);
  simopts = simset('FixedStep',1/Fs,'SrcWorkspace','current');
  if strcmp(musc{rr},'SOL') == 1
    mfft = 0.20;
  elseif strcmp(musc{rr},'GAS') == 1
    mfft = 0.50;
  elseif strcmp(musc{rr},'TA') == 1
    mfft = 0.35;
  elseif strcmp(musc{rr},'BFSH') == 1    
    mfft = 0.35;
  elseif strcmp(musc{rr},'BFLH') == 1
    mfft = 0.35;
  elseif strcmp(musc{rr},'SMEM') == 1
    mfft = 0.35;
  elseif strcmp(musc{rr},'VAS') == 1
    mfft = 0.50;
  elseif strcmp(musc{rr},'RF') == 1
    mfft = 0.65;
  elseif strcmp(musc{rr},'ILL') == 1    
    mfft = 0.50;
  elseif strcmp(musc{rr},'GLU') == 1
    mfft = 0.45;
  elseif strcmp(musc{rr},'GMED') == 1
    mfft = 0.50;
  elseif strcmp(musc{rr},'GLUU') == 1
    mfft = 0.50;  
  elseif strcmp(musc{rr},'ADDL') == 1
    mfft = 0.35;
  elseif strcmp(musc{rr},'ADDM') == 1
    mfft = 0.45;
  else
    fprintf('Erroneous muscle label.\n');
    return;
  end;      
  if strcmp(proc(end-1:end),'sh') == 1
    k1 = 1/(.01*8.4/(4.8*(1+1.5*mfft)))-1/(.05*8.4/(4.8*(1+1.5*mfft)));
    k2 = 1/(.05*8.4/(4.8*(1+1.5*mfft)));
  elseif strcmp(proc(end-1:end),'lo') == 1
    k1 = 1/(.08-.047*mfft)-1/(.09-.056*mfft);
    k2 = 1/(.09-.056*mfft);
  else fprintf('Invalid proc.\n');
    return;
  end;
  eval(['x_in = [tt thresh_stim_' musc{rr} '];']);
  sim('raasch.mdl',tt(end),simopts);
  eval(['a' musc{rr} '_thresh = {};']);
  eval(['s' musc{rr} '_thresh = {};']);
  eval(['n_a' musc{rr} '_thresh = [];']);
  eval(['n_s' musc{rr} '_thresh = [];']);
  eval(['num_gc = size(n_a' musc{rr} ',1);']);
  ind1 = 1;
  for qqq = 1:num_gc
    eval(['ind2 = ind1-1+length(s' musc{rr} '{qqq});']);
    eval(['a' musc{rr} '_thresh = [a' musc{rr} '_thresh; {a_out(ind1:ind2)}];']);
    eval(['s' musc{rr} '_thresh = [s' musc{rr} '_thresh; {thresh_stim_' musc{rr} '(ind1:ind2)}];']);
    eval(['x1 = 0:dt:(length(a' musc{rr} '_thresh{qqq})-1)*dt;']);
    x2 = linspace(0,x1(end),numpts);
    eval(['i_a' musc{rr} '_thresh=interp1(x1,a' musc{rr} '_thresh{qqq},x2,''spline'');']);
    %eval(['i_a' musc{rr} '_thresh=i_a' musc{rr} '_thresh-min(i_a' musc{rr} '_thresh);']);
    eval(['n_a' musc{rr} '_thresh=[n_a' musc{rr} '_thresh; i_a' musc{rr} '_thresh];']);
    eval(['i_s' musc{rr} '_thresh=interp1(x1,s' musc{rr} '_thresh{qqq},x2,''spline'');']);
    %eval(['i_s' musc{rr} '_thresh=i_s' musc{rr} '_thresh-min(i_s' musc{rr} '_thresh);']);
    eval(['n_s' musc{rr} '_thresh=[n_s' musc{rr} '_thresh; i_s' musc{rr} '_thresh];']);  
    ind1 = ind2 + 1;
  end;
  eval(['m_a' musc{rr} '_thresh = mean(n_a' musc{rr} '_thresh,1);']);  
  eval(['m_s' musc{rr} '_thresh = mean(n_s' musc{rr} '_thresh,1);']);  
  eval(['sd_a' musc{rr} '_thresh = std(n_a' musc{rr} '_thresh,0,1);']);  
  eval(['sd_s' musc{rr} '_thresh = std(n_s' musc{rr} '_thresh,0,1);']);  
  % Plot:
  if length(find(plotit==1)) > 0
    eval(['figure(' num2str(rr+50) ');']);
    hold on; 
    eval(['niceBars(0:(numpts-1),m_a' musc{rr} '_thresh,sd_a' musc{rr} '_thresh,''g'',0.5);']);
    axis tight;
  end;
  if length(find(plotit==2)) > 0
    eval(['figure(' num2str(rr+30) ');']);    
    hold on;
    eval(['niceBars(0:(numpts-1),m_s' musc{rr} '_thresh,sd_s' musc{rr} '_thresh,''g'',0.5);']);
    axis tight;
  end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Choose muscles and save:

if isfield(muscles,'SMEM')==1 && isfield(muscles,'BFLH')==1
  checkit1 = 0;
  while checkit1 < 1
    if exist('HAMin1')==0
      HAMin = input('Enter 0 if SMEM to be used for HAM, 1 if BFLH to be used:\n');
    else
      HAMin = HAMin1;
    end;
    if HAMin == 0
      m_aHAM = m_aSMEM; n_aHAM = n_aSMEM; sd_aHAM = sd_aSMEM; aHAM = aSMEM;
      m_sHAM = m_sSMEM; n_sHAM = n_sSMEM; sd_sHAM = sd_sSMEM; sHAM = sSMEM;   
      m_aHAM_thresh = m_aSMEM_thresh; n_aHAM_thresh = n_aSMEM_thresh; sd_aHAM_thresh = sd_aSMEM_thresh; aHAM_thresh = aSMEM_thresh;
      m_sHAM_thresh = m_sSMEM_thresh; n_sHAM_thresh = n_sSMEM_thresh; sd_sHAM_thresh = sd_sSMEM_thresh; sHAM_thresh = sSMEM_thresh;   
      checkit1 = 1;
    elseif HAMin == 1
      m_aHAM = m_aBFLH; n_aHAM = n_aBFLH; sd_aHAM = sd_aBFLH; aHAM = aBFLH;
      m_sHAM = m_sBFLH; n_sHAM = n_sBFLH; sd_sHAM = sd_sBFLH; sHAM = sBFLH;
      m_aHAM_thresh = m_aBFLH_thresh; n_aHAM_thresh = n_aBFLH_thresh; sd_aHAM_thresh = sd_aBFLH_thresh; aHAM_thresh = aBFLH_thresh;
      m_sHAM_thresh = m_sBFLH_thresh; n_sHAM_thresh = n_sBFLH_thresh; sd_sHAM_thresh = sd_sBFLH_thresh; sHAM_thresh = sBFLH_thresh;
      checkit1 = 1;
    else
    end;
  end;  
elseif isfield(muscles,'SMEM')==1 && isfield(muscles,'BFLH')==0 
  m_aHAM = m_aSMEM; n_aHAM = n_aSMEM; sd_aHAM = sd_aSMEM; aHAM = aSMEM;
  m_sHAM = m_sSMEM; n_sHAM = n_sSMEM; sd_sHAM = sd_sSMEM; sHAM = sSMEM;   
  m_aHAM_thresh = m_aSMEM_thresh; n_aHAM_thresh = n_aSMEM_thresh; sd_aHAM_thresh = sd_aSMEM_thresh; aHAM_thresh = aSMEM_thresh;
  m_sHAM_thresh = m_sSMEM_thresh; n_sHAM_thresh = n_sSMEM_thresh; sd_sHAM_thresh = sd_sSMEM_thresh; sHAM_thresh = sSMEM_thresh;   
  HAMin  = 0;  
elseif isfield(muscles,'SMEM')==0 && isfield(muscles,'BFLH')==1 
  m_aHAM = m_aBFLH; n_aHAM = n_aBFLH; sd_aHAM = sd_aBFLH; aHAM = aBFLH;
  m_sHAM = m_sBFLH; n_sHAM = n_sBFLH; sd_sHAM = sd_sBFLH; sHAM = sBFLH;
  m_aHAM_thresh = m_aBFLH_thresh; n_aHAM_thresh = n_aBFLH_thresh; sd_aHAM_thresh = sd_aBFLH_thresh; aHAM_thresh = aBFLH_thresh;
  m_sHAM_thresh = m_sBFLH_thresh; n_sHAM_thresh = n_sBFLH_thresh; sd_sHAM_thresh = sd_sBFLH_thresh; sHAM_thresh = sBFLH_thresh;
  HAMin  = 1;  
else  
  fprintf('No HAM saved.\n');
end;

if isfield(muscles,'GLUU')==1 && isfield(muscles,'GMED')==1
  checkit2 = 0;
  while checkit2 < 1
    if exist('HABin1')==0
      HABin = input('Enter 0 if GLUU to be used for HAB, 1 if GMED to be used:\n');
    else
      HABin = HABin1;
    end;
    if HABin == 0
      m_aHAB = m_aGLUU; n_aHAB = n_aGLUU; sd_aHAB = sd_aGLUU; aHAB = aGLUU;
      m_sHAB = m_sGLUU; n_sHAB = n_sGLUU; sd_sHAB = sd_sGLUU; sHAB = sGLUU;   
      m_aHAB_thresh = m_aGLUU_thresh; n_aHAB_thresh = n_aGLUU_thresh; sd_aHAB_thresh = sd_aGLUU_thresh; aHAB_thresh = aGLUU_thresh;
      m_sHAB_thresh = m_sGLUU_thresh; n_sHAB_thresh = n_sGLUU_thresh; sd_sHAB_thresh = sd_sGLUU_thresh; sHAB_thresh = sGLUU_thresh;   
      checkit2 = 1;
    elseif HABin == 1
      m_aHAB = m_aGMED; n_aHAB = n_aGMED; sd_aHAB = sd_aGMED; aHAB = aGMED;
      m_sHAB = m_sGMED; n_sHAB = n_sGMED; sd_sHAB = sd_sGMED; sHAB = sGMED;   
      m_aHAB_thresh = m_aGMED_thresh; n_aHAB_thresh = n_aGMED_thresh; sd_aHAB_thresh = sd_aGMED_thresh; aHAB_thresh = aGMED_thresh;
      m_sHAB_thresh = m_sGMED_thresh; n_sHAB_thresh = n_sGMED_thresh; sd_sHAB_thresh = sd_sGMED_thresh; sHAB_thresh = sGMED_thresh;   
      checkit2 = 1;
    else
    end;
  end;  
elseif isfield(muscles,'GLUU')==1 && isfield(muscles,'GMED')==0 
  m_aHAB = m_aGLUU; n_aHAB = n_aGLUU; sd_aHAB = sd_aGLUU; aHAB = aGLUU;
  m_sHAB = m_sGLUU; n_sHAB = n_sGLUU; sd_sHAB = sd_sGLUU; sHAB = sGLUU;   
  m_aHAB_thresh = m_aGLUU_thresh; n_aHAB_thresh = n_aGLUU_thresh; sd_aHAB_thresh = sd_aGLUU_thresh; aHAB_thresh = aGLUU_thresh;
  m_sHAB_thresh = m_sGLUU_thresh; n_sHAB_thresh = n_sGLUU_thresh; sd_sHAB_thresh = sd_sGLUU_thresh; sHAB_thresh = sGLUU_thresh;   
  HABin  = 0;  
elseif isfield(muscles,'GLUU')==0 && isfield(muscles,'GMED')==1 
  m_aHAB = m_aGMED; n_aHAB = n_aGMED; sd_aHAB = sd_aGMED; aHAB = aGMED;
  m_sHAB = m_sGMED; n_sHAB = n_sGMED; sd_sHAB = sd_sGMED; sHAB = sGMED;   
  m_aHAB_thresh = m_aGMED_thresh; n_aHAB_thresh = n_aGMED_thresh; sd_aHAB_thresh = sd_aGMED_thresh; aHAB_thresh = aGMED_thresh;
  m_sHAB_thresh = m_sGMED_thresh; n_sHAB_thresh = n_sGMED_thresh; sd_sHAB_thresh = sd_sGMED_thresh; sHAB_thresh = sGMED_thresh;   
  HABin  = 1;  
else  
  fprintf('No HAB saved.\n');
end;

clear a_*;
clear act*
clear i_*

path_out1_thr = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/avg_act_thresh/' proc ];
if exist(path_out1_thr) == 0
  eval(['mkdir ' path_out1_thr]);
end;
outfile1_thr = [path_out1_thr '/' listfile '.mat'];

save(outfile1_thr,'*_thresh','HA*in','thresh_values');

clear *_thresh simopts start1 stim stim1 stop1; 

path_out1 = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/avg_act/' proc ];
if exist(path_out1) == 0
  eval(['mkdir ' path_out1]);
end;
outfile1 = [path_out1 '/' listfile '.mat'];
save(outfile1,'m_*','n_*','sd_*','HA*in','a*','s*');

clear rr;
if exist('save_rem')>0 && save_rem == 1 
  path_out2 = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/removed/' procr ];
  if exist(path_out2) == 0
    eval(['mkdir ' path_out2]);
  end;
  outfile2 = [path_out2 '/' listfile '.mat'];
  save(outfile2,'r*','HA*in');
end;

if exist('save_thr')>0 && save_thr == 1 
  path_out3 = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/thresh_values/' procr ];
  if exist(path_out3) == 0
    eval(['mkdir ' path_out3]);
  end;
  outfile3 = [path_out3 '/' listfile '.mat'];
  save(outfile3,'thresh_values');
end;

load(outfile1_thr)
varlist = who('*_thresh');
for jjj = 1:length(varlist)
  eval([varlist{jjj}(1:end-7) '=' varlist{jjj} ';']);
end;
clear *_thresh;
save(outfile1_thr,'m_*','n_*','sd_*','HA*in','a*','s*','thresh_values');
