Skip to content

Commit

Permalink
Carry timeseries mean from runPSD.m and revamp plots (simplify the co…
Browse files Browse the repository at this point in the history
…de) in viewPSD.m
  • Loading branch information
Proulx-S committed May 24, 2023
1 parent 832472f commit 86aa618
Show file tree
Hide file tree
Showing 3 changed files with 243 additions and 94 deletions.
16 changes: 9 additions & 7 deletions runPSD.m
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
function [funPsd,psdStruct] = runPSD(funTs,W)
function funPsd = runPSD(funTs,W)
%Wrapper of the Chronux's mtspectrumc function for multitaper estimation of
%pds spectra, compatible with MRI data imported by MRIread.m.
% Parameterization is simplified to the halfbandwidth parameter W only.
% Includes minimal detrending to avoid distortions from the low frquency
% edge.
% See funPsd.psd for other useful parameters.
% funPsd.tr reflects the frequency resolution in Hz*1000

tMean = mean(funTs.vol,4);
funTs = vol2vec(funTs);
tr = funTs.tr/1000;

Expand Down Expand Up @@ -52,9 +52,11 @@
funPsd.psd.maskLabel = funPsd.vol2vecFlag;
funPsd.psd.param = param;

if nargout>1
psdStruct = funPsd.psd;
psdStruct.psd = funPsd.vec';
psdStruct = setNiceFieldOrder(psdStruct,{'dim' 'psd' 'f'});
end
funPsd.tMean = tMean;

% if nargout>1
% psdStruct = funPsd.psd;
% psdStruct.psd = funPsd.vec';
% psdStruct = setNiceFieldOrder(psdStruct,{'dim' 'psd' 'f'});
% end

174 changes: 87 additions & 87 deletions viewPSD.m
Original file line number Diff line number Diff line change
Expand Up @@ -46,89 +46,89 @@
f = funPsd.psd.f;
w = funPsd.psd.w;

%% Perform spatial normalization
switch normMethod
case 'psdNoise'
if isempty(funPsd.vec)
funPsd = vol2vec(funPsd);
end
noiseInd = funPsd.psd.f>4 & funPsd.psd.f<funPsd.psd.f(round(0.95*end));
normFact = exp( mean(log(funPsd.vec(noiseInd,:)),1) );
normFact0 = exp(mean(log(normFact),2));
case 'psdAv'
if isempty(funPsd.vec)
funPsd = vol2vec(funPsd);
end
normFact = exp( mean(log(funPsd.vec),1) );
normFact0 = exp(mean(log(normFact),2));
case 'smoothPsdAv'
error('this is probably broken')
%%% Compute average power and write to disk
funNorm = funPsd;
if isempty(funNorm.vol)
funNorm.vec = mean(log(funNorm.vec),1);
funNorm.nframes = 1;
normFact0 = exp(funNorm.vec);
funNorm = vec2vol(funNorm);
else
funNorm.vol = mean(log(funNorm.vol),4);
funNorm.nframes = 1;
tmp = vol2vec(funNorm);
normFact0 = exp(tmp.vec); clear tmp
end
funNormTmpName = [tempname '.nii.gz'];
MRIwrite(funNorm,funNormTmpName);


%%% Define mask for smoothing and write to disk
funMask = funNorm;
funMask.vol = ~isnan(funMask.vol);
funMaskTmpName = [tempname '.nii.gz'];
MRIwrite(funMask,funMaskTmpName);

%%% Apply smoothing (mri_fwhm) to normalization map
cmd = {'source /usr/local/freesurfer/fs-dev-env-autoselect'};
cmd{end+1} = ['mri_fwhm'...
' --i ' funNormTmpName...
' --o ' funNormTmpName...
' --smooth-only'...
' --mask ' funMaskTmpName...
' --fwhm ' num2str(fwhm)];
[status,cmdout] = system(strjoin(cmd,'; '));

% funNormSm = MRIread(funNormTmpName);
% figure('WindowStyle','docked');
% imagesc(exp(funNorm.vol(:,:,46)));
% ax = gca; ax.ColorScale = 'log';
% cLim = clim;
% figure('WindowStyle','docked');
% imagesc(exp(funNormSm.vol(:,:,46)));
% ax = gca; ax.ColorScale = 'log';
% clim(cLim);

funNorm = MRIread(funNormTmpName);
funNorm = vol2vec(funNorm,funPsd.vol2vec);
normFact = exp(funNorm.vec); clear funNorm
case 'none'
sz = size(funPsd.vec);
normFact = ones([1 sz(2)]);
normFact0 = ones([1 sz(2)]);
otherwise
error('')
end

