function norm_struct = define_norm(path1,chans,EMGleg,fpcut,save_it)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% define_norm.m - computes and returns normalization constants
% derived from both MVC trials and the fastest walking trials.
% These normalization constants are then applied to the 
% estimated muscle stimulations before they are filtered
% to estimate activation.
%
% Inputs:
% path1   - subject initials and date
% chans   - vector of EMG channels to be processed
% EMGleg  - L or R; refers to leg with EMG electrodes
% fpcut   - upper bounds for force plate noise floor (fastest trials)
% save_it - 1 to save norm_struct, 0 otherwise.
%
% Outputs:
% norm_struct - structure containing normalization constants for
% all EMG channels.  Negative values indicate normalization came
% from the fastest walking trials while positive values indicate
% they came from the MVC trials.
%
% Jared Markowitz
% Fall 2012
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute normalization factor based on fastest walking speed.
% This factor is the standard deviation of the data set (truncated
% to include an integer number of gait cycles).  Also process
% (i.e. recover envelope) to compare with MVC estimate:

% addpath '~/Dropbox/harvard/processEMG/preprocess'

infiles = {'walk200_1';'walk200_2';'walk200_3';'walk200_4';'walk200_5';'walk200_6';'walk200_7';'walk200_8'};

for k = 1:length(chans)
  eval(['totalEMG.chan' num2str(chans(k)) ' = [];']);
  eval(['n_chan' num2str(chans(k)) ' = [];']); 
end;

for i = 1:length(infiles)
  try
    infile1 = infiles{i};
    EMGdat  = loadEMG(infile1,path1,1,chans,EMGleg,fpcut);
    Fs      = EMGdat.Fs; 
    
    eval(['ind1 = EMGdat.' upper(EMGleg) 'hs(1);']);
    eval(['ind2 = EMGdat.' upper(EMGleg) 'hs(end);']);
    eval(['heels = EMGdat.' upper(EMGleg) 'hs;']);
  
    for j = 1:length(chans)      
      % Low pass signal and break into gait cycles for averaging
      eval(['dat1 = EMGdat.chan' num2str(chans(j)) ';']);
      
      % Remove mean
      dat2 = dat1 - mean(dat1);
      
      % Remove 60 Hz noise (notch filter)
      f0    = 59/(Fs/2);
      f1    = 61/(Fs/2);
      [b,a] = butter(2,[f0 f1],'stop');  
      dat3  = filtfilt(b,a,dat2);
    
      % Subtract out moving average to remove motion artifact:
      smooth1 = round(Fs/20);
      if mod(smooth1,2) > 0
        smooth1 = smooth1 - 1;
      end;  
      remdat = smooth(dat3,smooth1);
      dat4  = dat3 - remdat'; 

      % Preserve full gait cycles for sigma comparison
      eval(['totalEMG.chan' num2str(chans(j)) '=cat(2,totalEMG.chan' num2str(chans(j)) ',dat4(ind1:ind2));']);
      
      % Rectify:
      dat5  = abs(dat4); 
  
      % Smooth:
      smooth2 = round(0.1*Fs);
      if mod(smooth2,2) > 0
        smooth2 = smooth2-1;
      end;
      dat6 = smooth(dat5,smooth2);          
      %dat7 = dat6-min(dat6(100:end));
      dat8 = max(dat6,0);
      
      for m = 1:(length(heels)-1)
        indstart = heels(m);
        indstop  = heels(m+1);
        x1 = 0:1/Fs:(indstop-indstart)/Fs;
        x2 = linspace(0,x1(end),101);
        eval(['n_chan' num2str(chans(j)) '=[n_chan' num2str(chans(j)) '; interp1(x1,dat8(indstart:indstop),x2,''spline'',''extrap'')];']);     
      end;
    end;    
  catch
    fprintf(['No ' infiles{i} '\n']);
  end; 
end;  

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

% Compare with MVC trials and choose normalization method based on
% which is stronger.

eval(['load ''~/Documents/MATLAB/data/' path1 '/processed/EMG/' path1 '_calib_new.mat'';']);

for l = 1:length(chans)
  eval(['fast_norm.chan' num2str(chans(l)) '=std(totalEMG.chan' num2str(chans(l)) ',0);']);       
  eval(['m_chan' num2str(chans(l)) ' = mean(n_chan' num2str(chans(l)) ',1);']);
  eval(['sd_chan' num2str(chans(l)) ' = std(n_chan' num2str(chans(l)) ',0,1);']);
  eval(['fast_max.chan' num2str(chans(l)) '=max(m_chan' num2str(chans(l)) '+2*sd_chan' num2str(chans(l)) ');']);
  eval(['quant1 = fast_max.chan' num2str(chans(l)) ';']);
  eval(['quant2 = calib.chan' num2str(chans(l)) '.FMCval;']);

  %if quant1 >= quant2 
    %% Use fastest walking trials for normalization
    %% "fast_norm" for sigma normalization, fast_max for smoothed normalization
    %  eval(['norm_struct.mvc.chan' num2str(chans(l)) ' = -fast_max.chan' num2str(chans(l)) ';']);
  %else
    % Use MVC for normalization
    eval(['norm_struct.mvc.chan' num2str(chans(l)) ' = quant2;']);
  %end; 
  
  eval(['norm_struct.fw.chan' num2str(chans(l)) ' = fast_max.chan' num2str(chans(l)) ';']);
    
  clear quant1 quant2;
end;

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

% Save norm_struct:
if exist('save_it') > 0 && save_it == 1 
  outfile = ['~/Documents/MATLAB/data/' path1 '/processed/EMG/' path1 '_norm.mat'];
  save(outfile,'norm_struct');
end;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  