diff --git a/addWin.m b/addWin.m
index 936a204..86baa62 100644
--- a/addWin.m
+++ b/addWin.m
@@ -13,22 +13,34 @@ function addWin(H,gram,x,y)
 Fs = gram.param.Fs;
 winSz = gram.win(1) + 0.5/Fs;
 if contains(H.Title.String,'timeseries')
-    yLim = ylim;
-    xLim = xlim;
-    if isempty(x)
-        x = xLim(1); x = x.*[1 1];
-        xx = x + [0 1].*winSz;
+    if size(gram.t,2)>1
+        y = ylim;
+        xx = gram.t(:,:,:,:,:,:,1);
+        yy = repmat(y(1),[size(xx,1) 1]);
+        line(xx,yy,'Color','m','linewidth',3)
+
+        xx = gram.t(:,:,:,:,:,:,end);
+        yy = repmat(y(2),[size(xx,1) 1]);
+        line(xx,yy,'Color','m','linewidth',3)
     else
-        x = x.*[1 1];
-        xx = x + [-0.5 0.5].*winSz;
-    end
-    if isempty(y) || y==-inf
-        y = yLim(1);
-    elseif y==inf
-        y = yLim(2);
+        yLim = ylim;
+        xLim = xlim;
+        if isempty(x)
+            x = xLim(1); x = x.*[1 1];
+            xx = x + [0 1].*winSz;
+        else
+            x = x.*[1 1];
+            xx = x + [-0.5 0.5].*winSz;
+        end
+        if isempty(y) || y==-inf
+            y = yLim(1);
+        elseif y==inf
+            y = yLim(2);
+        end
+        y = y.*[1 1];
+        line(xx,y,'Color','m','linewidth',3)
     end
-    y = y.*[1 1];
-    line(xx,y,'Color','m','linewidth',3)
+
 elseif contains(H.Title.String,'spectrogram') || contains(H.Title.String,'coherogram')
     im = findobj(H.Children,'type','image');
     if isempty(x)
diff --git a/getJ2.m b/getJ2.m
index 320656b..764bc44 100644
--- a/getJ2.m
+++ b/getJ2.m
@@ -12,10 +12,16 @@
 
 
 tp2 = reshape(  tp .* exp(-f.*t*2*pi*1i)  ,[N Et*R*K*F*1*Wt]);
-d2  = reshape(  d - mean(d,1)             ,[N  E*R*1*1*V*W ]); % removing the mean here
+d2  = reshape(  d - mean(d,1)             ,[N E*R*1*1*V*W]); % removing the mean here
 
 d2 = permute(d2,[2 1]);
 J = d2*tp2; % [V E*R*K*F*1*W]
+% [E*R*1*1*V*W Et*R*K*F*1*Wt]
 
-J = reshape(J,[V E R K F 1 W]); % [V E R K F 1 W]
-J = permute(J,[6 2 3 4 5 1 7]); % [N E R K F V W]
\ No newline at end of file
+J = reshape(J,[E*R*1*1*V*W Et R K F 1 Wt]);
+J = permute(J,[2 3 4 5 6 7 1]);
+J = reshape(J,[Et R K F 1 Wt E R 1 1 V W]);
+J = permute(J,[13 7 2 3 4 11 6 1 5 8 9 10 12]); %[N ,E ,R ,K,F,V,W]
+if any(size(J,9:15)~=1); dbstack; error('X'); end
+% J = reshape(J,[V E R K F 1 W]); % [V E R K F 1 W]
+% J = permute(J,[6 2 3 4 5 1 7]); % [N E R K F V W]
\ No newline at end of file
diff --git a/plotPsd2.m b/plotPsd2.m
index 8f544bd..24fc748 100644
--- a/plotPsd2.m
+++ b/plotPsd2.m
@@ -1,7 +1,9 @@
 function [ax,F] = plotPsd2(volPsd,H,tWin)
 if ~exist('H','var'); H = [];                             end
 if isempty(H);        H = figure('WindowStyle','docked'); end
