estimate_source_current_meeg

 Estimate source currents from MEG data

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function estimate_source_current_meeg(p)
0002 % Estimate source currents from MEG data
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 disp(mfilename);
0008 
0009 % Set parameters
0010 time_signal = [0 0.3]; % Period for stimulus-evoked activity [s]
0011 time_noise = [-0.3 0]; % Period for baseline [s]
0012 
0013 % Make directory to save current
0014 current_dir = fullfile(p.proj_root, p.meeg_dirname, p.current_dirname);
0015 if exist(current_dir, 'dir') ~= 7
0016     mkdir(current_dir);
0017 end
0018 
0019 for ta = 1:length(p.task_list)
0020     
0021     %% Estimate current variance
0022     clear bp
0023     bayes_parm = vb_set_bayes_default_parameters;
0024     
0025     % Set input files (relative path from proj_root)
0026     bayes_parm.brainfile = fullfile(p.brain_dirname, [p.struct_name, '.brain.mat']);
0027     bayes_parm.areafile = fullfile(p.brain_dirname, [p.struct_name, '.area.mat']);
0028     bayes_parm.actfile = fullfile(p.brain_dirname, [p.struct_name, '.act.mat']);
0029     bayes_parm.eegfile{1} = fullfile(p.eeg_dirname, p.trial_dirname, ['crri_' p.task_list{ta} '.info.mat']);
0030     bayes_parm.megfile{1} = fullfile(p.meg_dirname, p.trial_dirname, ['rri_' p.task_list{ta} '.info.mat']);
0031     bayes_parm.basisfile_eeg = fullfile(p.eeg_dirname, p.leadfield_dirname, [p.struct_name '.basis.mat']);
0032     bayes_parm.basisfile_meg = fullfile(p.meg_dirname, p.leadfield_dirname, [p.task_list{ta} '.basis.mat']);
0033     bayes_parm.eegfile_baseline = bayes_parm.eegfile;
0034     bayes_parm.megfile_baseline = bayes_parm.megfile;
0035     
0036     % Specify current area and fMRI activity
0037     bayes_parm.area_key  = 'Cortex';% Current area
0038     bayes_parm.act_key = p.prior{ta};% Prior for current variance
0039     bayes_parm.prior_weight = p.prior_weight(ta);% Relative influence of prior information (0-1)
0040     
0041     % Load time information
0042     [~, ~, time_info] = vb_load_meg_data(fullfile(p.proj_root, bayes_parm.megfile{1}));
0043     time = time_info.time;
0044     
0045     % Set time window for stimulus-evoked activity
0046     [~, from] = min(abs(time-time_signal(1)));
0047     [~, to] = min(abs(time-time_signal(2)));
0048     bayes_parm.twin_meg = [from to]; % Time window for analysis
0049     bayes_parm.Tperiod = to-from+1; % Time period for current estimation
0050     bayes_parm.Tnext = to-from+1; % Time step for next period
0051     
0052     % Set time window for baseline
0053     [~, from] = min(abs(time-time_noise(1)));
0054     [~, to] = min(abs(time-time_noise(2)));
0055     bayes_parm.twin_noise = [from to];
0056     bayes_parm.twin_baseline = [from to];
0057 
0058     bayes_parm.patch_norm = ON;% Patch area normalization is applied
0059     
0060     % Set output file
0061     bayes_parm.bayesfile = fullfile(p.meeg_dirname, p.current_dirname, [p.task_list{ta} '.bayes.mat']);
0062     
0063     % Set parameters for noise covariance estimation
0064     bayes_parm = vb_set_noise_estimation_model(bayes_parm, p.proj_root);
0065     
0066     % Estimate current variance
0067     vb_job_vb_meeg(p.proj_root, bayes_parm)
0068     
0069     
0070     %% Estimate source current
0071     
0072     % Set input files (relative path from proj_root)
0073     current_parm.eegfile = bayes_parm.eegfile;
0074     current_parm.megfile = bayes_parm.megfile;
0075     current_parm.bayesfile = bayes_parm.bayesfile;
0076     
0077     % Trial-averaged current will be estimated.
0078     current_parm.trial_average = ON;
0079     
0080     % Set time window for estimating current
0081     a = vb_load_meg_data(fullfile(p.proj_root, current_parm.megfile{1}));
0082     current_parm.twin_meg = [1 size(a, 2)];
0083     current_parm.Tperiod = size(a, 2);
0084     
0085     % Set output file (relative path from proj_root)
0086     current_parm.currfile = fullfile(p.meeg_dirname, p.current_dirname, [p.task_list{ta} '.curr.mat']);
0087     
0088     % Estimate source current
0089     vb_job_current_meeg(p.proj_root, current_parm);
0090 end
0091 
0092 
0093 
0094 
0095