function z=Imem_noisemaker_bonham(sampLen,sampFreq,noiseSD)


nyquistFreq = sampFreq/2;  % half the sampling rate is the Nyquist freq

taxis = (0:(sampLen-1))/sampFreq;  % time axis (for plotting)

N = sampLen/2;

fInc = sampFreq/sampLen;  % frequency increment in fft frequency axis
faxis = (1:N)*fInc;       % frequency axis (for fft)

halfSpect = (2*rand(1,N)-1) + sqrt(-1)*(2*rand(1,N)-1);
halfSpect = halfSpect .* (1 + 300./faxis);
halfSpect=halfSpect./length(halfSpect);
% halfSpect does not include 0 Hz component

% to do ifft, need to complete the spectrum by copying it -- since for a real
%   time waveform, the spectrum for negative frequencies is the complex conjugate
%   of that for positive frequencies -- and set 0 Hz component to 0 so there's no
%   DC in the time waveform (average is 0)
if mod(sampLen,2)==0,
  faxis = [0 faxis -faxis(end-1:-1:1)];
  halfSpect(end) = real(halfSpect(end));
  randSpect = [0 halfSpect conj(halfSpect((end-1):-1:1))];
else
  faxis = [0 faxis -faxis(end:-1:1)];
  randSpect = [0 halfSpect conj(halfSpect(end:-1:1))];
end

z = ifft(randSpect);

if max(abs(imag(z)))>0,
  % should never get here -- this is just a sanity check
  error('non-zero imaginary half of signal');
end

% normalize so time waveform has desired standard deviation
z = noiseSD/std(z) * z;

% plot
subplot(2,1,1);
loglog(faxis(1:N),abs(randSpect(1:N)));
xlabel('frequency (Hz)');
ylabel('spectral amplitude');
grid on;
subplot(2,1,2);
plot(taxis,z);
xlabel('time (s)');
ylabel('signal amplitude');
grid on;