-if ~exist('tWin','var'); tWin = [];                             end
+if ~exist('tWin','var'); tWin = [];                       end
+if isfield(volPsd,'dsgn') && isfield(volPsd.dsgn,'f0'); f0 = volPsd.dsgn.f0; else; f0 = []; end
+
 
 switch class(H)
     case 'matlab.graphics.layout.TiledChartLayout'
@@ -52,6 +54,10 @@
     addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList)
 end
 
+if ~isempty(f0)
+    xline(f0,'--b')
+end
+
 if isfield(volPsd,'psdTrialGramMD') && ~isempty(volPsd.psdTrialGramMD) && ~isempty(tWin)
     t   = volPsd.psdTrialGramMD.t(1,1,1,1,1,1,:,1);
     f   = volPsd.psdTrialGramMD.f(1,1,1,1,:,1,1,1);
diff --git a/plotSpec.m b/plotSpec.m
new file mode 100644
index 0000000..b31dc25
--- /dev/null
+++ b/plotSpec.m
@@ -0,0 +1,104 @@
+function [ax,F] = plotSpec(volPsd,metricLabel,H,tWin)
+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'
+
+
+switch class(H)
+    case 'matlab.graphics.layout.TiledChartLayout'
+        F = H.Parent;
+    case 'matlab.ui.Figure'
+        F = H;
+    otherwise
+end
+
+figure(F);
+ax = {};
+ax{end+1} = nexttile;
+
+
+%% Select approrpiate data
+switch metricLabel
+    case 'psd'
+        mt    = volPsd.psd;
+        vec   = mt.PSD;
+        label = 'spatially averaged spectrum';
+    case 'coh'
+        mt    = volPsd.svd;
+        vec   = mt.COH;
+        label = 'coherence spectrum';
+    otherwise
+        dbstack; error('code trhat')
+end
+
+
+
+
+%%% average across space
+vec = squeeze(mean(vec,6));
+f = squeeze(mt.f);
+% f   = squeeze(volPsd.f);
+% psd = squeeze(mean(volPsd.vec,2));
+plot(f,vec,'k')
+grid on
+axis tight
+xlabel('f (Hz)')
+switch metricLabel
+    case 'psd'
+        ylabel('psd')
+        ax{end}.YScale = 'log';
+    case 'coh'
+        ylabel('coherence')
+    otherwise
+        dbstack; error('code trhat')
+end
+
+
+K = mt.K;
+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])
+
+yLim = vec(f>0.01);
+yLim = [min(yLim) max(yLim)];
+ylim(yLim)
+xlim([0 f(end)])
+
+addW([],volPsd.psd)
+
+
+
+if isfield(volPsd,'psdTrialGram') && ~isempty(volPsd.psdTrialGram)
+    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList)
+end
+
+if ~isempty(f0)
+    xline(f0,'--b')
+end
+
+if isfield(volPsd,'psdTrialGramMD') && ~isempty(volPsd.psdTrialGramMD) && ~isempty(tWin)
+    dbstack; error('double-check that');
+    t   = volPsd.psdTrialGramMD.t(1,1,1,1,1,1,:,1);
+    f   = volPsd.psdTrialGramMD.f(1,1,1,1,:,1,1,1);
+    if tWin==inf
+        psd = mean(mean(volPsd.psdTrialGramMD.vec.psdPC(:,:,:,:,:,:,:,:),6),7);
+    else
+        [~,wInd] = min(abs(t-tWin));
+        psd = mean(volPsd.psdTrialGramMD.vec.psdPC(:,:,:,:,:,:,wInd,:),6);
+    end
+    plot(squeeze(f),squeeze(psd),'--k')
+end
+
+
+ax = [ax{:}];
+
+
+
+
diff --git a/plotSpecGram.m b/plotSpecGram.m
new file mode 100644
index 0000000..4b19587
--- /dev/null
+++ b/plotSpecGram.m
@@ -0,0 +1,202 @@
+function [ax,F] = plotSpecGram(volPsd,timeLabel,metricLabel,H)
+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('timeLabel','var');     timeLabel = [];                             end
+if ~exist('metricLabel','var'); metricLabel = [];                             end
+if isempty(timeLabel);            timeLabel = 'gram'; end % 'gram' 'trialGram' 'trialGramMD'
+if isempty(metricLabel);        metricLabel = 'psd'; end % 'psd' 'psdEPC' 'coh' 'cohEPC' 'cohEK'
+
+
+switch class(H)
+    case 'matlab.graphics.layout.TiledChartLayout'
+        F = H.Parent;
+    case 'matlab.ui.Figure'
+        F = H;
+    otherwise
+end
+
+figure(F);
+ax = {};
+ax{end+1} = nexttile;
+
+%% Select approrpiate data
+switch timeLabel
+    case 'gram'
+        switch metricLabel
+            case 'psd'
+                mt    = volPsd.psdGram;
+                vec   = mt.PSD;
+                label = 'spectrogram';
+            case 'coh'
+                mt    = volPsd.svdGram;
+                vec   = mt.COH;
+                label = 'coherogram';
+            otherwise
+                dbstack; error('code trhat')
+        end
+    case 'trialGram'
+    %     case 'av'
+    %     title(['evoked spectrogram ' paramStr])
+    % case 'pc'
+    %     title(['phase-locked evoked spectrogram ' paramStr])
+    % case 'MD'
+    %     title(['discontinuous-taper evoked spectrogram ' paramStr])
+    % case 'MDpc'
+    %     title(['discontinuous-taper phase-locked evoked spectrogram ' paramStr])
+        switch metricLabel
+            case 'psd'
+                mt    = volPsd.psdTrialGram;
+                vec   = mt.vec.psd;
+                label = 'evoked spectrogram';
+            case 'psdEPC'
+                mt    = volPsd.psdTrialGram;
+                vec   = mt.vec.psdPC;
+                label = 'phase-locked evoked spectrogram';
+            case 'coh'
+                mt    = volPsd.svdTrialGram;
+                vec   = mt.vec.coh;
+                label = 'evoked coherogram';
+            case 'cohEPC'
+                mt    = volPsd.svdTrialGram;
+                vec   = mt.vec.cohEPC;
+                label = 'phase-locked evoked coherogram';
+            case 'cohEK'
+                mt    = volPsd.svdTrialGram;
+                vec   = mt.vec.cohEK;
+                label = 'evoked cross-trial coherogram';
+            otherwise
+                dbstack; error('code trhat')
+        end
+    case 'trialGramMD'
+        switch metricLabel
+            case 'psd'
+                mt    = volPsd.psdTrialGramMD;
+                vec   = mt.vec.psd;
+                label = 'discontinuous-taper evoked spectrogram';
+            case 'psdEPC'
+                mt    = volPsd.psdTrialGramMD;
+                vec   = mt.vec.psdPC;
+                label = 'discontinuous-taper phase-locked evoked spectrogram';
+            case 'coh'
+                mt    = volPsd.svdTrialGramMD;
+                vec   = mt.vec.coh;
+                label = 'discontinuous-taper evoked coherogram';
+            case 'cohEPC'
+                mt    = volPsd.svdTrialGramMD;
+                vec   = mt.vec.cohEPC;
+                label = 'discontinuous-taper phase-locked evoked coherogram';                
+            otherwise
+                dbstack; error('code trhat')
+        end
+    otherwise
+        dbstack; error('code trhat')
+end
+
+
+
+%% spatial average
+vec = permute(mean(vec(:,:,:,:,:,:,:,1),6),[5 7 2 8 1 3 4 6]);
+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]);
+K   = mt.K;
+T   = mt.T;
+[TW,W] = K2W(T,K,0);
+
+imagesc(t,f,vec)
+
+xlabel('time (s)')
+ylabel('freq (Hz)')
+switch metricLabel
+    case {'psd' 'psdEPC'}
+        ylabel(colorbar,'psd');
+        ax{end}.ColorScale = 'log';
+    case {'coh' 'cohEPC' 'cohEK'}
+        ylabel(colorbar,'coherence');
+    otherwise
+        dbstack; error('code trhat')
+end
+
+
+
+
+% if contains(timeLabel,'trial')
+%     mt.t(:,1,:,:,:,:,:)
+% else
+runDur = volPsd.nframes/mt.param.Fs;
+% 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])
+
+
+cLim = vec(f>0.01,:);
+cLim = [min(cLim(:)) max(cLim(:))];
+clim(cLim)
+
+
+
+
+addWin([],mt)
+addW([],mt)
+
+
+
+
+if isfield(volPsd,'psdTrialGram') && ~isempty(volPsd.psdTrialGram)
+    addFreq([],volPsd.psdTrialGram.onsetList,volPsd.psdTrialGram.durList)
+end
+
+if ~isempty(f0)
+    yline(f0,'--b')
+end
+
+
+
+
+
+
+ax = [ax{:}];
+
+
+%% %%%%%%%%%%%%%%%%%%
+% Add to timeseries %
+%%%%%%%%%%%%%%%%%% %%
+
+axTs = findobj(allchild(F.Children),'type','axes'); ttl = get(axTs,'Title'); ttl = get([ttl{:}],'String');
+axTs = axTs(contains(ttl,'timeseries'));
+
+%%% delete previous window size visual elements (magentat lines)
+hLine = findobj(axTs.Children,'type','Line'); %hLine = {hLine(:)};
+mLine = get(hLine,'Color'); if iscell(mLine); mLine = cat(1,mLine{:}); end
+mLine = all(cat(1,mLine)==[1 0 1],2);
+delete(hLine(mLine));
+
+% % if length(hLine)==1
+% %     hLine = {hLine};
+% % end
+% % 
+% %     mLine = get(hLine,'Color');
+% % else
+% %     mLine = {get([axTs.Children(:)],'Color')};
+% % end
+% mLine = get(hLine,'Color'); if ~iscell(mLine); mLine = {mLine}; end
+% ind = all(cat(1,mLine{:})==[1 0 1],2);
+% mLine(ind)
+% 
+% mLine = hLine(all(cat(1,mLine{:})==[1 0 1],2));
+% delete(mLine);
+% % if length(axTs.Children)>1
+% %     mLine = get(hLine,'Color');
+% % else
+% %     mLine = {get([axTs.Children(:)],'Color')};
+% % end
+% % mLine = axTs.Children(all(cat(1,mLine{:})==[1 0 1],2));
+% % delete(mLine);
+
+%%% add window size
+addWin(axTs,mt)
diff --git a/runFullMT2.m b/runFullMT2.m
index af5ec50..7a6d193 100644
--- a/runFullMT2.m
+++ b/runFullMT2.m
@@ -84,7 +84,13 @@
 param.durList = durList;
 
 if param.win(1)==inf; skipGram = true; skipTrialGram = true; else skipGram = false; skipTrialGram = false; end
