SVMCrossVal.git

1) %function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
2) function pst = calculatePST(timeline,pstopts,data)
3)     des         = pstopts.des;
4)     eventList   = pstopts.eventList;
5)     sessionList = pstopts.sessionList;
6) 
7)     bstart          = timeline.baselineStart;
8)     bend            = timeline.baselineEnd;

9)     edur            = 12;

10)     pre             = timeline.psthStart;
11)     post            = timeline.psthEnd;

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

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

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

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

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

130)     

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

132)     col_bias_removal=1;

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

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

164) 

165)     norm4SVM='mean'; %Normalization method for SVM

166)     disp(['normalization: ' norm4SVM]);   
167)     % none - no normalization
168)     % mean - mean normalization (meanPSTH max at .5, baseline at -.5]
169)     % minmax - all PSTHs between [0 1]
170)     
171)     for zr=1:evntrow
172)           tmp_maxmean(zr)=max(nanmean(pst{zr}));
173)           tmp_max(zr)=max(max(pst{zr}));
174)           tmp_min(zr)=min(min(pst{zr}));
175)     end;

176)     

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

190)     
191)     
192)     
193)