function [z,fAxis]=Imem_noisemaker(sampLen,sampFreq,rollOff,noiseSD,normalization,plotting)
%description: 
% From Clements and Bekkers Biophys J 1997: "Two arrays, representing real andimaginary components of the complex noise spectrum,were filled with random numbers between -1 and 1. 
% Each point of the real and imaginary spectrum was scaled by (1+ 300/frequency), and then the total power of the spectrum was normalized back to 1.
% A complex inverse fast Fourier transform was performed on the spectrum, and the resulting real array was scaled by
% Noise SD/(Nd/rT)1/2
% where Noise_SD is the desired standard deviation of simulated noise and Nd is the number of data points in the array. 
% The resulting noise spectrum is approximately l/f at low frequencies, then rolls off to white noise around 300 Hz. 
% The value of 300 Hz is typical for whole-cell recording from cultured hippocampal neurons, but may vary with different preparations and recording conditions. 
% The simulated recording noise had a standard deviation of 1 pA and a peak-to-peak range of approximately +5 pA."

%Tried to recapitulate this, with tremendous help from B. Bonham, though in the end the algorithm had problems with normalization, but appears quite reliable as long as sample length is 10000.

%called by:
%calls:
%notes: 
%code confirmed date:


tAxis=(0:(sampLen-1))/sampFreq;  % time axis (for plotting); if sampFreq is samples/second, and sampLen is total samples, units of tAxis will be seconds

N=sampLen/2;  %the created noise will be sampLen samples long;  we will work initially with a vector half that long and scale it to let that vector represent frequency
fInc=sampFreq/sampLen;  %frequency increment, (samples per second) per sample, i.e., seconds^-1
fAxis=(1:N)*fInc;       %frequency axis (for fft and for plotting); units are ?samples/second

halfSpect=(2*rand(1,N)-1) + sqrt(-1)*(2*rand(1,N)-1); %random real and random imaginary parts; note that halfSpect does not include 0 Hz component
if rollOff
    halfSpect=halfSpect.*(1+rollOff./fAxis);
else
    halfSpect=halfSpect.*(1+1./fAxis);
end

if normalization %normalization of total power -- this doesn't work
    a=real(halfSpect);
    b=a./sum(a);
    c=imag(halfSpect)*sqrt(-1);
    halfSpect=b+c;
end


%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;

if plotting
    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;
end