prt_creatertc.m

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

function hfile2 = prt_CreateRTC(hfile, params)
% PRT::CreateRTC  - create an RTC from a PRT
%
% FORMAT:       rtc = prt.CreateRTC(params);
%
% Input fields:
%
%       params      1x1 struct with optional fields
%        .hshape    HRF shape to use ({'twogamma'}, 'boynton')
%        .hpttp     HRF time to peak for positive response (5)
%        .hnttp     HRF time to peak for negative response (15)
%        .hpnr      HRF positive/negative ratio (6)
%        .hons      HRF onset (0)
%        .hpdsp     HRF positive response dispersion (1)
%        .hndsp     HRF negative response dispersion (1)
%        .ndcreg    deconvolution number of lags (regressors, 8)
%        .nderiv    add derivatives to HRF (1xN double, [])
%        .nvol      number of volumes (data points, 200)
%        .ortho     if given and true, orthogonalize derivatives
%        .params    1xN struct for parametrical designs (empty)
%          .cond    condition name or number
%          .name    parameter name
%          .opts    optional settings for parameter
%            .comp  compress parameter, 'log', 'sqrt'
%            .norm  z-normalization of parameter (otherwise mean removed)
%            .ortho normalize convolved regressor against design
%          .pval    parameter values (must match number of onsets)
%        .prtr      TR (in ms, 2000)
%        .rcond     Conditions to remove (rest, 1)
%        .regi      VxR regressors of interest or RTC (not orthogonalized)
%        .regni     VxR regressors of no interest or RTC (orthogonalized)
%        .rnorm     Normalization (default: normalize HRF plateau to 1)
%        .tshift    Temporally shift convoluted regressor (in ms, 0)
%        .type      HRF or FIR model ({'hrf'}, 'fir')
%
% Output fields:
%
%       rtc         RTC object (design matrix)

% Version:  v0.7b
% Build:    7091215
% Date:     Sep-12 2007, 3:44 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

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

% default parameters
defpar = { ...
    'hshape',  'char',   'nonempty', 'twogamma'; ...  % HRF type
    'hpttp',   'double', 'nonempty', 5; ...           % pos. time-to-peak
    'hnttp',   'double', 'nonempty', 15; ...          % neg. time-to-peak
    'hpnr',    'double', 'nonempty', 6; ...           % pos/neg ratio
    'hons',    'double', 'nonempty', 0; ...           % HRF onset
    'hpdsp',   'double', 'nonempty', 1; ...           % pos. dispersion
    'hndsp',   'double', 'nonempty', 1; ...           % neg. dispersion
    'moco',    'double', 'noinfnan', []; ...          % motion correction
    'ndcreg',  'double', 'nonempty', 8; ...           % # of deconv. lags
    'nderiv',  'double', 'noinfnan', []; ...          % add derivatives
    'nvol',    'double', 'nonempty', 200; ...         % number of volumes
    'ortho',   'logical','nonempty', false; ...       % orthogonalization
    'prtr',    'double', 'nonempty', 2000; ...        % TR (ms)
    'rcond',   'double', '',         1; ...           % remove which cond.
    'rnorm',   'double', 'nonempty', 0; ...           % normalization
    'tshift',  'double', 'nonempty', 0; ...           % temporal shift
    'type',    'char',   'nonempty', 'hrf' ...        % model type hrf/fir
};

% normalization factor (prior)
fnorm = 833.65780517;

% check optional paramters
if nargin > 1 && ...
    isstruct(params) && ...
    numel(params) == 1
    params = checkstruct(params, defpar);
else
    params = checkstruct(struct, defpar);
end
if isinf(params.ndcreg(1)) || ...
    isnan(params.ndcreg(1)) || ...
    params.ndcreg(1) < 1 || ...
    params.ndcreg(1) > 30
    ndcreg = 8;
else
    ndcreg = round(params.ndcreg(1));
