runPSD2.m

function funPsd = runPSD2(funTs,W,K,cFlag,mask)
%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
if ~exist('W','var'); W = []; end
if ~exist('K','var'); K = []; end
if ~exist('cFlag','var') || isempty(cFlag); cFlag = 0; end
if ~exist('mask','var'); mask = []; end

%% Detrend time series (detrend up to order-2 polynomial, since this is the highest order not fitting a sinwave)
tic
funTs2 = cell(size(funTs));
for i = 1:length(funTs)
    funTs2{i} = dtrnd4psd(funTs(i));
end
funTs = [funTs2{:}];
toc

%% Put runs in same structure
if length(funTs)>1
    funTs2 = funTs(1);
    funTs2.vol = cat(6,funTs.vol);
    funTs2.volMean = cat(6,funTs.volMean);
    funTs2.nruns = length(funTs);
    funTs = funTs2;
end

%% Compute PSD
funPsd = doIt(funTs,W,K,cFlag,mask);

%% Add image mean
funPsd.imMean = funTs.volMean;

%% Split back in different structures
nruns = size(funPsd.vec,4);
if nruns>1
    funPsd2 = funPsd;
    funPsd2.vec = [];
    funPsd2.nruns = 1;
    funPsd2 = repmat(funPsd2,[1 nruns]);
    for i = 1:nruns
        funPsd2(i).vec = funPsd.vec(:,:,:,i);
        funPsd2(i).vecMean = funPsd.vecMean(:,:,:,i);
        funPsd2(i).tr = funPsd.tr(:,:,:,i);
    end
end
funPsd = funPsd2; clear funPsd2

function funPsd = doIt(funTs,W,K,cFlag,mask)
if exist('mask','var') && ~isempty(mask)
    funTs = vol2vec(funTs,mask,1);
else
    funTs = vol2vec(funTs);
end
tr = funTs.tr/1000;

if any(all(funTs.vec==0,1))
    warning('Some voxels are all 0s. Adjust your mask to avoid later problems')
end

%% Set parameters
Wflag = exist('W','var') && ~isempty(W);
Kflag = exist('K','var') && ~isempty(K);
if Wflag && Kflag
    error('Cannot specify both W and K');
elseif Wflag
    T = tr.*funTs.nframes;
    TW = T*W;
    K = round(TW*2-1);
    TW = (K+1)/2;
    param.tapers = [TW K];
    Wreal = TW/T;
    display(['w  (halfwidth) requested  : ' num2str(W,'%0.5f ')])
    display(['w  (halfwidth) used       : ' num2str(Wreal,'%0.5f ')])
    display(['tw (time-halfwidth) used  : ' num2str(TW)])
    display(['k  (number of tapers) used: ' num2str(K)])
elseif Kflag
    TW = (K+1)/2;
    T = tr.*funTs.nframes;
    W = TW/T; Wreal = W;
    param.tapers = [TW K];
    display(['k  (number of tapers) requested  : ' num2str(K)])
    display(['w  (halfwidth) used       : ' num2str(W,'%0.5f ')])
    display(['tw (time-halfwidth) used  : ' num2str(TW)])
end


%% Perform the multitaper PSD estimation
funPsd = funTs; funPsd.vec = [];
param.Fs = 1/tr;
param.complex = cFlag;
% [funPsd.vec,f] = mtspectrumc(funTs.vec, param);
i=1;
[funPsd.vec,vecC,f] = mtspectrumc2(funTs.vec(:,:,:,i), param);
if size(funTs.vec,4)>1
    vec = repmat(funPsd.vec,[1 1 1 size(funTs.vec,4)]);
    vecC = repmat(vecC,[1 1 1 size(funTs.vec,4)]);
    f = repmat(f,[1 1 1 size(funTs.vec,4)]);
    parfor i = 2:size(funTs.vec,4)
        [vec(:,:,:,i),vecC(:,:,:,i),f(:,:,:,i)] = mtspectrumc2(funTs.vec(:,:,:,i), param);
    end
    funPsd.vec = vec;
end
funPsd.nframes = size(funPsd.vec,1);
funPsd.tr = mode(diff(f))*1000;

%% Output some stuff
if isfield(funPsd,'psd')
    funPsd = rmfield(funPsd,'psd');
end
funPsd.psd.dim = strjoin({'space' 'freq' 'taper' 'run'},' X ');
tmp = diff(f,[],4); if any(tmp(:)); dbstack; error('X'); end
funPsd.psd.f = f(:,:,:,1);
if param.complex
    funPsd.psd.spec = permute(vecC,[2 1 3 4]);
else
    funPsd.psd.spec = [];
end
funPsd.psd.w = Wreal;
funPsd.psd.tw = TW;
funPsd.psd.mask = funTs.vol2vec;
funPsd.psd.maskLabel = funPsd.vol2vecFlag;
funPsd.psd.param = param;