function [SCM] = spm_roi_fir_desmtx2(xX, Sess, FIRconf)

%  creates the design matrix for a finite impulse response model
%  - allows sub-TR effective resolution of design matrix
% 
%  Arguments:
%  xX, Sess - returned by spm results
%  FIRconf - structure containing the following fields:
%      nPosEst: number of positive peristimulus time elements to estimate
%      nNegEst: number of negative (i.e. prestimulus) PST elements to estimate
%      condition: the condition number to be analyzed
%      TER: the effective temporal resolution for averaging
%  
%  based on fmri_par2scm.m from the FS-FAST toolbox by Doug Greve,
%  MGH-NMR Center.  Modified by Russ Poldrack for use with the SPM
%  Toolbox.   This code is made available under the terms of the GNU 
%  Public License. This software is meant for research only and should 
%  not be used for clinical purposes.
%
% $Id: spm_roi_fir_desmtx2.m,v 1.1 2002/10/28 22:12:56 jcooper Exp $


if nargin<3, 
  fprintf('Error: too few arguments!\n');
  return;
end;

nHEst = FIRconf.NPosEst + FIRconf.NNegEst; % # of estimated points in PSTH
tPreStim = round(FIRconf.NNegEst*FIRconf.TER); % amount of prestimulus time
TR = xX.RT;
TER=FIRconf.TER;
nTP=(max(size(Sess.row))*TR)/TER; % # of estimated timepoints
nRuns = 1;  % only estimate one run in this version
NullCondId = 0;

% doug's code uses parfile notation - first we need to create this from
% trial onsets

Nconds = 2; % we will always do one condition at a time
            % Nconds=2 because null also counts as a condition
	    
% Par is a variable denoting the condition occuring at each timepoint
c=FIRconf.condition;
Par=zeros(nTP,2);
Par(:,1)=[0:TER:(nTP-1)*TER]';
for ons=1:size(Sess.ons{c},1),
    ons_int=round(Sess.ons{c}(ons)/TER)*TER; % find nearest TR for each
                                             % onset
                                             % 
    Par(find(Par(:,1)==ons_int),2)=1;
end;

%% Check for events beyond the last sample %%
tMax = (nTP-1)*TER;

%% Compute the number of StimTypes %%
FirstCondId = 1; 
LastCondId  = Nconds-1;
nNNCond = Nconds -1;

% Get nPerCond %
nPerCond = [];
for n = 1:nNNCond + 1,
  cid = (n-1) + FirstCondId;
  nPerCond(n) = length(squeeze(find(Par(:,2)==cid)));
end

nPreStim = floor(tPreStim/TER);
tTR = [0:TER:tMax];


X =  [];
c1 = zeros(1,nHEst);

% get the indicies of all the excitations for StimId %
iStimIndex = find( Par(:,2) == 1 );

% convert the stimulus times to sample indicies
tPres = Par(iStimIndex,1);
iok = find(tPres >= 0 & tPres <= tMax);
tPres = tPres(iok);
nPres = hist(tPres,tTR);

% construct the conv mtx for this Stimulus Type %
Pulses = [nPres zeros(1,nPreStim) ];
c1(1) = Pulses(1);
E = toeplitz(Pulses,c1);

% add to global conv mtx %
X = cat(2,X,E);

if(nPreStim ~= 0)
    X = X(1+nPreStim:nTP+nPreStim,:);
end

SCM(:,:) = X;

% remove the mean from each column
% testing suggests that this makes no difference whatsoever in the
% results (probably because we use pinv() for estimation), 
% but it seems like the right thing to do

 for col=1:size(SCM,2),
  SCM(:,col)=SCM(:,col)-mean(SCM(:,col));
 end;


return;














