2 % UCP(C) tests if the polynomial C is a Unit-Circle-Polynomial.
3 % It tests if all roots lie inside the unit circle like
4 % B=ucp(C) does the same as
5 % B=all(abs(roots(C))<1) but much faster.
6 % The Algorithm is based on the Jury-Scheme.
7 % C are the elements of the Polynomial
8 % C(1)*X^N + ... + C(N)*X + C(N+1).
11 % O. Foellinger "Lineare Abtastsysteme", Oldenburg Verlag, Muenchen, 1986.
12 % F. Gausch "Systemtechnik", Textbook, University of Technology Graz, 1993.
15 % $Id: ucp.m 5090 2008-06-05 08:12:04Z schloegl $
16 % Copyright (C) 1996-1999,2008 by Alois Schloegl <a.schloegl@ieee.org>
18 % This program is free software: you can redistribute it and/or modify
19 % it under the terms of the GNU General Public License as published by
20 % the Free Software Foundation, either version 3 of the License, or
21 % (at your option) any later version.
23 % This program is distributed in the hope that it will be useful,
24 % but WITHOUT ANY WARRANTY; without even the implied warranty of
25 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 % GNU General Public License for more details.
28 % You should have received a copy of the GNU General Public License
29 % along with this program. If not, see <http://www.gnu.org/licenses/>.
37 lambda = c(:,lc)./c(:,1);
38 % disp([lc,size(lambda), sum(b),toc]);
39 % ratio must be less then 1
40 b = b & (abs(lambda) < 1);
41 % and reduced polynomial must be a UCP, too.
42 c(:,1:lc-1) = c(:,1:lc-1) - lambda(:,ones(1,lc-1)).*c(:,lc:-1:2);