x=linspace(-4,16,1e4);
p=exp(-(x-6).^2/(2*3^2))/(3*sqrt(2*pi));
z=normrnd(6*ones(1,1e4),3);

[n50,x50]=hist(z(1:50),10);
dx50=x50(2)-x50(1);
p50=n50/50/dx50;


[n500,x500]=hist(z(1:500),100);
dx500=x500(2)-x500(1);
p500=n500/500/dx500;

[n1e4,x1e4]=hist(z,2000);
dx1e4=x1e4(2)-x1e4(1);
p1e4=n1e4/1e4/dx1e4;

[n1e4,x1e4]=hist(z,2000);
dx1e4=x1e4(2)-x1e4(1);
p1e4=n1e4/1e4/dx1e4;

figure(1)
clf

subplot(131)
stairs(x50,p50,'b')
axis([-4 16 0 .15])
hold on
plot(x,p,'k')
xlabel('x','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('p(x)','FontAngle','italic','FontSize',24,'FontName','Times')
title('50 points, 10 bins')
hold off

subplot(132)
stairs(x500,p500,'b')
axis([-4 16 0 .15])
hold on
plot(x,p,'k')
xlabel('x','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('p(x)','FontAngle','italic','FontSize',24,'FontName','Times')
title('500 points, 100 bins')
hold off

subplot(131)
stairs(x50,p50,'b')
axis([-4 16 0 .15])
hold on
plot(x,p,'k')
xlabel('x','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('p(x)','FontAngle','italic','FontSize',24,'FontName','Times')
title('50 points, 10 bins')
hold off

subplot(133)
stairs(x1e4,p1e4,'b')
axis([-4 16 0 .15])
hold on
plot(x,p,'k')
xlabel('x','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('p(x)','FontAngle','italic','FontSize',24,'FontName','Times')
title('10000 points, 2000 bins')
hold off

shg

clear mu sigma
for kk=1:20;
    mu(kk)=mean(z((kk-1)*500+1:kk*500));
    sigma(kk)=std(z((kk-1)*500+1:kk*500));
    
end

clear r
for kk=0:25;
    z1=z(1:50-kk);
    z2=z(kk+1:50);
    r25(kk+1)=mean(z1.*z2);
end

for kk=0:200;
    z1=z(1:400-kk);
    z2=z(kk+1:400);
    r200(kk+1)=mean(z1.*z2);
end


for kk=0:2000;
    z1=z(1:4000-kk);
    z2=z(kk+1:4000);
    r2000(kk+1)=mean(z1.*z2);
end

figure(2)
clf
subplot(211)
plot([0:25],r25,'*-',[0:200],r200,'*-',[0:2000],r2000,'*-',[0 100],36*[1 1],'k-.')
xlim([0 100])
xlabel('\tau (sec)','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('R(\tau)','FontAngle','italic','FontSize',24,'FontName','Times')
title('Autocorrelation function')
legend('T=50s','T=400s','T=4000s')


subplot(212)
semilogx([0:2000],r2000,'r*-')
xlabel('\tau (sec)','FontAngle','italic','FontSize',24,'FontName','Times')
ylabel('R(\tau)','FontAngle','italic','FontSize',24,'FontName','Times')
title('Autocorrelation function')
xlim([1 2e3])
legend('T=4000s')
grid


[p12,w12]=pwelch(z,ones(1,2^12),0);

[p10,w10]=pwelch(z,ones(1,2^10),0);

[p9,w9]=pwelch(z,ones(1,2^9),0);

figure(3)
clf

subplot(211)
semilogy(w12/(2*pi),p12,'*-',w9/(2*pi),p9,'*-')
grid
ylim([.2 30])
xlabel('freq (Hz)','FontSize',24,'FontName','Times')
ylabel('S(f)','FontSize',24,'FontName','Times')
legend('N=2, T/N = 512s, \Delta f= 2 mHz','N=19, T/N = 4096s, \Delta f= 0.2 mHz')

subplot(212)
plot(w12/(2*pi),p12,'*-',w9/(2*pi),p9,'*-')
xlim([0 0.01])
grid
xlabel('freq (Hz)','FontSize',24,'FontName','Times')
ylabel('S(f)','FontSize',24,'FontName','Times')
legend('N=2, T/N = 512s, \Delta f= 2 mHz','N=19, T/N = 4096s, \Delta f= 0.2 mHz')
shg