-if isempty(param.onsetList); skipTrialGram = true; end
+if isempty(param.onsetList)
+    skipTrialGram = true;
+else
+    if all((size(param.onsetList)==1)==[1 0])
+        dbstack; error('onsetList must be a column vector')
+    end
+end
 
 
 % if length(param.win)>2
@@ -522,7 +528,7 @@
 if isempty(nShuf);         nShuf = 0; end
 if isempty(nRun);           nRun = 1; end
 if isempty(cohFrange); cohFrange = [0 inf]; end
-orderWin = -1;
+orderWin = 0;
 % skip.gram = any(ismember(param.win(2),[0 inf nan]));
 for runInd = 1:nRun
     if verbose && nRun>1; disp(['---Run ' num2str(runInd) '/' num2str(nRun) '---']); end
@@ -823,7 +829,7 @@
         NFFT=max(2^(nextpow2(N*E)+pad),N*E);
         [f,fInd]=getfgrid(Fs,NFFT,[0 Fs/2]);
         f = permute(f,[1 3 4 5 2 6 7 8]); % frequencies[time x trial x run x taper x freq x vox x window x mode]
-        [Nf,Ef,Rf,Kf,Ff,Vf,Wf] = size(f);
+        [Nf,Ef,Rf,Kf,Ff,Vf,Wf,Mf] = size(f);
         F = Ff;
 
 