end
nderiv = unique(params.nderiv(:)');
nderiv(nderiv ~= fix(nderiv) | nderiv < 1 | nderiv > 8) = [];
ndnum = numel(nderiv);
if isinf(params.nvol(1)) || ...
    isnan(params.nvol(1)) || ...
    params.nvol(1) < 1 || ...
    params.nvol(1) > 3000
    nvol = 200;
else
    nvol = round(params.nvol(1));
end
ortho = params.ortho(1);
olist = struct;
if isinf(params.prtr(1)) || ...
    isnan(params.prtr(1)) || ...
    params.prtr(1) < 10 || ...
    params.prtr(1) > 30000
    prtr = 2000;
else
    prtr = round(params.prtr(1));
end
ptype = params.type(:)';
if ~strcmpi(ptype, 'hrf')
    ptype = false;
else
    ptype = true;
end
rcond = params.rcond(:)';
rcond(isinf(rcond) | isnan(rcond) | rcond < 1 | rcond ~= fix(rcond)) = [];
if isinf(params.rnorm(1)) || ...
    isnan(params.rnorm(1))
    params.rnorm = 1;
else
    params.rnorm = params.rnorm(1);
end
if isinf(params.tshift(1)) || ...
    isnan(params.tshift(1))
    tshift = 0;
else
    tshift = round(params.tshift(1));
end

% create empty RTC
hfile2 = BVQXfile('new:rtc');
bc2 = bvqxfile_getcont(hfile2.L);

% get shape of HRF (and derivatives)
if ptype
    hf = hrf(params.hshape, 1e-3, params.hpttp, params.hnttp, ...
        params.hpnr, params.hons, params.hpdsp, params.hndsp);
end

% check protocol for volumes
usevol = 0;
if ~isempty(bc.ResolutionOfTime) && ...
    lower(bc.ResolutionOfTime(1)) == 'v'
    usevol = prtr;
end

% get conditions
cond = bc.Cond;
ncond = numel(cond);

% remove indices
rcond(rcond > ncond) = [];

% check condition parameters
if ~isfield(params, 'params') || ...
   ~isstruct(params.params) || ...
    isempty(params.params) || ...
   ~isfield(params.params, 'cond') || ...
   ~isfield(params.params, 'name') || ...
   ~isfield(params.params, 'pval')
    params.params = ...
        cell2struct(cell(0, 0, 4), {'cond', 'name', 'opts', 'pval'}, 3);
else
    params.params = params.params(:)';
    if ~isfield(params.params, 'opts')
        params.params(1).opts = [];
    end
    cnames = cell(numel(bc.Cond), 1);
    for cc = 1:numel(cnames)
        cnames{cc} = bc.Cond(cc).ConditionName{1};
    end
    for pc = numel(params.params):-1:1
        p = params.params(pc);
        if isempty(p.cond) || ...
           (~ischar(p.cond) && ...
            (~isa(p.cond, 'double') || ...
             numel(p.cond) ~= 1 || ...
             isinf(p.cond) || ...
             isnan(p.cond) || ...
             p.cond < 1 || ...
             p.cond > numel(bc.Cond) || ...
             p.cond ~= fix(p.cond))) || ...
            isempty(p.name) || ...
           ~ischar(p.name) || ...
            isempty(p.pval) || ...
           ~isa(p.pval, 'double') || ...
            any(isinf(p.pval(:)) | isnan(p.pval(:)))
            warning( ...
                'BVQXfile:BadArgument', ...
                'Invalid parameter %d.', ...
                pc ...
            );
            params.params(pc) = [];
            continue;
        end
        if ischar(p.cond)
            cf = find(strcmpi(p.cond, cnames));
            if isempty(cf)
                warning( ...
                    'BVQXfile:BadArgument', ...
                    'Invalid condition name (%s) in parameter %d.', ...
                    p.cond, pc ...
                );
                params.params(pc) = [];
                continue;
            end
            p.cond = cf(1);
        end
        if any(rcond == p.cond)
            params.params(pc) = [];
            continue;
        end
        pvs = size(p.pval);
        npv = size(bc.Cond(p.cond).OnOffsets, 1);
        svf = find(pvs == npv);
        if isempty(svf)
            warning( ...
                'BVQXfile:BadArgument', ...
                'Invalid number of parameter values for condition %s.', ...
                cnames(p.cond) ...
            );
            params.params(pc) = [];
        end
        if svf(1) > 1
            p.pval = reshape(permute( ...
                p.pval, [svf(1), setdiff(1:ndims(p.pval), svf(1))]), ...
                [npv, numel(p.pval) / npv]);
        end
        if size(p.pval, 2) > 3
            warning( ...
                'BVQXfile:BadArgument', ...
                'Only up to three sub-parameters with interaction supported.' ...
            );
            params.params(pc) = [];
            continue;
        elseif size(p.pval, 2) == 3
            p.pval(:, 7) = p.pval(:, 1) .* p.pval(:, 2) .* p.pval(:, 3);
            p.pval(:, 4) = p.pval(:, 1) .* p.pval(:, 2);
            p.pval(:, 5) = p.pval(:, 1) .* p.pval(:, 3);
            p.pval(:, 6) = p.pval(:, 2) .* p.pval(:, 3);
        elseif size(p.pval, 2) == 2
            p.pval(:, 3) = p.pval(:, 1) .* p.pval(:, 2);
        end
        if ~isstruct(p.opts)
            p.opts = struct( ...
                'comp' , 'no', ...
                'norm' , true, ... ...
                'ortho', false);
        end
        if ~isfield(p.opts, 'comp') || ...
           ~ischar(p.opts.comp) || ...
           ~any(strcmpi(p.opts.comp, {'log', 'no', 'sqrt'}))
            p.opts.comp = 'no';
        else
            p.opts.comp = lower(p.opts.comp(:)');
        end
        if ~isfield(p.opts, 'norm') || ...
           ~islogical(p.opts.norm) || ...
            numel(p.opts.norm) ~= 1
            p.opts.norm = true;
        end
        if ~isfield(p.opts, 'ortho') || ...
           ~islogical(p.opts.ortho) || ...
            numel(p.opts.ortho) ~= 1
            p.opts.ortho = false;
        end
        if ~strcmp(p.opts.comp, 'no')
            pvs = sign(p.pval);
            p.pval = abs(p.pval);
            if strcmp(p.opts.comp, 'log')
                if any(p.pval(:) < 1)
                    p.pval = p.pval + 1;
                end
                p.pval = pvs .* log(p.pval);
            else
                p.pval = pvs .* sqrt(p.pval);
            end
        end
        p.pval = p.pval - repmat(mean(p.pval), [size(p.pval, 1), 1]);
        if p.opts.norm
            p.pval = p.pval ./ repmat(std(p.pval), [size(p.pval, 1), 1]);
        end
        if size(p.pval, 2) > 1
            rep = repmat(p, [1, size(p.pval, 2)]);
            for rpc = 1:numel(rep)
                rep(rpc).pval = rep.(rpc).pval(:, rpc);
                if numel(rep) == 3
                    if rpc < 3
                        rep(rpc).name = ...
                            sprintf('%s - p%d', rep(rpc).name, rpc);
                    else
                        rep(rpc).name = sprintf('%s - p1Xp2', rep(rpc).name);
                    end
                else
                    switch (rpc)
                        case {1, 2, 3}
                            rep(rpc).name = ...
                                sprintf('%s - p%d', rep(rpc).name, rpc);
                        case {4}
                            rep(rpc).name = sprintf('%s - p1Xp2', rep(rpc).name);
                        case {5}
                            rep(rpc).name = sprintf('%s - p1Xp3', rep(rpc).name);
                        case {6}
                            rep(rpc).name = sprintf('%s - p2Xp3', rep(rpc).name);
                        case {7}