
function [sqrtrho,theta]=env(Ahat,bhat,y,noiselevel)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function [sqrtho,theta]=env(y,Ahat,bhat,noiselevel)
%
% Computes the envelope of spectra consistent with the data
% by using state covariance statistics for the system SYS.
%
% If 0 < NOISELEVEL <1 is assigned a value, then a level of
% white noise at the input is assumed which correspond to the
% specified ``noiselevel'' fraction of the smallest sv of the
% the state covariance.
%
% System matrices need to be normalized so that:
% Ahat*Ahat'+bhat*bhat' = Identity
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ny=length(y);
dim=length(Ahat); ell=1;

y=y(:); y=y-mean(y);
   [x ,t]=rlsim(Ahat,bhat,y);              x=x(ell:Ny,:); 
   Pf=cov(conj(x));
   [xr,t]=rlsim(Ahat,bhat,conj(flipud(y))); xr=xr(ell:Ny,:);
   Pr=cov(conj(xr));
   Paverage=(Pf+Pr)/2;

E=eye(size(Ahat));
W=commutant(Paverage,Ahat,E);
P=W*E+E*W';

if nargin==4, P=P-min(eig(P))*noiselevel; end
 
Pinvhalf=inv(sqrtm(P));

theta=linspace(0,2*pi,1024);
exptheta=exp(j*theta);
I=eye(size(Ahat));
sv=[];
for i=1:length(theta),
sv=[sv; norm( Pinvhalf*inv( I-conj(exptheta(i))*Ahat ) *bhat ) ];
end

sqrtrho=1./sv;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tryphon Georgiou  (revised May 8, 2000)
