function data = caoSimulateData(neural,TR,HIRFamp,SNR,crop)

% Usage:  data = caoSimulateData(neural,TR,HIRFamp,SNR,[crop])
%  Expectation is that "neural" represents the normalized neural response
%  at each TR in the experiment.  HIRFamp is the peak BOLD response to
%  neural activity of amplitude 1 lasting 1 TR.  SNR is used to scale the
%  noise that is added to the simulated BOLD.  Normally distributed (white) 
%  noise is added to the data with standard deviation (1/SNR).  If crop is
%  equal to 0 or not given, the full convolution result is returned.  If
%  crop is provided and equal to 1, then the returned data are the same
%  length as the provided neural response.

HIRF = spm_hrf(TR);

HIRF = HIRF*HIRFamp/max(HIRF(:));

data = conv(neural,HIRF);

data = data + (1/SNR)*randn(size(data));

if exist('crop')
    if crop == 1
        data = data(1:length(neural));
    end
end