X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fgrpdelay.m;fp=octave_packages%2Fsignal-1.1.3%2Fgrpdelay.m;h=5bddc5aa0fd51973e80e0be0109bbde6a1ae65f6;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/signal-1.1.3/grpdelay.m b/octave_packages/signal-1.1.3/grpdelay.m new file mode 100644 index 0000000..5bddc5a --- /dev/null +++ b/octave_packages/signal-1.1.3/grpdelay.m @@ -0,0 +1,321 @@ +## Copyright (C) 2000 Paul Kienzle +## Copyright (C) 2004 Julius O. Smith III +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + +## 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 + +function [gd,w] = grpdelay(b,a=1,nfft=512,whole,Fs) + + if (nargin<1 || nargin>5) + print_usage; + end + HzFlag=0; + if length(nfft)>1 + if nargin>4 + print_usage(); + elseif nargin>3 % grpdelay(B,A,F,Fs) + Fs = whole; + HzFlag=1; + else % grpdelay(B,A,W) + Fs = 1; + end + w = 2*pi*nfft/Fs; + nfft = length(w)*2; + whole = ''; + else + if nargin<5 + Fs=1; % return w in radians per sample + if nargin<4, whole=''; + elseif ~ischar(whole) + Fs = whole; + HzFlag=1; + whole = ''; + end + if nargin<3, nfft=512; end + if nargin<2, a=1; end + else + HzFlag=1; + end + + if isempty(nfft), nfft = 512; end + if ~strcmp(whole,'whole'), nfft = 2*nfft; end + w = Fs*[0:nfft-1]/nfft; + end + + if ~HzFlag, w = w * 2*pi; end + + oa = length(a)-1; % order of a(z) + if oa<0, a=1; oa=0; end % a can be [] + ob = length(b)-1; % order of b(z) + if ob<0, b=1; ob=0; end % b can be [] as well + oc = oa + ob; % order of c(z) + + c = conv(b,fliplr(conj(a))); % c(z) = b(z)*conj(a)(1/z)*z^(-oa) + cr = c.*[0:oc]; % cr(z) = derivative of c wrt 1/z + num = fft(cr,nfft); + den = fft(c,nfft); + minmag = 10*eps; + polebins = find(abs(den)