diff --git a/addFreq.m b/addFreq.m
index 545509d..3d71a56 100644
--- a/addFreq.m
+++ b/addFreq.m
@@ -1,5 +1,7 @@
-function addFreq(H,onsets,ondurs)
-if ~exist('ondurs','var'); ondurs = []; end
+function addFreq(H,onsets,ondurs,freqFlag)
+if ~exist('ondurs','var');     ondurs = []; end
+if ~exist('freqFlag','var'); freqFlag = []; end
+if isempty(freqFlag);        freqFlag = 1; end
 if isempty(H)
     H = gca;
 else
@@ -17,32 +19,39 @@ function addFreq(H,onsets,ondurs)
     yLim = ylim;
     text(fStim,yLim(2),'fStim','HorizontalAlignment','right','VerticalAlignment','top','Color','r')
 end
-if ~isempty(ondurs) && length(unique(ondurs))==1
+if freqFlag && ~isempty(ondurs) && length(unique(ondurs))==1
     fOn = 1/ondurs(1);
     if any(contains(H.Title.String,{'spectrogram' 'coherogram'}))
         yline(fOn,'Color','g','linestyle',':')
         text(xLim(2),fOn,'fOn','HorizontalAlignment','right','VerticalAlignment','top','Color','g')
-        yline(fOn/2,'Color','g','linestyle',':')
-        text(xLim(2),fOn/2,'fOn','HorizontalAlignment','right','VerticalAlignment','top','Color','g')
+        if freqFlag>1
+            yline(fOn/2,'Color','g','linestyle',':')
+            text(xLim(2),fOn/2,'fOn/2','HorizontalAlignment','right','VerticalAlignment','top','Color','g')
+        end
     elseif any(contains(H.Title.String,{'spectrum'}))
         xline(fOn,'Color','g','linestyle',':')
         text(fOn,yLim(2),'fOn','HorizontalAlignment','left','VerticalAlignment','top','Color','g')
-        xline(fOn/2,'Color','g','linestyle',':')
-        text(fOn/2,yLim(2),'fOn','HorizontalAlignment','left','VerticalAlignment','top','Color','g')
+        if freqFlag>1
+            xline(fOn/2,'Color','g','linestyle',':')
+            text(fOn/2,yLim(2),'fOn/2','HorizontalAlignment','left','VerticalAlignment','top','Color','g')
+        end
     end
 end
-if ~isempty(ondurs)
+
+if freqFlag>1 && ~isempty(ondurs)
     offsetList = onsets + ondurs;
     offdurList = onsets(2:end) - offsetList(1:end-1);
-    if length(unique(offdurList))==1
-        fOff = 1/offdurList(1);
-        if any(contains(H.Title.String,{'spectrogram' 'coherogram'}))
-            yline(fOff,'Color','b','linestyle',':')
-            text(xLim(2),fOff,'fOff','HorizontalAlignment','right','VerticalAlignment','top','Color','b')
-        elseif any(contains(H.Title.String,{'spectrum'}))
-            xline(fOff,'Color','b','linestyle',':')
-            text(fOff,yLim(2),'fOff','HorizontalAlignment','left','VerticalAlignment','top','Color','b')
-        end
+    fOff = 1/mean(offdurList);
+    if any(contains(H.Title.String,{'spectrogram' 'coherogram'}))
+        yline(fOff,'Color','b','linestyle',':')
+        text(xLim(2),fOff,'fOff','HorizontalAlignment','right','VerticalAlignment','top','Color','b')
+        yline(fOff,'Color','b','linestyle',':')
+        text(xLim(2),fOff/2,'fOff/2','HorizontalAlignment','right','VerticalAlignment','top','Color','b')
+    elseif any(contains(H.Title.String,{'spectrum'}))
+        xline(fOff,'Color','b','linestyle',':')
+        text(fOff,yLim(2),'fOff','HorizontalAlignment','left','VerticalAlignment','top','Color','b')
+        xline(fOff/2,'Color','b','linestyle',':')
+        text(fOff/2,yLim(2),'fOff/2','HorizontalAlignment','left','VerticalAlignment','top','Color','b')
     end
 end
 
