From b9c62d543f8d9f02fbc7fd3d1ce43f963fe9f65f Mon Sep 17 00:00:00 2001
From: Proulx-S <sebastien.proulx2@mail.mcgill.ca>
Date: Fri, 2 Aug 2024 01:16:18 -0400
Subject: [PATCH] -refactor space-time-frequency analysis -start refactoring
 ploting

---
 dtrnd2.m        |   1 +
 getResp.m       |  54 +++++++++++------
 getTapers.m     |   2 +-
 plotSpecAll.m   |  29 ++++++---
 plotTs2.m       |   8 ++-
 runFullMT2.m    | 157 ++++++++++++++++++++++++++++++++++++------------
 tryL2svd.m      |  26 ++++++--
 volPsdFullMt2.m | 141 +++++++++++++++++--------------------------
 volTsGetResp.m  | 125 +++++++++++++++++++++++++++++---------
 9 files changed, 348 insertions(+), 195 deletions(-)

diff --git a/dtrnd2.m b/dtrnd2.m
index 017e963..795bc4e 100644
--- a/dtrnd2.m
+++ b/dtrnd2.m
@@ -32,6 +32,7 @@
         [funTs(I).vec,funTs(I).poly] = dtrnd2(funTs(I).vec,funTs(I).t,prcBOLDflag,order);
         funTs(I).poly.vol2vec = funTs(I).vol2vec;
     end
+    order = size(funTs(I).poly.beta,2)-1;
     funTs = setNiceFieldOrder(funTs,{'vol' 'vol2vec' 'vec' 't' 'poly' 'volInfo' 'vecInfo'});
     poly = [];
 
diff --git a/getResp.m b/getResp.m
index 9a2f2c9..c220291 100644
--- a/getResp.m
+++ b/getResp.m
@@ -20,24 +20,32 @@
 
 
 %% Mask files
-if isfield(volAnat.fun.mask,'head') && ~isempty(volAnat.fun.mask.head)
-    %%%head
-    fMask = volAnat.fun.mask.head.mri.fspec;
-elseif isfield(volAnat.fun.mask,'brain') && ~isempty(volAnat.fun.mask.brain)
-    %%%brain
-    fMask = volAnat.fun.mask.brain.mri.fspec;
-    % mask = mask & any(volAnat.fun.mask.brain.mri.vol,4);
-else
+if isempty(volAnat)
     fMask = [];
+else
+    if isfield(volAnat.fun.mask,'head') && ~isempty(volAnat.fun.mask.head)
+        %%%head
+        fMask = volAnat.fun.mask.head.mri.fspec;
+    elseif isfield(volAnat.fun.mask,'brain') && ~isempty(volAnat.fun.mask.brain)
+        %%%brain
+        fMask = volAnat.fun.mask.brain.mri.fspec;
+        % mask = mask & any(volAnat.fun.mask.brain.mri.vol,4);
+    else
+        fMask = [];
+    end
 end
 
-
 %% Functional design
 if isfield(volTs,'dsgn')
     param.funDsgn.trStim   = volTs.dsgn.dt;
     param.funDsgn.k        = 1;
-    param.funDsgn.startSeq = volTs.dsgn.onsets;
-    param.funDsgn.durSeq   = volTs.dsgn.ondurs;
+    if isfield(volTs.dsgn,'onsets')
+        param.funDsgn.startSeq = volTs.dsgn.onsets;
+        param.funDsgn.durSeq   = volTs.dsgn.ondurs;
+    else
+        param.funDsgn.startSeq = volTs.dsgn.onsetList;
+        param.funDsgn.durSeq   = volTs.dsgn.ondurList;
+    end
     param.funDsgn.condSeq  = ones(size(param.funDsgn.startSeq));
     param.funDsgn.label    = volTs.dsgn.label;
     param.funDsgn.trDecon  = param.trDecon;
@@ -903,14 +911,14 @@
         TENTzeroFlag = 1;
 
 
-        %% Contruct afni command
+        %% Construct afni command
         cmdTmp = {srcAfni};
         if size(fList,2)>1
             cmdTmp{end+1} = ['echo ''  ''run' num2str(I) '/' num2str(size(fList,1)) ' -- echo' num2str(E) '/' num2str(size(fList,2))];
         else
             cmdTmp{end+1} = ['echo ''  ''run' num2str(I) '/' num2str(size(fList,1))];
         end
-        if ~exist(fStat,'file') || force
+        % if ~exist(fStat,'file') || force
             if ~isfield(param.funDsgn,'trDecon')
                 if isfield(param.funDsgn,'trStim')
                     param.funDsgn.trDecon = param.funDsgn.trStim;
@@ -925,14 +933,20 @@
             cmdTmp{end+1} = ['echo ''   ''' fStat];
             cmdTmp{end+1} = ['echo ''   ''' fMat];
             cmdTmp{end+1} = 'echo ''   ''done';
-        else
-            cmdTmp{end+1} = ['echo ''   ''' fResp];
-            cmdTmp{end+1} = ['echo ''   ''' fStat];
-            cmdTmp{end+1} = ['echo ''   ''' fMat];
-            cmdTmp{end+1} = 'echo ''   ''already done, skipping';
-        end
+        % else
+        %     cmdTmp{end+1} = ['echo ''   ''' fResp];
+        %     cmdTmp{end+1} = ['echo ''   ''' fStat];
+        %     cmdTmp{end+1} = ['echo ''   ''' fMat];
+        %     cmdTmp{end+1} = 'echo ''   ''already done, skipping';
+        % end
         fRun(I,E).cmd =  strjoin(cmdTmp,newline);
 
+        if exist(fStat,'file') && ~force
+            cmdTmp(startsWith(cmdTmp,'-') | startsWith(cmdTmp,'3dDeconvolve')) = [];
+            cmdTmp{ismember(cmdTmp,{'echo ''   ''done'})} = 'echo ''   ''already done, skipping';
+            cmdTmp(endsWith(cmdTmp,'  ''done')) = [];
+        end
+
         cmd = [cmd cmdTmp];
         % [status,cmdout] = system(strjoin(cmdTmp,newline),'-echo');
     end
