viewPSD2.m

function [hF,tl,axSpec,maskC] = viewPSD2(funPsd,f0,mask,fpass,id,threshFlag)
%% Prepare some stuff
if ~exist("threshFlag",'var') || isempty(threshFlag)
    threshFlag = false;
end
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
    roiFlag = 0;
    maskC = [];
end
if ~exist('f0','var')
    f0=0.1;
else
    f0 = sort(f0);
end
if ~exist('id','var')
    id = {};
end
funPsd = vol2vec(funPsd);


%% 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])
if isfield(funPsd,'imMean') && ~all(size(funPsd.imMean)==1)
    im = squeeze(mean(funPsd.imMean,6));
    imagesc(im)
    ax = gca;
    ax.Colormap = gray;
    ax.YTick = []; ax.XTick = [];
    ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
    ax.CLim = prctile(im(:),[0 99]);
    title('timeseries mean')
else
    ax = gca;
    ax.Visible = 'off';
end

%% Mask
nexttile(2,[2 1])
if isempty(maskC)
    ax = gca;
    ax.Visible = 'off';
else
    im = squeeze(funPsd.imMean);
    imagesc(im)
    ax = gca;
    ax.Colormap = gray;
    ax.YTick = []; ax.XTick = [];
    ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
    ax.CLim = prctile(im(:),[0 99]);
    hold on
    hMask1 = plot(maskC);
    hMask1.FaceColor = 'r';
    hMask1.FaceAlpha = 0.1;
    hMask1.EdgeColor = 'none';
    title('mask')
end

%% Normalization
% drawnow
nexttile(3,[2 1])
if isfield(funPsd.psd,'norm') && ~isfield(funPsd.psd.norm,'fact2')
    normFlag = 1;
    
    normFact = nan([size(funPsd.vec,[3 4]) size(funPsd.vol2vec)]);
    normFact(:,:,funPsd.vol2vec) = permute(funPsd.psd.norm.fact,[3 4 2 1]);
    normFact = permute(normFact,[3 4 5 1 2]);
    if ~isreal(normFact)
        normFact = conj(normFact).*normFact;
    end
    imagesc(squeeze(mean(normFact,5)))
    hold on
    if ~isempty(maskC)
        hMask2 = plot(maskC);
        hMask2.FaceColor = 'none';
        hMask2.EdgeColor = 'k';
    end
    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
f = funPsd.psd.f;
allF0 = f0;
for fInd = 1:length(allF0)
    f0 = allF0(fInd);
    if fInd==1
        nexttile(4,[2 1])
    elseif fInd>5
        warning('too many f0 to plot, skipping')
        break
    else
        nexttile(8+fInd-1,[2 1])
    end
    if size(funPsd.vec,2)>1
        [~,f0Ind] = min(abs(f-f0'),[],2);
        f0 = f(f0Ind);
        tmpPsd = vec2vol(funPsd);
        tmpIm = squeeze(tmpPsd.vol(:,:,:,f0Ind));
        if ~isreal(tmpIm)
            tmpIm = conj(tmpIm).*tmpIm;
        end
        hIm = imagesc(tmpIm);
        hold on
        if ~isempty(maskC)
            hMask3 = plot(maskC);
            hMask3.FaceColor = 'none';
            hMask3.EdgeColor = 'k';
        end
        ax = gca;
        ax.YTick = []; ax.XTick = [];
        ax.ColorScale = 'log'; ax.Colormap = jet;
        ax.PlotBoxAspectRatio = [1 1 1]; ax.DataAspectRatio = [1 1 1];
        if normFlag
            ylabel(colorbar,'normalized psd')
        else
            ylabel(colorbar,'raw psd')
        end

        if threshFlag && isfield(funPsd.psd,'aboveNoiseInd')
            tmp = zeros(size(mask));
            tmp(:) = funPsd.psd.aboveNoiseInd(:,f0Ind);
            hIm.AlphaData = tmp;
            ax.CLim = prctile(tmpIm(logical(tmp)&logical(mask)),[0 95]);
            ax.Color = 'k';
            hMask3.EdgeColor = 'w';
            title(['f0=' num2str(f0,'%0.3f') 'Hz (thresholded)'])
        else
            title(['f0=' num2str(f0,'%0.3f') 'Hz'])
        end
        if normFlag
            ax.CLim(1) = 1;
            if exist('mask','var') && ~isempty(mask)
                ax.CLim(2) = max(tmpIm(logical(mask)));
            end
        end
    else
        ax = gca;
        ax.Visible = 'off';
    end
%     %%% scale every image to the first
%     if fInd==1
%         cLim  = ax.CLim;
%     else
%         ax.CLim = cLim;
%     end
end
f0 = allF0;



%% Spectrum
% nexttile(9,[1 4])
axSpec = 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
%     yyaxis right
%     h{end+1} = plot(f,diff(psdErr,[],2)./psd);
    if ~isempty(psdErr)
        h{end+1} = plot(f,psdErr,'r'); legLabel{end+1} = '95%CI';
        h{end} = h{end}(1);
    end
else
    tmpPsd = vec2vol(funPsd);
    if ~exist('mask','var')
        mask = true(size(tmpPsd.vol,1:3));
    end
    tmpPsd = vol2vec(tmpPsd,mask,1);
    if ~isreal(tmpPsd.vec)
        tmpPsd.vec = conj(tmpPsd.vec).*tmpPsd.vec;
    end
    psd = mean(tmpPsd.vec,2,'omitnan');
    psd = mean(psd,4);
    W = funPsd.psd.w;
    h{end+1} = plot(f,squeeze(psd),'k'); legLabel{end+1} = 'mean spectrum';
    hold on
end
ax = gca;
ax.YScale = 'log';
ax.YAxisLocation = 'right';
axis tight
if exist('fpass','var') && ~isempty(fpass)
    ax.XLim = fpass;
end
yLim = [min(psd) max(psd)];
if max(f)>3.5
    yLim(1) = mean(psd(f>3.5));
end
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
allF0 = f0;
for fInd = 1:length(allF0)
    f0 = allF0(fInd);
    [~,f0Ind] = min(abs(f-f0'),[],2);
    f0 = f(f0Ind);
    if fInd==1
        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';
    else
        plot(f0.*[1 1],yLim,':r');
        plot(f0+W.*[-1 1],yLim(1).*[1 1],'g','LineWidth',5);
    end
end
f0 = allF0;
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
if iscell(id)
    titleStr1 = {strjoin(id,'; ')};
else
    titleStr1 = {id};
end
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