vtc_crosscorrelate.m

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

function ccvmp = vtc_CrossCorrelate(hfile, hfile2, opts)
% VTC::CrossCorrelate  - create CC map of two VTCs
%
% FORMAT:       ccvmp = vtc.CrossCorrelate(vtc2 [, opts])
%
% Input fields:
%
%       vtc2        second VTC (must match in dims and layout)
%       opts        options settings
%        .reverse   reverse time courses of second VTC
%
% Output fields:
%
%       ccvmp       cross-correlation r-VMP
%
% Note: the toolbox internal cov_nd function is used which gives
%       slightly different r values than corrcoef.

% Version:  v0.7b
% Build:    7091116
% Date:     Sep-11 2007, 4:05 PM 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, 'vtc') || ...
    numel(hfile2) ~= 1 || ...
   ~isBVQXfile(hfile2, 'vtc')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bb1 = aft_BoundingBox(hfile);
bb2 = aft_BoundingBox(hfile2);
if any(bb1.BBox(:) ~= bb2.BBox(:)) || ...
    any(bb1.DimXYZ ~= bb2.DimXYZ)
    error( ...
        'BVQXfile:BadArgument', ...
        'Dimension/Layout mismatch.' ...
    );
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'reverse') || ...
    isempty(opts.reverse) || ...
   ~islogical(opts.reverse)
    opts.reverse = false;
else
    opts.reverse = opts.reverse(1);
end
if opts.reverse
    revstr = ' TC-reversed';
else
    revstr = '';
end

% get contents
sc1 = bvqxfile_getscont(hfile.L);
bc1 = sc1.C;
sc2 = bvqxfile_getscont(hfile2.L);
bc2 = sc2.C;

% create vmp
ccvmp = newnatresvmp(bb1.BBox, bc1.Resolution, 2);
bcvmp = bvqxfile_getcont(ccvmp.L);
bcvmp.OriginatingVTC = sc1.F;
bcvmp.LinkedPRT = bc1.NameOfLinkedPRT;
bcvmp.Map.Name = sprintf('CC %s <-> %s%s', sc1.F, sc2.F, revstr);
bcvmp.Map.DF1 = size(bc1.VTCData, 1) - 1;

% iterate over last spatial dim
for z = 1:size(bc1.VTCData, 4);
    
    % get components for cov_nd
    r1 = double(permute(bc1.VTCData(:, :, :, z), [2, 3, 1, 4]));
    if ~opts.reverse
        r2 = double(permute(bc2.VTCData(:, :, :, z), [2, 3, 1, 4]));
    else
        r2 = double(permute(bc2.VTCData(end:-1:1, :, :, z), [2, 3, 1, 4]));
    end
    
    % compute r value
    [cc, cr] = cov_nd(r1, r2);
    cr(isnan(cr)) = 0;
    bcvmp.Map.VMPData(:, :, z) = single(cr);
end

% set data to VMP
bcvmp.Map.BonferroniValue = sum(bcvmp.Map.VMPData(:) ~= 0);
fdrl = [0.1, 0.05, 0.04, 0.03, 0.02, 0.01, 0.005, 0.001];
fdrt2 = applyfdr(double(bcvmp.Map.VMPData), 'r', fdrl(:), bcvmp.Map.DF1, 0, true);
bcvmp.Map.FDRThresholds = [fdrl(:), fdrt2];
bcvmp.Map.NrOfFDRThresholds = numel(fdrl);
bcvmp.Map.LowerThreshold = fdrt2(end);
bcvmp.Map.UpperThreshold = min(1, 2 * fdrt2(end));
bvqxfile_setcont(ccvmp.L, bcvmp);