@@ -1117,7 +1131,7 @@
 % trMri = 3;
 cmd{end+1} = ['-stim_times_subtract ' num2str(trMri*nDummy,'%f') ' \'];
 cmd{end+1} = '-num_stimts 1 \';
-cmd{end+1} = ['-stim_label 1 ' label{1} ' \'];
+cmd{end+1} = ['-stim_label 1 ' char(label) ' \'];
 
 % write design to file
 fido = fopen(fStim, 'w');
diff --git a/getTapers.m b/getTapers.m
index b78a8e0..bf9e38c 100644
--- a/getTapers.m
+++ b/getTapers.m
@@ -3,7 +3,7 @@
 
 [TW,W,K] = K2W(N*tr,K,0);
 
-if isempty(t)
+if isempty(t) || max(abs(diff(t,2)))<1e-10
     [tp,eigs] = dpsschk([TW K],N,1/tr); % check tapers
     eigs = permute(eigs,[2 1]);
     % t = permute(0:tr:((N-1)*tr),[2 1]);
diff --git a/plotSpecAll.m b/plotSpecAll.m
index 9b04a0f..2a5aebc 100644
--- a/plotSpecAll.m
+++ b/plotSpecAll.m
@@ -1,7 +1,8 @@
-function plotSpecAll(volPsd,volTs,respQthresh)
+function plotSpecAll(volPsd,volTs,volResp,respQthresh)
+saveFlag = 1;
 
 if ~exist('respQthresh','var'); respQthresh = []; end
-if isempty(respQthresh);        respQthresh = 0.05; end % 'psd' 'coh'
+if isempty(respQthresh);        respQthresh = 1; end % respQthresh = 1 means don't threshold
 
 Fpsd = figure('WindowStyle','docked');
 HtPsd = tiledlayout(8,1); HtPsd.TileSpacing = 'tight'; HtPsd.Padding = 'tight';
@@ -14,14 +15,15 @@ function plotSpecAll(volPsd,volTs,respQthresh)
 axPsd{end+1} = nexttile;
 axPsd{end+1} = plotSpecGram(volPsd,'trialGramMD','psd'   ,HtPsd,volTs,respQthresh);
 axPsd{end+1} = plotSpecGram(volPsd,'trialGramMD','psdEPC',HtPsd,volTs,respQthresh);
-cLim = get([axPsd{3:5}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
-set([axPsd{3:5}],'CLim',cLim);
-cLim = get([axPsd{7:8}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
-set([axPsd{7:8}],'CLim',cLim);
+% cLim = get([axPsd{3:5}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
+% set([axPsd{3:5}],'CLim',cLim);
+% cLim = get([axPsd{7:8}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
+% set([axPsd{7:8}],'CLim',cLim);
 
-[~,b,~] = fileparts(replace(volPsd.fspec,'.nii.gz',''));
+[~,b,~] = fileparts(fileparts(volPsd.fspec));
+% [~,b,~] = fileparts(replace(volPsd.fspec,'.nii.gz',''));
 b = strsplit(b,'_');
-title(HtPsd,b{1})
+title(HtPsd,strjoin(b(1:3)))
 
 
 Fcoh = figure('WindowStyle','docked');
@@ -73,4 +75,13 @@ function plotSpecAll(volPsd,volTs,respQthresh)
 cLim = get([axCoh{7:8}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
 set([axCoh{7:8}],'CLim',cLim);
 
-title(HtCoh,b{1})
+% title(HtCoh,b{1})
+title(HtCoh,strjoin(b(1:3)))
+
+
+%% save
+if saveFlag
+    [a,b,~] = fileparts(fileparts(volPsd.fspec));
+    saveas(Fpsd,fullfile(a,[b '_tPsd.fig']))
+    saveas(Fcoh,fullfile(a,[b '_tCoh.fig']))
+end
\ No newline at end of file
diff --git a/plotTs2.m b/plotTs2.m
index f99a7c2..4d06ca0 100644
--- a/plotTs2.m
+++ b/plotTs2.m
@@ -4,9 +4,9 @@
 if ~exist('onsets','var');            onsets = []; end
 if ~exist('ondurs','var');            ondurs = []; end
 if ~exist('roiInd','var');            roiInd = []; end
-if ~exist('respQthresh','var'); respQthresh = []; end
+if ~exist('respQthresh','var');  respQthresh = []; end
 if ~exist('volAnat','var');          volAnat = []; end
-if isempty(respQthresh);         respQthresh = 0.05; end
+if isempty(respQthresh);         respQthresh = 1; end % respQthresh = 1 means don't threshold
 if isempty(volAnat);                  volRoi = [];
                                      volMask = [];
 elseif ~isempty(roiInd)
@@ -135,9 +135,13 @@
 
 if isempty(onsets) && isfield(volTs,'dsgn') && isfield(volTs.dsgn,'onsets')
     onsets = volTs.dsgn.onsets;
+elseif isempty(onsets) && isfield(volTs,'dsgn') && isfield(volTs.dsgn,'onsetList')
+    onsets = volTs.dsgn.onsetList';
 end
 if isempty(ondurs) && isfield(volTs,'dsgn') && isfield(volTs.dsgn,'ondurs')
     ondurs = volTs.dsgn.ondurs;
+elseif isempty(ondurs) && isfield(volTs,'dsgn') && isfield(volTs.dsgn,'ondurList')
+    ondurs = volTs.dsgn.ondurList';
 end
 
 
diff --git a/runFullMT2.m b/runFullMT2.m
index 7a6d193..a9265c9 100644
--- a/runFullMT2.m
+++ b/runFullMT2.m
@@ -24,8 +24,20 @@
 if ~exist('onsets','var'); onsets = []; end
 if ~exist('ondurs','var'); ondurs = []; end
 
-if isempty(onsets); onsets = funTs.dsgn.onsets; end
-if isempty(ondurs); ondurs = funTs.dsgn.ondurs; end
+if isempty(onsets)
+    if isfield(funTs.dsgn,'onsets')
+        onsets = funTs.dsgn.onsets;
+    else
+        onsets = funTs.dsgn.onsetList;
+    end
+end
+if isempty(ondurs)
+    if isfield(funTs.dsgn,'ondurs')
+        ondurs = funTs.dsgn.ondurs;
+    else
+        ondurs = funTs.dsgn.ondurList;
+    end
+end
 
 % if K==1; skipSVD = true; end
 
@@ -74,7 +86,7 @@
 if ~exist('durList','var');     durList = []; end
 
 if isempty(win); win = [inf 0]; end
-
+windFlag = 1;
 
 
 if funTs.nvoxels==1; if verbose; disp('only one timeseries, skipping SVD'); end; skipSVD = true; end
@@ -154,7 +166,7 @@
     allWin = repmat(1:param.win(1),[param.win(3) 1]);
     allWin = allWin + (((1:param.win(3))-1)*param.win(2))'; % win x t
     allWin(any(allWin>funTs.nframes,2),:) = [];
-    allWin(end+1,:) = (funTs.nframes-param.win(1)+1:funTs.nframes)';
+    % allWin(end+1,:) = (funTs.nframes-param.win(1)+1:funTs.nframes)';
     param.win(3) = [];
     param.win = param.win*tr;
     if verbose && param.win(2)~=inf
@@ -172,8 +184,10 @@
 if ~skipTrialGram
     allWin; % [win X timeIndex]
     onsetList = param.onsetList;
+    if windFlag; onsetList(1) = []; end
     winSz = param.win(1)./tr;
-    n = max(allWin(:));
+    n = funTs.nframes;
+    % n = max(allWin(:));
     nWin = size(allWin,1);
     nTrial = size(onsetList,1);
     allWin2 = repmat({zeros(n,nWin)},[nTrial 1]);
@@ -182,28 +196,48 @@
             if trialInd == 1
                 allWin2{trialInd}(allWin(winInd,:),winInd) = 1;
             else
-                offsetInd = floor((onsetList(trialInd) - onsetList(1)) ./ tr);
+                % offsetInd = floor((onsetList(trialInd) - onsetList(1)) ./ tr);
+                offsetInd = floor((onsetList(trialInd) - onsetList(1)) ./ tr) + 1;
                 tInd = allWin(winInd,:) + offsetInd;
                 tInd(tInd>n) = [];
                 allWin2{trialInd}(tInd,winInd) = 1;
             end
         end
     end
-    allWin3 = any(cat(3,allWin2{:}),3); % [win X time]
 
+    %%% allign windows to the end of the timeseries
+    timeShiftVol = n-find(allWin2{1}(:,end),1,'last');
+    for trialInd = 1:nTrial
+        allWin2{trialInd}(end-timeShiftVol+1:end,:) = false;
+        allWin2{trialInd} = circshift(allWin2{trialInd},timeShiftVol,1);
+    end
+    
     %%% remove windows exceeding timeseries
-    endInd = find(allWin3(end,:)==1,1)+1;
+    allWin3 = any(cat(3,allWin2{:}),3); % [win X time]
+    endInd = find(allWin3(end,:),1)+1;
+    % endInd = find(allWin3(end,:),1);
     for trialInd = 1:nTrial
         allWin2{trialInd}(:,endInd:end) = [];
     end
     % allWin3(:,endInd:end) = [];
 
+    % %%% remove completely overlapping windows
+    % endInd = find(sum((allWin2{1}+allWin2{2}(:,1))==2,1)==winSz);
+    % for trialInd = 1:nTrial
+    %     allWin2{trialInd}(:,endInd:end) = [];
+    % end
+    % % allWin3(:,endInd:end) = [];
+
+
     %%% remove completely overlapping windows
-    endInd = find(sum((allWin2{1}+allWin2{2}(:,1))==2,1)==winSz);
+    [a,endInd] = max(sum(allWin2{end}==allWin2{end-1}(:,end),1));
+    if a<n; endInd = endInd - 1; end
     for trialInd = 1:nTrial
-        allWin2{trialInd}(:,endInd:end) = [];
+        allWin2{trialInd}(:,1:endInd) = [];
     end
-    % allWin3(:,endInd:end) = [];
+    % imagesc(any(cat(3,allWin2{:}),3))
+    
+    
 
     % %%% visualize
     % figure('WindowStyle','docked');
@@ -231,7 +265,7 @@
         allWin4(:,winInd) = find(allWin3(:,winInd));
     end
     allWinTrialLock = permute(allWin4,[2 1 3]); % [win X timeIndex]
-    % reshape(allWinTrialLock,size(allWinTrialLock,1),prod(size(allWinTrialLock,[2 3])))
+    % allWinTrialLock = reshape(allWinTrialLock,size(allWinTrialLock,1),prod(size(allWinTrialLock,[2 3])));
     clear allWin2 allWin3 allWin4 winSz2
 else
     allWinTrialLock = [];
@@ -360,7 +394,7 @@
 
     % tic
     [funPsd.psd,funPsd.psdGram,funPsd.psdTrialGram,funPsd.psdTrialGramMD,funPsd.svd,funPsd.svdGram,funPsd.svdTrialGram,funPsd.svdTrialGramMD]...
-        = computeAll(funTs,TP,param,skip,verbose);
+        = computeAll(funTs,TP,param,windFlag,skip,verbose);
     % disp('+++++')
     % toc
     % disp('+++++')
@@ -521,7 +555,7 @@
 end
 
 
-function [psd,psdGram,psdTrialGram,psdTrialGramMD,svd,svdGram,svdTrialGram,svdTrialGramMD] = computeAll(funTs,TP,param,skip,verbose)
+function [psd,psdGram,psdTrialGram,psdTrialGramMD,svd,svdGram,svdTrialGram,svdTrialGramMD] = computeAll(funTs,TP,param,windFlag,skip,verbose)
 if ~exist('nShuf','var');         nShuf = []; end
 if ~exist('nRun','var');           nRun = []; end
 if ~exist('cohFrange','var'); cohFrange = []; end
@@ -540,7 +574,8 @@
     % Over full timeseires %
     %%%%%%%%%%%%%%%%%%%%% %%
     %[time x trial x run x taper x freq x vox x window x mode]
-
+    tic
+    disp('full timeseries analysis')
     %%% tapers
     tp = permute(TP.full.tp,[1 3 4 2 5 6 7 8]); % tapers[time x trial x run x taper x freq x vox x window x mode]
     [Nk,Ek,Rk,Kk,Fk,Vk,Wk,M] = size(tp);
@@ -597,9 +632,10 @@
         dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
         prm = [ 6   4   2   1   3   5   7   8  ];
         dim = strjoin(dim(prm),' ');
-        j = permute(J,prm); %[V K E N R F W Mk]
-        j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*Mk]
-        [u,s,v] = pagesvd(j,'econ','vector'); % s[Mk V E*N*R*F*W]
+        % j = permute(J,prm); %[V K E N R F W Mk]
+        % j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*Mk]
+        % [u,s,v] = pagesvd(j,'econ','vector'); % s[Mk V E*N*R*F*W]
+        [u,s,~] = pagesvd(reshape(permute(J,prm),[V K E*1*R*F*1*1]),'econ','vector'); % s[Mk V E*N*R*F*W]
         coh     = s.^2./sum(s.^2,1); % coherence[Mk V E*N*R*F*W]
         spSVmag = abs(u);
         % coh = reshape(coh,[M 1 E 1 R F 1]); % [Mk V E N R F W]
@@ -700,11 +736,15 @@
     res.full.dim     = [N E R K F V W M];
     res.full.dimInfo = '[N E R K F V W M]';
 
+    toc
+
 
 
     %% %%%%%%%%%%%%%%%%%%%%%
     % Over each timewindow %
     %%%%%%%%%%%%%%%%%%%%% %%
+    tic
+    disp('time-resolved analysis')
     if ~skip.gram
         %[time x trial x run x taper x freq x vox x window x mode]
         %[   7       2     1       3      5     8       20]
@@ -743,6 +783,7 @@
         res.gram.COH    = zeros(1,E,1,1,F,1,W,M  ); % coherence [time x trial x run x taper x freq x vox x window x mode] at each trial
 
         %%% loop over windows
+        fprintf([repmat('|',1,W) '\n\n']);
         for wInd = 1:W
             %%% Compute J
             ind  = w(:,:,:,:,:,:,wInd);
@@ -767,9 +808,10 @@
                 dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
                 prm = [ 6   4   2   1   3   5   7   8  ];
                 dim = strjoin(dim(prm),' ');
-                j = permute(J,prm); %[V K E N R F W Mk]
-                j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*Mk]
-                [u,s,v] = pagesvd(j,'econ','vector'); % s[Mk V E*N*R*F*W]
+                % j = permute(J,prm); %[V K E N R F W Mk]
+                % j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*Mk]
+                % [u,s,v] = pagesvd(j,'econ','vector'); % s[Mk V E*N*R*F*W]
+                [~,s,~] = pagesvd(reshape(permute(J,prm),[V K E*1*R*F*1*1]),'econ','vector'); % s[Mk V E*N*R*F*W]
                 coh = s.^2./sum(s.^2,1); % coherence[Mk V E*N*R*F*W]
                 coh = reshape(coh,[M 1 E 1 R F 1]); % [Mk V E N R F W]
                 dim = {'Mk' 'V' 'E' 'N' 'R' 'F' 'W' 'K'};
@@ -781,6 +823,7 @@
 
             %%% Output window time
             res.gram.t(:,:,:,:,:,:,wInd,:) = tWin;
+            fprintf('\b''\n');
         end
         if K == 1
             res.gram.COH = [];
@@ -797,11 +840,19 @@
         res.gram = [];
     end
 
+    toc
 
 
     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % Over each event-related timewindow %
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
+    % Much time could be saved by reusing the output of the time-resolved
+    % analysis, as long as the step sizes allow to construct every
+    % event-related window from available time-resolved windows. However,
+    % the missing-data event-related analysis seems superior to the
+    % event-related analysis so the latter might not be needed.
+    tic
+    disp('trial-locked time-resolved analysis')
     %[time x trial x run x taper x freq x vox x window x mode]
     %[   7       2     1       3      5     8       20]
     if ~skip.trialGram
@@ -851,6 +902,7 @@
 
 
         %%% loop over windows
+        fprintf([repmat('|',1,W) '\n\n']);
         for wInd = 1:W
             %%% Compute J
             ind  = w(:,:,:,:,:,:,wInd,:);
@@ -884,9 +936,10 @@
                 dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'M'};
                 prm = [ 6   4   2   1   3   5   7   8  ];
                 dim = strjoin(dim(prm),' ');
-                j = permute(J,prm); %[V K E N R F W M]
-                j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*M]
-                [u,s,v] = pagesvd(j,'econ','vector'); % s[M V E*N*R*F*W]
+                % j = permute(J,prm); %[V K E N R F W M]
+                % j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*M]
+                % [u,s,v] = pagesvd(j,'econ','vector'); % s[M V E*N*R*F*W]
+                [~,s,~] = pagesvd(reshape(permute(J,prm),[V K E*1*R*F*1*1]),'econ','vector'); % s[M V E*N*R*F*W]
                 coh = s.^2./sum(s.^2,1); % coherence[M V E*N*R*F*W]
                 coh = reshape(coh,[M 1 E 1 R F 1]); % [M V E N R F W]
                 %       1   2   3   4   5   6   7   8
@@ -904,9 +957,10 @@
                 dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
                 prm = [ 6   4   2   1   3   5   7   8  ];
                 dim = strjoin(dim(prm),' ');
-                j = permute(mean(J,2),prm); %[V K E N R F W Mk]
-                j = reshape(j,[V K 1*1*R*F*1*1]); %[V K E*N*R*F*W*M]
-                [u,s,v] = pagesvd(j,'econ','vector'); % s[M V E*N*R*F*W]
+                % j = permute(mean(J,2),prm); %[V K E N R F W Mk]
+                % j = reshape(j,[V K 1*1*R*F*1*1]); %[V K E*N*R*F*W*M]
+                % [u,s,v] = pagesvd(j,'econ','vector'); % s[M V E*N*R*F*W]
+                [~,s,~] = pagesvd(reshape(permute(mean(J,2),prm),[V K 1*1*R*F*1*1]),'econ','vector'); % s[M V E*N*R*F*W]
                 coh = s.^2./sum(s.^2,1); % coherence[M V E*N*R*F*W]
                 coh = reshape(coh,[M 1 1 1 R F 1]); % [M V E N R F W]
                 dim = {'Mk' 'V' 'E' 'N' 'R' 'F' 'W' 'K'};
@@ -917,18 +971,24 @@
             end
 
             %%% Compute coherence with trials concatenated as extra sets of tapers (equivalent to averaging across trials)
-            dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
-            prm = [ 6   4   2   1   3   5   7   8  ];
-            dim = strjoin(dim(prm),' ');
-            j = permute(J,prm); %[V K E N R F W Mk]
-            j = reshape(j,[V K E 1*R*F*1*1]); %[V K E N*R*F*W*Mk]
-            dim = {'V' 'K' 'E' 'N*R*F*W*Mk'};
-            prm = [ 1   4   2   3          ];
-            dim = strjoin(dim(prm),' ');
-            j = permute(j,prm); %[V N*R*F*W*Mk K E]
-            j = reshape(j,[V 1*R*F*1*1 K*E]); %[V N*R*F*W*Mk K*E]
-            j = permute(j,[1 3 2]); %[V K*E N*R*F*W*Mk];
-            [u,s,v] = pagesvd(j,'econ','vector'); % s[Mke V N*R*F*W K E]
+            % dim1 = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
+            % prm1 = [ 6   4   2   1   3   5   7   8  ];
+            % dim1 = strjoin(dim1(prm1),' ');
+            % j = permute(J,prm1); %[V K E N R F W Mk]
+            % j = reshape(j,[V K E 1*R*F*1*1]); %[V K E N*R*F*W*Mk]
+            % dim2 = {'V' 'K' 'E' 'N*R*F*W*Mk'};
+            % prm2 = [ 1   4   2   3          ];
+            % dim2 = strjoin(dim2(prm2),' ');
+            % j = permute(j,prm2); %[V N*R*F*W*Mk K E]
+            % j = reshape(j,[V 1*R*F*1*1 K*E]); %[V N*R*F*W*Mk K*E]
+            % j = permute(j,[1 3 2]); %[V K*E N*R*F*W*Mk];
+            dim1 = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
+            prm1 = [ 6   4   2   1   3   5   7   8  ];
+            dim1 = strjoin(dim1(prm1),' ');
+            dim2 = {'V' 'K' 'E' 'N*R*F*W*Mk'};
+            prm2 = [ 1   4   2   3          ];
+            dim2 = strjoin(dim2(prm2),' ');
+            [~,s,~] = pagesvd(permute(reshape(permute(reshape(permute(J,prm1),[V K E 1*R*F*1*1]),prm2),[V 1*R*F*1*1 K*E]),[1 3 2]),'econ','vector'); % s[Mke V N*R*F*W K E]
             coh = s.^2./sum(s.^2,1); % coherence[Mke V N*R*F*W K E]
             coh = reshape(coh,[Mek 1 1 1 F 1 1 1]); % [Mke V N R F W K E]
             %       1   2   3   4   5   6   7   8
@@ -941,6 +1001,7 @@
 
             %%% Output window time
             res.trialGram.t(:,:,:,:,:,:,wInd,:) = tWin;
+            fprintf('\b''\n');
         end
         if K == 1
             res.trialGram.COH    = [];
@@ -960,10 +1021,15 @@
         res.trialGram = [];
     end
 
+    toc
+
+
     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % Over each event-related timewindow   %
     % (with missing data seperating tials) %
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
+    tic
+    disp('trial-locked time-resolved analysis (using missing data tapers)')
     %[time x trial x run x taper x freq x vox x window x mode]
     %[   7       2     1       3      5     8       20]
     if ~skip.trialGram
@@ -985,7 +1051,15 @@
         t = TP.trialGramMD.t;
         [Nt,Et,Rt,Kt,Ft,Vt,Wt,Mt] = size(t);
         %%%% shifted to 0 at stim onset (phase coherent cross-trial avg)
-        tPC = reshape(reshape(t,Nw,E) - funTs.dsgn.onsets',Nw*E,1);
+        if isfield(funTs.dsgn,'onsets')
+            onsets = funTs.dsgn.onsets;
+        else
+            onsets = funTs.dsgn.onsetList;
+        end
+        if windFlag
+            onsets(1) = [];
+        end
+        tPC = reshape(reshape(t,Nw,E) - onsets,Nw*E,1);
         
 
         %%% freq
@@ -1009,6 +1083,7 @@
         res.trialGramMD.COHepc = zeros(1,1,1,1,F,1,W,M); % coherence [time x trial x run x taper x freq x vox x window x mode] trials concatenated as extra sets of tapers (eVENT AS TAPERS k   )
 
         %%% loop over windows
+        fprintf([repmat('|',1,W) '\n\n']);
         for wInd = 1:W
             %%% Compute J
             ind  = w(:,:,:,:,:,:,wInd);
@@ -1082,6 +1157,7 @@
 
             %%% Output window time
             res.trialGramMD.t(:,:,:,:,:,:,wInd,:) = tWin;
+            fprintf('\b''\n');
         end
         if K == 1
             res.trialGramMD.COH    = [];
@@ -1100,6 +1176,7 @@
         res.trialGramMD = [];
     end
 
+    toc
 
 
 
diff --git a/tryL2svd.m b/tryL2svd.m
index 95ff92b..bfc9cda 100644
--- a/tryL2svd.m
+++ b/tryL2svd.m
@@ -1,6 +1,7 @@
 function tryL2svd(volPsd)
 method = 1;
 adjustPhase = 1;
+saveFlag = 1;
 
 N = volPsd.svd.dim(1);
 W = volPsd.svd.dim(2);
@@ -11,7 +12,8 @@ function tryL2svd(volPsd)
 W = volPsd.svd.dim(7);
 M = volPsd.svd.dim(8);
 
-subStr = strsplit(replace(volPsd.fspec,'.nii.gz',''),filesep); subStr = strsplit(subStr{end},'_');
+% subStr = strsplit(replace(volPsd.fspec,'.nii.gz',''),filesep); subStr = strsplit(subStr{end},'_');
+subStr = strsplit(replace(volPsd.fspec,'.nii.gz',''),filesep); subStr = strsplit(subStr{end-1},'_');
 
 spSV = permute(volPsd.svd.spSV,[5 6 8 1 2 3 4 7]); % [F V M]
 switch method
@@ -64,7 +66,7 @@ function tryL2svd(volPsd)
 cLim = get(axIm,'CLim');
 cLim = [min([cLim{:}]) max([cLim{:}])];
 set(axIm,'CLim',cLim);
-title(ht,subStr{1})
+title(ht,strjoin(subStr(1:3)))
 
 
 % add underlay
@@ -129,7 +131,7 @@ function tryL2svd(volPsd)
 axIm = [ax{2:end}];
 cLim = [-pi pi];
 set(axIm,'CLim',cLim);
-title(ht,subStr{1})
+title(ht,strjoin(subStr(1:3)))
 
 
 % add underlay
@@ -170,7 +172,7 @@ function tryL2svd(volPsd)
 
 
 %% Frequency singular vector
-figure('WindowStyle','docked');
+fSV = figure('WindowStyle','docked');
 ht = tiledlayout(GridSize(1),GridSize(2)); ax = {};
 ht.TileSpacing = 'tight'; ht.Padding = 'tight';
 ax{end+1} = nexttile;
@@ -190,12 +192,24 @@ function tryL2svd(volPsd)
     title(['Level-2 mode ' num2str(ml2) '; method' num2str(method)])
     grid on
 
-    xline(1/unique(diff(volPsd.dsgn.onsets)),'--r')
+    if isfield(volPsd.dsgn,'onsets')
+        xline(1/mean(diff(volPsd.dsgn.onsets)),'--r')
+    else
+        xline(1/mean(diff(volPsd.dsgn.onsetList)),'--r')
+    end
 end
 axIm = [ax{2:end}];
 yLim = get(axIm,'YLim');
 yLim = [min([yLim{:}]) max([yLim{:}])];
 set(axIm,'YLim',yLim);
 
-title(ht,subStr{1})
+title(ht,strjoin(subStr(1:3)))
+
+
 
+if saveFlag
+    [a,b,~] = fileparts(fileparts(volPsd.fspec));
+    saveas(fMag,  fullfile(a,[b '_L2SVD_spSvMag.fig']))
+    saveas(fPhase,fullfile(a,[b '_L2SVD_spSvPhase.fig']))
+    saveas(fSV,   fullfile(a,[b '_L2SVD_tSvMag.fig']))
+end
\ No newline at end of file
diff --git a/volPsdFullMt2.m b/volPsdFullMt2.m
index f33edad..855725a 100644
--- a/volPsdFullMt2.m
+++ b/volPsdFullMt2.m
@@ -68,95 +68,56 @@
 %% Do the processing here
 %%%%%%%%%%%%%%%%%%%%%%%%%
 
+
+    
+
 %% Mask
 if ~isempty(volAnat)
-    volTs = vec2vol(vol2vec(volTs));
-    for I = 1:length(volAnat)
-        if isfield(volAnat(I),'fun') && isfield(volAnat(I).fun,'mask')
-
-            mask = volAnat(I).fun.mask.crop.vol;
-            
-            if isfield(volAnat(I).fun.mask,'head')
-                mask = mask & volAnat(I).fun.mask.head.mri.vol;
-            elseif isfield(volAnat(I).fun.mask,'brain')
-                mask = mask & volAnat(I).fun.mask.brain.mri.vol;
-            end
+    if length(volAnat)==1 && ...
+            isfield(volAnat,'mask') && isfield(volAnat.mask,'crop') && isfield(volAnat.mask.crop,'mri') && isfield(volAnat.mask.crop.mri,'vol') && ~isempty(volAnat.mask.crop.mri.vol)
+        %%%crop
+        mask = volAnat.mask.crop.mri.vol;
+        %%%head
+        mask = mask & any(volAnat.mask.head.mri.vol,4);
+        % %%%brain
+        % mask = mask & any(volAnat.mask.brain.mri.vol,4);
+        %%%apply
+        volTs = applyMask(volTs,mask);
+    else
+        volTs = vec2vol(vol2vec(volTs));
+        for I = 1:length(volAnat)
+            if isfield(volAnat(I),'fun') && isfield(volAnat(I).fun,'mask')
+
+                mask = volAnat(I).fun.mask.crop.vol;
+
+                if isfield(volAnat(I).fun.mask,'head')
+                    mask = mask & volAnat(I).fun.mask.head.mri.vol;
+                elseif isfield(volAnat(I).fun.mask,'brain')
+                    mask = mask & volAnat(I).fun.mask.brain.mri.vol;
+                end
+
+                if ~isempty(volPsd)
+                    threshPrctl = 99;
+                    Mi = 1;
+                    [~,Fi] = min(abs(volPsd.svd.f - 1/mean(diff(volPsd.dsgn.onsets))));
+                    threshMask = false(size(volPsd.vol2vec));
+                    threshMask(volPsd.vol2vec) = abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi))>prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl);
+                    % figure('WindowStyle','docked');
+                    % imagesc(mask>prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl))
+                    % histogram(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi))); xline(prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl))
+                    mask = mask & threshMask;
+                end
+
+                %%%apply
+                volTs(I) = applyMask(volTs(I),mask);
 
