%displays the 1/f noise of intermittency left of pd3
%hold on
r=3.8284;
samp=2^13;
b=zeros(1,samp);
b(1)=.1;
for i=2:5000
    b(1)=r*b(1)*(1-b(1));
end
for i=2:samp
    b(i)=r*b(i-1)*(1-b(i-1));
end
%window=samp;
window=512;
sound(b);
Y=fft(b,window);
Pyy = Y.* conj(Y) / window;
hold on
%Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: 
f = 1000*(0:window/2)/window;
loglog(f,Pyy(1:window/2+1))
tail=85;
fl=log(f(2:tail));
pl=log(Pyy(2:tail));
p=polyfit(fl,pl,1)
myp=p(1)*fl+p(2);
loglog(f(2:tail),exp(myp),'k')
r=3.8285;
samp=2^13;
b=zeros(1,samp);
b(1)=.1;
for i=2:5000
    b(1)=r*b(1)*(1-b(1));
end
for i=2:samp
    b(i)=r*b(i-1)*(1-b(i-1));
end
sound(b);
Y=fft(b,window);
Pyy = Y.* conj(Y) / window;

%Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: 
f = 1000*(0:window/2)/window;
loglog(f,Pyy(1:window/2+1),'r')
title('Frequency content of y')
xlabel('frequency (Hz)')
hold off