function halfint(opt)
% M-file to illustrate white and various shades of pink and brown noise
%by using fractional integration of a white noise signal
%alpha 1/2 corresponds to pink noise
%If your version of Matlab doesnt make a sound, use instead
%player=audioplayer(f,mylen);
%play(player);
%opt=0; %either opt=1 plot waves and fourier spec or 0 loglog as in pink
mylen=8000;
freq=mylen;
alpha=1/2;
gm=gamma(alpha);
c=ones(1,mylen);
f=2*rand(1,mylen)-c;
g=zeros(1,mylen);
%Plot the white noise signal
if opt==1
    plot(f);
    hold on
end
%sound(f,freq);
%Now make a fractional integral g of the white noise signal f
for i=1:mylen
    x=i/mylen;
    for j=1:i-1
        t=j/mylen;
        g(i)=g(i)+(x-t)^(alpha-1)*f(j)/gm;
    end
end
%Normalize g to have a max +/-1 like f did
g=g./max(abs(g));
if opt==1;
    %Plot g in red and make its sound
    plot(g,'r')
end
sound(g,mylen);
if opt==1
    %Now make a Fourier transform of g so we can see its power spectrum
    L=mylen; 
    Fs=mylen;
    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    Y = fft(g,NFFT);
    %plot(real(Y));
    Y=Y/max(Y);
    f = Fs*linspace(0,1,NFFT/2);
    % Plot single-sided amplitude spectrum in green
    fy=2*abs(Y(1:NFFT/2));
    plot(f,fy,'g')
    hold off
else
    fftsamp=1024;
    Y=fft(g,fftsamp);
    Pyy = Y.* conj(Y) / fftsamp;
    head=32;
    %Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: 
    f = 1000*(0:fftsamp/2)/fftsamp;
    loglog(f(head:fftsamp/2+1),Pyy(head: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
end