/usr/share/octave/packages/signal-1.3.0/private/inv_residue.m is in octave-signal 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 | ## Copyright (C) 2007 R.G.H. Eschauzier <reschauzier@yahoo.com>
##
## 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 <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{b_out}, @var{a_out}] =} inv_residue (@var{r_in}, @var{p_in}, @var{k_in}, @var{tol})
## Undocumented internal function. This function is used by the impinvar
## and invimpinvar functions in the signal package.
## @end deftypefn
## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com>
## Inverse of Octave residue function
function [b_out, a_out] = inv_residue(r_in, p_in, k_in, tol)
n = length(r_in); # Number of poles/residues
k = 0; # Capture contstant term
if (length(k_in)==1) # A single direct term (order N = order D)
k = k_in(1); # Capture constant term
elseif (length(k_in)>1) # Greater than one means non-physical system
error("Order numerator > order denominator");
endif
a_out = poly(p_in);
b_out = zeros(1,n+1);
b_out += k*a_out; # Constant term: add k times denominator to numerator
i=1;
while (i<=n)
term = [1 -p_in(i)]; # Term to be factored out
p = r_in(i)*deconv(a_out,term); # Residue times resulting polynomial
p = prepad(p, n+1, 0, 2); # Pad for proper length
b_out += p;
m = 1;
mterm = term;
first_pole = p_in(i);
while (i<n && abs(first_pole-p_in(i+1))<tol) # Multiple poles at p(i)
i++; # Next residue
m++;
mterm = conv(mterm, term); # Next multiplicity to be factored out
p = r_in(i) * deconv(a_out, mterm); # Resulting polynomial
p = prepad(p, n+1, 0, 2); # Pad for proper length
b_out += p;
endwhile
i++;
endwhile
endfunction
|