function [tg,F] = iirgroupdelay(b,a)
% usage
%        [tg,F] = iirgroupdelay(b,a)
% computes the frequency response of the IIR filter
% 
%   H(z) = b(1) + b(2) z^{-1}  + ... + b(M) z^{-M}
%          ---------------------------------------
%          a(1) + a(2)z^{-1} + .... + a(N)z^{-N}
% 
% The output is the group delay for the frequency points
% F (in cycles/sample).  Group delay is defined as
%  
%    tg(F) = -1/(2*pi) * d phase(H(F))/dF 
%
% If no output arguments are used then the response is plotted.
%
% R Kakarala
% UC Berkeley Extension

if (nargin < 2) % no denominator
    a=1;
end;

N = max(length(b),length(a));
N2 = 2^(4+ceil(log2(N)));  % zero-pad by a factor of 16.  
N2 = max(N2,1024);         % make sure its large enough
numerator = fft(b,N2);  % compute response of numerator (moving avg)
denominator  = fft(a,N2); % response of denominator (feedback part)

denominator(abs(denominator)==0) = eps;   % avoid divide by zero by setting zero to eps, machine precision

H = numerator./denominator;
H=H(2:N2/2-1);  % take only half the freq response, due to conjugate symmetry
% get phase angle
phaH=angle(H);

% since phase is limited to -pi, pi,
% unwrap phase angle into a continuous function
unwrap_phaH = unwrap(phaH);

% compute group delay: this is the negative of the 
% derivative of unwrapped phase with respect to frequency
grd = -1*diff(unwrap_phaH);  
% diff is MATLAB's difference approximation to derivative
grd = grd/(2*pi);  %divide by 2pi

% since the fft in this case evaluates phase at frequences separated by
% 1/N2 cycles/sample, convert to units of samples by multiplying by N2;
grd_in_samples = grd * N2;

% now plot
F=linspace(0,0.5,length(grd_in_samples));

if (nargout == 0)
    figure;
    plot( F, grd_in_samples); 
    xlabel('freq (Cycles/sample)'); 
    ylabel('group delay (samples)');
    grid on;
else
    tg=grd_in_samples;
end;