@@ -847,7 +853,7 @@
         %%% loop over windows
         for wInd = 1:W
             %%% Compute J
-            ind  = w(:,:,:,:,:,:,wInd);
+            ind  = w(:,:,:,:,:,:,wInd,:);
             tWin = reshape(funTs.t(ind,:,:,:),size(ind)); tWin = tWin([1 end],:);
             d = funTs.vec(ind,:,:,:);   % [time  x vox x taper x run               ]
             d = permute(d,[3 2 4 1]);   % [taper x vox x run   x time*trial        ]
@@ -857,27 +863,31 @@
                 d = d - mean(d,1);
             end
             tp = tp;
-            t = t; % adjust t here for sub-tr stimulus onsets
+            t = t; %!!!!!!!!!!!!!!! adjust t here for sub-tr stimulus onsets
+
+
+
+            % dbstack; error('the error is somewhere here: in tp,t or f, or within getJ2')
             J = getJ2(d,tp,t,f)/Fs; % [time x trial x run x taper x freq x vox x window]
             %                         [N      E       R     K       F      V     W]
 
             %%% Compute psd at each trial
-            res.trialGram.PSD(:,:,:,:,:,:,wInd)    = mean(  conj(J).*J  ,4);
+            res.trialGram.PSD(:,:,:,:,:,:,wInd,:)    = mean(  conj(J).*J  ,4); % Compute power then average across tapers...
             %%% Compute psd averaged across trials
