function h = gaussian_noise(Sh, sampleRate, duration, T, N)
% function h = gaussian_noise(Sh, sampleRate, duration, T, N)
% E. Thrane: (borrows from code in readTimeSeries2) creates a NxT array of h(t)
% given a noise power spectrum Sh
%   T = the number of trials
%   sampleRate = sample rate
%   duration = signal duration
%   Sh = noise power spectral density
%   N is the number of data points

% default number of trials is one
try
  T;
catch
  T=1;
end

% calculate N or use user-specified value
try
  N;
catch
  N = duration * sampleRate;
end

% prepare for FFT
if ( mod(N,2)== 0 )
  numFreqs = N/2 - 1;
else
  numFreqs = (N-1)/2;
end
deltaF = 1/duration;
f = deltaF*[1:1:numFreqs]';
flow = deltaF;

% Next power of 2 from length of y
NFFT = 2^nextpow2(N);
amp_values = Sh(:,2);
f_transfer1 = Sh(:,1);
Pf1 = interp1(f_transfer1, amp_values, f, 'spline');
%
deltaT = 1/sampleRate;
norm1 = sqrt(N/(2*deltaT)) * sqrt(Pf1);
norm1 = repmat(norm1, 1, T);
%
re1 = norm1*sqrt(1/2).* randn(numFreqs, T);
im1 = norm1*sqrt(1/2).* randn(numFreqs, T);
z1  = re1 + 1i*im1;

% freq domain solution for htilde1, htilde2 in terms of z1, z2
htilde1 = z1;
% convolve data with instrument transfer function
otilde1 = htilde1*1;
% set DC and Nyquist = 0, then add negative freq parts in proper order
if ( mod(N,2)==0 )
% note that most negative frequency is -f_Nyquist when N=even
  otilde1 = [ zeros(1,T); otilde1; zeros(1,T); flipud(conj(otilde1)) ];
else
% no Nyquist frequency when N=odd
  otilde1 = [ zeros(1,T); otilde1; flipud(conj(otilde1)) ];
end
% fourier transform back to time domain
o1_data = ifft(otilde1);
% take real part (imag part = 0 to round-off)
h = real(o1_data);
