function [gd,w] = grpdelay(b,a,n,dum)
%GRPDELAY Group delay of a digital filter.
%	[Gd,W] = GRPDELAY(B,A,N) returns length N vectors Gd and W
%	containing the group delay and the frequencies at which it 
%	is evaluated. Group delay is -d{angle(w)}/dw.  The frequency
%	response is evaluated at N points equally spaced around the
%	upper half of the unit circle.   When N is a power of two,
%	the computation is done faster using FFTs.
%	GRPDELAY(B,A,N,'whole') uses N points around the whole unit
%	circle.  See also FREQZ.

%	Unpublished algorithm from J. O. Smith 5-9-88
%	L. Shure, 5-9-88
%	(c) Copyright 1988, by The MathWorks, Inc.
	
nb = max(size(b));
na = max(size(a));
nn = max(size(n));
s = 2;
if nargin == 4
	s = 1;
end
if nn == 1
	w = (2*pi/s*(0:n-1)/n)';
else
	w = n;
	n = nn;
end
c = conv(b, conj(a(na:-1:1)));
c = c(:).';	% make a row vector
nc = max(size(c));
cr = c.*(0:(nc-1));
if s*n >= nc	% pad with zeros to get the n values needed
	gd = real(fft([cr zeros(1,s*n-nc)]) ./ fft([c zeros(1,s*n-nc)]));
	gd = gd(1:n) - ones(1,n)*(na-1);
else	% find multiple of s*n points greater than nc
	nfact = s*ceil(nc/(s*n));
	mmax = n*nfact;
	gd = real(fft(cr,mmax) ./ fft(c,mmax));
	gd = gd(1:nfact:mmax) - ones(1,n)*(na-1);
end
gd = gd(:);
