SVMCrossVal.git

1) %function pst = calculatePST(timeline,pstopts,data)

2) function pst = calculatePST(timeline,pstopts,data)

3) %function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)

4)     des         = pstopts.des;
5)     eventList   = pstopts.eventList;
6)     sessionList = pstopts.sessionList;

7)     col_bias_removal = pstopts.colBias;
8)     norm4SVM         = pstopts.psthNorm; %Normalization method for SVM

9) 
10)     bstart          = timeline.baselineStart;
11)     bend            = timeline.baselineEnd;

12)     edur            = 12;

13)     pre             = timeline.psthStart;
14)     post            = timeline.psthEnd;

15)     tr              = timeline.tr;
16)     res             = tr * timeline.trFactor;

17) 
18)     normz           = 'file';
19)     pm              = 0;
20) 
21)     lsess           = getNumberOfScans(des);
22)     nSessions       = getNumberOfSessions(des);
23) 
24)     [evntrow evntcol]=size(eventList);
25)     
26) 
27)     hsec=str2num(des.xsDes.High_pass_Filter(8:end-3)); % Highpass filter [sec] from SPM.mat
28) 
29)     if strcmp(des.xBF.UNITS,'secs')
30)         unitsecs=1;
31)     end;
32) 
33)     nScansPerSession=getNumberOfScans(des);
34)     %stime=[0:tr:max(nScansPerSession)*tr+post-tr]; % Stimulus time for raw data plot
35)     stime=0:tr:max(nScansPerSession)*tr+round(post/tr)*tr-tr; % Stimulus time for raw data plot
36) 
37) 
38) 
39)     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40)     % RUN
41)     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42) 
43) 
44)     % Digital Highpass
45)     Rp=0.5;
46)     Rs=20;
47)     NO=1;
48)     Wp=1/((1/2/tr)/(1/hsec));
49)     [B, A] = ellip(NO,Rp,Rs,Wp,'high');
50) 
51)     sdata(1:max(nScansPerSession)+round(post/tr),1:nSessions)=nan; % Open Data Matrix
52)     for z=1:nSessions % Fill Data Matrix sessionwise
53)         sdata(1:nScansPerSession(z),z)=data(sum(nScansPerSession(1:z))-nScansPerSession(z)+1:sum(nScansPerSession(1:z)))';
54)     end;
55) %         usdata=sdata; % Keep unfiltered data
56) 
57)     sdatamean=nanmean(nanmean(sdata(:,:)));
58)     for z=1:nSessions
59) %             X(:,z)=[1:1:max(nScansPerSession)]'; % #Volume
60)         sdata(1:nScansPerSession(z),z)=filtfilt(B,A,sdata(1:nScansPerSession(z),z)); %Filter Data (Highpass)
61)     end;
62)     sdata=sdata+sdatamean;
63) 
64) 
65)     %%%%Parametric Modulation Modus%%%%
66)     if pm %Find Parameters for Event of Interest
67)         [imods modss mods erow evntrow eventList] = getParametricMappingEvents(eventList,evntrow,des,pmf);
68)     end;
69)     %%%%PM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70) 
71) 
72)     for zr=1:evntrow
73)         n{zr}=0;
74)         nn{zr}=0; 
75)         nnn{zr}=0;
76)         sstart{zr}=1;
77)     end;
78) 
79) 
80)     sesst0=0; 
81)     for sessionID=sessionList
82)         if sessionID>1
83)             sesst0(sessionID)=sum(lsess(1:sessionID-1))*tr;  
84)         end;
85)         for zr=1:evntrow  %LABEL NUMBER, EVENT GROUP
86)             sstart{zr}=n{zr}+1;
87)             for ze=1:evntcol % EVENT INDEX in EventList
88)                 if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
89)                     for zz=1:length(des.Sess(sessionID).U(eventList(zr,ze)).ons) % EVENT REPETITION NUMBER
90)                         if ~unitsecs
91)                             des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)-1)*tr;
92)                             des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)=(des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)-1)*tr;
93)                         end;
94) 
95)                         nnn{zr}=nnn{zr}+1; % INFO for rawdataplot start
96)                         if des.Sess(sessionID).U(eventList(zr,ze)).dur(zz)<edur
97)                             mev{zr}(nnn{zr},1:2)=[des.Sess(sessionID).U(eventList(zr,ze)).ons(zz)+sesst0(sessionID) edur]; % modeled event [onset length]
98)                         else
99)                             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)];
100)                         end; % INFO for rawdataplot end
101) 
102)                         n{zr}=n{zr}+1;
103)                         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');

