gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\grpdelay.m
function [gd_out,w_out] = grpdelay(b,a,n,dum,Fs) %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 (in radians) 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. If you don't specify % N, it defaults to 512. % % GRPDELAY(B,A,N,'whole') uses N points around the whole unit circle. % % [Gd,F] = GRPDELAY(B,A,N,Fs) and [Gd,F] = GRPDELAY(B,A,N,'whole',Fs) % given sampling frequency Fs in Hz return a vector F in Hz. % % Gd = GRPDELAY(B,A,W) and Gd = GRPDELAY(B,A,F,Fs) return the group delay % evaluated at the points in W (in radians) or F (in Hz), where Fs is the % sampling frequency in Hz. % % GRPDELAY(B,A,...) with no output arguments plots the group delay versus % normalized frequency (Nyquist == 1) in the current figure window. % % See also FREQZ. % Unpublished algorithm from J. O. Smith, 5-9-88 % Author(s): L. Shure, 5-9-88 % T. Krauss, 4-19-93, revised % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.2 $ $Date: 1998/08/18 18:09:04 $ error(nargchk(1,5,nargin)) if nargin == 1, a = 1; n = 512; whole = 'no'; samprateflag = 'no'; elseif nargin == 2, n = 512; whole = 'no'; samprateflag = 'no'; elseif nargin == 3, whole = 'no'; samprateflag = 'no'; elseif nargin == 4, if isstr(dum), whole = 'yes'; samprateflag = 'no'; else whole = 'no'; samprateflag = 'yes'; Fs = dum; end elseif nargin == 5, whole = 'yes'; samprateflag = 'yes'; end nb = length(b); na = length(a); nn = length(n); if strcmp(whole,'yes'), s = 1; else s = 2; end c = conv(b, conj(a(na:-1:1))); c = c(:).'; % make a row vector nc = length(c); cr = c.*(0:(nc-1)); if length(n)==1, w = (2*pi/s*(0:n-1)/n)'; if s*n >= nc % pad with zeros to get the n values needed % dividenowarn temporarily supresses warnings to avoid "Divide by zero" gd = dividenowarn(fft([cr zeros(1,s*n-nc)]),... fft([c zeros(1,s*n-nc)])); gd = real(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; % dividenowarn temporarily supresses warnings to avoid "Divide by zero" gd = dividenowarn(fft(cr,mmax), fft(c,mmax)); gd = real(gd(1:nfact:mmax)) - ones(1,n)*(na-1); end gd = gd(:); else if strcmp(samprateflag,'no'), w = n; else w = 2*pi*n/Fs; end s = exp(j*w); gd = real(polyval(cr,s)./polyval(c,s)); gd = gd - ones(size(gd))*(na-1); end if strcmp(samprateflag,'yes') f = w*Fs/2/pi; else f = w; end if nargout == 0 % do plots newplot; if strcmp(samprateflag,'no'), plot(f/pi,gd) xlabel('Normalized frequency (Nyquist == 1)') else plot(f,gd) xlabel('Frequency (Hertz)') end ylabel('Group delay (in samples)') set(gca,'xgrid','on','ygrid','on') elseif nargout == 1 gd_out = gd; elseif nargout == 2 gd_out = gd; w_out = f; end