Skip to content

Commit

Permalink
commit before adapting for vsmDriven analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
@Proulx-S
Proulx-S committed Jul 30, 2024
1 parent 5f462aa commit b2bb93c
Show file tree
Hide file tree
Showing 8 changed files with 293 additions and 106 deletions.
43 changes: 26 additions & 17 deletions addFreq.m
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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

140 changes: 85 additions & 55 deletions getResp.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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];
Expand All @@ -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
Expand All @@ -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');
Expand Down Expand Up @@ -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 ' \'];
Expand All @@ -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;
Expand Down
34 changes: 25 additions & 9 deletions plotSpec.m
View file
Original file line number Diff line number Diff line change
@@ -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'
Expand Down Expand Up @@ -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));
Expand All @@ -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)];
Expand All @@ -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)
Expand Down
Loading

0 comments on commit b2bb93c

Please sign in to comment.