% 20060414 Mitch Withers
% script to generate various colors of random noise
clear all
close all

N=2^14;
fs=40;
dt=1/fs;
df=fs/N;
T=N*dt;

f=df:df:fs;
f=f';
t=0:dt:T-dt;
t=t';

% get some white noise
% purely theoretical, though here we bandlimit to the nyquist interval
% contain equal power at all frequencies
rand('seed',1);
x=rand(N,1)-0.5;
X=fft(x);
phase=angle(X);
y=(abs(X).^2);

% plot white noise flat psd
subplot(2,1,1)
plot(t,x)
xlabel('time (s)')
ylabel('x(t)')
title('white noise')
subplot(2,1,2)
loglog(f,y)
xlabel('frequency (Hz)')
ylabel('|X(f)|^2')
orient('tall')
%print -dps whitenoise.ps
disp(' ')
disp('White noise has a flat power spectrum')
pause

% get the autocorrelation of the white noise.
tt=(-T+dt:dt:T-dt)';
C = xcorr(x)/sum(abs(x.*x));
plot(tt,C)
xlabel('lag (s)')
title('Normalized autocorrelation of whitenoise')
disp(' ')
disp('White noise has a delta function autocorrelation')
pause

% make pink noise psd proportional to 1/f
% often used in audio engineering because the human ear
% is sensitive to equal power on octaves (i.e. this has
% a flat spectrum in logspace
y1=(f.^(-0.5)).*y;

y2=y1.*exp(-i*(phase+2*pi*f));

% this cheating because if we forced hermetian
% it would be approximately real, but this cheat
% will also show what we're trying to show
xpink = real(ifft(y2)); 

subplot(2,1,1)
plot(t,xpink)
xlabel('time (s)')
ylabel('x(t)')
title('pink noise')
subplot(2,1,2)
loglog(f,y1)
xlabel('frequency (Hz)')
ylabel('|X(f)|^2')
orient('tall')
%print -dps pinknoise.ps
disp(' ')
disp('pink noise has a psd proportional to 1/f')
disp('often used in audio engineering because the human ear')
disp('is sensitive to equal power on octaves')
pause

% find the autocorrelation
Cpink = xcorr(xpink)/sum(abs(xpink.*xpink));
plot(tt,Cpink)
xlabel('lag (s)')
title('Normalized autocorrelation of pink noise')
disp(' ')
disp('Pink noise has a wider impulse function autocorrelation')
pause

% make brown noise psd proportional to 1/f^2
% often modeled by algorithms for brownian motion, the physics
% to show this is left as an exercise
y1=(f.^(-1.0)).*y;
y2=y1.*exp(-i*(phase+2*pi*f));
xbrown = real(ifft(y2));      

subplot(2,1,1)
plot(t,xbrown)
xlabel('time (s)')
ylabel('x(t)')
title('brown noise')
subplot(2,1,2)
loglog(f,y1)
xlabel('frequency (Hz)')
ylabel('|X(f)|^2')
orient('tall')
%print -dps brownnoise.ps
disp(' ')
disp('brown noise has a psd proportional to 1/f^2')
disp('often modeled by algorithms for brownian motion')
pause

% find the autocorrelation
Cbrown = xcorr(xbrown)/sum(abs(xbrown.*xbrown));
plot(tt,Cbrown)
xlabel('lag (s)')
title('Normalized autocorrelation of brown noise')
disp(' ')
disp('Brown noise autocorrelation has an even wider central peak')
pause

% make blue noise psd proportional to f
% useful in the electronic graphics community 
y1=(f.^(0.5)).*y;
y2=y1.*exp(-i*(phase+2*pi*f));
xblue = real(ifft(y2));      

subplot(2,1,1)
plot(t,xblue)
xlabel('time (s)')
ylabel('x(t)')
title('blue noise')
subplot(2,1,2)
loglog(f,y1)
xlabel('frequency (Hz)')
ylabel('|X(f)|^2')
orient('tall')
%print -dps bluenoise.ps
disp(' ')
disp('blue noise has a psd proportional to f')
disp('useful in the electronic graphics community')
pause

% find the autocorrelation
Cblue = xcorr(xblue)/sum(abs(xblue.*xblue));
plot(tt,Cblue)
xlabel('lag (s)')
title('Normalized autocorrelation of blue noise')
disp(' ')
disp('Blue noise autocorrelation is still highly localized')
pause

% make violet noise psd proportional to f^2
% not sure where this is used.  Its differentiated white noise in
% psd space.
y1=(f.^(1.0)).*y;
y2=y1.*exp(-i*(phase+2*pi*f));
xviolet = real(ifft(y2));      

subplot(2,1,1)
plot(t,xviolet)
xlabel('time (s)')
ylabel('x(t)')
title('violet noise')
subplot(2,1,2)
loglog(f,y1)
xlabel('frequency (Hz)')
ylabel('|X(f)|^2')
orient('tall')
%print -dps violetnoise.ps
disp(' ')
disp('violet noise has a psd proportional to f^2')
pause

% find the autocorrelation
Cviolet = xcorr(xviolet)/sum(abs(xviolet.*xviolet));
plot(tt,Cviolet)
xlabel('lag (s)')
title('Normalized autocorrelation of violet noise')
disp(' ')
disp('Violet noise autocorrelation is also highly localized')
pause

disp('Can you explain why the autocorrelation of noise on the blue end of')
disp('the spectrum is highly localized while the red end is not?')

% grey noise, band limited by a complex xfer function not explored here
% arguably the white noise above is psuedo-grey

