% 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

Generated on Sat 16-May-2009 00:04:49 by