%%% Apply normalization
if isempty(funPsd.vec)
funPsd = vol2vec(funPsd);
end
funPsdNorm = funPsd;
switch globalNormMethod
case 'toImageAverage'
funPsdNorm.vec = funPsd.vec./normFact .* normFact0;
case 'to1'
funPsdNorm.vec = funPsd.vec./normFact;
end
funPsdNorm.normFact = normFact;
% %% Perform spatial normalization
% switch normMethod
% case 'psdNoise'
% if isempty(funPsd.vec)
% funPsd = vol2vec(funPsd);
% end
% noiseInd = funPsd.psd.f>4 & funPsd.psd.f<funPsd.psd.f(round(0.95*end));
% normFact = exp( mean(log(funPsd.vec(noiseInd,:)),1) );
% normFact0 = exp(mean(log(normFact),2));
% case 'psdAv'
% if isempty(funPsd.vec)
% funPsd = vol2vec(funPsd);
% end
% normFact = exp( mean(log(funPsd.vec),1) );
% normFact0 = exp(mean(log(normFact),2));
% case 'smoothPsdAv'
% error('this is probably broken')
% %%% Compute average power and write to disk
% funNorm = funPsd;
% if isempty(funNorm.vol)
% funNorm.vec = mean(log(funNorm.vec),1);
% funNorm.nframes = 1;
% normFact0 = exp(funNorm.vec);
% funNorm = vec2vol(funNorm);
% else
% funNorm.vol = mean(log(funNorm.vol),4);
% funNorm.nframes = 1;
% tmp = vol2vec(funNorm);
% normFact0 = exp(tmp.vec); clear tmp
% end
% funNormTmpName = [tempname '.nii.gz'];
% MRIwrite(funNorm,funNormTmpName);
%
%
% %%% Define mask for smoothing and write to disk
% funMask = funNorm;
% funMask.vol = ~isnan(funMask.vol);
% funMaskTmpName = [tempname '.nii.gz'];
% MRIwrite(funMask,funMaskTmpName);
%
% %%% Apply smoothing (mri_fwhm) to normalization map
% cmd = {'source /usr/local/freesurfer/fs-dev-env-autoselect'};
% cmd{end+1} = ['mri_fwhm'...
% ' --i ' funNormTmpName...
% ' --o ' funNormTmpName...
% ' --smooth-only'...
% ' --mask ' funMaskTmpName...
% ' --fwhm ' num2str(fwhm)];
% [status,cmdout] = system(strjoin(cmd,'; '));
%
% % funNormSm = MRIread(funNormTmpName);
% % figure('WindowStyle','docked');
% % imagesc(exp(funNorm.vol(:,:,46)));
% % ax = gca; ax.ColorScale = 'log';
% % cLim = clim;
% % figure('WindowStyle','docked');
% % imagesc(exp(funNormSm.vol(:,:,46)));
% % ax = gca; ax.ColorScale = 'log';
% % clim(cLim);
%
% funNorm = MRIread(funNormTmpName);
% funNorm = vol2vec(funNorm,funPsd.vol2vec);
% normFact = exp(funNorm.vec); clear funNorm
% case 'none'
% sz = size(funPsd.vec);
% normFact = ones([1 sz(2)]);
% normFact0 = ones([1 sz(2)]);
% otherwise
% error('')
% end
%
% %%% Apply normalization
% if isempty(funPsd.vec)
% funPsd = vol2vec(funPsd);
% end
% funPsdNorm = funPsd;
% switch globalNormMethod
% case 'toImageAverage'
% funPsdNorm.vec = funPsd.vec./normFact .* normFact0;
% case 'to1'
% funPsdNorm.vec = funPsd.vec./normFact;
% end
% funPsdNorm.normFact = normFact;


%% Compute summary statistics for later
Expand All @@ -154,10 +154,10 @@
end
titleStr{end+1} = ['tw=' num2str(funPsd.psd.tw) ', ' 'w=' num2str(w) 'Hz' ', ' 'k=' num2str(funPsd.psd.param.tapers(2))];
% titleStr{end+1} = ['w=' num2str(w) 'Hz'];
switch normMethod
case 'smoothPsdAv'
titleStr{end+1} = ['fwhm=' num2str(fwhm) 'mm'];
end
% switch normMethod
% case 'smoothPsdAv'
% titleStr{end+1} = ['fwhm=' num2str(fwhm) 'mm'];
% end
[~,b,c] = fileparts(funPsd.fspec);
titleStr{end+1} = [b c];
% titleStr = [strjoin(titleStr(1:end-1),'; ') titleStr(end)];
Expand Down
147 changes: 147 additions & 0 deletions viewPSD2.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
function hF = viewPSD2(funPsd,f0,mask,fpass,id)


