n = 8192;
nb = 64;
t = 1:n;
waiting = 'please hit any key to continue';

%---------------------------------------------------------------
%     White noise
%---------------------------------------------------------------

white = randn(n,1);
close all
plot(t,white)
title('white noise = normal random values')
disp(waiting)
pause

disp('now playing a soundtrack of white noise')
sound(white)
disp(waiting)
pause

uniform = fft(white);
close all
semilogy(1:n/2, abs(uniform(1:n/2,1)) )
title('the power spectrum of white noise is uniform over frequency')
disp(waiting)
pause

%---------------------------------------------------------------
%     Brown noise
%---------------------------------------------------------------

brown = cumsum(white);
close all
plot(t,brown)
title('brown noise = the cumulative sum (integral) of white noise')
xlabel('Brownian motion is the cumulative sum of many random normal steps')
disp(waiting)
pause

close all
plot(cumsum(randn(n,1)), cumsum(randn(n,1)))
title('two-dimensional Brownian motion (also known as a drunken walk)')
xlabel(' plot(cumsum(randn(n,1)), cumsum(randn(n,1))) ')
disp(waiting)
pause

disp('now playing a soundtrack of brown noise')
sound(brown * 4 / max(abs(brown)))
disp(waiting)
pause

inversesquare = fft(brown);
close all
loglog(1:n/2, abs(inversesquare(1:n/2,1)),  1:n/2, (1:n/2).^(-2) )
title('the power spectrum of brown noise drops off as 1/f^2')
disp(waiting)
pause

%---------------------------------------------------------------
%     Pink noise
%---------------------------------------------------------------

f = (1:n/2)';
randphase = exp(-2i * pi .* rand(n/2,1));
noise = randn(n/2,1) .* randphase;
pinkspectrum = zeros(n,1);
pinkspectrum(1:n/2) = (1./f) .* noise;
for i=2:n/2; pinkspectrum(n+2-i) = conj(pinkspectrum(i)); end

close all
loglog(f,abs(pinkspectrum(1:n/2)), f, 1./f, 'r')
title('the power spectrum of pink noise drops off as 1/f')
disp(waiting)
pause

pink = real(ifft(pinkspectrum));

close all
plot(pink)
title('the power spectrum of pink noise drops off as 1/f')
xlabel('pink noise tends to look like "nasty Brownian mountains"')
disp(waiting)
pause

disp('now playing a soundtrack of pink noise')
sound(pink * 4 / max(abs(pink)))

