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;
   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);
dl=[dl2 dld];
ol1=old(1:1,1:128);
ol2=-1.*fliplr(ol1);
ol=[ol2 old];
s=0:-3:-60;
s1=.05.*s;
ss=(4096)*10.^s1;
x0=[0 16 32 48 56 64];
for u=1:6
nmin=x0(u)-ss;
nmax=x0(u)+ss;
for j=1:length(ss)
   i(j)=1;
while(i(j)<=256&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)>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);

for l=1:q
   for k=1:256
      ii(l,k)=i0sin(dl(1,k),dl(1,k+1),x0(u) ss(1,l));
      hh(l,k)=i1sin(dl(1,k),dl(1,k+1),x0(u),ss(1,l));
      gg(l,k)=i2sin(dl(1,k),dl(1,k+1),x0(u),ss(1,l));	
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);

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*i0sin(nmin(1,j),dl(1,imin(1,j)+1),x0(u),ss(j))+...
ol(1,imax(1,j))^2*i0sin(dl(1,imax(1,j)),nmax(1,j),x0(u),ss(j))+rr(j);
ey(j)=ol(imin(j))*i0sin(nmin(j),dl(imin(j)+1),x0(u),ss(j))+ol(imax(j))*i0sin(dl(imax(j)),nmax(j),x0(u),ss(j))+aa(j);
disp(j)=ey2(j)-ey(j)^2;
exy(j)=ol(imin(j))*i1sin(nmin(j),dl(imin(j)+1),x0(u),ss(j))+ol(imax(j))*i1sin(dl(imax(j)),nmax(j),x0(u),ss(j))+bb(j);
ex(j)=i1sin(nmin(j),dl(imin(j)+1),x0(u),ss(j))+i1sin(dl(imax(j)),nmax(j),x0(u),ss(j))+cc(j);
ex2(j)=i2sin(nmin(j),dl(imin(j)+1),x0(u),ss(j))+i2sin(dl(imax(j)),nmax(j),x0(u),ss(j))+dd(j);

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;
end
 plot(s,sdr)
hold on
end