diff --git a/getResp.m b/getResp.m
index 2cc84c7..9a2f2c9 100644
--- a/getResp.m
+++ b/getResp.m
@@ -885,12 +885,13 @@
         %% Define files
         fIn = fList{I,E};
         fOut = fIn;
-        fStat = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.nii.gz');
-        fMat = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.xmat.1D');
-        fStim = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_startTime.1D');
-        fResp = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_resp.nii.gz');
-        fFit = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_fit.nii.gz');
-        fResid = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_resid.nii.gz');
+        fStat   = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.nii.gz');
+        fMat    = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.xmat.1D');
+        fMatFig = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.xmat.fig');
+        fStim   = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_startTime.1D');
+        fResp   = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_resp.nii.gz');
+        fFit    = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_fit.nii.gz');
+        fResid  = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_resid.nii.gz');
 
         fRun(I,E).fResp = fResp;
         fRun(I,E).fFit = fFit;
@@ -917,7 +918,7 @@
                     dbstack; error('please specify deconvolution time resolution (param.funDsgn.trDecon)');
                 end
             end
-            [tmp,n] = afniCmd(fIn,fMask,fStim,param.nDummy,startSeq,condSeq,dt,param.funDsgn.label,TENTzeroFlag,fResp,fFit,fResid,fMat,fStat,verbose);
+            [tmp,n] = afniCmd(fIn,fMask,fStim,param.nDummy,startSeq,condSeq,dt,param.funDsgn.label,TENTzeroFlag,fResp,fFit,fResid,fMat,fStat,verbose,param.nDummyRemoved);
             cmdTmp = [cmdTmp tmp];
 
             cmdTmp{end+1} = ['echo ''   ''' fResp];
@@ -937,65 +938,77 @@
     end
 end
 