-            if ~isempty(volPsd)
-                threshPrctl = 99;
-                Mi = 1;
-                [~,Fi] = min(abs(volPsd.svd.f - 1/mean(diff(volPsd.dsgn.onsets))));
-                threshMask = false(size(volPsd.vol2vec));
-                threshMask(volPsd.vol2vec) = abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi))>prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl);
-                % figure('WindowStyle','docked');
-                % imagesc(mask>prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl))
-                % histogram(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi))); xline(prctile(abs(volPsd.svd.spSV(:,:,:,:,Fi,:,:,Mi)),threshPrctl))
-                mask = mask & threshMask;
-            end
 
-            %%%apply
-            volTs(I) = applyMask(volTs(I),mask);
-
-
-
-
-            % % volTs(I) = applyMask(volTs(I),volAnat(I).fun.mask.crop.vol);
-            % 
-            % % %%%arteries
-            % % volRoi = volAnat(I).fun.roi.vesselCenter;
-            % % ind  = squeeze(contains(volRoi.label,'a'));
-            % % volRoi.mri.vol(:,:,:,~ind) = [];
-            % % volRoi.label(:,:,:,~ind)   = [];
-            % % 
-            % % combineFlag = 1;
-            % % vol2vec(volTs(I))
-            % % volTs(I) = applyRoi(volTs(I),volRoi,[],combineFlag);
-            % 
-            % 
-            % % 
-            % % 
-            % % 
-            % %%%crop
-            % mask = volAnat(I).fun.mask.crop.vol;
-            % % 
-            % % % %%%vessel and surroundings
-            % % % mask = mask & any(volAnat(I).fun.mask.vessel.mri.vol,4);
-            % % 
-            % % 
-            % % 
-            % % 
-            % % volMask.mri.vol   = any(volRoi.mri.vol,4);
-            % % volMask.mri.label = strjoin(volRoi.label,'+');
-            % % 
-            % % % %%%brain
-            % % % mask = mask & any(volAnat(I).fun.mask.brain.mri.vol,4);
-            % % 
-            % % % %%%head
-            % % % mask = mask & any(volAnat(I).fun.mask.head.mri.vol,4);
-            % % 
-            % % 
-            % % if exist('volMask','var')
-            % %     volMask.mri.vol
-            % %     roi = permute(arteries,[4 1 2 3]);
-            % %     volTs(I).vecRoi = permute(roi(:,any(roi,1))~=0,[3 2 4 5 1]);
-            % %     volTs(I).vecInfo = strjoin({volTs(I).vecInfo 'roi'},' x ');
-            % % end
-            % % 
-            % % 
-            % % mask = mask & arteries(:,:,:,1)~=0;
-            % % if exist('arteries','var')
-            % %     roi = permute(arteries,[4 1 2 3]);
-            % %     volTs(I).vecRoi = permute(roi(:,any(roi,1))~=0,[3 2 4 5 1]);
-            % %     volTs(I).vecInfo = strjoin({volTs(I).vecInfo 'roi'},' x ');
-            % % end
-            % % 
-            % % mask = mask & arteries(:,:,:,1)~=0;
-            % % 
+            end
         end
