function [ITL,ITU,IZCT,energy,zc,N1,N2]=epdStats(signal,fs)
%[ITL,ITU,IZCT,ENERGY,ZC]=EPDSTATS(SIGNAL,FS)  
%  Calculates the energy, zero crossing rate and statistics for determining  
%  the thresholds used for isolated word endpoint detection as described 
%  by Rabiner and Sambur [1].
%
%  [1] Rabiner,L.R. and Sambur,M.R., "An Algorithm for Determining the Endpoints of
%      Isolated Utterances". The Bell System Technical Journal, 1975.
%
%  Author: Stuart N Wrigley       
%  MAD - Matlab Auditory Demonstrations
%  (c) University of Sheffield 1998
%  Revision 0.02: 27 July 1998

winSizeMs=10;   % should be 10ms
winShiftMs=10;  % should be 10ms
silenceEndMs=100;
backoffMs=250;

winSize=millitosamples(winSizeMs,fs);
winShift=millitosamples(winShiftMs,fs);
silenceEnd=millitosamples(silenceEndMs,fs);
mywindow=window(winSize,'hamming');

signedSignal=sgn(signal);
j=1;
for i=1:winShift:length(signal)-winSize
   energy(j)=(sum((abs(signal(i:i+winSize-1))).*mywindow));%/winSize;
   zc(j)=sum(abs(signedSignal(i+1:i+winSize)-signedSignal(i:i+winSize-1)));
   j=j+1;
end
silenceRange=1:length(1:winShift:silenceEnd-winSize);

% Zero Crossing Threshold
% fixed threshold is 25 zero crossings per 10ms if signal has fs of 10kHz
% 1. assume can linearly extrapolate this to given fs.
% 2. assume IZC is average zc for initial 100ms of signal. [1] vague on this issue.
IF=(((25*winSizeMs)/10)/10000)*fs;
IZC=mean(zc(silenceRange));
zcstd=std(zc(silenceRange));
IZCT=min(IF,IZC+2*zcstd);

% Energy Thresholds
% 1. assume IMN is average energy for initial 100ms of signal. [1] vague on this issue.
IMX=max(energy);
IMN=mean(energy(silenceRange));
I1=0.03*(IMX-IMN)+IMN;
I2=4*IMN;
ITL=min(I1,I2);
ITU=5*ITL;


N1=0;
N2=0;
%function [N1,N2]=getEndPoints(ITL,ITU,IZCT,energy,zc)
duration=length(energy);
backoffLength=length(1:winShift:millitosamples(backoffMs,fs)-winSize);

done=0;
% start estimation
for m=1:duration
   if and(energy(m)>=ITL,~done)
      for i=m:duration
         if energy(i)<ITL
            break
         else   
            if energy(i)>=ITU
               if ~done
                  N1=i-(i==m);
                  done=1;
               end
               break
            end
         end
      end
   end   
end
startID=max(N1-backoffLength,1);
endID=N1;
M1=sum(zc(startID:endID)>=IZCT);
if M1>=3
   for i=startID:endID
      if zc(i)>=IZCT
         N1=i;
         break;
      end      
   end
end

done=0;
% end estimation
for m=duration:-1:1
   if and(energy(m)>=ITL,~done)
      for i=m:-1:1
         if energy(i)<ITL
            break;
         else
            if energy(i)>=ITU
               if ~done
                  N2=i+(i==m);
                  done=1;
               end
               break
            end
         end
      end   
   end
end
startID=N2;
endID=min(N2+backoffLength,length(zc));
M2=sum(zc(startID:endID)>=IZCT);
if M2>=3
   for i=max([1 startID]):endID   % max added by MPC to prevent neg indices
      if zc(i)>=IZCT
         N2=i;
         break;
      end      
   end
end

warpRatio=round(length(signal)/length(energy));
N1=N1*warpRatio;
N2=N2*warpRatio;