%function pst = calculatePST(timeline,pstopts,data)
function pst = calculatePST(timeline,pstopts,data)
%function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
des = pstopts.des;
eventList = pstopts.eventList;
sessionList = pstopts.sessionList;
col_bias_removal = pstopts.colBias;
norm4SVM = pstopts.psthNorm; %Normalization method for SVM
bstart = timeline.baselineStart;
bend = timeline.baselineEnd;
edur = 12;
pre = timeline.psthStart;
post = timeline.psthEnd;
tr = timeline.tr;
res = tr * timeline.trFactor;
normz = 'file';
pm = 0;
lsess = getNumberOfScans(des);
nSessions = getNumberOfSessions(des);
[evntrow evntcol]=size(eventList);
hsec=str2num(des.xsDes.High_pass_Filter(8:end-3)); % Highpass filter [sec] from SPM.mat
if strcmp(des.xBF.UNITS,'secs')
unitsecs=1;
end;
nScansPerSession=getNumberOfScans(des);
%stime=[0:tr:max(nScansPerSession)*tr+post-tr]; % Stimulus time for raw data plot
stime=0:tr:max(nScansPerSession)*tr+round(post/tr)*tr-tr; % Stimulus time for raw data plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Digital Highpass
Rp=0.5;
Rs=20;
NO=1;
Wp=1/((1/2/tr)/(1/hsec));
[B, A] = ellip(NO,Rp,Rs,Wp,'high');
sdata(1:max(nScansPerSession)+round(post/tr),1:nSessions)=nan; % Open Data Matrix
for z=1:nSessions % Fill Data Matrix sessionwise
sdata(1:nScansPerSession(z),z)=data(sum(nScansPerSession(1:z))-nScansPerSession(z)+1:sum(nScansPerSession(1:z)))';
end;
% usdata=sdata; % Keep unfiltered data
sdatamean=nanmean(nanmean(sdata(:,:)));
for z=1:nSessions
% X(:,z)=[1:1:max(nScansPerSession)]'; % #Volume
sdata(1:nScansPerSession(z),z)=filtfilt(B,A,sdata(1:nScansPerSession(z),z)); %Filter Data (Highpass)
end;
sdata=sdata+sdatamean;
%%%%Parametric Modulation Modus%%%%
if pm %Find Parameters for Event of Interest
[imods modss mods erow evntrow eventList] = getParametricMappingEvents(eventList,evntrow,des,pmf);
end;
%%%%PM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for zr=1:evntrow
n{zr}=0;
nn{zr}=0;
nnn{zr}=0;
sstart{zr}=1;
end;
sesst0=0;
for sessionID=sessionList
if sessionID>1
sesst0(sessionID)=sum(lsess(1:sessionID-1))*tr;
end;
for zr=1:evntrow %LABEL NUMBER, EVENT GROUP
sstart{zr}=n{zr}+1;
for ze=1:evntcol % EVENT INDEX in EventList
if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
for zz=1:length(des.Sess(sessionID).U(eventList(zr,ze)).ons) % EVENT REPETITION NUMBER
if ~unitsecs
des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)-1)*tr;
des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)-1)*tr;
end;
nnn{zr}=nnn{zr}+1; % INFO for rawdataplot start
if des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)<edur
mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) edur]; % modeled event [onset length]
else
mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)];
end; % INFO for rawdataplot end
n{zr}=n{zr}+1;
pst{zr}(n{zr},:)=interp1(stime,sdata(:,sessionID),[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+pre:res:des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+post],'linear');
pstevnt{zr}(n{zr})=eventList(zr,ze); % NEW 090110 Specify #Event/Condition for each PST
if strcmp(normz,'epoc')
bline=nanmean(pst{zr}(n{zr},round(-pre/res+(bstart)/res+1):round(-pre/res+(bend)/res+1)));
if isnan(bline)
pst{zr}(n{zr},1:-pre/res+post/res+1)=nan;
else
nn{zr}=nn{zr}+1;
pst{zr}(n{zr},:)=(pst{zr}(n{zr},:)-bline)/bline*100; % 'epoch-based' normalization
end;
end;
end;
end;
end;
if ~strcmp(normz,'epoc')
bline(zr)=nanmean(nanmean(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
bstd(zr)=nanmean(nanstd(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
nn{zr}=n{zr};
end;
end;
if strcmp(normz,'filz')
for zr=1:evntrow
pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bstd); % session-based z-score normalization
end;
elseif strcmp(normz,'file')
for zr=1:evntrow
pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bline)*100; % session-based normalization
end;
end;
end;
%%%%COL BIAS REMOVAL, added 090110 %%%%%
if col_bias_removal
for zr=1:evntrow % Event ROW
pst4row(zr,:)=nanmean(pst{zr}); % Average PST for overall ROW
n_col=0;
for ze=1:evntcol % Even COL
if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
n_col=n_col+1;
tmp_index=find(pstevnt{zr}(:)==eventList(zr,ze));
% Calculate Average PST for each COL within ROW
pst4col(n_col,:)=nanmean(pst{zr}(tmp_index,:));
% COL BIAS: Calculate Difference between COL Ave and Overall Row Ave
pst_delta(n_col,:)=pst4col(n_col,:)-pst4row(zr,:);
% Remove COL BIAS from single PSTs
for z=1:length(tmp_index)
pst{zr}(tmp_index(z),:)=pst{zr}(tmp_index(z),:)-pst_delta(n_col,:);
end;
end;
end;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% new 090109 Axel: "Normalization" for SVM
% disp(['normalization: ' norm4SVM]);
% none - no normalization
% mean - mean normalization (meanPSTH max at .5, baseline at -.5]
% minmax - all PSTHs between [0 1]
for zr=1:evntrow
tmp_maxmean(zr)=max(nanmean(pst{zr}));
tmp_max(zr)=max(max(pst{zr}));
tmp_min(zr)=min(min(pst{zr}));
end;
pstmaxmean=max(tmp_maxmean);
pstmax=max(tmp_max);
pstmin=min(tmp_min);
for zr=1:evntrow
switch norm4SVM
case {'none'}
case {'mean'}
pst{zr}=pst{zr}./pstmaxmean-.5;
case {'minmax'}
pst{zr}=(pst{zr}-pstmin)./(pstmax-pstmin);
end;
end;
end