+        % if isfield(volAnat,'fun') && isfield(volAnat.fun,'mask') && isfield(volAnat.fun.mask,'crop') && ~isempty(volAnat.fun.mask.crop.vol)
+        %     volTs = applyMask(volTs,volAnat.fun.mask.crop.vol);
+        % end
     end
-    % if isfield(volAnat,'fun') && isfield(volAnat.fun,'mask') && isfield(volAnat.fun.mask,'crop') && ~isempty(volAnat.fun.mask.crop.vol)
-    %     volTs = applyMask(volTs,volAnat.fun.mask.crop.vol);
-    % end
 end
 
 % %% Add time
@@ -188,7 +149,10 @@
 
 phaseRand = 0;
 taperPerm = info.perm;
-% info.onsetList
+if isempty(info.onsetList)
+    info.onsetList = volTs.dsgn.onsetList';
+    info.ondurList = volTs.dsgn.ondurList';
+end
 volPsd = runFullMT2(volTs,W,K,[winSz winStep],info.onsetList,info.ondurList,[],1,info.skipSvd,[],[],taperPerm,phaseRand);
 
 % phaseRand = 2^2;
@@ -210,7 +174,8 @@
 volPsd.tr = mean(diff(volPsd.psd.f))*1000;
 tmp = vec2vol(volPsd);
 tmp.vol(1:10,1:10,1,:) = repmat(permute(tmpMean,[2 3 4 1]),[10 10 1 1]);
