| 2095645b |
rawdata=cell2mat({extr.mean}); % Raw Data
pst{voxel} = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,rawdata,sessionList);
end
decodePerformance = [];
for timeShift = timeLineStart:1:timeLineEnd
frameStart = floor(-globalStart+1+timeShift - 0.5*decodeDuration);
frameEnd = min(ceil(frameStart+decodeDuration + 0.5*decodeDuration),-globalStart+globalEnd);
tmp =[];
anyvoxel = 1;
for label = 1:size(pst{1,anyvoxel},2)
for dp = 1:size(pst{1,anyvoxel}{1,label},1) % data point
row = label;
for voxel = 1:size(pst,2)
row = [row, pst{1,voxel}{1,label}(dp,frameStart:frameEnd)]; % label,value,value
end
tmp = [tmp; row];
end
end
svmdata = tmp(:,2:size(tmp,2));
svmlabel = tmp(:,1);
performance = svmtrain(svmlabel, svmdata, svmargs);
minPerformance = min(minPerformance,performance);
maxPerformance = max(maxPerformance,performance);
decodePerformance = [decodePerformance; performance];
end
outputStruct.decodePerformance = decodePerformance;
outputStruct.svmdata = svmdata;
outputStruct.svmlabel = svmlabel;
outputStruct.rawTimeCourse = pst;
outputStruct.minPerformance = minPerformance;
outputStruct.maxPerformance = maxPerformance;
% display(sprintf('Min CrossVal Accuracy: %g%% \t Max CrossVal Accuracy: %g%%',minPerformance,maxPerformance));
end
|
| 2095645b |
% rad = 0; % one voxel
% opt = 1; % xyz coordinates [mm]
vox = voxelList;
nRoi = size(vox,1);
nImg = numel(V);
for k=1:nImg
extr(k) = struct(...
'val', repmat(NaN, [1 nRoi]),...
'mean', repmat(NaN, [1 nRoi]),...
'sum', repmat(NaN, [1 nRoi]),...
'nvx', repmat(NaN, [1 nRoi]),...
'posmm', repmat(NaN, [3 nRoi]),...
'posvx', repmat(NaN, [3 nRoi]));
roicenter = round(inv(V(k).mat)*[vox, ones(nRoi,1)]');
for l = 1:nRoi
% if rad==0
x = roicenter(1,l);
y = roicenter(2,l);
z = roicenter(3,l);
% else
% tmp = spm_imatrix(V(k).mat);
% vdim = tmp(7:9);
% vxrad = ceil((rad*ones(1,3))./(ones(nRoi,1)*vdim))';
% [x y z] = ndgrid(-vxrad(1,l):sign(vdim(1)):vxrad(1,l), ...
% -vxrad(2,l):sign(vdim(2)):vxrad(2,l), ...
% -vxrad(3,l):sign(vdim(3)):vxrad(3,l));
% sel = (x./vxrad(1,l)).^2 + (y./vxrad(2,l)).^2 + ...
% (z./vxrad(3,l)).^2 <= 1;
% x = roicenter(1,l)+x(sel(:));
% y = roicenter(2,l)+y(sel(:));
% z = roicenter(3,l)+z(sel(:));
% end;
dat = spm_sample_vol(V(k), x, y, z,0);
[maxv maxi] = max(dat);
tmp = V(k).mat*[x(maxi); y(maxi); z(maxi);1]; % Max Pos
extr(k).val(l) = maxv;
extr(k).sum(l) = sum(dat);
extr(k).mean(l) = nanmean(dat);
extr(k).nvx(l) = numel(dat);
extr(k).posmm(:,l) = tmp(1:3);
extr(k).posvx(:,l) = [x(maxi); y(maxi); z(maxi)]; % Max Pos
end;
end;
end
% disp(sprintf('Extracted at %.1f %.1f %.1f [xyz(mm)], average of %i voxel(s) [%.1fmm radius Sphere]',vox,length(x),rad));
function pst = calculatePST(des,globalStart,baselineStart,baselineEnd,globalEnd,eventList,data,sessionList)
bstart = baselineStart;
bend = baselineEnd;
edur = 12;
pre = globalStart;
post = globalEnd;
res = 1;
normz = 'file';
pm = 0;
lsess = getNumberOfScans(des);
nSessions = getNumberOfSessions(des);
tr = 2;
[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');
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;
end
|