%samp=2^13;
samp=8000;
b=zeros(1,samp);
%for i=(2:samp)
%	b(i)=b(i)+b(i-1);
%end
c=zeros(1,samp);
x=0;
r=0.37;
s=sqrt(1-r^2);
for i=(1:samp)
	x=r*x+s*randn;
	b(i)=x;
end
r=0.9;
s=sqrt(1-r^2);
for i=(1:samp)
	x=r*x+s*randn;
	b(i)=b(i)+x;
end
r=0.99;
s=sqrt(1-r^2);
for i=(1:samp)
	x=r*x+s*randn;
	b(i)=b(i)+x;
end
mx=max(b);
mn=min(b);
b=b/(mx-mn)/4;
%plot(b);
sound(b);
%player=audioplayer(b,samp);
%play(player);
%c=fft(b);
% plot(real(c));
% hold on
% plot(imag(c),'r');
% hold off
fftsamp=1024;
Y=fft(b,fftsamp);
Pyy = Y.* conj(Y) / fftsamp;


head=36;
%Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: 
f = 1000*(0:fftsamp/2)/fftsamp;
loglog(f,Pyy(1:fftsamp/2+1))
title('Frequency content of y')
xlabel('frequency (Hz)')
hold on
fl=log(f(head:fftsamp/2+1));
pl=log(Pyy(head:fftsamp/2+1));
p=polyfit(fl,pl,1)
myp=p(1)*fl+p(2);
loglog(f(head:fftsamp/2+1),exp(myp),'k')
hold off