show_ic_meg

 Show ICs

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

0001 function show_ic_meg(p)
0002 % Show ICs
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 disp(mfilename);
0008 
0009 % Classification program
0010 correct_ic_prog='correct_classify_ic_meg';
0011 
0012 % Make directory to save figures
0013 save_dir = fullfile(p.fig_root, mfilename, p.sub);
0014 if exist(save_dir, 'dir') ~= 7
0015     vb_mkdir(save_dir);
0016 end
0017 
0018 % Load MEG data
0019 data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, 'all.info.mat');
0020 [meg, ~, time_info] = vb_load_meg_data(data_file);
0021 [Nch, Nt, Ntr] = size(meg);
0022 samplingrate = time_info.sample_frequency;
0023 time = time_info.time;
0024 
0025 % Load positions of channels
0026 [pos, channel_info] = vb_load_channel(data_file);
0027 ix_axial = find(channel_info.Type == 2);% axial sensor
0028 pos = pos(ix_axial, :);
0029 
0030 % Load ICA result
0031 ica_dir = fullfile(p.proj_root, p.meg_dirname, p.ica_dirname);
0032 ica_result = fullfile(ica_dir, [correct_ic_prog '.mat']);
0033 load(ica_result,'w','s','v','unmix','mix','brain_ic','noise_ic')
0034 
0035 % Calculate ICs
0036 ic = zeros(Nch, Nt, Ntr);
0037 for tr = 1:Ntr
0038     ic(:,:,tr) = unmix*meg(:,:,tr);  % [component x sample x trial]
0039 end
0040 
0041 % Calculate stimulu-triggered average and power spectrum
0042 mic = mean(ic, 3);% Average each IC across trials
0043 [pow, f] = vb_fftpower(ic, samplingrate);% Calculate power spectrum
0044 mp = mean(log(pow), 3);% Average power spectrum across trials
0045 
0046 % Obtain classification
0047 classification = -ones(1, Nch);
0048 classification(brain_ic) = 1;% 1: brain_ic, -1: noise_ic
0049 
0050 % Show and save ICs
0051 h = figure;
0052 cc = 1;
0053 while cc <= Nch
0054     clf(h)
0055     set(0, 'CurrentFigure', h)
0056     k = 1;% Position of figure
0057     for c = cc:min([cc+4,Nch]);
0058         if classification(c) == 1
0059             color = 'b';
0060         else
0061             color = 'r';
0062         end
0063         
0064         % Plot spatial pattern of mixing vector
0065         subplot(5, 3, 1+3*(k-1))
0066         vb_plot_sensor_2d(pos, mix(ix_axial,c));
0067         vb_plot_sensor_2d_head_plot_add(gca);
0068         axis square off
0069         title(['IC ' num2str(c)])
0070         
0071         % Plot stimulus-triggered average of components
0072         subplot(5, 3, 2+3*(k-1))
0073         plot(time, mic(c, :), color);
0074         xlim([-0.5 1])
0075         if k == 1
0076             title('Stimulus-triggered average')
0077         elseif k == 5
0078             xlabel('Time [s]')
0079         end
0080         
0081         % Plot power spectrum
0082         subplot(5, 3, 3+3*(k-1))
0083         plot(f, mp(c, :), color);
0084         xlim([0 50])
0085         if k == 1
0086             title('Power spectrum (log)')
0087         elseif k == 5
0088             xlabel('Frequency [Hz]')
0089         end
0090         k = k+1;
0091     end
0092     ax1 = axes('Position', [0 0 1 1], 'Visible', 'off');
0093     text1 = {'Blue : Brain'; 'Red : Noise'};
0094     axes(ax1)
0095     text(0.93, 0.87, text1)
0096     
0097     % Save figure
0098     vb_savefig_as_shown(h, fullfile(save_dir, [num2str(cc) '-' num2str(cc+4)]))
0099     cc = cc+5;
0100 end
0101