-volPsd.psd.fspec = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['logPsd']},'_')) '.nii.gz'];
+volPsd.psd.fspec = replace(volTs.fspec,'_volTs.nii.gz','_volPsd.nii.gz');
+% volPsd.psd.fspec = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['logPsd']},'_')) '.nii.gz'];
 MRIwrite(tmp,volPsd.psd.fspec);
 
 %%% Coherence (first singular vector)
@@ -224,7 +189,8 @@
     if info.perm
         tmp.vol(end-9:end,end-9:end,1,:) = repmat(permute(volPsd.svd.COH_permMean(:,:,:,:,:,:,:,1),[1 2 3 5 4 6 7 8]),[10 10 1 1]);
     end
-    volPsd.svd.fspec.spSVmag = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['part-mag'] ['spSV']},'_')) '.nii.gz'];
+    volPsd.svd.fspec.spSVmag = replace(volTs.fspec,'_volTs.nii.gz','_volCohMag.nii.gz');
+    % volPsd.svd.fspec.spSVmag = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['part-mag'] ['spSV']},'_')) '.nii.gz'];
     MRIwrite(tmp,volPsd.svd.fspec.spSVmag);
     tmp = vec2vol(volPsd);
     tmp.vol = angle(tmp.vol);
@@ -232,7 +198,8 @@
     if info.perm
         tmp.vol(end-9:end,end-9:end,1,:) = repmat(permute(volPsd.svd.COH_permMean(:,:,:,:,:,:,:,1),[1 2 3 5 4 6 7 8]),[10 10 1 1]);
     end
