function [eventResponse, time] = cuixuNIRSretrieve(signal, eventOnset, duration, firstTimePoint, tr)
% [eventResponse, time] = cuixuNIRSretrieve(signal, eventOnset, duration,
% firstTimePoint, tr)
% This function is to retrive signal for events
% 
% signal: a column vector, the original signal
% eventOnset: structure, the onset time (in seconds) of event.
% duration: (a single number) the duration (in seconds) of response you are
%   interested.
% firstTimePoint: the relative start time (s) of signal you are interested in (relative the
%   eventOnset). 0 means the retrieved signal starts from the onset of event.
%   -1 means the retrived BOLD starts 1s before the onset of 
%   event. Of couse, firstTimePoint can be any integer value (such as -10).
%   -10 by default
% tr: dt (time between two adjacent sampling, in second) 0.1 by default
%
% eventResponse: A structure. signal for the events.
% time: the time (in seconds) of the response. It's a vector.
%
% Xu Cui
% 2007-09-12
% revised 2008-5-7

if ~exist('firstTimePoint')
    firstTimePoint = -5;
end
if ~exist('tr')
    tr = 0.1;
end


wholeAbsoluteBOLD = signal;

if isempty(eventOnset)
    eventResponse = {};
    time = {};
    warning('nothing has been done b/c event timing data not available.');
    return
end

% convert data format: from structure to cell array if the original
% input is structure
eventisstructure = 0;
if isstruct(eventOnset)
    eventisstructure = 1;
    names = fields(eventOnset);
    eventOnset = struct2cell(eventOnset);
end

% Duration is of unit 'data points', duration is of unit 'second'
Duration = floor(duration/tr);

% adjust the eventOnset time relative to firstTimePoint
for ii=1:length(eventOnset)
    eventOnset{ii} = eventOnset{ii} + firstTimePoint;
end

eventResponse = cell(1,length(eventOnset));
for eventindex=1:length(eventOnset)
    eventResponse{eventindex} = nan*zeros(length(eventOnset{eventindex}), Duration);
end
    
for eventindex = 1:length(eventOnset)
    for j=1:length(eventOnset{eventindex})
        for i=1:Duration
		    s0 = eventOnset{eventindex}(j)/tr + 1 + i - 1;
            if ceil(s0)<=length(wholeAbsoluteBOLD) & floor(s0)>=1
    			eventResponse{eventindex}(j,i) = wholeAbsoluteBOLD(floor(s0)) + (wholeAbsoluteBOLD(ceil(s0))-wholeAbsoluteBOLD(floor(s0))).*(s0 - floor(s0));
            end
        end
    end
end

time = [firstTimePoint:tr:duration+firstTimePoint];
try
    time = time(1:size(eventResponse{1},2));
end

% change the result format (from cell array to structure)
if eventisstructure == 1
    for ii=1:length(names)
        eval(['result2.' names{ii} '=eventResponse{ii};']);
    end
    eventResponse = result2;

end

return