-            res.trialGram.PSDeav(:,:,:,:,:,:,wInd)   = mean(  res.trialGram.PSD(:,:,:,:,:,:,wInd)  ,2);
+            res.trialGram.PSDeav(:,:,:,:,:,:,wInd,:)   = mean(  res.trialGram.PSD(:,:,:,:,:,:,wInd,:)  ,2); % then average across trials.
             %%% Compute psd phase-coherently averaged across trials
-            res.trialGram.PSDepc(:,:,:,:,:,:,wInd) = mean(  conj(mean(J,2)).*mean(J,2)  ,4);
+            res.trialGram.PSDepc(:,:,:,:,:,:,wInd,:) = mean(  conj(mean(J,2)).*mean(J,2)  ,4); % Average across trials then compute power then average across tapers.
 
 
             if K > 1
                 %%% Compute coherence at each trial
-                dim = {'N' 'E' 'R' 'K' 'F' 'V' 'W' 'Mk'};
+                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 Mk]
-                j = reshape(j,[V K E*1*R*F*1*1]); %[V K E*N*R*F*W*Mk]
+                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]
-                coh = s.^2./sum(s.^2,1); % coherence[Mk 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
                 dim = {'M' 'V' 'E' 'N' 'R' 'F' 'W' 'K'};
@@ -895,14 +905,14 @@
                 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*Mk]
-                [u,s,v] = pagesvd(j,'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 1 1 R F 1]); % [Mk V E N R F W]
+                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]
+                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'};
                 prm = [ 2    4   3   6   5   7   8   1 ];
                 dim = strjoin(dim(prm),' ');
-                coh = permute(coh,prm); %[V N E F R W K Mk]
+                coh = permute(coh,prm); %[V N E F R W K M]
                 res.trialGram.COHepc(:,:,:,:,:,:,wInd,:) = coh;
             end
 
diff --git a/simPsd3.m b/simPsd3.m
index f754da9..f8111b7 100644
--- a/simPsd3.m
+++ b/simPsd3.m
@@ -1,71 +1,109 @@
 function [volPsdSim, volTsSim] = simPsd3(dsgn,info,plotFlag)
-% volTs = vec2vol(volTs); volTs.vol = [];
-% if isfield(volTs,'resp'); volTs = rmfield(volTs,'resp'); end
-% if isfield(volTs,'dsgn'); volTs = rmfield(volTs,'dsgn'); end
 if ~exist('plotFlag','var'); plotFlag = []; end
 if isempty(plotFlag);        plotFlag = 0; end