-    volPsd.svd.fspec.spSVphase = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['part-phase'] ['spSV']},'_')) '.nii.gz'];
+    volPsd.svd.fspec.spSVphase = replace(volTs.fspec,'_volTs.nii.gz','_volCohPhase.nii.gz');
+    % volPsd.svd.fspec.spSVphase = [fullfile(info.preprocDir,strjoin({['sub-' info.sub] ['ses-' info.ses] ['part-phase'] ['spSV']},'_')) '.nii.gz'];
     MRIwrite(tmp,volPsd.svd.fspec.spSVphase);
 end
 if isfield(volPsd.svd,'spSV_pVal')
diff --git a/volTsGetResp.m b/volTsGetResp.m
index edff5b2..f819aea 100644
--- a/volTsGetResp.m
+++ b/volTsGetResp.m
@@ -1,11 +1,13 @@
-function [out, info] = volTsGetResp(do,info,volTs,volAnat)
+function [out, info] = volTsGetResp(do,info,volTs,volAnat,force)
 if isempty(do)
     do.loadIt = 0;
     do.doIt   = 1;
     do.saveIt = 0;
 end
 
+if ~exist('force','var');     force = []; end
 if ~exist('volAnat','var'); volAnat = []; end
+if isempty(force);            force = 0; end
 
 % if ~isfield(info,'K');                   info.K = []        ; end
 % if ~isfield(info,'win');               info.win = zeros(0,2); end
