function EMGsmooth = e_smooth(EMGdet, Window, varargin)
%E_SMOOTH Compute moving average/smoothing filter of detected EMG channels.
%   EMGsmooth = e_smooth(EMGdet, Window, Options)
%             or
%   EMGsmooth = e_smooth(EMGdet, Order, Wn, LinFilt [, Causality])
%
%   EMGsmooth = e_smooth(EMGdet, Window, Options) applies a mean-absolute-
%   value (MAV) filter to multiple-channel detected EMG.  Matrix EMGdet
%   contains vectors of detected EMG (the longer matrix dimension denotes
%   row vs. column data).  Window specifies the windowing values.  If
%   window is a scalor, that window value is applied as a
%   fixed window.  Else, Window must have the same length as
%   an EMG vector, and specifies the window length to use at
%   each sample location.
%
%   Processing option combinations are:
%   1) 'MAV': Specifies that MAV processing is to be performed by
%      e_smooth.  This method is the default method and the default
%      interpretation of the function arguments.
%   2) 'Noncausal' or 'Causal' or h:  h specifies the number
%      of samples to shift the smoothing window forward in
%      time.  h is specified either as a scalor or vector the
%      same length as an EMG vector, similar to Window.  Note
%      that 'Noncausal' is equivolent to h=(Window-1)/2, for
%      Window an odd value.  'Causal' processing is equivolent
%      to h=0.
%   3) 'EdgesBest' or 'EdgesFixed' or 'EdgesZero' or 'EdgesNaN'
%   4) 'SumFast' or 'SumFastLong' or 'SumFilter' or 'SumFull'
%      (Note: 'SumFast', 'SumFastLong' and 'SumFilter' must use
%      a fixed-length smoothing window with h also fixed.)
%      (Note: 'SumFastLong' is not yet implimented.)
%
%   EMGsmooth = e_smooth(EMGdet, Order, Wn, LinFilt [, Causality]) filters
%   the input EMG data with a low-pass linear filter.  EMGdet is specified as
%   described above.  Order is the scalor filter order.  LinFilt
%   specifies the type of linear filter and must be one of (see the
%   Signal Processing Toolbox) "butter" (Butterworth), "ellip"
%   (Elliptical, with 0.5 decibels of ripple in the passband and a
%   stopband 20 decibels down), "cheby1" (Chebyshev type I, with 0.5
%   decibels of ripple in the passband), or "cheby2" (Chebyshev type II,
%   with the stopband ripple 20 decibels down).  Optional argument
%   Causality must be one of "Causal" or "Noncausal".  If "Causal" is
%   selected, then the data are filtered with "filter()".  If
%   "Noncausal" is selected, then the data are filtered with
%   "filtfilt()", thus a zero-phase filter of twice the filter order
%   results.  Noncausal filtering is the default.
%
%   Null-valued options are ignored.  Null specification of the processing
%   method results in 'MAV' processing and interpretation of the
%   function arguments.
%
%   Output EMGsmooth will have the dimensions of a
%   detected EMG waveform vector.
%
%   EMG Amplitude Estimation Toolbox - Ted Clancy - WPI

% Copyright (c) Edward A. Clancy, 2007.
% This work is licensed under the Aladdin free public license.
% For copying permissions see license.txt.
% email: ted@wpi.edu

% 5 February 2007.

%%%%%%%%%%%%%%%%%%%%% Process First Argument %%%%%%%%%%%%%%%%%%%%%

% Check number of input arguments
if (nargin < 2  |  nargin > 6)
  error(['Bogus argument count (' int2str(nargin) '), expected 2-6.']);
end

% Get detected EMG waveform vectors and form composite COLUMN vector.
if (length(size(EMGdet)) > 2), error('EMGdet is not a vector or 2D-matrix.'); end
if (min(min(EMGdet)) < 0.0),   error('EMGdet has a negative valued entry.');  end
[Erow, Ecol] = size(EMGdet);
if (Erow == Ecol)
  error(['EMGdet has same number of rows as columns (' int2str(Erow) ').']);
end
if ( ((Erow-1)*(Ecol-1))  ==  0 )  % Vector.
  if (Erow == 1), EMGcomp = EMGdet';
  else,           EMGcomp = EMGdet;  end
