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  <lang@dsp.rice.edu>%%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 angleind_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 circleif isodd(n_orig_u);     disp('There are roots on unit circle with odd multiplicity!!');   returnend% check whether all roots could be sortedif n_orig_u+n_orig_i+n_orig_o ~= n;   disp('Not all roots could be sorted');   returnend% find roots on the unit circle (differentiated pol.) and sort by angleind_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 togetherif 2*n_dif_u ~= n_orig_u   disp('The numbers of roots on unit circle do not fit together!!');   returnend% check whether maximum or minimum phase polynomial is desiredr_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];   endend   % compute  output polynomial[dum,ind] = sort(angle(r_out)); r_out = r_out(ind);  p_out = poly(r_out(ransort(length(r_out))));