function [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) %function [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) % % The polynomial Beta(z)=z^n + b_1*z^(n-1) + ... + b_n % and the real-valued parameter r are computed such that % % r*Beta(z)*Beta(1/z) = rhs % % An iterative Newton algorithm by V. Kucera, % Discrete Linear Control (1979) is used. % Convergence problems may occur for root clusters close % to the unit circle. % % Inputs: % 1. rhs - "the right hand side" % 2. limit - max number of iterations (default 300) % 3. epsilon - accuracy limit for termination (default 1e-10) % % Outputs: % 1. Beta - spectral factor, with zeros inside the unit circle % 2. r - real-valued scalar % 3. ierr - Error codes: % a) ierr=0, convergence has occured. % b) ierr=1, no convergence within limited iterations. % c) ierr=2, right hand side (rhs) not positive definite % % 4. L1_norm - the L1-norm of the error (r*Beta(z)*Beta(1/z) - rhs) % % % [Beta,r,ierr,L1_norm]=spefac(rhs,limit,epsilon) % % % M. Sternad (1995-08) if nargin < 2 ,limit=300;end if nargin < 3 ,epsilon=1E-10;end stop=0; ierr=0; n=(length(rhs)+1)/2; rr = rhs(n:length(rhs)); if rr(1) < max(abs(rr)),ierr=2;stop=1;end if rhs(n) > 0 A=rhs(n:2*n-1)/sqrt(rhs(n)); stop=0; else ierr=1; stop=1; end % start iteration i=0; while ~stop i=i+1; B=rhs(n:2*n-1); W=A; k=0; while k~=n-1 for j=1:n-k,rhs1(j)=A(n-k+1-j);end Al(k+1)=A(n-k)/rhs1(n-k); if abs(Al(k+1)) < 0 ierr=2; stop=1; else A(1:n-k-1)=A(1:n-k-1)-Al(k+1)*rhs1(1:n-k-1); E(n-k)=2*B(n-k)/A(1); if k0 k=k-1; for j=1:n-k,P(j)=E(n-k+1-j);end E(1:n-k)=E(1:n-k)-Al(k+1)*P(1:n-k); end A=(E+W)/2; if sum(A.^2)-B(1) < epsilon*B(1) ierr=0; stop=1; elseif i==limit ierr=1; stop=1; end end Beta = A/A(1); r = A(1)^2; L1_norm = norm(r*xcorr(Beta)-rhs,1);