runPSD.m

function funPsd = runPSD(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('cFlag','var') || isempty(cFlag)
    cFlag = 0;
end

tMean = mean(funTs.vol,4);
if exist('mask','var') && ~isempty(mask)
    funTs = vol2vec(funTs,mask);
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

%% Detrend time series (detrend up to order-2 polynomial, since this is the highest order not fitting a sinwave)
funTs = dtrnd4psd(funTs);

%% Perform the multitaper PSD estimation
funPsd = funTs; funPsd.vec = [];
param.Fs = 1/tr;
param.complex = cFlag;
% [funPsd.vec,f] = mtspectrumc(funTs.vec, param);
[funPsd.vec,vecC,f] = mtspectrumc2(funTs.vec, param);
funPsd.nframes = size(funPsd.vec,1);
funPsd.tr = mode(diff(f))*1000;

% plot(f,mean(funPsd.vec,2))
% funPsd = vec2vol(funPsd);
% [~,b] = min(abs(f-0.931351));
% imagesc(funPsd.vol(:,:,1,b));
% ax = gca; ax.ColorScale = 'log';

%% Output some stuff
if isfield(funPsd,'psd')
    funPsd = rmfield(funPsd,'psd');
end
funPsd.psd.dim = strjoin({'space' 'freq' 'taper'},' X ');
funPsd.psd.f = f;
if param.complex
    funPsd.psd.spec = permute(vecC,[2 1 3]);
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;

funPsd.tMean = tMean;

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