%% Prepare some stuff
if exist('mask','var') && ~isempty(mask)
roiFlag = 0;
maskC = getMaskOutline(mask,5);
elseif isfield(funPsd,'roi')
roiFlag = 1;
mask = funPsd.roi.mask;
maskC = getMaskOutline(mask,5);
else
error('code that')
end

%% Prepare figure
hF = figure('WindowStyle','docked');
% tl = tiledlayout(3,4);
tl = tiledlayout(5,4);
tl.TileSpacing = "tight"; tl.Padding = "tight";
tl.TileIndexing = 'rowmajor';

%% Brain
nexttile(1,[2 1])
imagesc(funPsd.tMean)
ax = gca;
ax.Colormap = gray;
ax.YTick = []; ax.XTick = [];
ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
title('timeseries mean')

%% Mask
nexttile(2,[2 1])
imagesc(funPsd.tMean)
ax = gca;
ax.Colormap = gray;
ax.YTick = []; ax.XTick = [];
ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
hold on
hMask1 = plot(maskC);
hMask1.FaceColor = 'r';
hMask1.FaceAlpha = 0.1;
hMask1.EdgeColor = 'none';
title('mask')

%% Normalization
nexttile(3,[2 1])
if isfield(funPsd.psd,'norm')
normFlag = 1;
normFact = nan(size(funPsd.vol2vec));
normFact(funPsd.vol2vec) = funPsd.psd.norm.fact;
imagesc(normFact)
hold on
hMask2 = plot(maskC);
hMask2.FaceColor = 'none';
hMask2.EdgeColor = 'k';
ax = gca;
ax.YTick = []; ax.XTick = [];
ax.ColorScale = 'log'; ax.Colormap = jet;
ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
ylabel(colorbar,'noise floor psd')
title('normalization factor')
else
normFlag = 0;
ax = gca;
ax.Visible = 'off';
end

%% f0
nexttile(4,[2 1])
f = funPsd.psd.f;
[~,f0Ind] = min(abs(f-f0'),[],2);
tmpPsd = vec2vol(funPsd);
imagesc(tmpPsd.vol(:,:,:,f0Ind));
hold on
hMask3 = plot(maskC);
hMask3.FaceColor = 'none';
hMask3.EdgeColor = 'k';
ax = gca;
ax.YTick = []; ax.XTick = [];
ax.ColorScale = 'log'; ax.Colormap = jet;
ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
ylabel(colorbar,'raw psd')
title(['f0=' num2str(f0,'%0.3f') 'Hz'])

%% Spectrum
% nexttile(9,[1 4])
nexttile(17,[1 4])
legLabel = {};
h = {};
if roiFlag
psd = funPsd.roi.psd;
psdErr = funPsd.roi.psdErr;
W = funPsd.roi.w;
h{end+1} = plot(f,psd,'k'); legLabel{end+1} = 'mean spectrum';
hold on
h{end+1} = plot(f,psdErr,'r'); legLabel{end+1} = '95%CI';
h{end} = h{end}(1);
else
tmpPsd = vec2vol(funPsd);
tmpPsd = vol2vec(tmpPsd,mask,1);
psd = mean(tmpPsd.vec,2);
W = funPsd.psd.w;
h{end+1} = plot(f,psd,'k'); legLabel{end+1} = 'mean spectrum';
hold on
end
ax = gca;
ax.YScale = 'log';
ax.YAxisLocation = 'right';
axis tight
ax.XLim = fpass;
yLim = [min(psd) max(psd)];
yLim(1) = mean(psd(f>3.5));
yLim(1) = exp(log(yLim(1)) - range(log(yLim))*0.05);
ylim(yLim);
grid on
ax.XMinorGrid = 'on';
xlabel('Hz')
if normFlag
ylabel('normalized psd')
else
ylabel('raw psd')
end
h{end+1} = plot(f0.*[1 1],yLim,':r'); legLabel{end+1} = 'f0';
h{end+1} = plot(f0+W.*[-1 1],yLim(1).*[1 1],'g','LineWidth',5); legLabel{end+1} = 'BW';
legend([h{:}],legLabel,'box','off')
if normFlag
title('spectrum average within mask, normalized voxel-wise to noise=1')
else
title('spectrum average within mask')
end

%% Title
titleStr1 = {strjoin(id,'; ')};
if normFlag
titleStr2 = {'normalized voxel-wise'};
else
titleStr2 = {'raw'};
end
titleStr2{end+1} = ['W=' num2str(funPsd.psd.w,'%0.4f')];
titleStr2{end+1} = ['K=' num2str(funPsd.psd.param.tapers(2))];
titleStr2 = {strjoin(titleStr2,'; ')};

titleStr = [titleStr1 titleStr2];
tlTitle = title(tl,titleStr,'interpreter','none');
drawnow

0 comments on commit 86aa618

Please sign in to comment.