% Function media2
% Calcula media e mediana, com time stamp no INÍCIO do período
% Admite NaNs na matriz de entrada
% Written by Luciana Rizzo
%
% input: data_ (time | values to be averaged)        
% output: stat (time | var1_avg| var1_std| var1_median| var1_p25| var1_p75| var1_max| var1_min| var1_#| ... same for var2, etc )
%         statorder (time | avgs(var1,var2,...) | std(var1,var2,...) |  medians(var1,var2,...) | p25(var1,var2,...) | p75(var1,var2,...) | max(var1,var2,...) | min(var1,var2,...) | #(var1,var2,...) |        
% tolerancia para o final de cada período em dias (ex: 5seg = 5/60/60/24) 
% (se tol<0, diminui o período de média; se tol>0, aumenta o período de média)  
%
function [stat,statorder] = fmedia2(data_,first,last,delta,tol)
tempo_=first:delta:last;    % average time line
tempo_=tempo_';
nvar=size(data_,2)-1;
stat=nan(length(tempo_),8*nvar+1);
stat(:,1)=tempo_;
%
for i=1:length(tempo_)
    if(mod(i,1000)==0)
        disp([num2str(100*i/length(tempo_)),'%'])
    end
    a=find(data_(:,1)>=(tempo_(i)+tol) & data_(:,1)<(tempo_(i)+delta+tol));
    if (~isempty(a))
        for j=1:nvar
            stat(i,8*j-6)=nanmean(data_(a,j+1)); %avg
            stat(i,8*j-5)=nanstd(data_(a,j+1)); %std
            stat(i,8*j-4)=prctile(data_(a,j+1),50); %median
            stat(i,8*j-3)=prctile(data_(a,j+1),25); %p25
            stat(i,8*j-2)=prctile(data_(a,j+1),75); %p75
            stat(i,8*j-1)=max(data_(a,j+1)); %max
            stat(i,8*j)=min(data_(a,j+1)); %min
            stat(i,8*j+1)=sum(~isnan(data_(a,j+1))); %num
        end
    end
end
%
statorder(:,1)=stat(:,1); % time
jj=2;
for j=2:8:nvar*8-6 % averages
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
%jj=nvar+2;
for j=3:8:nvar*8-5 % std
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
%jj=2*nvar+2;
for j=4:8:nvar*8-4 % medians
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
for j=5:8:nvar*8-3 % p25
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
for j=6:8:nvar*8-2 % p75
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
for j=7:8:nvar*8-1 % max
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
for j=8:8:nvar*8-0 % min
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
for j=9:8:nvar*8+1 % num
    statorder(:,jj)=stat(:,j);
    jj=jj+1;
end
%
clear tempo_ last first delta a i j data_ nvar jj