-%% Extract and plot design matrix
-if verbose>1
-    mri = MRIread(fIn,1); trMri = mri.tr/1000; nFrame = mri.nframes-param.nDummy; clear mri
-    cmdX = cmdTmp;
-    cmdX{contains(cmdTmp,'-input ')} = ['-nodata ' num2str(nFrame) ' ' num2str(trMri,'%f') ' \'];
-    ind = find(contains(cmdX,' \'),1,'last');
-    cmdX = [cmdX(1:ind) {'-x1D_stop \'} cmdX(ind+1:end)];
-    % cmdX = [cmdX(1:end-2) {'-x1D_stop \'} cmdX(end-1:end)];
-    cmdX{end} = [cmdX{end}];
-    % [status,cmdout] = system(strjoin(cmdX,newline));
-    [status,cmdout] = system(strjoin(cmdX,newline),'-echo'); if status || isempty(cmdout) || contains(cmdout,'ERROR'); dbstack; error(cmdout); error('x'); end
-    cmdX = {srcAfni};
-    cmdX{end+1} = ['1dcat ' fMat];
-    % [status,cmdout] = system(strjoin(cmdX,newline));
-    [status,cmdout] = system(strjoin(cmdX,newline),'-echo'); if status || isempty(cmdout) || contains(cmdout,'ERROR'); dbstack; error(cmdout); error('x'); end
-    param.funDsgn.mat = str2num(cmdout);
-    figure('WindowStyle','docked');
-    pInd = false([1 size(param.funDsgn.mat,2)]);
-    if TENTzeroFlag
-        pInd(1:(size(param.funDsgn.mat,2) - n + 2)) = true;
-    else
-        pInd(1:(size(param.funDsgn.mat,2) - n)) = true;
-    end
-    % imagesc(param.funDsgn.mat(:,pInd:end)); colormap gray
-    imagesc(param.funDsgn.mat); colormap gray
-    ax = gca;
-    ax.XTick = 1:size(param.funDsgn.mat,2);
-    ax.XTickLabel(pInd) = cell(1,nnz(pInd));
-    ax.XTickLabel(~pInd) = cellstr(num2str(((0:nnz(~pInd)-1)*dt)','%0.1f'));
-    clim([-1 1])
-end
+% %% Extract and plot design matrix
+% if verbose>1
+%     mri = MRIread(fIn,1); trMri = mri.tr/1000; nFrame = mri.nframes-(param.nDummy - param.nDummyRemoved); clear mri
+%     cmdX = cmdTmp;
+%     cmdX{contains(cmdTmp,'-input ')} = ['-nodata ' num2str(nFrame) ' ' num2str(trMri,'%f') ' \'];
+%     ind = find(contains(cmdX,' \'),1,'last');
+%     cmdX = [cmdX(1:ind) {'-x1D_stop \'} cmdX(ind+1:end)];
+%     % cmdX = [cmdX(1:end-2) {'-x1D_stop \'} cmdX(end-1:end)];
+%     cmdX{end} = [cmdX{end}];
+%     % [status,cmdout] = system(strjoin(cmdX,newline));
+%     [status,cmdout] = system(strjoin(cmdX,newline),'-echo'); if status || isempty(cmdout) || contains(cmdout,'ERROR'); dbstack; error(cmdout); error('x'); end
+%     cmdX = {srcAfni};
+%     cmdX{end+1} = ['1dcat ' fMat];
+%     % [status,cmdout] = system(strjoin(cmdX,newline));
+%     [status,cmdout] = system(strjoin(cmdX,newline),'-echo'); if status || isempty(cmdout) || contains(cmdout,'ERROR'); dbstack; error(cmdout); error('x'); end
+%     param.funDsgn.mat = str2num(cmdout);
+%     figure('WindowStyle','docked');
+%     pInd = false([1 size(param.funDsgn.mat,2)]);
+%     if TENTzeroFlag
+%         pInd(1:(size(param.funDsgn.mat,2) - n + 2)) = true;
+%     else
+%         pInd(1:(size(param.funDsgn.mat,2) - n)) = true;
+%     end
+%     % imagesc(param.funDsgn.mat(:,pInd:end)); colormap gray
+%     imagesc(param.funDsgn.mat); colormap gray
+%     ax = gca;
+%     ax.XTick = 1:size(param.funDsgn.mat,2);
+%     ax.XTickLabel(pInd) = cell(1,nnz(pInd));
+%     ax.XTickLabel(~pInd) = cellstr(num2str(((0:nnz(~pInd)-1)*dt)','%0.3f'));
+%     clim([-1 1])
+% end
 
 %% Run afni command
 cmd = strjoin(cmd,newline); % disp(strjoin(cmd,newline))
 [status,cmdout] = system(cmd,'-echo'); if status || isempty(cmdout) || contains(cmdout,'ERROR'); dbstack; error(cmdout); error('x'); end
 
 % Extract and plot design matrix
+% dt = trStim;
+hMat = figure('Visible','off');
+cmdX = {srcAfni};
+cmdX{end+1} = ['1dcat ' fMat];
+[status,cmdout] = system(strjoin(cmdX,newline));
+param.funDsgn.mat = str2num(cmdout);
+pInd = false([1 size(param.funDsgn.mat,2)]);
+if TENTzeroFlag
+    pInd(1:(size(param.funDsgn.mat,2) - n + 2)) = true;
+else
+    pInd(1:(size(param.funDsgn.mat,2) - n)) = true;
+end
+imagesc(param.funDsgn.mat); colormap gray
+ax = gca;
+ax.XTick = 1:size(param.funDsgn.mat,2);
+ax.XTickLabel(pInd) = cell(1,nnz(pInd));
+ax.XTickLabel(~pInd) = cellstr(num2str(((0:nnz(~pInd)-1)*dt)','%0.3f'));
+clim([-1 1])
+set(hMat, 'CreateFcn', 'set(gcbo,''Visible'',''on'')'); 
+savefig(hMat,fMatFig,'compact')
 if verbose>1
-    % dt = trStim;
-    cmdX = {srcAfni};
-    cmdX{end+1} = ['1dcat ' fMat];
-    [status,cmdout] = system(strjoin(cmdX,newline));
-    param.funDsgn.mat = str2num(cmdout);
-    figure('WindowStyle','docked');
-    imagesc(param.funDsgn.mat); colormap gray
-    ax = gca;
-    ax.XTick = 1:size(param.funDsgn.mat,2);
-    ax.XTickLabel(pInd) = cell(1,nnz(pInd));
-    ax.XTickLabel(~pInd) = cellstr(num2str(((0:nnz(~pInd)-1)*dt)','%0.1f'));
-    clim([-1 1])
+    hMat.Visible = 'on';
+    hMat.WindowStyle = 'docked';
 end
 
 
 
+
 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Fixed-effect model on concatenated runs (separate baselines) %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
 cmd = {srcAfni};
 if size(fList,1)>1
+    dbstack; error('double-check that')
     for E = 1:size(fList,2)
 
         %% Define files
@@ -1005,6 +1018,7 @@
 
         fStat = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.nii.gz');
         fMat = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.xmat.1D');
+        fMatFig = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_stats.xmat.fig');
         fStim = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_startTime.1D');
         fResp = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_resp.nii.gz');
         fFit = fullfile(fileparts(replace(fOut,'.nii.gz','')),'cond-visOn_fit.nii.gz');
@@ -1068,13 +1082,24 @@
 end
 
 
-function [cmd,n] = afniCmd(fIn,fMask,fStim,nDummy,startSeq,condSeq,dt,label,TENTzeroFlag,fResp,fFit,fResid,fMat,fStat,verbose)
+function [cmd,n] = afniCmd(fIn,fMask,fStim,nDummy,startSeq,condSeq,dt,label,TENTzeroFlag,fResp,fFit,fResid,fMat,fStat,verbose,nDummyRemoved)
+if ~exist('nDummyRemoved','var'); nDummyRemoved = []; end
+if isempty(nDummyRemoved);        nDummyRemoved = nDummy; warning('param.nDummyRemoved not specified, assuming param.nDummyRemoved = param.nDummy'); end
+if nDummyRemoved && nDummyRemoved~=nDummy; dbstack; error('param.nDummyRemoved specified, but does not match param.nDummy'); end
 cmd = {'3dDeconvolve -overwrite \'};
 % cmdTmp{end+1} = ['-force_TR ' num2str(trStim) ' \'];
 if iscell(fIn)
-    cmd{end+1} = ['-input ' strjoin(fIn,' ') '[' num2str(nDummy) '..$] \'];
+    if nDummyRemoved
+        cmd{end+1} = ['-input ' strjoin(fIn,' ') ' \'];
+    else
+        cmd{end+1} = ['-input ' strjoin(fIn,' ') '[' num2str(nDummy) '..$] \'];
+    end
 else
-    cmd{end+1} = ['-input ' fIn '[' num2str(nDummy) '..$] \'];
+    if nDummyRemoved
+        cmd{end+1} = ['-input ' fIn ' \'];
+    else
+        cmd{end+1} = ['-input ' fIn '[' num2str(nDummy) '..$] \'];
+    end
 end
 if ~isempty(fMask)
     cmd{end+1} = ['-mask ' fMask ' \'];
@@ -1096,12 +1121,17 @@
 
 % write design to file
 fido = fopen(fStim, 'w');
-for i = 1:length(fIn)
+if iscell(fIn)
+    for i = 1:length(fIn)
+        fprintf(fido,'%.3f ',startSeq(condSeq==1));
+        fprintf(fido,'\n');
+    end
+else
     fprintf(fido,'%.3f ',startSeq(condSeq==1));
     fprintf(fido,'\n');
 end
 fclose(fido);
-deconWin = mode(diff(startSeq(condSeq==1)));
+deconWin = min(diff(startSeq(condSeq==1)));
 deconWin = floor(deconWin/dt)*dt;
 b = 0;
 c = deconWin-dt;
diff --git a/plotSpec.m b/plotSpec.m
index b31dc25..de16118 100644
--- a/plotSpec.m
+++ b/plotSpec.m
@@ -1,13 +1,16 @@
-function [ax,F] = plotSpec(volPsd,metricLabel,H,tWin)
+function [ax,F] = plotSpec(volPsd,metricLabel,H,tWin,volTs,respQthresh)
 if ~exist('H','var'); H = [];                             end
 if isempty(H);        H = figure('WindowStyle','docked'); end
 
 if isfield(volPsd,'dsgn') && isfield(volPsd.dsgn,'f0'); f0 = volPsd.dsgn.f0; else; f0 = []; end
 
-if ~exist('tWin','var'); tWin = [];                       end
-if ~exist('metricLabel','var'); metricLabel = [];                             end
-if isempty(metricLabel);        metricLabel = 'psd'; end % 'psd' 'coh'
-
+if ~exist('tWin','var');                                     tWin = []; end
+if ~exist('metricLabel','var');                       metricLabel = []; end
+if ~exist('volTs','var');                                   volTs = []; end
+if ~exist('respQthresh','var');                       respQthresh = []; end
+if isempty(metricLabel);                              metricLabel = 'psd'; end % 'psd' 'coh'
+if isempty(respQthresh) && strcmp(metricLabel,'psd'); respQthresh = 0.05; end % 'psd' 'coh'
+threshAvFlag = ismember(metricLabel,{'psd' 'psdEPC'}) && ~isempty(volTs) && isfield(volTs,'resp') && ~isempty(respQthresh) && respQthresh~=inf && respQthresh~=1;
 
 switch class(H)
     case 'matlab.graphics.layout.TiledChartLayout'
@@ -40,7 +43,11 @@
 
 
 %%% average across space
-vec = squeeze(mean(vec,6));
+if threshAvFlag
+    vec = squeeze(mean(vec(:,:,:,:,:,volTs.resp.Fq.vol(volPsd.vol2vec)<=respQthresh,:,1),6));
+else
+    vec = squeeze(mean(vec(:,:,:,:,:,:,:,1),6));
+end
 f = squeeze(mt.f);
 % f   = squeeze(volPsd.f);
 % psd = squeeze(mean(volPsd.vec,2));
@@ -63,8 +70,17 @@
 T = mt.T;
 [TW,W] = K2W(T,K,0);
 
-paramStr = ['(K=' num2str(K) '; 2W=' num2str(W*2,'%0.4f') 'Hz; T=' num2str(T,'%0.2f') 'sec); TW=' num2str(TW)];
-title(['spatially averaged spectrum ' paramStr])
+paramStr = ['(K=' num2str(K) '; 2W=' num2str(W*2,'%0.4f') 'Hz; T=' num2str(T,'%0.2f') 'sec; TW=' num2str(TW) ')'];
+if threshAvFlag
+    title(['spectrum ' paramStr ' averaged across voxels with Q<=' num2str(respQthresh,'%0.2f')])
+else
+    switch metricLabel
+        case 'psd'
+            title(['spectrum ' paramStr ' averaged across voxels'])
+        case 'coh'
+            title(['coherence spectrum ' paramStr])
+    end
+end
 
 yLim = vec(f>0.01);
 yLim = [min(yLim) max(yLim)];
@@ -76,7 +92,7 @@
 
 
 if isfield(volPsd,'psdTrialGram') && ~isempty(volPsd.psdTrialGram)
-    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList)
+    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList,2)
 end
 
 if ~isempty(f0)
diff --git a/plotSpecAll.m b/plotSpecAll.m
new file mode 100644
index 0000000..9b04a0f
--- /dev/null
+++ b/plotSpecAll.m
@@ -0,0 +1,76 @@
+function plotSpecAll(volPsd,volTs,respQthresh)
+
+if ~exist('respQthresh','var'); respQthresh = []; end
+if isempty(respQthresh);        respQthresh = 0.05; end % 'psd' 'coh'
+
+Fpsd = figure('WindowStyle','docked');
+HtPsd = tiledlayout(8,1); HtPsd.TileSpacing = 'tight'; HtPsd.Padding = 'tight';
+axPsd = {};
+axPsd{end+1} = plotTs2(volTs,HtPsd);
+axPsd{end+1} = plotSpec(volPsd,'psd',HtPsd,[],volTs);
+axPsd{end+1} = plotSpecGram(volPsd,'gram'       ,'psd'   ,HtPsd,volTs,respQthresh);
+axPsd{end+1} = plotSpecGram(volPsd,'trialGram'  ,'psd'   ,HtPsd,volTs,respQthresh);
+axPsd{end+1} = plotSpecGram(volPsd,'trialGram'  ,'psdEPC',HtPsd,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);
+
+[~,b,~] = fileparts(replace(volPsd.fspec,'.nii.gz',''));
+b = strsplit(b,'_');
+title(HtPsd,b{1})
+
+
+Fcoh = figure('WindowStyle','docked');
+HtCoh = tiledlayout(8,1); HtCoh.TileSpacing = 'tight'; HtCoh.Padding = 'tight';
+axCoh = {};
+axCoh{end+1} = plotTs2(volTs,HtCoh,[],[],[],[],respQthresh);
+if isfield(volPsd,'svd') && isfield(volPsd.svd,'COH') && ~isempty(volPsd.svd.COH)
+    axCoh{end+1} = plotSpec(volPsd,'coh',HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdGram')        && isfield(volPsd.svdGram,'COH') && ~isempty(volPsd.svdGram.COH)
+    axCoh{end+1} = plotSpecGram(volPsd,'gram'       ,'coh'   ,HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdTrialGram')   && isfield(volPsd.svdTrialGram,'vec')   && isfield(volPsd.svdTrialGram.vec,'coh')    && ~isempty(volPsd.svdTrialGram.vec.coh)
+    axCoh{end+1} = plotSpecGram(volPsd,'trialGram'  ,'coh'   ,HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdTrialGram')   && isfield(volPsd.svdTrialGram,'vec')   && isfield(volPsd.svdTrialGram.vec,'cohEPC') && ~isempty(volPsd.svdTrialGram.vec.cohEPC)
+    axCoh{end+1} = plotSpecGram(volPsd,'trialGram'  ,'cohEPC',HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdTrialGram')   && isfield(volPsd.svdTrialGram,'vec')   && isfield(volPsd.svdTrialGram.vec,'cohEK')  && ~isempty(volPsd.svdTrialGram.vec.cohEK)
+    axCoh{end+1} = plotSpecGram(volPsd,'trialGram'  ,'cohEK' ,HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdTrialGramMD') && isfield(volPsd.svdTrialGramMD,'vec') && isfield(volPsd.svdTrialGramMD.vec,'coh')     && ~isempty(volPsd.svdTrialGram.vec.coh)
+    axCoh{end+1} = plotSpecGram(volPsd,'trialGramMD','coh'   ,HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+if isfield(volPsd,'svdTrialGramMD') && isfield(volPsd.svdTrialGramMD,'vec') && isfield(volPsd.svdTrialGramMD.vec,'cohEPC')  && ~isempty(volPsd.svdTrialGram.vec.cohEPC)
+    axCoh{end+1} = plotSpecGram(volPsd,'trialGramMD','cohEPC',HtCoh);
+else
+    axCoh{end+1} = nexttile;
+end
+
+ind = 3:6; indX = false(size(ind)); for i = 1:length(ind); indX(i) = ~isempty(findobj(axCoh{ind(i)}.Children,'Type','Image')); end; ind = ind(indX);
+if length(ind)>1
+    cLim = get([axCoh{ind}],'CLim'); cLim = cat(1,cLim{:}); cLim = [min(cLim(:)) max(cLim(:))];
+    set([axCoh{ind}],'CLim',cLim);
+end
+
+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})
diff --git a/plotSpecGram.m b/plotSpecGram.m
index 4b19587..8487d37 100644
--- a/plotSpecGram.m
+++ b/plotSpecGram.m
@@ -1,4 +1,4 @@
-function [ax,F] = plotSpecGram(volPsd,timeLabel,metricLabel,H)
+function [ax,F] = plotSpecGram(volPsd,timeLabel,metricLabel,H,volTs,respQthresh)
 if ~exist('H','var'); H = [];                             end
 if isempty(H);        H = figure('WindowStyle','docked'); end
 
@@ -6,9 +6,12 @@
 
 if ~exist('timeLabel','var');     timeLabel = [];                             end
 if ~exist('metricLabel','var'); metricLabel = [];                             end
+if ~exist('volTs','var');             volTs = []; end
+if ~exist('respQthresh','var'); respQthresh = []; end
 if isempty(timeLabel);            timeLabel = 'gram'; end % 'gram' 'trialGram' 'trialGramMD'
 if isempty(metricLabel);        metricLabel = 'psd'; end % 'psd' 'psdEPC' 'coh' 'cohEPC' 'cohEK'
-
+if isempty(respQthresh);        respQthresh = 0.05; end % 'psd' 'coh'
+threshAvFlag = ismember(metricLabel,{'psd' 'psdEPC'}) && ~isempty(volTs) && isfield(volTs,'resp') && ~isempty(respQthresh) && respQthresh~=inf && respQthresh~=1;
 
 switch class(H)
     case 'matlab.graphics.layout.TiledChartLayout'
@@ -98,7 +101,11 @@
 
 
 %% spatial average
-vec = permute(mean(vec(:,:,:,:,:,:,:,1),6),[5 7 2 8 1 3 4 6]);
+if threshAvFlag
+    vec = permute(mean(vec(:,:,:,:,:,volTs.resp.Fq.vol(volPsd.vol2vec)<=respQthresh,:,1),6),[5 7 2 8 1 3 4 6]);
+else
+    vec = permute(mean(vec(:,:,:,:,:,:,:,1),6),[5 7 2 8 1 3 4 6]);
+end
 f   = permute(mt.f       ,[5 7 2 8 1 3 4 6]);
 Fs  = mt.param.Fs;
 t   = permute(mean(mt.t(:,1,:,:,:,:,:,1),1),[5 7 2 8 1 3 4 6]);
@@ -126,12 +133,16 @@
 % if contains(timeLabel,'trial')
 %     mt.t(:,1,:,:,:,:,:)
 % else
-runDur = volPsd.nframes/mt.param.Fs;
+runDur = volPsd.t(end);
 % end
 xlim([0 runDur])
 
 paramStr = ['(K=' num2str(K) '; 2W=' num2str(W*2,'%0.4f') 'Hz; T=' num2str(T,'%0.2f') 'sec; TW=' num2str(TW) ')'];
-title([label ' ' paramStr])
+if threshAvFlag
+    title([label ' ' paramStr ' averaged across voxels with Q<=' num2str(respQthresh,'%0.2f')])
+else
+    title([label ' ' paramStr])
+end
 
 
 cLim = vec(f>0.01,:);
@@ -148,7 +159,7 @@
 
 
 if isfield(volPsd,'psdTrialGram') && ~isempty(volPsd.psdTrialGram)
-    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList)
+    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList,1)
 end
 
 if ~isempty(f0)
diff --git a/plotTs2.m b/plotTs2.m
index b0d6f44..f99a7c2 100644
--- a/plotTs2.m
+++ b/plotTs2.m
@@ -1,18 +1,21 @@
-function [ax,F] = plotTs2(volTs,H,onsets,ondurs,volAnat,roiInd)
-if ~exist('H','var');            H = []; end
-if isempty(H);                   H = figure('WindowStyle','docked'); end
-if ~exist('onsets','var');  onsets = []; end
-if ~exist('ondurs','var');  ondurs = []; end
-if ~exist('roiInd','var');  roiInd = []; end
-if ~exist('volAnat','var'); volRoi = [];
-                           volMask = [];
+function [ax,F] = plotTs2(volTs,H,onsets,ondurs,volAnat,roiInd,respQthresh)
+if ~exist('H','var');                      H = []; end
+if isempty(H);                             H = figure('WindowStyle','docked'); end
+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('volAnat','var');          volAnat = []; end
+if isempty(respQthresh);         respQthresh = 0.05; end
+if isempty(volAnat);                  volRoi = [];
+                                     volMask = [];
 elseif ~isempty(roiInd)
-                            volRoi = volAnat; clear volAnat
-                           volMask = [];
+                                      volRoi = volAnat; clear volAnat
+                                     volMask = [];
 else
-                            volRoi = [];
-                           volMask = volAnat; clear volAnat
-                   volMask.mri.vol = any(volMask.mri.vol,4);
+                                      volRoi = [];
+                                     volMask = volAnat; clear volAnat
+                             volMask.mri.vol = any(volMask.mri.vol,4);
 end
 
 
@@ -54,7 +57,11 @@
 %     ts = volTs.vec;
 %     ts = mean(ts(:,q<0.001),2);
 % else
+if isfield(volTs,'resp') && (respQthresh~=inf || respQthresh~=1)
+    ts = squeeze(mean(volTs.vec(:,volTs.resp.Fq.vol(volTs.vol2vec)<=respQthresh),2));
+else
     ts = squeeze(mean(volTs.vec,2));
+end
 % end
 plot(squeeze(volTs.t),ts,'k')
 grid on
@@ -64,10 +71,15 @@
 
 
 %% Add response if available
-if isfield(volTs,'resp')
+
+if ~isempty(volTs) && isfield(volTs,'resp')
     hold on
     ts = permute(volTs.resp.ts.vol,[4 1 2 3]);
-    ts = mean(ts(:,volTs.vol2vec),2);
+    if respQthresh~=inf && respQthresh~=1
+        ts = mean(ts(:,volTs.resp.Fq.vol<=respQthresh),2);
+    else
+        ts = mean(ts(:,volTs.vol2vec),2);
+    end
     % %%%% weighted average
     % weights = volTs.resp.F.vol;
     % tmp = weights(weights~=0); tmp = tmp - min(tmp); tmp = tmp./mean(tmp);
@@ -112,7 +124,12 @@
 elseif ~isempty(volRoi)
     dbstack; error('code that');
 end
-title(['timeseries (' strjoin(paramStr,'; ') ')'])
+
+if ~isempty(volTs) && isfield(volTs,'resp') && (respQthresh~=inf && respQthresh~=1)
+    title(['timeseries (' strjoin(paramStr,'; ') ') averaged across voxels with Q<=' num2str(respQthresh,'%0.2f')])
+else
+    title(['timeseries (' strjoin(paramStr,'; ') ')'])
+end
 
 
 
diff --git a/volPsdFullMt2.m b/volPsdFullMt2.m
index 91725ed..f33edad 100644
--- a/volPsdFullMt2.m
+++ b/volPsdFullMt2.m
@@ -1,4 +1,4 @@
-function [out, info] = volPsdFullMt2(do,info,volTs,volAnat)
+function [out, info] = volPsdFullMt2(do,info,volTs,volAnat,volPsd)
 if isempty(do)
     do.loadIt = 0;
     do.doIt   = 1;
@@ -6,6 +6,7 @@
 end
 
 if ~exist('volAnat','var'); volAnat = []; end
+if ~exist('volPsd','var');   volPsd = []; end % for purpose of a second-step analysis using a mask derived from the first-step analysis
 
 if ~isfield(info,'K');                   info.K = []        ; end
 if ~isfield(info,'win');               info.win = zeros(0,2); end
@@ -81,6 +82,18 @@
                 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);
 
@@ -201,7 +214,7 @@
 MRIwrite(tmp,volPsd.psd.fspec);
 
 %%% Coherence (first singular vector)
-if isfield(volPsd.svd,'spSV')
+if isfield(volPsd.svd,'spSV') && ~isempty(volPsd.svd.spSV)
     volPsd.vec = permute(abs(volPsd.svd.spSV(:,:,:,:,:,:,:,1)),[5 6 8 1 2 3 4 7]);
     volPsd.nframes = size(volPsd.vec,1);
     volPsd.tr = mean(diff(volPsd.svd.f))*1000;
diff --git a/volTsGetResp.m b/volTsGetResp.m
index 1d1bd60..edff5b2 100644
--- a/volTsGetResp.m
+++ b/volTsGetResp.m
@@ -114,7 +114,22 @@
 
 
 param.nDummy = info.dummy;
-param.trDecon = 1;
+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({''...
+        ['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;
+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);