@@ -17,8 +19,20 @@
 
 
 if ~isfield(volTs,'dsgn');   volTs.dsgn = []        ; end
-if ~isempty(volTs.dsgn); onsets = volTs.dsgn.onsets; end
-if ~isempty(volTs.dsgn); ondurs = volTs.dsgn.ondurs; end
+if ~isempty(volTs.dsgn)
+    if isfield(volTs.dsgn,'onsets')
+        onsets = volTs.dsgn.onsets;
+    else
+        onsets = volTs.dsgn.onsetList;
+    end
+end
+if ~isempty(volTs.dsgn)
+    if isfield(volTs.dsgn,'ondurs')
+        ondurs = volTs.dsgn.ondurs;
+    else
+        ondurs = volTs.dsgn.ondurList;
+    end
+end
 
 
 
@@ -73,27 +87,41 @@
 
 %% Mask
 if ~isempty(volAnat)
-    volTs = vol2vec(volTs);
-    for I = 1:length(volAnat)
-        if isfield(volAnat(I),'fun') && isfield(volAnat(I).fun,'mask') && isfield(volAnat(I).fun.mask,'crop') && ~isempty(volAnat(I).fun.mask.crop.vol)
-            % volTs(I) = applyMask(volTs(I),volAnat(I).fun.mask.crop.vol);
-
-            %%%crop
-            mask = volAnat(I).fun.mask.crop.vol;
-            if isfield(volAnat(I).fun.mask,'head') && ~isempty(volAnat(I).fun.mask.head)
-                %%%head
-                mask = mask & any(volAnat(I).fun.mask.head.mri.vol,4);
-            elseif isfield(volAnat(I).fun.mask,'brain') && ~isempty(volAnat(I).fun.mask.brain)
-                %%%brain
-                mask = mask & any(volAnat(I).fun.mask.brain.mri.vol,4);
+    if length(volAnat)==1 && ...
+            isfield(volAnat,'mask') && isfield(volAnat.mask,'crop') && isfield(volAnat.mask.crop,'mri') && isfield(volAnat.mask.crop.mri,'vol') && ~isempty(volAnat.mask.crop.mri.vol)
+        %%%crop
+        mask = volAnat.mask.crop.mri.vol;
+        %%%head
+        mask = mask & any(volAnat.mask.head.mri.vol,4);
+        % %%%brain
+        % mask = mask & any(volAnat.mask.brain.mri.vol,4);
+        %%%apply
+        volTs = applyMask(volTs,mask);
+    else
+
+        volTs = vol2vec(volTs);
+        for I = 1:length(volAnat)
+            if isfield(volAnat(I),'fun') && isfield(volAnat(I).fun,'mask') && isfield(volAnat(I).fun.mask,'crop') && ~isempty(volAnat(I).fun.mask.crop.vol)
+                % volTs(I) = applyMask(volTs(I),volAnat(I).fun.mask.crop.vol);
+
+                %%%crop
+                mask = volAnat(I).fun.mask.crop.vol;
+                if isfield(volAnat(I).fun.mask,'head') && ~isempty(volAnat(I).fun.mask.head)
+                    %%%head
+                    mask = mask & any(volAnat(I).fun.mask.head.mri.vol,4);
+                elseif isfield(volAnat(I).fun.mask,'brain') && ~isempty(volAnat(I).fun.mask.brain)
+                    %%%brain
+                    mask = mask & any(volAnat(I).fun.mask.brain.mri.vol,4);
+                end
+                %%%apply
+                volTs(I) = applyMask(volTs(I),mask);
             end
