Home > vbmeg > demo > tutorial_for_vbmeg2 > advanced > show_trial_average_meg.m

show_trial_average_meg

PURPOSE ^

Show trial-averaged MEG data

SYNOPSIS ^

function show_trial_average_meg(p)

DESCRIPTION ^

 Show trial-averaged MEG data

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function show_trial_average_meg(p)
0002 % Show trial-averaged 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_to_show = [-0.1 0.4];% sec
0011 time_of_interest = [0 0.24];% sec
0012 Npeak = 2;% Number of peaks at which spatial patterns is shown
0013 
0014 % Make directory to save figures
0015 save_dir = fullfile(p.fig_root, mfilename, p.sub);
0016 if exist(save_dir, 'dir') ~= 7
0017     vb_mkdir(save_dir);
0018 end
0019 
0020 % Load time information
0021 data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, 'rri_all.info.mat');
0022 [~, ~, time_info] = vb_load_meg_data(data_file);
0023 [~, from_show] = min(abs(time_info.time-time_to_show(1)));
0024 [~, to_show] = min(abs(time_info.time-time_to_show(2)));
0025 time = time_info.time(from_show: to_show);
0026 
0027 % Load positions of channels
0028 [pos, channel_info] = vb_load_channel(data_file);
0029 ix_axial = find(channel_info.Type == 2);% axial sensor
0030 
0031 % Show trial-averaged MEG data
0032 for ta = 1:length(p.task_list)
0033     
0034     % Set input file
0035     data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, ...
0036         ['rri_' p.task_list{ta} '.info.mat']);
0037     
0038     % Load MEG data
0039     data = vb_load_meg_data(data_file);
0040     
0041     % Average MEG data across trials
0042     y = mean(data(:, from_show:to_show, :), 3);
0043     
0044     % Detect peaks of power
0045     [~, from_peak] = min(abs(time-time_of_interest(1)));
0046     [~, to_peak] = min(abs(time-time_of_interest(2)));
0047     Nsample = to_peak-from_peak+1;
0048     py = mean(y(:,from_peak:to_peak).^2, 1);
0049     dpy = diff(py);
0050     ix = find(dpy(1:Nsample-2) > 0 & dpy(2:Nsample-1) < 0)+1;
0051     ix2 = sortrows([ix(:) py(ix)'], 2);
0052     ix_peak = sort(ix2(end-Npeak+1:end, 1))+from_peak-1;% index of peak time
0053     
0054     % Plot time series of averaged MEG data
0055     h = figure;
0056     max_y = max(y(:));
0057     min_y = min(y(:));
0058     subplot(2, 1, 1)
0059     plot(time, y)
0060     hold on
0061     for peak = 1:Npeak
0062         ix = ix_peak(peak);
0063         plot([time(ix) time(ix)], [min_y max_y])
0064     end
0065     xlabel('Time [s]')
0066     ylabel('Magnetic field [T]')
0067     axis([time(1) time(end) min_y max_y])
0068     
0069     % Plot spatial patterns at peaks
0070     for peak = 1:Npeak
0071         ix = ix_peak(peak);
0072         subplot(2, Npeak, Npeak+peak)
0073         vb_plot_sensor_2d(pos(ix_axial, :), y(ix_axial, ix), [min_y max_y]);
0074         vb_plot_sensor_2d_head_plot_add(gca);
0075         axis square off
0076         colorbar
0077         title([num2str(time(ix)) ' s'])
0078     end
0079     
0080     % Save figure
0081     vb_savefig_as_shown(h, fullfile(save_dir, p.task_list{ta}))
0082 end
0083

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005