-refactor space-time-frequency analysis

-start refactoring ploting
Proulx-S · Aug 2, 2024 · b9c62d5 · b9c62d5
1 parent b2bb93c
commit b9c62d5
Show file tree

Hide file tree

Showing 9 changed files with 348 additions and 195 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
diff --git 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