dl=ones(1,256);
ol=ones(1,255);


for l=0:7,
   for v=0:15,
   m=l*16+v+1;
   dld1(1,m)=2^l*(v+16.0)-16.5;
   old(1,m)=2^l*(v+16.5)-16.5;
end
end
dv=[4079.5];
dld=[dld1 dv];
dl1=dld(1:1,3:129);
dl2=-1.*fliplr(dl1);
dl=[dl2 dld];
ol1=old(1:1,2:128);
ol2=-1.*fliplr(ol1);
ol=[ol2 old];
s=-87.18:1-57.18;
s1=.05.*s;
ss=(2831.31/sqrt(2))*10.^s1;
%x0=[0 7.5 15.5 31.5 39.5 47.5];
x0=.3;
%for u=1:1
nmin=x0-6.*ss;
nmax=x0+6.*ss;
for j=1:length(ss)
   i(j)=1;
while(i(j)<=255&dl(1,i(j))<nmin(j))
   i(j)=i(j)+1;
end
imin(j)=i(j)-1;
if imin(j)<1
   imin(j)=1;
end
if imin(j)>255
   imin(j)=255;
end
end
for j=1:length(ss)
   i(j)=256;
   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)>255
      imax(j)=255;
   end
end
q=length(ss);
ii=ones(q,255);

for l=1:q
   for k=1:255
      ii(l,k)=gau(dl(1,k),dl(1,k+1),ss(1,l),x0,101);
      hh(l,k)=tau(dl(1,k),dl(1,k+1),ss(1,l),x0,101);
      gg(l,k)=var(dl(1,k),dl(1,k+1),ss(1,l),x0,101);	
end
end
e2=ol.^2;

for l=1:q
   for k=1:255
      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,nn)+...
ol(1,imax(1,j))^2*gau(dl(1,imax(1,j)),nmax(1,j),ss(1,j),x0,nn)+rr(j);
ey(j)=ol(imin(j))*gau(nmin(j),dl(imin(j)+1),ss(j),x0,nn)+ol(imax(j))*gau(dl(imax(j)),nmax(j),ss(j),x0, 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,nn)+ol(imax(j))*tau(dl(imax(j)),nmax(j),ss(j),x0, nn)+bb(j);
ex(j)=tau(nmin(j),dl(imin(j)+1),ss(j),x0,nn)+tau(dl(imax(j)),nmax(j),ss(j),x0, nn)+cc(j);
ex2(j)=var(nmin(j),dl(imin(j)+1),ss(j),x0,nn)+var(dl(imax(j)),nmax(j),ss(j),x0, nn)+dd(j);

covxy(j)=exy(j)-x0*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;
b(j)=20*log10(a(j));
z(j)=-10*log10(disp(j));
end
 plot(s,z)
hold on
%end