V_ROOTSTAB determines the number of polynomial roots outside, inside and on the unit circle [NO,NI,NC]=v_rootstab(P)
Inputs: p Polynomial with real or complex coefficients
Outputs: no Number of roots outside the unit circle
ni Number of roots inside the unit circle
nc Number of roots lying on the unit circle
This routine uses the algorithm given in [1]. Rounding errors may cause roots that lie on the unit circle to migrate to either inside or outside.
Refs:
[1] Messaoud Benidir. On the root distribution of general polynomials with respect to the unit circle.
Signal Processing, 53 (1): 75–82, August 1996. ISSN 0165-1684. doi: 10.1016/0165-1684(96)00077-1.
[ Note error in equation (52) which should read k^ = +k - Q(0)/(conj(c0)*(1/c - c)) ]0001 function [no,ni,nc]=v_rootstab(p) 0002 %V_ROOTSTAB determines the number of polynomial roots outside, inside and on the unit circle [NO,NI,NC]=v_rootstab(P) 0003 % 0004 % Inputs: p Polynomial with real or complex coefficients 0005 % 0006 % Outputs: no Number of roots outside the unit circle 0007 % ni Number of roots inside the unit circle 0008 % nc Number of roots lying on the unit circle 0009 % 0010 % This routine uses the algorithm given in [1]. Rounding errors may cause roots that lie on the unit circle to migrate to either inside or outside. 0011 % 0012 % Refs: 0013 % [1] Messaoud Benidir. On the root distribution of general polynomials with respect to the unit circle. 0014 % Signal Processing, 53 (1): 75–82, August 1996. ISSN 0165-1684. doi: 10.1016/0165-1684(96)00077-1. 0015 % [ Note error in equation (52) which should read k^ = +k - Q(0)/(conj(c0)*(1/c - c)) ] 0016 0017 % Copyright (C) Mike Brookes 2025 0018 % 0019 % VOICEBOX is a MATLAB toolbox for speech processing. 0020 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0021 % 0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 % This program is free software; you can redistribute it and/or modify 0024 % it under the terms of the GNU General Public License as published by 0025 % the Free Software Foundation; either version 2 of the License, or 0026 % (at your option) any later version. 0027 % 0028 % This program is distributed in the hope that it will be useful, 0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0031 % GNU General Public License for more details. 0032 % 0033 % You can obtain a copy of the GNU General Public License from 0034 % http://www.gnu.org/copyleft/gpl.html or by writing to 0035 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0037 0038 no=0; % initialize count of roots outside unit circle 0039 nc=0; % initialize count of roots on the unit circle 0040 if all(p==0) % p is all zero 0041 ni=0; 0042 else 0043 p=p(:).'; % force to be a row vector 0044 p=p(find(p~=0,1):end); % trim leading zeros 0045 np=length(p); % initial 1+order(p) 0046 np0=np; % save initial order for calculating ni later 0047 npd=0; % initialize saved order for calculating nc 0048 while length(p)>1 0049 p=p/sqrt(p*p'); % normalize p each loop 0050 pf=conj(p(np:-1:1)); % flipped version of p 0051 k=-p(np)/pf(np); 0052 q=p+k*pf; % null out the constant coefficient 0053 q(np)=0; % force exact zero in case of rounding errors 0054 if all(q==0) 0055 p=p(1:end-1).*(np-1:-1:1); % take derivative 0056 if npd==0 0057 npd=np; % save current 1+order(p) 0058 nod=no; % save current count of roots outside unit circle 0059 end 0060 elseif q(1)==0 % if |k|=1 and q~=0 0061 q=q(1:find(q~=0,1,'last')); % trim trailing zeros from q 0062 dr=-q(end)/(pf(end)*k); % direction we will move normalized k in 0063 if abs(real(dr))>abs(imag(dr)) 0064 cf=abs(real(dr)); % increase or decrease real part by 0.5 0065 else 0066 cf=0.25*abs(imag(dr)); % increase or decrease imag part by 2 0067 end 0068 c=sqrt(1+cf^2)-cf; % choose c in range (0,1) 0069 p=p/c+k*c*pf+[zeros(1,np-length(q)) q]; 0070 elseif abs(k)>1 0071 q=q(1:find(q~=0,1,'last')); % trim trailing zeros from q 0072 p=conj(q(end:-1:1)); % flip q 0073 no=no+np-length(p); % increement count of zeros outside unit circle 0074 else 0075 p=q(1:find(q~=0,1,'last')); % trim trailing zeros from q and take as p 0076 end 0077 np=length(p); % 1+order(p) 0078 end 0079 if npd>0 0080 nc=npd-1-2*(no-nod); 0081 end 0082 ni=np0-1-nc-no; 0083 end