104)                         pstevnt{zr}(n{zr})=eventList(zr,ze); % NEW 090110 Specify #Event/Condition for each PST

105)                         if strcmp(normz,'epoc')
106)                             bline=nanmean(pst{zr}(n{zr},round(-pre/res+(bstart)/res+1):round(-pre/res+(bend)/res+1)));
107)                             if isnan(bline)
108)                                 pst{zr}(n{zr},1:-pre/res+post/res+1)=nan;
109)                             else

110)                                 nn{zr}=nn{zr}+1;

111)                                 pst{zr}(n{zr},:)=(pst{zr}(n{zr},:)-bline)/bline*100; % 'epoch-based' normalization
112)                             end;
113)                         end;
114)                     end;
115)                 end;
116)             end;
117)             if ~strcmp(normz,'epoc')
118)                 bline(zr)=nanmean(nanmean(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
119)                 bstd(zr)=nanmean(nanstd(pst{zr}(sstart{zr}:n{zr},-pre/res+(bstart)/res+1:-pre/res+(bend)/res+1)));
120)                 nn{zr}=n{zr};
121)             end;
122)         end;
123)         if strcmp(normz,'filz')
124)             for zr=1:evntrow
125)                 pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bstd); % session-based z-score normalization
126)             end;
127)         elseif strcmp(normz,'file')
128)             for zr=1:evntrow
129)                 pst{zr}(sstart{zr}:n{zr},:)=(pst{zr}(sstart{zr}:n{zr},:)-mean(bline))/mean(bline)*100; % session-based normalization
130)             end;
131)         end;
132)     end;

133)     

134)     %%%%COL BIAS REMOVAL, added 090110 %%%%%

135) 

136)     if col_bias_removal
137) 
138)         for zr=1:evntrow % Event ROW
139)             pst4row(zr,:)=nanmean(pst{zr}); % Average PST for overall ROW
140)             n_col=0;
141)             for ze=1:evntcol % Even COL
142)                 if ze==1 || (ze>1 && eventList(zr,ze)~=eventList(zr,ze-1))
143)                     n_col=n_col+1;
144)                     tmp_index=find(pstevnt{zr}(:)==eventList(zr,ze));
145)                     % Calculate Average PST for each COL within ROW
146)                     pst4col(n_col,:)=nanmean(pst{zr}(tmp_index,:));
147)                     % COL BIAS: Calculate Difference between COL Ave and Overall Row Ave
148)                     pst_delta(n_col,:)=pst4col(n_col,:)-pst4row(zr,:);
149)                     % Remove COL BIAS from single PSTs
150)                     for z=1:length(tmp_index)
151)                         pst{zr}(tmp_index(z),:)=pst{zr}(tmp_index(z),:)-pst_delta(n_col,:);
152)                     end;
153)                 end;
154)             end;
155)         end;
156)     end;
157)     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
158)     

159)     %%%%%%%%%%% new 090109 Axel: "Normalization" for SVM

160) 

161) %     disp(['normalization: ' norm4SVM]);   

162)     % none - no normalization
163)     % mean - mean normalization (meanPSTH max at .5, baseline at -.5]
164)     % minmax - all PSTHs between [0 1]
165)     
166)     for zr=1:evntrow
167)           tmp_maxmean(zr)=max(nanmean(pst{zr}));
168)           tmp_max(zr)=max(max(pst{zr}));
169)           tmp_min(zr)=min(min(pst{zr}));
170)     end;

171)     

172)     pstmaxmean=max(tmp_maxmean);
173)     pstmax=max(tmp_max);
174)     pstmin=min(tmp_min);
175)     
176)     for zr=1:evntrow
177)         switch norm4SVM
178)             case {'none'}
179)             case {'mean'}
180)                 pst{zr}=pst{zr}./pstmaxmean-.5;
181)             case {'minmax'}
182)                 pst{zr}=(pst{zr}-pstmin)./(pstmax-pstmin);
183)         end;
184)     end;

185)     
186)     
187)     
188)