%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% process_perry.m takes as input the digitized version of the 
% Perry stimulation profiles and runs them through the Raasch
% bilinear filter to obtain activation estimates.
%
% Jared Markowitz
% Fall 2012
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load '~/Documents/MATLAB/data/Perry_muscles/stim.mat'

% Pass through Raasch bilinear dynamics
Fs = 100/1.08;
tt = 0:1.08/100:3*1.08;

muscles = {'AL';'AM';'BFSH';'GAS';'GLU_low';'GLU_upp';'GMED';'ILL';'RF';'SMEM';'SOL';'TA';'VLAT';'BFLH'};

for i=1:length(muscles)
  
  simopts = simset('FixedStep',1/Fs,'SrcWorkspace','current');

  muscles1 = muscles{i};
  
  if strcmp(muscles{i},'SOL') == 1
    mfft = 0.20;
  elseif strcmp(muscles{i},'GAS') == 1
    mfft = 0.50;
  elseif strcmp(muscles{i},'TA') == 1
    mfft = 0.35;
  elseif strcmp(muscles{i},'BFSH') == 1    
    mfft = 0.35;
  elseif strcmp(muscles{i},'SMEM') == 1
    mfft = 0.35;
    muscles1 = 'HAM';
  elseif strcmp(muscles{i},'BFLH') == 1
    mfft = 0.35;
  elseif strcmp(muscles{i},'VLAT') == 1
    mfft = 0.50;
    muscles1 = 'VAS';
  elseif strcmp(muscles{i},'RF') == 1
    mfft = 0.65;
  elseif strcmp(muscles{i},'ILL') == 1    
    mfft = 0.50;
  elseif strcmp(muscles{i},'GLU_low') == 1
    mfft = 0.45;
    muscles1 = 'GLU';
  elseif strcmp(muscles{i},'GLU_upp') == 1
    mfft = 0.50;
    muscles1 = 'HAB1';
  elseif strcmp(muscles{i},'GMED') == 1
    mfft = 0.50;
    muscles1 = 'HAB2';
  elseif strcmp(muscles{i},'AL') == 1
    mfft = 0.35;
    muscles1 = 'ADDL';
  elseif strcmp(muscles{i},'AM') == 1
    mfft = 0.45;
    muscles1 = 'ADDM';
  else
    fprintf('Erroneous muscle label.\n');
    return;
  end;
    
  % Short
  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)));

  eval(['stim = cat(2,p_s' muscles{i} ',p_s' muscles{i} '(2:end),p_s' muscles{i} '(2:end));']);
  stim = max(stim,0.001);
  
  eval(['x_in = [tt'' stim''];']);
  sim('raasch.mdl',tt(end),simopts);
  eval(['act_sh.' muscles1 ' = a_out(end-100:end);']);
  
  figure
  plot(stim(end-100:end),'linewidth',2);
  hold on;
  plot(a_out(end-100:end),'r','linewidth',2);  
  
  % Long
  k1 = 1/(.08-.047*mfft)-1/(.09-.056*mfft);
  k2 = 1/(.09-.056*mfft);
  
  eval(['x_in = [tt'' stim''];']);
  sim('raasch.mdl',tt(end),simopts);
  eval(['act_lo.' muscles1 ' = a_out(end-100:end);']);
  
  plot(a_out(end-100:end),'g','linewidth',2);  
  grid on; axis tight;
  legend('stim','sh','lo','location','best');
  eval(['title(''' muscles1 ''');']);
  xlabel('Percent GC');
  ylabel('Stim/Act Estimate');
end;

save('~/Documents/MATLAB/data/Perry_muscles/stim_act.mat','p_*','act*');

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