home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / GRPDELAY.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.3 KB  |  43 lines
  1. function [gd,w] = grpdelay(b,a,n,dum)
  2. %GRPDELAY Group delay of a digital filter.
  3. %    [Gd,W] = GRPDELAY(B,A,N) returns length N vectors Gd and W
  4. %    containing the group delay and the frequencies at which it 
  5. %    is evaluated. Group delay is -d{angle(w)}/dw.  The frequency
  6. %    response is evaluated at N points equally spaced around the
  7. %    upper half of the unit circle.   When N is a power of two,
  8. %    the computation is done faster using FFTs.
  9. %    GRPDELAY(B,A,N,'whole') uses N points around the whole unit
  10. %    circle.  See also FREQZ.
  11.  
  12. %    Unpublished algorithm from J. O. Smith 5-9-88
  13. %    L. Shure, 5-9-88
  14. %    (c) Copyright 1988, by The MathWorks, Inc.
  15.     
  16. nb = max(size(b));
  17. na = max(size(a));
  18. nn = max(size(n));
  19. s = 2;
  20. if nargin == 4
  21.     s = 1;
  22. end
  23. if nn == 1
  24.     w = (2*pi/s*(0:n-1)/n)';
  25. else
  26.     w = n;
  27.     n = nn;
  28. end
  29. c = conv(b, conj(a(na:-1:1)));
  30. c = c(:).';    % make a row vector
  31. nc = max(size(c));
  32. cr = c.*(0:(nc-1));
  33. if s*n >= nc    % pad with zeros to get the n values needed
  34.     gd = real(fft([cr zeros(1,s*n-nc)]) ./ fft([c zeros(1,s*n-nc)]));
  35.     gd = gd(1:n) - ones(1,n)*(na-1);
  36. else    % find multiple of s*n points greater than nc
  37.     nfact = s*ceil(nc/(s*n));
  38.     mmax = n*nfact;
  39.     gd = real(fft(cr,mmax) ./ fft(c,mmax));
  40.     gd = gd(1:nfact:mmax) - ones(1,n)*(na-1);
  41. end
  42. gd = gd(:);
  43.