-            %%%apply
-            volTs(I) = applyMask(volTs(I),mask);
         end
+        % if isfield(volAnat,'fun') && isfield(volAnat.fun,'mask') && isfield(volAnat.fun.mask,'crop') && ~isempty(volAnat.fun.mask.crop.vol)
+        %     volTs = applyMask(volTs,volAnat.fun.mask.crop.vol);
+        % end
     end
-    % if isfield(volAnat,'fun') && isfield(volAnat.fun,'mask') && isfield(volAnat.fun.mask,'crop') && ~isempty(volAnat.fun.mask.crop.vol)
-    %     volTs = applyMask(volTs,volAnat.fun.mask.crop.vol);
-    % end
+
 end
 
 % % %% Add time
@@ -112,27 +140,64 @@
 % 
 % % volTs = vol2vec(volTs);
 
-
-param.nDummy = info.dummy;
+if length(volTs)==1 && isfield(volTs,'nDummy') && ~isempty(volTs.nDummy)
+    param.nDummy = volTs.nDummy;
+else
+    param.nDummy = info.dummy;
+end
 param.nDummyRemoved = param.nDummy;
 if volTs.tr/1000 ~= volTs.dsgn.dt
     volTs.tr/1000 
-    diff(volTs.dsgn.onsets)
-    volTs.dsgn.onsets / (volTs.tr/1000)
-    volTs.dsgn.onsets / 4
-    warning(strjoin({''...
+    if isfield(volTs.dsgn,'onsets')
+        diff(volTs.dsgn.onsets)
+        volTs.dsgn.onsets / (volTs.tr/1000)
+        volTs.dsgn.onsets / 4
+
+        warning(strjoin({''...
         ['volume TR   =  ' sprintf('%7.3f ',volTs.tr/1000) 'sec']...
         ['stim dt     =  ' sprintf('%7.3f ',volTs.dsgn.dt) 'sec']...
         ['stim onsets = [' sprintf('%7.3f ',volTs.dsgn.onsets) ']sec']...
         ['            = [' sprintf('%7.3f ',(volTs.dsgn.onsets / (volTs.tr/1000))) ']vol']...
         'Defaulting to deconvolution TR = 1sec'},newline))
-    param.trDecon = 1;
+        param.trDecon = 1;
+    else
+        diff(volTs.dsgn.onsetList)
+        volTs.dsgn.onsetList / (volTs.tr/1000)
+        volTs.dsgn.onsetList / volTs.dsgn.dt
+        % volTs.dsgn.onsetList / volTs.dsgn.dt - round(volTs.dsgn.onsetList / volTs.dsgn.dt)
+
+        tError = (volTs.dsgn.onsetList(end) / (volTs.tr/1000) - volTs.dsgn.onsetList(end) / (volTs.dsgn.dt)) * 1000;
+
+        warning(strjoin({''...
+        ['volume TR   =  ' sprintf('%7.3f ',volTs.tr) 'ms']...
+        ['stim dt     =  ' sprintf('%7.3f ',volTs.dsgn.dt*1000) 'ms']...
+        ['stim onsets = [' sprintf('%7.0f ',volTs.dsgn.onsetList*1000) ']ms']...
+        ['            = [' sprintf('%7.3f ',(volTs.dsgn.onsetList / (volTs.tr/1000))) ']vol']...
+        ''...
+        [num2str(tError) 'ms error by the last onset']},newline))
+
+        if abs(tError)<10
+            param.trDecon = volTs.dsgn.dt;
+        else
+            dbstack; error('timing problem here')
+        end
+        
+
+        
+    end
+    
 else
     param.trDecon = volTs.tr/1000; %volTs.tr/1000/4; %1
 end
 param.verbose = 1;
-force = 1;
-[files,fRun,fSes,fSes_echoCat,param] = getResp(volTs,volAnat,param,force);
+forceThis = force;
+switch info.dataSetLabel
+    case 'vsmDriven'
+        [files,fRun,fSes,fSes_echoCat,param] = getResp(volTs,[],param,forceThis);
+    otherwise
+        dbstack; error('double-check that');
+        [files,fRun,fSes,fSes_echoCat,param] = getResp(volTs,volAnat,param,forceThis);
+end
 volResp.ts = MRIread(files.resp.f{1});
 volResp.tsOnBase = MRIread(files.respOnBase.f{1});
 volResp.base = MRIread(files.base.f{1});