dl=ones(1,257);
ol=ones(1,256);
l=0:7;


for i=1:8,
if l(i)<=0
eta=0;
else
eta=1;
end
   for v=0:15,
   m=l(i)*16+v+1;                              % se descrie legea A
   dld1(1,m)=2^(l(i)-eta+1)*(v+16*eta); 
   old(1,m)=2^(l(i)-eta+1)*(v+16*eta+.5);
end
end
dv=[4096];
dld=[dld1 dv];
dl1=dld(1:1,2:129);
dl2=-1.*fliplr(dl1);    %         nivele de decizie
dl=[dl2 dld];
ol1=old(1:1,1:128);    % niv. de reconstructie
ol2=-1.*fliplr(ol1);
ol=[ol2 old];
s=-87.18:1:-57.18;     %nivele de intrare
s1=.05.*s;
ss=(2852.9/sqrt(2))*10.^s1;    %(ss) --tensiunea efectiva a zgomotului   2852-amplit. fundam.
x0=[0 .1]%x0=[0 .1 .2 .3 .4 .5];       % decalajul
for u=1:2
nmin=x0(u)-2.*ss;         %limitele de integrare
nmax=x0(u)+2.*ss;     
for j=1:length(ss)
   i(j)=1;
while(i(j)<=256&dl(1,i(j))<nmin(j))        %calc. ce trepte vor interveni (de la imin la imax
   i(j)=i(j)+1;
end
imin(j)=i(j)-1;
if imin(j)<1
   imin(j)=1;
end
if imin(j)>256
   imin(j)=256;
end
end
for j=1:length(ss)
   i(j)=257;
   while (i(j)>=1&dl(1,i(j))>nmax(j))
      i(j)=i(j)-1;
   end
   imax(j)=i(j);
   if imax(j)<1
      imax(j)=1;
   end
   if imax(j)>256
      imax(j)=256;
   end
end
q=length(ss);
ii=ones(q,256);      % un vector cu toate elementele 1 de lungime 256

for l=1:q
   for k=1:256
      ii(l,k)=gau(dl(1,k),dl(1,k+1),ss(1,l),x0(u),101); % calc. integrale pg. 1963.
      hh(l,k)=tau(dl(1,k),dl(1,k+1),ss(1,l),x0(u),101);
      gg(l,k)=var(dl(1,k),dl(1,k+1),ss(1,l),x0(u),101);	
end
end
e2=ol.^2;

for l=1:q
   for k=1:256
      jj(l,k)=ii(l,k)* e2(1,k);
      ff(l,k)=ii(l,k)*ol(1,k);
      tt(l,k)=hh(l,k)*ol(1,k);
   end
end
ey2=ones(1,q);ey=ones(1,q);
exy=ones(1,q);disp=ones(1,q);
nn=101;
for j=1:q
   rl=jj(j:j,imin(j)+1:imax(j)-1); rr(j)=sum(rl);
   al=ff(j:j,imin(j)+1:imax(j)-1); aa(j)=sum(al);
   bl=tt(j:j,imin(j)+1:imax(j)-1); bb(j)=sum(bl);
cl=hh(j:j,imin(j)+1:imax(j)-1); cc(j)=sum(cl);
xl=gg(j:j,imin(j)+1:imax(j)-1); dd(j)=sum(xl);
ey2(1,j)=ol(1,imin(1,j))^2*gau(nmin(1,j),dl(1,imin(1,j)+1),ss(1,j),x0(u),nn)+...
ol(1,imax(1,j))^2*gau(dl(1,imax(1,j)),nmax(1,j),ss(1,j),x0(u),nn)+rr(j);%calc. momentu' de ord2
ey(j)=ol(imin(j))*gau(nmin(j),dl(imin(j)+1),ss(j),x0(u),nn)+ol(imax(j))*gau(dl(imax(j)),nmax(j),ss(j),x0(u), nn)+aa(j);
disp(j)=ey2(j)-ey(j)^2;
%exy(j)=ol(imin(j))*tau(nmin(j),dl(imin(j)+1),ss(j),x0(u),nn)+ol(imax(j))*tau(dl(imax(j)),nmax(j),ss(j),x0(u), nn)+bb(j);
ex(j)=tau(nmin(j),dl(imin(j)+1),ss(j),x0(u),nn)+tau(dl(imax(j)),nmax(j),ss(j),x0(u), nn)+cc(j);
ex2(j)=var(nmin(j),dl(imin(j)+1),ss(j),x0(u),nn)+var(dl(imax(j)),nmax(j),ss(j),x0(u), nn)+dd(j);
dispx(j)=ex2(j)-ex(j)^2;


%covxy(j)=exy(j)-x0(u)*ey(j);
%r0(j)=covxy(j)/(ss(j)*sqrt(disp(j)));
%sdr(j)=-10*log10(1-r0(j)^2);
%a(j)=covxy(j)/ss(j)^2;

y1=['zgo',num2str(u),'(j)'];
y2=['=10*log10(2*disp(j)/2852.9^2);'];
y=[y1 y2];
eval(y)
end
z=['plot(s,zgo',num2str(u), ') '];
eval(z)
% plot(s,sdr)
hold on
end
%z=['plot(s,zgo',num2str(u), ') '];
%eval(z)