factorize.m

Polynomial Root Finder is a reliable and fast C program (+ Matlab gateway) for finding all roots of

function [p_out] = factorize(p_in,string)
% function [p_out] = factorize(p_in,string)
%
%    Input:   p_in, string
%
%    Output:  p_out
%   
%    Description: Program yields (spectral) factorization of the polynomial 
%                 p_in regarding to the the unit circle. For reliable
%                 results p_in should have single roots off the unit circle
%                 and double roots on the unit circle. If the input variable
%                 string contains at least one 'a' the maximum phase portion
%                 is determined (taking the double roots of p_in on the
%                 unit circle once + all roots of p_in outside the
%                 unit circle). In all other cases the function yields the
%        	  minimum phase portion. 
%
%    subroutines: ransort.m bitrev.m rootsl.mex 
%
%

%File Name: factorize.m
%Last Modification Date: %G%	%U%
%Current Version: %M%	%I%
%File Creation Date: Tue Sep  7 09:29:06 1993
%Author: Markus Lang  <[e-mail]>
%
%Copyright: All software, documentation, and related files in this distribution
%           are Copyright (c) 1993  Rice University
%
%Permission is granted for use and non-profit distribution providing that this
%notice be clearly maintained. The right to distribute any portion for profit
%or as part of any commercial product is specifically reserved for the author.
%
%Change History:
%

p_in = p_in(:)';
n = length(p_in) - 1;

% find roots of original and differentiated polynomial
[r_orig,e] = rootsl(p_in);  r_orig = r_orig(:).';  error_orig = e;
[r_dif, error_dif] = rootsl((n:-1:1).*p_in(n+1:-1:2));  r_dif = r_dif(:).';
  
% find roots on, inside and outside the unit circle (original pol.)
% and sort by angle
ind_orig_u = find(abs(1-abs(r_orig))<10*error_orig);
ind_orig_i = find(abs(r_orig)<abs(1-10*error_orig));
ind_orig_o = find(abs(r_orig)>abs(1+10*error_orig));
r_orig_u = r_orig(ind_orig_u);  [dum,ind] = sort(angle(r_orig_u)); 
r_orig_u = r_orig_u(ind);  n_orig_u = length(r_orig_u);
r_orig_i = r_orig(ind_orig_i);  [dum,ind] = sort(angle(r_orig_i)); 
r_orig_i = r_orig_i(ind);  n_orig_i = length(r_orig_i);
r_orig_o = r_orig(ind_orig_o);  [dum,ind] = sort(angle(r_orig_o)); 
r_orig_o = r_orig_o(ind);  n_orig_o = length(r_orig_o);


% check multiplicity of roots on unit circle
if isodd(n_orig_u);  
   disp('There are roots on unit circle with odd multiplicity!!');
   return
end

% check whether all roots could be sorted
if n_orig_u+n_orig_i+n_orig_o ~= n;
   disp('Not all roots could be sorted');
   return
end

% find roots on the unit circle (differentiated pol.) and sort by angle
ind_dif_u = find(abs(1-abs(r_dif))<10*error_dif);
r_dif_u = r_dif(ind_dif_u);  [dum,ind] = sort(angle(r_dif_u)); 
r_dif_u = r_dif_u(ind);  n_dif_u = length(r_dif_u);

% check whether the numbers of roots on the unit circle of both polynomials
% fit together
if 2*n_dif_u ~= n_orig_u
   disp('The numbers of roots on unit circle do not fit together!!');
   return
end


% check whether maximum or minimum phase polynomial is desired
r_out = [r_dif_u r_orig_i];
if exist('string') == 1;
   if length(find(string=='a')) > 0;  
     r_out = [r_dif_u r_orig_o];
   end
end   

% compute  output polynomial
[dum,ind] = sort(angle(r_out)); r_out = r_out(ind);  
p_out = poly(r_out(ransort(length(r_out))));