function [TTA,TTAsem,count] = caoComputeTTA(design,data,TTAlength);
% Usage:  [TTA,TTAsem,count] = caoComputeTTA(design,data,TTAlength);
%
%  design is a matrix either 1 x nTpts (if all the data were acquired with
%  the same design) or nScans x nTpts (if each scan had a different
%  stimulus presentation order).  Data is a matrix that is nScans x nTpts.
%  The code loops through the design matrix, and for each stimulus type
%  (different stimulus types should be represented by integers from
%  1:nStimTypes), extracts TTAlength timepoints (starting at the stimulus
%  presentation).  These are then averaged to estimate a mean response to
%  each stimulus type, and the standard deviation across the repe

if size(design,2) ~= size(data,2)
    error('Number of timepoints in design and number of timepoints in data don''t match.');
end

% if only 1 design matrix is supplied, then replicate it to match size of
% data
if size(design,1) == 1
    design = repmat(design,[size(data,1) 1]);
end

nStimTypes = max(design(:)); % assumed that each stimulus type is represented 
  % an integer, starting at 1, so the max value in the matrix is the number
  % of stimulus types.
  
% This is a bit of a hack.  On the last couple trials in each scan, when TTAlength 
% is larger than the minimum ISI, you'll run out the end of the scan while
% trying to collect points for the TTA.  So pad the matrix by replicating
% the last timepoint for TTAlength.
data = cat(2,data,repmat(data(:,end),[1 TTAlength]));
  
count = zeros([nStimTypes 1]); % matrix for keeping track of number of stim presentations
for nScan = 1:size(design,1)
    for nTpt = 1:size(design,2)
        if design(nScan,nTpt) % need zeros everywhere there's no stim
            stimType = design(nScan,nTpt); % crucial that design matrix have integers matched to stimType
            count(stimType) = count(stimType) + 1; % increment appropriate count
            exerpts(count(stimType),:,stimType) = data(nScan,nTpt:(nTpt+TTAlength-1));
        end
    end
end
  
if nStimTypes > 1
    TTA = squeeze(mean(exerpts)); % take the mean across repetitions to form the TTA
    TTAsem = squeeze(std(exerpts))./repmat(sqrt(count'),[TTAlength 1]); % stdev works only across the 1st dimension
else
    TTA = mean(exerpts)'; % take the mean across repetitions to form the TTA
    TTAsem = std(exerpts)'./repmat(sqrt(count'),[TTAlength 1]); % stdev works only across the 1st dimension
end
