private/calculatePST.m
62678829
 %function pst = calculatePST(timeline,pstopts,data)
c54095b9
 function pst = calculatePST(timeline,pstopts,data)
62678829
 %function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
c54095b9
     des         = pstopts.des;
     eventList   = pstopts.eventList;
     sessionList = pstopts.sessionList;
fb860baf
     col_bias_removal = pstopts.colBias;
     norm4SVM         = pstopts.psthNorm; %Normalization method for SVM
c54095b9
 
     bstart          = timeline.baselineStart;
     bend            = timeline.baselineEnd;
9478fa59
     edur            = 12;
c54095b9
     pre             = timeline.psthStart;
     post            = timeline.psthEnd;
f625db4f
     tr              = timeline.tr;
     res             = tr * timeline.trFactor;
9478fa59
 
     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');
6157c8a0
                         pstevnt{zr}(n{zr})=eventList(zr,ze); % NEW 090110 Specify #Event/Condition for each PST
9478fa59
                         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
6157c8a0
                                 nn{zr}=nn{zr}+1;
9478fa59
                                 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;
90f42025
     
6157c8a0
     %%%%COL BIAS REMOVAL, added 090110 %%%%%
fb860baf
 
6157c8a0
     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;
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     
90f42025
     %%%%%%%%%%% new 090109 Axel: "Normalization" for SVM
c343847b
 
fb860baf
 %     disp(['normalization: ' norm4SVM]);   
c343847b
     % 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;
90f42025
     
c343847b
     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;
90f42025
     
     
     
     
9478fa59
 end