octave_packages/signal-1.1.3/residuez.m

   1 %% Copyright (C) 2005 Julius O. Smith III <jos@ccrma.stanford.edu>
   2 %%
   3 %% This program is free software; you can redistribute it and/or modify it under
   4 %% the terms of the GNU General Public License as published by the Free Software
   5 %% Foundation; either version 3 of the License, or (at your option) any later
   6 %% version.
   7 %%
   8 %% This program is distributed in the hope that it will be useful, but WITHOUT
   9 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 %% details.
  12 %%
  13 %% You should have received a copy of the GNU General Public License along with
  14 %% this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {[@var{r}, @var{p}, @var{f}, @var{m}] =} residuez (@var{B}, @var{A})
  18 %% Compute the partial fraction expansion of filter @math{H(z) = B(z)/A(z)}.
  19 %%
  20 %% INPUTS:
  21 %% @var{B} and @var{A} are vectors specifying the digital filter @math{H(z) = B(z)/A(z)}.
  22 %% Say @code{help filter} for documentation of the @var{B} and @var{A}
  23 %% filter coefficients.
  24 %%
  25 %% RETURNED:
  26 %%   @itemize
  27 %%   @item @var{r} = column vector containing the filter-pole residues@*
  28 %%   @item @var{p} = column vector containing the filter poles@*
  29 %%   @item @var{f} = row vector containing the FIR part, if any@*
  30 %%   @item @var{m} = column vector of pole multiplicities
  31 %%   @end itemize
  32 %%
  33 %% EXAMPLES:
  34 %% @example
  35 %%   Say @code{test residuez verbose} to see a number of examples.
  36 %% @end example
  37 %%
  38 %% For the theory of operation, see
  39 %% @indicateurl{http://ccrma.stanford.edu/~jos/filters/residuez.html}
  40 %%
  41 %% @seealso{residue residued}
  42 %% @end deftypefn
  43
  44 function [r, p, f, m] = residuez(B, A, tol)
  45   % RESIDUEZ - return residues, poles, and FIR part of B(z)/A(z)
  46   %
  47   % Let nb = length(b), na = length(a), and N=na-1 = no. of poles.
  48   % If nb<na, then f will be empty, and the returned filter is
  49   %
  50   %             r(1)                      r(N)
  51   % H(z) = ----------------  + ... + ----------------- = R(z)
  52   %        [ 1-p(1)/z ]^m(1)         [ 1-p(N)/z ]^m(N)
  53   %
  54   % If, on the other hand, nb >= na, the FIR part f will not be empty.
  55   % Let M = nb-na+1 = order of f = length(f)-1). Then the returned filter is
  56   %
  57   % H(z) = f(1) + f(2)/z + f(3)/z^2 + ... + f(M+1)/z^M + R(z)
  58   %
  59   % where R(z) is the parallel one-pole filter bank defined above.
  60   % Note, in particular, that the impulse-response of the one-pole
  61   % filter bank is in parallel with that of the the FIR part.  This can
  62   % be wasteful when matching the initial impulse response is important,
  63   % since F(z) can already match the first N terms of the impulse
  64   % response. To obtain a decomposition in which the impulse response of
  65   % the IIR part R(z) starts after that of the FIR part F(z), use RESIDUED.
  66   %
  67   % J.O. Smith, 9/19/05
  68
  69   if nargin==3
  70     warning("tolerance ignored");
  71   end
  72   NUM = B(:)'; DEN = A(:)';
  73   % Matlab's residue does not return m (since it is implied by p):
  74   [r,p,f,m]=residue(conj(fliplr(NUM)),conj(fliplr(DEN)));
  75   p = 1 ./ p;
  76   r = r .* ((-p) .^m);
  77   if f, f = conj(fliplr(f)); end
  78 end
  79
  80 %!test
  81 %! B=[1 -2 1]; A=[1 -1];
  82 %! [r,p,f,m] = residuez(B,A);
  83 %! assert(r,0,100*eps);
  84 %! assert(p,1,100*eps);
  85 %! assert(f,[1 -1],100*eps);
  86 %! assert(m,1,100*eps);
  87
  88 %!test
  89 %! B=1; A=[1 -1j];
  90 %! [r,p,f,m] = residuez(B,A);
  91 %! assert(r,1,100*eps);
  92 %! assert(p,1j,100*eps);
  93 %! assert(f,[],100*eps);
  94 %! assert(m,1,100*eps);
  95
  96 %!test
  97 %! B=1; A=[1 -1 .25];
  98 %! [r,p,f,m] = residuez(B,A);
  99 %! [rs,is] = sort(r);
 100 %! assert(rs,[0;1],1e-7);
 101 %! assert(p(is),[0.5;0.5],1e-8);
 102 %! assert(f,[],100*eps);
 103 %! assert(m(is),[1;2],100*eps);
 104
 105 %!test
 106 %! B=1; A=[1 -0.75 .125];
 107 %! [r,p,f,m] = residuez(B,A);
 108 %! [rs,is] = sort(r);
 109 %! assert(rs,[-1;2],100*eps);
 110 %! assert(p(is),[0.25;0.5],100*eps);
 111 %! assert(f,[],100*eps);
 112 %! assert(m(is),[1;1],100*eps);
 113
 114 %!test
 115 %! B=[1,6,2]; A=[1,-2,1];
 116 %! [r,p,f,m] = residuez(B,A);
 117 %! [rs,is] = sort(r);
 118 %! assert(rs,[-10;9],1e-7);
 119 %! assert(p(is),[1;1],1e-8);
 120 %! assert(f,[2],100*eps);
 121 %! assert(m(is),[1;2],100*eps);
 122
 123 %!test
 124 %! B=[6,2]; A=[1,-2,1];
 125 %! [r,p,f,m] = residuez(B,A);
 126 %! [rs,is] = sort(r);
 127 %! assert(rs,[-2;8],1e-7);
 128 %! assert(p(is),[1;1],1e-8);
 129 %! assert(f,[],100*eps);
 130 %! assert(m(is),[1;2],100*eps);
 131
 132 %!test
 133 %! B=[1,6,6,2]; A=[1,-2,1];
 134 %! [r,p,f,m] = residuez(B,A);
 135 %! [rs,is] = sort(r);
 136 %! assert(rs,[-24;15],2e-7);
 137 %! assert(p(is),[1;1],1e-8);
 138 %! assert(f,[10,2],100*eps);
 139 %! assert(m(is),[1;2],100*eps);
 140
 141 %!test
 142 %! B=[1,6,6,2]; A=[1,-(2+j),(1+2j),-j];
 143 %! [r,p,f,m] = residuez(B,A);
 144 %! [rs,is] = sort(r);
 145 %! assert(rs,[-2+2.5j;7.5+7.5j;-4.5-12j],1E-6);
 146 %! assert(p(is),[1j;1;1],1E-6);
 147 %! assert(f,-2j,1E-6);
 148 %! assert(m(is),[1;2;1],1E-6);
 149
 150 %!test
 151 %! B=[1,0,1]; A=[1,0,0,0,0,-1];
 152 %! [r,p,f,m] = residuez(B,A);
 153 %! [as,is] = sort(angle(p));
 154 %! rise = [ ...
 155 %!  0.26180339887499 - 0.19021130325903i; ...
 156 %!  0.03819660112501 + 0.11755705045849i; ...
 157 %!  0.4; ...
 158 %!  0.03819660112501 - 0.11755705045849i; ...
 159 %!  0.26180339887499 + 0.19021130325903i;];
 160 %! pise = [ ...
 161 %! -0.80901699437495 - 0.58778525229247i; ...
 162 %!  0.30901699437495 - 0.95105651629515i; ...
 163 %!  1; ...
 164 %!  0.30901699437495 + 0.95105651629515i; ...
 165 %! -0.80901699437495 + 0.58778525229247i];
 166 %! assert(r(is),rise,100*eps);
 167 %! assert(p(is),pise,100*eps);
 168 %! assert(f,[],100*eps);
 169 %! assert(m,[1;1;1;1;1],100*eps);