function emg = naive_process_emg( raw_emg, thresh, baseline_interval, suppression )
% function emg = naive_process_emg( raw_emg, thresh, baseline_interval, suppression )
%
% Processes raw emg data into on/off activation data.
%
% This version uses the naive threshold approach.  In computer science terms,
% naive isn't really derogatory, it just means the first thing you might have thought of.
%
% Parameters:
%  raw_emg           - The raw emg data
%  thresh            - Threshold to use, 0 means calculate one based on the data
%  baseline_interval - How many samples at the beginning to consider for
%                      determining "baseline activation". Defualt: 400
%  suppression       - Size of gaps in activation to fill. Default: 0
%
% Nigel Stepp <stepp@atistar.net>
%


if size(raw_emg, 1) == 1
	raw_emg = raw_emg';
end;

if baseline_interval == 0
	baseline_interval = 400;
end;

do_plots = false;

if do_plots
	figure;
	plot(raw_emg/max(raw_emg), 'b');
	hold on;
end;

% recenter, rectify, and normalize
emg = raw_emg - mean(raw_emg);
emg = abs(emg);
emg = emg / max(emg);

if do_plots
	plot(emg, 'm');
end;

if thresh == 0
	% examine initial noise to determine baseline activation
	basestd = std(emg(1:baseline_interval));
	basemean = mean(emg(1:baseline_interval));

	thresh = basemean + 5*basestd
end;

% Filter
[B,A] = butter(2, 50/500);
emg = filtfilt(B, A, emg);

if do_plots
	plot(emg, 'k');

	plot(repmat([thresh],1,length(emg)), 'r');
end;

emg(emg >= thresh) = 1;
emg(emg < thresh) = 0;

% coalesce activation blocks based on "suppression"
%
% Ok, this may look a little mysterious.  Here is what is going on.
% First we find all the points where emg is equal to 0, but we remove
% the first 'suppression' number of points at the beginning.
zero_inds = find(emg==0);
zero_inds = zero_inds(zero_inds > suppression);

% Now we traverse the zero points from the end to the beginning.
for i = length(zero_inds)-1:-1:1
	% If there are 1 points in this supppression sized chunk, set this point to 1
	% As we count backwards, all points in suppression sized chunks will be set to 1.
	if sum( emg( zero_inds(i)-suppression : zero_inds(i)-1 ) ) > 0 && emg(zero_inds(i)+1) == 1
		emg(zero_inds(i)) = 1;
	end;
end;

if do_plots
	plot(emg, 'g');
	hold off;
end;