-SNR = inf; %10;
-nVoxNull = 0; %10;
+if ~isfield(dsgn,'f0');       dsgn.f0 = []; end
+SNR = 5; %10;
+nVoxNull = 12; %10;
 
 
+if ~isempty(dsgn.onsets) && all((size(dsgn.onsets)==1)==[1 0])
+    dbstack; error('onsets must be a column vector')
+end
+
 
 % ndummy  = dsgn.dummy;
 nframes = dsgn.nframes;
-dt = dsgn.dt;
-
-onsets = (0:1/dsgn.onsetFreq:nframes*dt)';
-ondurs = repmat(1/dsgn.onFreq,size(onsets));
-nframes = ceil((onsets(end) + mode(diff(onsets))) ./ dt);
-onsets = onsets + 1/dsgn.onsetFreq/2;
-
-rndSeed = 0;
-[t, d, c] = simTs3(onsets,ondurs,nframes,dt,dt/100,SNR,plotFlag,rndSeed);
-% t(1:ndummy) = [];
-% d(1:ndummy) = [];
-% c(1:ndummy) = [];
-rndSeed = 1;
-[tShift, dShift, cShift] = simTs3(onsets+1,ondurs,nframes,dt,dt,SNR,plotFlag,rndSeed);
-% tShift(1:ndummy) = [];
-% dShift(1:ndummy) = [];
-% cShift(1:ndummy) = [];
-
-
-% T = ceil((onsetList(end)+mean(diff(onsetList)))./tr);
-% t = 0:tr:(T-1)*tr; % 5 seconds of data (time)
-% f = 0.05;
-% rng(rndSeed)
-% phi = 0;
-% c = sin(2*pi*f*t+phi).*SNR + randn(size(t));
-% phi = 2*pi/16;
-% cShift = sin(2*pi*f*t+phi).*SNR + randn(size(t));
-
-
-
-
-
-% volTsSim = vec2vol(volTs);
-
-
-volTsSim.vol = permute(cat(1,c,[]),[1 3 4 2]);
-% volTsSim.vol = permute(cat(1,c,cShift),[1 3 4 2]);
-volTsSim.vol = cat(1,volTsSim.vol,randn([nVoxNull 1 1 size(volTsSim.vol,4)]));
-
-
-
-volTsSim.dsgn.onsets = onsets;
-volTsSim.dsgn.ondurs = ondurs;
-volTsSim.t = t';
-volTsSim.nframes = length(t);
-volTsSim.nvoxels = prod(size(volTsSim.vol,1:3));
-volTsSim.tr = dt*1000;
-info.win(1) = round(dsgn.winLength/info.tr);
-
-% volTsSim.dsgn.onsets([1 end]) = [];
-% volTsSim.dsgn.ondurs([1 end]) = [];
-volTsSim.dsgn.onsets(end) = [];
-volTsSim.dsgn.ondurs(end) = [];
-
-volPsdSim = volPsdFullMt2([],info,volTsSim);
+dt = dsgn.dt; % [sec]
+
+if isempty(dsgn.f0)
+    onsets = (0:1/dsgn.onsetFreq:nframes*dt)';
+    ondurs = repmat(1/dsgn.onFreq,size(onsets));
+    nframes = ceil((onsets(end) + mode(diff(onsets))) ./ dt);
+    onsets = onsets + 1/dsgn.onsetFreq/2;
+
+    rndSeed = 0;
+    [t, d, c] = simTs3(onsets,ondurs,nframes,dt,dt/100,SNR,plotFlag,rndSeed);
+    % t(1:ndummy) = [];
+    % d(1:ndummy) = [];
+    % c(1:ndummy) = [];
+    rndSeed = 1;
+    [tShift, dShift, cShift] = simTs3(onsets+1,ondurs,nframes,dt,dt,SNR,plotFlag,rndSeed);
+    % tShift(1:ndummy) = [];
+    % dShift(1:ndummy) = [];
+    % cShift(1:ndummy) = [];
+
+
+    % T = ceil((onsetList(end)+mean(diff(onsetList)))./tr);
+    % t = 0:tr:(T-1)*tr; % 5 seconds of data (time)
+    % f = 0.05;
+    % rng(rndSeed)
+    % phi = 0;
+    % c = sin(2*pi*f*t+phi).*SNR + randn(size(t));
+    % phi = 2*pi/16;
+    % cShift = sin(2*pi*f*t+phi).*SNR + randn(size(t));
+
+
+
+
+
+    % volTsSim = vec2vol(volTs);
+
+
+    volTsSim.vol = permute(cat(1,c,[]),[1 3 4 2]);
+    % volTsSim.vol = permute(cat(1,c,cShift),[1 3 4 2]);
+    volTsSim.vol = cat(1,volTsSim.vol,randn([nVoxNull 1 1 size(volTsSim.vol,4)]));
+
+
+
+    volTsSim.dsgn.onsets = onsets;
+    volTsSim.dsgn.ondurs = ondurs;
+    volTsSim.t = t';
+    volTsSim.nframes = length(t);
+    volTsSim.nvoxels = prod(size(volTsSim.vol,1:3));
+    volTsSim.tr = dt*1000;
+    info.win(1) = round(dsgn.winLength/info.tr);
+
+    % volTsSim.dsgn.onsets([1 end]) = [];
+    % volTsSim.dsgn.ondurs([1 end]) = [];
+    volTsSim.dsgn.onsets(end) = [];
+    volTsSim.dsgn.ondurs(end) = [];
+
+    volPsdSim = volPsdFullMt2([],info,volTsSim);
+
+
+else
+    tS = 0; tE = (dsgn.nframes+dsgn.dummy-1)*dsgn.dt;
+    t = tS:dt:tE;
+    theta  = 0;
+    amp    = 1;
+    c      = real(   amp .* exp(1i.* (2*pi*dsgn.f0.*t - theta)  )   )  +  1;
+    theta  = 1/8 .* (2*pi);
+    amp    = 1;
+    cShift = real(   amp .* exp(1i.* (2*pi*dsgn.f0.*t - theta)  )   )  +  1;
+    % figure('WindowStyle','docked');
+    % plot(t,real(c),'.-'); hold on
+    % plot(t,real(cShift),'.-r');
+
+    %% Add noise and null voxels
+    volTsSim.vol = permute(cat(1,c,cShift),[1 3 4 2]);
+    volTsSim.vol = cat(1,volTsSim.vol,zeros([nVoxNull 1 1 size(volTsSim.vol,4)])   +  1);
+    if SNR~=inf
+        rng(0)
+        volTsSim.vol = volTsSim.vol.*SNR + randn(size(volTsSim.vol)).*std(c);
+    end
+    
+    volTsSim.t = t';
+    volTsSim.vol(:,:,:,1:dsgn.dummy) = [];
+    volTsSim.t(1:dsgn.dummy) = [];
+
+    volTsSim.dsgn = dsgn;
+    volTsSim.nframes = dsgn.nframes;
+    volTsSim.dummy = dsgn.dummy;
+    volTsSim.nvoxels = prod(size(volTsSim.vol,1:3));
+    volTsSim.tr = dt*1000;
+
+    volPsdSim = [];
+end
 
 
 
diff --git a/volPsdFullMt2.m b/volPsdFullMt2.m
index 583cf82..91725ed 100644
--- a/volPsdFullMt2.m
+++ b/volPsdFullMt2.m
@@ -13,7 +13,9 @@
 if ~isfield(info,'dtrndOrder'); info.dtrndOrder = []        ; end
 if ~isfield(info,'onsetList');   info.onsetList = []        ; end
 if ~isfield(info,'ondurList');   info.ondurList = []        ; end
+if ~isfield(info,'perm');             info.perm = []         ; end
 
+if isempty(info.perm); info.perm = 0; end
 
 
 %% User variables