clear;
N=100;
m=2*N;
r=rand(m,1);
rn=r(1:N);
rnsum=sum(rn);
rn2=rn'*rn;
for k=1:N
    c(k)=0;
    for i=1:N
        c(k)=c(k)+r(i)*r(i+k);
    end
    c(k)=(N*c(k)-rnsum^2)/(N*rn2-rnsum^2);
end
figure(1);
bar(c);