check_removing_eog_from_meg

 Check the result of removing EOG components from MEG data

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

0001 function check_removing_eog_from_meg(p)
0002 % Check the result of removing EOG components from MEG data
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 % Set parameters
0008 file = 3;% First run of visual task
0009 
0010 % Set input files
0011 original_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, [p.file_list{file} '.meg.mat']);
0012 removed_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, ['r_' p.file_list{file} '.meg.mat']);
0013 
0014 % Load MEG data
0015 [o, ~, time_info] = vb_load_meg_data(original_file);
0016 r = vb_load_meg_data(removed_file);
0017 
0018 % Load positions of channels
0019 [pos, channel_info] = vb_load_channel(original_file);
0020 ix_axial = find(channel_info.Type == 2);% axial sensor
0021 pos = pos(ix_axial, :);
0022 o = o(ix_axial, :, :);
0023 r = r(ix_axial, :, :);
0024 
0025 % Calculate power
0026 po=sum(sum(o.^2,2),3);
0027 pr=sum(sum(r.^2,2),3);
0028 dp=po-pr;
0029 
0030 % Get minimum and maximum values of power
0031 ma=max(po(:));
0032 mi=min(pr(:));
0033 
0034 % Show power of original MEG data
0035 h = figure;
0036 subplot(2, 2, 1)
0037 vb_plot_sensor_2d(pos, po, [mi ma]);
0038 vb_plot_sensor_2d_head_plot_add(gca);
0039 axis square off
0040 colorbar
0041 title('Power (before)')
0042 
0043 % Show power of EOG-removed MEG data
0044 subplot(2, 2, 2)
0045 vb_plot_sensor_2d(pos, pr, [mi ma]);
0046 vb_plot_sensor_2d_head_plot_add(gca);
0047 axis square off
0048 colorbar
0049 title('Power (after)')
0050 
0051 % Show an example of original and EOG-removed MEG data
0052 [~, ch] = max(dp);% Select channel
0053 d = squeeze(o(ch, :, :)-r(ch, :, :));
0054 pd = sum(d.^2, 1);
0055 [~, tr] = max(pd);% Select trial
0056 
0057 subplot(2,2,3)
0058 plot(time_info.time, [o(ch, :, tr)', r(ch, :, tr)'])
0059 xlim([time_info.time(1) time_info.time(end)])
0060 legend('Before','After')
0061 xlabel('Time [s]')
0062 ylabel('Magnetic field [T]')
0063 title(['Channel ' channel_info.Name{ch} ', Trial ' num2str(tr)])
0064 
0065 % Show location of selected channel
0066 subplot(2,2,4)
0067 x = zeros(length(ix_axial), 1);
0068 x(ch) = 1;
0069 vb_plot_sensor_2d(pos, x);
0070 vb_plot_sensor_2d_head_plot_add(gca);
0071 axis square off
0072 title(['Channel ' channel_info.Name{ch}])
0073 
0074 % Make directory to save figures
0075 save_dir = fullfile(p.fig_root, mfilename);
0076 if exist(save_dir, 'dir') ~= 7
0077     vb_mkdir(save_dir);
0078 end
0079 
0080 % Save figure
0081 vb_savefig_as_shown(h, fullfile(save_dir, p.sub))