function PN = PinkNoise(n, exponent)

if (nargin < 2)
  exponent = 1; % pure Pink Noise
end

if (nargin < 1)
  n = 1024; % default length of vector of Pink Noise to be generated
end
half = (n/2+1);

randPhase = exp(2 * pi * 1i * rand(n,1));

PinkSpectrum = randn(n,1);
PinkSpectrum(1) = PinkSpectrum(1) * randPhase(1);
for f = 2:half
  PinkSpectrum(f) = (1/f)^exponent * PinkSpectrum(f) * randPhase(f);
  PinkSpectrum(n+2-f) = conj(PinkSpectrum(f));
end   

plot( 2:half, log10(abs(PinkSpectrum(2:half))) )
title('Pink Noise has a 1/f spectrum')

PN = abs(ifft(PinkSpectrum));
plot( 2:half, PN(2:half) )
title('Pink Noise is 1/f noise')

