function [sqrtrho]=envlp(P,Ahat,bhat,theta,noiselevel)

% function [sqrtho,theta]=envlp(P,A,B,theta,noiselevel) : spectral envelop
%
% Computes the envelope of spectra consistent with the data
% by using state covariance statistics P for the input-to-state system (A,B).
%
% 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 (default: NOISELEVEL=0).
%
% Values of NOISELEVEL > 0 tend to increase resolution between spectral lines.
%
% System matrices need to be normalized so that:
% A*A'+B*B' = I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=eye(size(Ahat));
W=commutant(P,Ahat,I);
P=W+W';

if nargin==5, P=P-min(eig(P))*noiselevel; end
if nargin==3, theta=linspace(0,2*pi,200); end
 
Pinvhalf=inv(sqrtm(P));

exptheta=exp(j*theta);
sv=[];
for i=1:length(theta),
sv=[sv; norm( Pinvhalf*inv( I-conj(exptheta(i))*Ahat ) *bhat ) ];
end

sqrtrho=1./sv;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%