else  % Matrix.
  if (Erow > Ecol), EMGcomp = mean(EMGdet')';
  else,             EMGcomp = mean(EMGdet)';
  end
end

%%%%%%%%%%%%%%% Determine Smoothing Method, Call Subfunction %%%%%%%%%%%%%%%

FlagMethod = 0;
for i=1:length(varargin)
  if ischar(varargin{i})
    switch varargin{i}
      case 'MAV'  % Default.
        if FlagMethod==1, error('Repeated method selection.'); end
        FlagMethod = 1;
      case {'butter', 'ellip', 'cheby1', 'cheby2'}
        if FlagMethod==1, error('Repeated method selection.'); end
          switch nargin
            case 4, EMGsmooth = MethodLinFilt(EMGcomp, Window, ...
                                varargin{1}, varargin{2}, 'Noncausal');
            case 5, EMGsmooth = MethodLinFilt(EMGcomp, Window, ...
                                varargin{1}, varargin{2}, varargin{3});
            otherwise, error('Bogus arg count for linear filter usage.');
          end
       % Return vector of desired dimension.
       if (Ecol > Erow), EMGsmooth = EMGsmooth'; end
       return;  %%%%%%%%%% NORMAL LinFilt RETURN.  %%%%%%%%%%
    end  % switch varargin{i}.
  end
end

% 'MAV' method selected (or defaulted).
EMGsmooth = MethodMAV(EMGcomp, Window, varargin);
% Return vector of desired dimension.
if (Ecol > Erow), EMGsmooth = EMGsmooth'; end
return;  %%%%%%%%%% NORMAL MAV RETURN.  %%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MAV Method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function EMGsmooth = MethodMAV(EMGcomp, Window, Option)

%%%%%%%%%%%%%%%%%%%% Continue Processing Arguments %%%%%%%%%%%%%%%%%%%%

% Prepare Window COLUMN vector.
N = length(EMGcomp);
if (length(size(Window)) > 2), error('Window is not a scalor or vector.'); end
if (min(Window) <= 0), error('Window has a zero or negative valued entry.');  end
if (length(Window) == 1)  % Scalor context.
  Window = Window * ones(N,1);
else  % Vector context.
  [Wrow, Wcol] = size(Window);
  if ~( (Wrow==N & Wcol==1) | (Wrow==1 & Wcol==N) )
    error(['Window has bogus dimensions (' int2str(Wrow) 'x' int2str(Wcol) ').']);
  end
  if (Wcol > Wrow) Window = Window'; end  % Cast to column vector.
end

% Optional arguments.
% Set defaults and flags.
h = fix(Window/2);  % Default for Noncausal processing.
DefCausal = 'Noncausal';
DefEdges  = 'EdgesBest';
DefSum    = 'SumFast';
FlagCausal = 0;    % Zero denotes option in this catagory
FlagEdges  = 0;    %   has not yet been specified.  One
FlagSum    = 0;    %   denotes otherwise.

% Process the options.
for i=1:length(Option)
  if isstr( Option{i} )
    switch Option{i}
      case 'MAV'  % Don't need to do anything.  (Checked redundancy above.)
      case 'Noncausal'
        if (FlagCausal > 0), error('Repeated causality selection.'); end
          FlagCausal = 1;  % h default set for this case above.
      case 'Causal'
        if (FlagCausal > 0), error('Repeated causality selection.'); end
        FlagCausal = 1;
        DefCausal = 'Causal';
        h = 0;
      case {'EdgesBest', 'EdgesFixed', 'EdgesZero', 'EdgesNaN'}
        if (FlagEdges > 0), error('Repeated edge selection.'); end
        FlagEdges = 1;
        DefEdges = Option{i};
      case {'SumFast', 'SumFastLong', 'SumFilter', 'SumFull'}
        if (FlagSum > 0 ), error('Repeated sum selection.'); end
        FlagSum = 1;
        DefSum = Option{i};
      otherwise
        error(['Bogus option "' Option{i} '".']);
    end
  elseif isempty(Option{i})  % Do nothing for null arguments.
  else  % Must be h.
    if (FlagCausal > 0 ), error('Repeated causality selection.'); end
    FlagCausal = 1;
    DefCausal = 'h';
    if (length(size(Option{i})) > 2)
      error('h is not a scalor or vector.');
    end
    if (length(Option{i}) == 1)  % Scalor context.
      h = Option{i} * ones(N,1);
    else  % Vector context.
      [hRow, hCol] = size(Option{i});
      if ~( (hRow==N & hCol==1) | (hRow==1 & hCol==N) )
        error(['h has bogus dimensions (' int2str(hRow) 'x' int2str(hCol) ').']);
      end
      h = Option{i};
      if (hCol > hRow) h = h'; end  % Cast to column vector.
    end
  end  % Must be h.
end  % for.

% Be sure that the window lengths are appropriate for the selected causality.
switch DefCausal
  case 'Noncausal',     Window = ( fix(Window/2) .* 2) + 1;  % Nearest odd integer.
  case {'Causal', 'h'}, Window = round(Window);              % Nearest integer.
end
  
% Be sure that all h values are whole numbers and in range.
h = round(h);
if( any(h<0) ),       error('All "h" values must be >= zero.');    end
if( any(h>=Window) ), error('All "h" values must be < "Window".'); end

% If fast summing is selected, be sure Window and h are fixed.
if ( strcmp(DefSum, 'SumFast') | strcmp(DefSum, 'SumFastLong') | strcmp(DefSum, 'SumFilter') )
  if ( max(h) ~= min(h) ), error('"h" must be fixed for fast summing options.'); end
  if ( max(Window) ~= min(Window) ), error('"Window" must be fixed for fast summing options.'); end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute %%%%%%%%%%%%%%%%%%%%%%%%%

% Create the summing indexes, then adjust Window for edge processing.
Istart = [1:N]' - Window + h + 1;
Iend   = Istart + Window - 1;
BadEdge = find( Istart<=0 | Istart>N | Iend<=0 | Iend>N );
Istart = max(  min(Istart, N), 1  );
Iend   = max(  min(Iend,   N), 1  );
FixWin = Window(1);  % If window is fixed, this saves its length.
FixH   = h(1);       % If h      is fixed, this saves its length.
clear h;  % Not used again.

switch DefEdges
  case 'EdgesBest',  Window = Iend - Istart + 1;
  case 'EdgesFixed'  % Do nothing.
  case 'EdgesZero',  Window(BadEdge) = Inf;
  case 'EdgesNaN',   Window(BadEdge) = NaN;
end
clear BadEdge;  % Not used again.

% Sum in "switch" statement, then divide by "Window" to smooth.
if (strcmp(DefSum, 'SumFull')==0), clear Iend Istart; end  % Only used by 'SumFull'.
EMGsmooth = ones(N,1);  % Pre-allocate.
switch DefSum
  case 'SumFast'
    if( exist('e__smo') ~= 3 ), error('Fast summing mex file not found.'); end
    e__smo(EMGsmooth, EMGcomp, FixWin, FixH, 1);  % Fast mex sum.
  case 'SumFastLong'
    error('SumFastLong not programmed yet.')
    if( exist('e__smo') ~= 3 ), error('Fast summing mex file not found.'); end
    e__smo(EMGsmooth, EMGcomp, FixWin, FixH, 2);  % Fast mex sum.
  case 'SumFilter'
    EMGcomp = [EMGcomp; zeros(FixH,1)];  % Pad for causality.
    EMGsmooth = filter(ones(FixWin,1), 1, EMGcomp);
    EMGsmooth = EMGsmooth(FixH+1:N+FixH);  % Extract with correct causality.
  case 'SumFull'
    for i = 1:N, EMGsmooth(i) = sum( EMGcomp(Istart(i) : Iend(i)) ); end
end
EMGsmooth = EMGsmooth ./ Window;

EMGsmooth = max(EMGsmooth, 0);  % Prevent round-off from causing a negative value.

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% LinFilt Methods %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function EMGsmooth = MethodLinFilt(EMGin, Order, Wn, LinFilt, Causality)

%%%%%%%%%%%%%%%%%%%% Continue Processing Arguments %%%%%%%%%%%%%%%%%%%%

% Order.
if ischar(Order), error('Order must be numeric.'); end
if length(Order)>1, error('Order must be a scalor.'); end
if Order <=0, error('Order must be a positive integer.'); end
% Wn.
if ischar(Wn) == 1, error('Argument 2 must be numeric.'); end
if length(Wn)>1, error('Argument 2 must be a scalor.'); end
if (Wn<=0 | Wn>=1), error('Argument 2 must be 0<Wn<1.'); end
% LinFilt.
if ischar(LinFilt)==0, error('LinFilt method must be a string.'); end
% Causality.
if ischar(Causality)==0, error('LinFilt causality must be a string.'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Filter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch LinFilt
  case 'butter', [b, a] = butter(Order, Wn);
  case 'ellip',  [b, a] = ellip(Order, 0.5, 20, Wn);
  case 'cheby1', [b, a] = cheby1(Order, 0.5, Wn);
  case 'cheby2', [b, a] = cheby2(Order, 0.5, Wn);
  otherwise, error(['Internal Error: Bogus method "' LinFilt '".']);
end

switch Causality
  case 'Causal',    EMGsmooth = filter(b, a, EMGin);
  case 'Noncausal', EMGsmooth = filtfilt(b, a, EMGin);
  otherwise, error(['Bogus causality "' Causality '".']);
end

EMGsmooth = max(EMGsmooth, 0);  % Prevent round-off from causing a negative value.

return;