tonicEMG(varargin)

function varargout = tonicEMG(varargin)

% returns the tonic EMG activity or level
%
% [EMGL]=tonicEMG(TIME,EMG,INTERVAL[,MOMENT][,statistic])
%
% TIME: column vector of time in milliseconds
% EMG: signal in any unit
% INTERVAL:  interval in milliseconds to compute EMG level
%            if no MOMENT is specified, the whole signal is divided in
%            INTERVAL blocks and for each the statistic is computed
% MOMENT (optional): colum vector of time points in milliseconds for tonic EMG measurement
%                    EMG level is computed for INTERVAL milliseconds after MOMENT  
% statistic (optional): statistical parameter (format: string) to compute for INTERVAL
%                       Default: mean
%                       has to be a Matlab-defined function like median etc.
%
%************output*********************
% [EMGL]=tonicEMG(...)
% [EMGL] contains the following columns:
%       [starting_time stop_time mean(EMG)*]
%       *or median, if specified in PARAMETER
%************output options*********************
% [EMGL,varnames]=tonicEDA(...)
% additionally returns variable names as a cell array

% written by Robert Schleicher
% r *DOT* schleicher *AT* uni *MINUS* koeln *DOT* de

%***************handle input parameters**************
if nargin<3
    disp('not enough input arguments')
    return
end
time=varargin{1};
EMG=varargin{2};
interval=varargin{3}; %applies to case 4 and 5 
moment=[1:interval:time(end)]'; 
statistic=@mean; %default statistic

if nargin>3
   switch nargin 
        case 4
            if isnumeric(varargin{4}) %interval given...
              moment=varargin{4};
            else
              temp=varargin{4}; %could be statisticial function
            end            
        case 5
            moment=varargin{4};
            temp=varargin{5}; %could be statisticial function
        otherwise        %more
   end
   %check whether temp is a valid statistical function
   if exist('temp') && isstr(temp) %statistical parameter given
        if temp(1)=='@' %anonymous function given, e.g. '@(x)mean(x)*0.7' 
            eval(['temp=' temp ';']);           
        else %name of a function given
            temp=str2func(temp);            
        end
        if isa(temp,'function_handle');
            statistic=temp;        
        end
   end
   
end

flag=~isnan(EMG);
EMG=EMG(flag); %remove NaNs 

%***************filtering**************
%comment this out if you don't want to filter
EMG_raw=EMG;
[b,a] = butter(8,1/5,'high'); %worked for me to remove EOG interferences from orbicularis EMG
EMG = filtfilt(b,a,EMG_raw); 
EMG=abs(EMG); %rectify

%%**************************************
%uncomment the following two lines if you want to use the envelope of your signal
%make sure you downloaded Envelope1.1 from Matlab file exchange before
%( see http://www.mathworks.com/matlabcentral/fileexchange/ )
%and put the file envelope.m in your working directory

% [up,down]=envelope(time,EMGrec); 
% EMGrec=up;



ret=[]; %Output vector
for i=1:size(moment,1)
    flag=((time(:,1)>=moment(i)) & (time(:,1)<=(moment(i)+interval)));
        if sum(flag)<1
            disp(['No data for moment: ' num2str(moment(i))]);             
        else 
            ret=[ret;moment(i) (moment(i)+interval) statistic(EMG(flag))];
        end    
    
end

varargout{1}=ret;
if nargout>1
    varargout{2}={'start'; 'stop'; 'EMG_level'};
end
%plot results:

f1=figure; %create new figure
h1=plot(time/1000, EMG,'b'); %filtered and rectified signal
leg1=legend('filtered & rectified signal','Location','NorthWest');
hold;
a1=gca;
set(gca,'NextPlot','add');
lh=[];
a2=axes;
set(a2,'NextPlot','add');


for i=1:size(ret,1)
    if (ret(i,1)+interval)>time(end) %in case last moment + interval exceed time
        interval=time(end)-ret(i,1);
    end    
    rh=fill([ret(i,1) (ret(i,1)+interval) (ret(i,1)+interval) ret(i,1)]/1000, ...
            [0 0 ret(i,3) ret(i,3)],[0.6 0.6 0],'LineStyle',':','FaceAlpha',0.5);
    
end

linkaxes([a1 a2])
set(a2,'Color','none')
legend2=[func2str(statistic) ' integral'];
leg2=legend(legend2);
set(leg2,'Color','w')