1######################################################################## 2## 3## Copyright (C) 2007-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {[@var{multp}, @var{idxp}] =} mpoles (@var{p}) 28## @deftypefnx {} {[@var{multp}, @var{idxp}] =} mpoles (@var{p}, @var{tol}) 29## @deftypefnx {} {[@var{multp}, @var{idxp}] =} mpoles (@var{p}, @var{tol}, @var{reorder}) 30## Identify unique poles in @var{p} and their associated multiplicity. 31## 32## The output is ordered from pole with largest magnitude to smallest 33## magnitude. 34## 35## If the relative difference of two poles is less than @var{tol} then they are 36## considered to be multiples. The default value for @var{tol} is 0.001. 37## 38## If the optional parameter @var{reorder} is zero, poles are not sorted. 39## 40## The output @var{multp} is a vector specifying the multiplicity of the poles. 41## @code{@var{multp}(n)} refers to the multiplicity of the Nth pole 42## @code{@var{p}(@var{idxp}(n))}. 43## 44## For example: 45## 46## @example 47## @group 48## p = [2 3 1 1 2]; 49## [m, n] = mpoles (p) 50## @result{} m = [1; 1; 2; 1; 2] 51## @result{} n = [2; 5; 1; 4; 3] 52## @result{} p(n) = [3, 2, 2, 1, 1] 53## @end group 54## @end example 55## 56## @seealso{residue, poly, roots, conv, deconv} 57## @end deftypefn 58 59function [multp, indx] = mpoles (p, tol, reorder) 60 61 if (nargin < 1 || nargin > 3) 62 print_usage (); 63 endif 64 65 if (nargin < 2 || isempty (tol)) 66 tol = 0.001; 67 endif 68 69 if (nargin < 3 || isempty (reorder)) 70 reorder = true; 71 endif 72 73 Np = numel (p); 74 75 ## force poles to be a column vector 76 77 p = p(:); 78 79 if (reorder) 80 ## sort with largest magnitude first 81 [~, ordr] = sort (abs (p), "descend"); 82 p = p(ordr); 83 else 84 ordr = (1:Np).'; 85 endif 86 87 ## find pole multiplicity by comparing relative difference of poles 88 89 multp = zeros (Np, 1); 90 indx = []; 91 n = find (multp == 0, 1); 92 while (n) 93 dp = abs (p-p(n)); 94 if (p(n) == 0.0) 95 if (any (abs (p) > 0 & isfinite (p))) 96 p0 = mean (abs (p(abs (p) > 0 & isfinite (p)))); 97 else 98 p0 = 1; 99 endif 100 else 101 p0 = abs (p(n)); 102 endif 103 k = find (dp < tol * p0); 104 ## Poles can only be members of one multiplicity group. 105 if (numel (indx)) 106 k = k(! ismember (k, indx)); 107 endif 108 m = 1:numel (k); 109 multp(k) = m; 110 indx = [indx; k]; 111 n = find (multp == 0, 1); 112 endwhile 113 multp = multp(indx); 114 indx = ordr(indx); 115 116endfunction 117 118 119%!test 120%! [mp, n] = mpoles ([0 0], 0.01); 121%! assert (mp, [1; 2]); 122 123%!test 124%! [mp, n] = mpoles ([-1e4, -0.1, 0]); 125%! assert (mp, ones (3, 1)); 126%! assert (n, [1; 2; 3]); 127