octave_packages/signal-1.1.3/private/inv_residue.m

   1 ## Copyright (C) 2007 R.G.H. Eschauzier <reschauzier@yahoo.com>
   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 ## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com>
  17
  18 ## This function is necessary for impinvar and invimpinvar of the signal package
  19
  20 ## Inverse of Octave residue function
  21 function [b_out, a_out] = inv_residue(r_in, p_in, k_in, tol)
  22
  23   n = length(r_in); % Number of poles/residues
  24
  25   k = 0; % Capture contstant term
  26   if (length(k_in)==1)    % A single direct term (order N = order D)
  27     k = k_in(1);          % Capture constant term
  28   elseif (length(k_in)>1) % Greater than one means non-physical system
  29     error("Order numerator > order denominator");
  30   endif
  31
  32   a_out = poly(p_in);
  33
  34   b_out  = zeros(1,n+1);
  35   b_out += k*a_out; % Constant term: add k times denominator to numerator
  36   i=1;
  37   while (i<=n)
  38     term   = [1 -p_in(i)];               % Term to be factored out
  39     p      = r_in(i)*deconv(a_out,term); % Residue times resulting polynomial
  40     p      = prepad(p, n+1, 0, 2);       % Pad for proper length
  41     b_out += p;
  42
  43     m          = 1;
  44     mterm      = term;
  45     first_pole = p_in(i);
  46     while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
  47        i++; %Next residue
  48        m++;
  49        mterm  = conv(mterm, term);              % Next multiplicity to be factored out
  50        p      = r_in(i) * deconv(a_out, mterm); % Resulting polynomial
  51        p      = prepad(p, n+1, 0, 2);           % Pad for proper length
  52        b_out += p;
  53     endwhile
  54   i++;
  55   endwhile
  56 endfunction