function hc = noise_gen(file, fs, segdur, nsegs, verbose)
% function hc = noise_gen(file, fs, segdur, nsegs, verbose)
%   E Thrane
% Generate colored Gaussian noise specified by the PSD describe in file.
% The first column of file should be frequency and the second column should be
% PSD.
% intput arguments:
%   fs = sampling frequency
%   segdur = duration of time series
%   nsegs = number of segments
%   verbose = print diagnostic information and plots
% output:
%   hc = colored noise

% default is verbose=false
try
  verbose;
catch
  verbose=false;
end

% the number of data points to sample in each segment
N = fs*segdur;

% begin with white noise
h = randn(N, nsegs);

% normalize so that mean(psd) = normfac2*mean(abs(hf).^2) = 1
h = h * sqrt(fs/2);

% fft white noise
[hf f] = fft_eht(h, fs);

% diagnostic print statements
if verbose
  % diagnostic with Welch's method
  [pp0 ff0] = psd(h(:,1), size(h,1), fs);
  normfac = 2/fs;
  fprintf('mean of white noise psd is = %1.3e\n', mean(normfac*pp0));
  % diagnostic with other method
  normfac2 = 2/(fs*N);
  fprintf('mean of white noise psd is = %1.3e\n', ...
    normfac2*mean(abs(hf(:,1)).^2));
end

% load psd file: the first column should be frequency and the second should be
% the psd in strain^2/Hz
g = load(file);

% parse psd file
ff = g(:,1);
Pf = g(:,2);

% interpolate
Pfi = interp1(ff, Pf, f);

% fill in out-of-band frequency bins using the nearest available bin
Pfi(f<ff(1)) = Pf(1);
Pfi(f>ff(end)) = Pf(end);

% color the noise so that the expectation value of mean(psd) = Pf
hfc = hf .* repmat(sqrt(Pfi), 1, nsegs);

% inverse Fourier transform to get time series
hc = ifft_eht(hfc);

% sanity check: convert back to frequency domain and make a psd
if verbose
  [ppc, ffc] = psd(hc(:,1), size(hc,1), fs);
  figure;
  loglog(ff, Pf, 'k', 'linewidth', 2);
  hold on;
  loglog(ffc, normfac*ppc, 'b', 'linewidth', 2);
  hold on;
  loglog(f, normfac2*abs(hfc).^2, 'r');
  grid on;
  xlabel('f (Hz)');
  ylabel('strain^2/Hz');
  pretty;
  print('-dpng', 'noise_gen');
end

return
