function H = unitarh(G,F,lim)
% H = unitarh(G,F,lim)
% Unitarity relation ; calculation of the numerator H of the reflexion
% coefficient of the s parameter matrix of an analog filter
% from the polynome G, the denominator of all s parameters,
% and the polynome F, the numerator of the transmission coefficient,
% using H*hurwitz(H) = G*hurwitz(G) - F*hurwitz(F) ;
% the roots with least (normally negative) real parts, in H*hurwitz(H),
% are retained for building H.
% The coefficients of H which are either less than sqrt(eps) or less
% than lim if 3 arguments are introduced, are set to zero.
% P. MURET 26-nov-2002
if prod(size(G)) ~= length(G) | prod(size(F)) ~= length(F)
error('the arguments must be vectors.');
end
if length(G) < length(F)
error('the degree of G must not be less than that of F');
end
% calculation of the square hurwitzian moduli:
Godd = hurwitz(G);
Fodd = hurwitz(F);
G2 = conv(G,Godd);
F2 = conv(F,Fodd);
% the length of the second polynome is fitted to that of the first one
% by adding zeros at the beginning, and the unitary relation is done
ng2=length(G2);
nf2=length(F2);
Z=zeros(1,ng2-nf2);
DH2=G2 - [Z F2];
% roots of the resulting polynome are extracted and sorted out,
% beginning with that with the most positive real part,
% in decreasing order. Then the first half is suppressed,
% leaving only the roots with negative real parts.
RH2=roots(DH2);
degh=length(RH2)/2;
RnH2=esort(RH2);
RnH2(1:degh)=[];
H=poly(RnH2);
% setting the too small coefficients of H to zero
for k = 1:degh
if (abs(H(k))-sqrt(eps)) < 0, H(k)=0; end
if nargin == 3
if (abs(H(k))- lim) <0, H(k)=0; end
end
end