function [proEMG, rectEMG,lowEMG] = processEMG_Outlier(rawEMG, AnalogRate, referenceEMG)
% Purpose:  From raw EMG data as input, the method performs the
%           filtering, rectifying, and normalizing operations needed to
%           obtain an EMG envelope.
% 
% USAGE:  [proEMG, rectEMG] = processEMG(rawEMG, AnalogRate)
%
% Input:   rawEMG are columns of EMG channels
%          AnalogRate is the sampling rate of the raw EMG data
%
% Output:  proEMG is column matrix of normalized EMG envelopes
%          rectEMG (optional) of just norm rectified EMG
%
% NOTE:    rawEMG signals must contain at least 12 data points
%          in order to apply a fourth order lowpass filter
%
% ASeth Nov-07, streamlined version of ASA, 9-05
% Stanford University

[nAnalogFrames, nc] = size(rawEMG);

%% Apply 4th order, 0-lag, Butterworth band-pass filter to raw EMG signal.
order = 4;
fs = AnalogRate;
if fs >= 1080
% 	cutoff = [20 400];                % default
    cutoff = [40 450];                % Grimmer et. al. 2022
% 	cutoff = [80 400];                % use when there is 60 Hz noise
elseif fs == 600                      % For Delaware EMG data
	cutoff = [11 222];                % default
	% cutoff = [44 222];              % use when there is 60 Hz noise
end
[b, a] = butter(order/2, cutoff/(0.5*fs));
bandEMG = filter(b, a, rawEMG, [], 1);  

%% Apply Second-order IIR notch filter
F0 = 50;   % Interference is at 60 Hz
BW = 6;    % Choose a bandwidth factor of 6Hz
[num1,den1] = iirnotch(F0/(fs/2),BW/(fs/2));
bandEMG = filter(num1, den1, bandEMG, [], 1); 


%% Rectify the filtered EMG signal.
rectEMG = abs(bandEMG);

%% Apply 4th order, 0-lag, Butterworth low-pass filter to rectified signal.
order = 4;
cutoff = 6;
[b, a] = butter(order, cutoff/(0.5*fs));
lowEMG = filter(b, a, rectEMG, [], 1);  


%% Normalize rectified and low-pass filtered EMG signals.
onesMat = ones(nAnalogFrames,1);
nMinSamples = round(0.01 * nAnalogFrames); % average 1% of samples to get "min".

% Output 1 (Not used in our work)
sortRect = sort(rectEMG);           
minRect = mean(sortRect(1:nMinSamples,:));
normRect = (rectEMG - onesMat*minRect)./(onesMat*(max(rectEMG) - minRect));

% Output 2 (We use this output)
sortLow = sort(lowEMG);
minLow = mean(sortLow(1:nMinSamples,:));
sortRef = sort(referenceEMG);
maxLow = mean(sortRef(end-20:end,:));
% minLow = min(lowEMG);
% normLow = (lowEMG - onesMat*minLow)./ (max(referenceEMG) - minLow); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
normLow = (lowEMG - minLow)./ (maxLow - minLow);
% normLow = (lowEMG - onesMat*minLow)./ (max(referenceEMG));

% Return processed EMG data.        
rectEMG = normRect;
proEMG = normLow;                

