rtc_calcbetas.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function [betas, irtc, ptc, se] = rtc_CalcBetas(hfile, tcd, dim)
% RTC::CalcBetas  - perform GLM calculcation
%
% FORMAT:       [betas, irtc, ptc, se] = rtc.CalcBetas(tc [, tdim])
%
% Input fields:
%
%       tc          time course data
%       tdim        required for numeric tc data, temporal dimension
%
% Output fields:
%
%       betas       betas, in time course data dimension
%       irtc        inverse design matrix
%       ptc         predicted tc
%       se          standard error
%
% Note: if no constant term is in the matrix, it will be added as the
%       last column

% Version:  v0.7b
% Build:    7083011
% Date:     Aug-30 2007, 11:30 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'rtc')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% further arguments
if nargin < 3 || ...
   ~isa(dim, 'double') || ...
    numel(dim) ~= 1 || ...
    isinf(dim) || ...
    isnan(dim) || ...
    fix(dim) ~= dim || ...
    dim < 1 || ...
    dim > 5
    rdim = 1;
else
    rdim = dim;
end

% get design matrix right
bc = bvqxfile_getcont(hfile.L);
if isfield(bc, 'Type') && ...
    strcmpi(bc.Type, 'designmatrix')
    rtc = bc.RTCMatrix;
else
    error( ...
        'BVQXfile:TypeNotSupported', ...
        'This RTC type is currently not supported for beta calculation.' ...
    );
end

% add all-1 column for confound
mrtc = mean(rtc);
vrtc = var(rtc);
if ~any(mrtc ~= 0 & vrtc == 0)
    rtc(:, end + 1) = 1;
end

% build rtc' and pinv(rtc' * rtc)
trtc = rtc';
irtc = pinv(trtc * rtc);
prtc = irtc * trtc;

% number of rows and betas
numrows = size(rtc, 1);
numbets = size(rtc, 2);

% check second argument -> BVQXfile object
if numel(tcd) == 1 && ...
    isBVQXfile(tcd, true)

    % get struct version
    sstc = bvqxfile_getscont(tcd.L);
    stcd = sstc.C;

    % get filetype
    try
        ftype = lower(sstc.S.Extensions{1});
    catch
        error( ...
            'BVQXfile:InternalError', ...
            'Invalid BVQXfile object given as time course data.' ...
        );
    end

    % switch over filetype
    switch lower(ftype)

        % FMR
        case {'fmr'}

            % STCs loaded?
            if isempty(stcd.Slice) || ...
               ~isstruct(stcd.Slice) || ...
               ~isfield(stcd.Slice, 'STCData')
                try
                    fmr_LoadSTCs(tcd);
                    stcd = bvqxfile_getcont(tcd.L);
                catch
                    error( ...
                        'BVQXfile:InternalError', ...
                        'Error loading STCs for FMR %s.', ...
                        sstc.F ...
                    );
                end
            end

            % check MTC
            nvol = numel(stcd.Slice(1).STCData(:, 1));
            if nvol ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of time points in FMR/STC.' ...
                );
            end

            % get STC data
            nslc = stcd.NrOfSlices;
            xres = stcd.ResolutionX;
            yres = stcd.ResolutionY;
            tcd = zeros(nvol, xres, yres, nslc);
            switch (stc.DataStorageFormat)
                case {1}
                    for sc = 1:nslc
                        tcd(:, :, :,sc) = permute( ...
                            stcd.Slice(sc).STCData(:, :, :), [3, 1, 2]);
                    end
                case {2}
                    for sc = 1:nslc
                        tcd(:, :, :,sc) = permute( ...
                            stcd.Slice.STCData(:, :, :, sc), [3, 1, 2]);
                    end
                otherwise
                    error( ...
                        'BVQXfile:InvalidObject', ...
                        'Unsupported FMR DataStorageFormat.' ...
                    );
            end

        % MTC
        case {'mtc'}

            % check MTC
            if size(stcd.MTCData, 1) ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of time points in MTC.' ...
                );
            end

            % get MTC data
            tcd = double(stcd.MTCData(:, :));

        % VTC
        case {'vtc'}

            % check VTC
            if size(stcd.VTCData, 1) ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of volumes in VTC.' ...
                );
            end

            % get VTC data
            tcd = double(stcd.VTCData(:, :, :, :));

        % otherwise
        otherwise

            error( ...
                'BVQXfile:TypeNotSupported', ...
                'BVQXfile object of type %s not supported.', ...
                ftype ...
            );
    end

% for other numerics
elseif isnumeric(tcd) && ...
   ~isa(tcd, 'double')

    % create double
    tcd = double(tcd);

end

% check dim
if rdim < 1 || ...
    rdim > length(size(tcd)) || ...
    size(tcd, rdim) ~= numrows
    error( ...
        'BVQXfile:InternalError', ...
        'Invalid data matrix vs. selected dimension.' ...
    );
end
if rdim > 1
    neword = [rdim, 1:length(size(tcd))];
    newodr = find(neword == rdim);
    neword(newodr(2)) = [];
    tcd = permute(tcd, neword);
end

% reshaping data to comply
tcds = size(tcd);
numvox = prod(tcds(2:end));
tcd = reshape(tcd, [tcds(1), numvox]);
tcds(1) = [];
tcdrs = tcds;
if length(tcdrs) < 2
    tcdrs(2) = 1;
end

% perform calculus
betas = zeros([numbets, numvox]);
for vc = 1:numvox
    betas(:, vc) = prtc * tcd(:, vc);
end

% more output?
if nargout > 2
    ptc = rtc * betas;
    ptc = shiftdim(reshape(shiftdim(ptc, 1), [tcds, numrows]), ...
        length(tcds));
    if nargout > 3
        se = squeeze(std(reshape(tcd, [numrows, tcdrs]) - ptc, 0, rdim)) ...
             * sqrt((numrows - 1) / (numrows - numbets));
    end
end

% reshape betas
betas = reshape(shiftdim(betas, 1), [tcds, numbets]);