msort.m

ADSP TOOLBOX: Version 2.0 and gui m-files

function [z, m] = msort(r,tol)
% [z, m] = msort(r)
% Sets up roots z (from given roots r) and their multiplicity m
% Real roots returned first (sorted)
% Imaf roots sorted by increasing imaginary part.
% tol =tolerence for repeated roots. Default: tol = 0.001
%
% Support file for lpsiggui

% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



if nargin<2,tol=0.001;end

% Remove real and imag parts less than tol
%rr=real(r);ri=imag(r);
%r=rr.*(abs(rr)>tol)+j*ri.*(abs(ri)>tol);

%Rounding.
rr=(rem(real(r),1));
rr=rr.*(abs(rr)>tol & (1-abs(rr))>tol)+sign(rr).*((1-abs(rr))<tol);
ri=(rem(imag(r),1));
ri=ri.*(abs(ri)>tol & (1-abs(ri))>tol)+sign(ri).*((1-abs(ri))<tol);
r=fix(real(r))+rr+j*(fix(imag(r))+ri);


%Rearrange by real roots first followed by only type a+jb if complex
r=flipud(cplxpair(r(:)));
l=length(r);

k=find(imag(r)==0);
r1=[];r2=[];r3=[];
if ~isempty(k),
rr1=r(k);r(k)=[];l=length(r);
[p,s]=sort(rr1);
r1=rr1(s);
end
if ~isempty(r)
r=r(1:2:l);
end
k=find(real(r)==0);
if ~isempty(k),
rr2=r(k);r(k)=[];l=length(r);
[p,s]=sort(imag(rr2));
r2=rr2(s);
end
if ~isempty(r)
[p,s]=sort(imag(r));
r3=r(s);
end
y=[r1(:);r2(:);r3(:)];


%find multiplicity
z(1)=y(1);m(1)=1;
i=1;y(1)=[];while ~isempty(y)
t=z(i)-y(1);
if abs(real(t))<tol & abs(imag(t))<tol
m(i)=m(i)+1;
else
i=i+1;z(i)=y(1);m(i)=1;
end
y(1)=[];
end
z=z(:);m=m(:);