Home > freetb4matlab > signal > grpdelay.m

grpdelay

PURPOSE ^

% Compute the group delay of a filter.

SYNOPSIS ^

function [gd,w] = grpdelay(b,a=1,nfft=512,whole,Fs)

DESCRIPTION ^

% Compute the group delay of a filter.
%
% [g, w] = grpdelay(b)
%   returns the group delay g of the FIR filter with coefficients b.
%   The response is evaluated at 512 angular frequencies between 0 and
%   pi. w is a vector containing the 512 frequencies.
%   The group delay is in units of samples.  It can be converted
%   to seconds by multiplying by the sampling period (or dividing by
%   the sampling rate fs).
%
% [g, w] = grpdelay(b,a)
%   returns the group delay of the rational IIR filter whose numerator
%   has coefficients b and denominator coefficients a.
%
% [g, w] = grpdelay(b,a,n)
%   returns the group delay evaluated at n angular frequencies.  For fastest
%   computation n should factor into a small number of small primes.
%
% [g, w] = grpdelay(b,a,n,'whole')
%   evaluates the group delay at n frequencies between 0 and 2*pi.
%
% [g, f] = grpdelay(b,a,n,Fs)
%   evaluates the group delay at n frequencies between 0 and Fs/2.
%
% [g, f] = grpdelay(b,a,n,'whole',Fs)
%   evaluates the group delay at n frequencies between 0 and Fs.
%
% [g, w] = grpdelay(b,a,w)
%   evaluates the group delay at frequencies w (radians per sample).
%
% [g, f] = grpdelay(b,a,f,Fs)
%   evaluates the group delay at frequencies f (in Hz).
%
% grpdelay(...)
%   plots the group delay vs. frequency.
%
% If the denominator of the computation becomes too small, the group delay
% is set to zero.  (The group delay approaches infinity when
% there are poles or zeros very close to the unit circle in the z plane.)
%
% Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}],  is the rate of change of
% phase with respect to frequency.  It can be computed as:
%
%               d/dw H(e^-jw)
%        g(w) = -------------
%                 H(e^-jw)
%
% where
%         H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
%
% By the quotient rule,
%                    A(z) d/dw B(z) - B(z) d/dw A(z)
%        d/dw H(z) = -------------------------------
%                               A(z) A(z)
% Substituting into the expression above yields:
%                A dB - B dA 
%        g(w) =  ----------- = dB/B - dA/A
%                    A B
%
% Note that,
%        d/dw B(e^-jw) = sum(k b_k e^-jwk)
%        d/dw A(e^-jw) = sum(k a_k e^-jwk)
% which is just the FFT of the coefficients multiplied by a ramp.
%
% As a further optimization when nfft>>length(a), the IIR filter (b,a)
% is converted to the FIR filter conv(b,fliplr(conj(a))).
% For further details, see 
% http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:
Generated on Sat 16-May-2009 00:04:49 by m2html © 2003