function delay = grdelays(b,a,w)
% The function delay = grdelays(b,a,w)
% evaluates the group delay of an analog filter provided
% the numerator and denominator coefficients in vectors b and a
% of the transfer function are given
%                            nb-1      nb-2
% 		 b(s)   b(1)s  +  b(2)s   + ... +  b(nb)
% 	  H(s) = ---- = --------------------------------
% 	                     na-1      na-2
% 		 a(s)   a(1)s  +  a(2)s   + ... +  a(na) 
% 
% w is the frequency range
% Miroslav Vlcek, December 1996 
na=length(a);
nb=length(b);
ac=a;
n=fix(na/2);
for i=1:n,
    ac(na-2*i+1)=-a(na-2*i+1);
end
bc=b;
m=fix(nb/2);
for i=1:m,
    bc(nb-2*i+1)=-b(nb-2*i+1);
end
bac=conv(b,ac);
abc=conv(a,bc);
odd=(abc-bac)/2 ;% pocita lichou cast citatele prenosove funkce
even=(abc+bac)/2; % pocita sudou cast citatele prenosove funkce

% vypocet faze prenosove funkce H(s)=b(s)/a(s) pro s=j*omega
s=sqrt(-1)*w;
re=polyval(even,s);
im=sqrt(-1)*polyval(odd,s);
ratio=im ./re;
phi=atan(ratio);
alpha=rozbal(phi,3.); % rozbaluje arctan do spojite funkce bez skoku 

ader=polyder(a);
acder=-polyder(ac);
bder=polyder(b);
bcder=-polyder(bc);


% vypocet skupinoveho zpozdeni prenosove funkce H(s)=b(s)/a(s)
cit1=polyval(ader,s);
jmen1=polyval(a,s);
cit2=polyval(acder,s);
jmen2=polyval(ac,s);
cit3=polyval(bder,s);
jmen3=polyval(b,s);
cit4=polyval(bcder,s);
jmen4=polyval(bc,s);

delay= (cit1 ./ jmen1 + cit2 ./ jmen2 - cit3 ./ jmen3 - cit4 ./ jmen4)/2;
