function varargout = grpdelay(f,varargin)
%GRPDELAY  DFILT object group delay.
%   [Gd,W] = GRPDELAY(G) computes the group delay (in samples) of the
%   DFILT object G, and returns it in vector Gd.  The vector of
%   frequencies at which Gd is evaluated is returned in vector W.
%
%   Group delay is given by -d(arg(H(w)))/dw where arg(H(w)) is the 
%   phase response of the filter G.
%
%   [Gd,W] = GRPDELAY(G,N) uses N frequency points to compute the
%   group delay.  If you don't specify N, it defaults to 512.  For
%   FIR filters, it is recommended to use a value of N that is a
%   power of two to perform the computation faster.
%
%   GRPDELAY(G,N,'whole') uses N points around the whole unit circle.  
%
%   [Gd,F] = GRPDELAY(G,N,Fs) and [Gd,F] = GRPDELAY(G,N,'whole',Fs)
%   uses the sampling frequency Fs (in Hz) for the computation and
%   returns a vector of frequencies F, also in Hz.
%
%   Gd = GRPDELAY(G,W) and Gd = GRPDELAY(G,F,Fs), where W or F are vectors,
%   return the group delay evaluated at the points in W (in radians/sample)
%   or F (in Hz).
%
%   GRPDELAY(...) with no output arguments plots the group delay.
%
%   See also DFILT/IMPZ, DFILT/FREQZ, DFILT/ZPLANE, DFILT/STEPZ.

%   Author(s): R. Losada
%   Copyright 1988-2000 The MathWorks, Inc.
%   $Revision: 1.4 $  $Date: 2000/06/13 21:21:00 $ 

[b,a] = tf(f);

[Gd,w] = grpdelay(b,a,varargin{:});
normalizedfreq = 1;

if (nargin == 4) | (nargin == 3 & isnumeric(varargin{2})),
    % Fs was specified,
    normalizedfreq = 0;
end

if nargout == 0,
    grpdelayplot(Gd,w,normalizedfreq)       
else
    varargout = {Gd,w};
end


%-------------------------------------------------------------------------
function grpdelayplot(Gd,w,normalizedfreq)

xlbl = 'Frequency (Hz)';
ylbl = 'Group delay (samples)';


if normalizedfreq,
    xlbl = 'Normalized Frequency (\times\pi rad/sample)';
    w = w./pi;    
end

newplot,
plot(w,Gd);
xlabel(xlbl);
ylabel(ylbl);
set(gca,'xgrid','on','ygrid','on');
