remove_noise_ic_meg

 Remove noise ICs from MEG data

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

0001 function remove_noise_ic_meg(p)
0002 % Remove noise ICs 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 % Classification program
0010 correct_ic_prog = 'correct_classify_ic_meg';
0011 
0012 % Set paremeters
0013 prefix = 'i';% Prefix of denoised MEG data
0014 
0015 % Load MEG data and MEGinfo
0016 data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, 'all.info.mat');
0017 load(data_file, 'fileinfo')
0018 Nfile = length(fileinfo.filename);
0019 [meg, ch_info] = vb_load_meg_data(data_file);
0020 [Nch, Nt, Ntr] = size(meg);
0021 
0022 % Load ICA result
0023 ica_dir = fullfile(p.proj_root, p.meg_dirname, p.ica_dirname);
0024 ica_result = fullfile(ica_dir, [correct_ic_prog '.mat']);
0025 load(ica_result,'unmix','mix','brain_ic')
0026 
0027 % Make denoising filter
0028 g = zeros(Nch);
0029 for ic = 1:length(brain_ic)
0030     g(brain_ic(ic), brain_ic(ic)) = 1;
0031 end
0032 denoise = mix*g*unmix;
0033 
0034 % Remove noise ICs
0035 for file=1:Nfile
0036 
0037     % Load MEG data
0038     old_file=fileinfo.filename{file};
0039     load(old_file)
0040     
0041     % Obtain indices of channels
0042     clear IA
0043     for ii = 1 : length(ch_info.ID)
0044         IA(ii) = find(MEGinfo.ChannelInfo.ID == ch_info.ID(ii));
0045     end
0046     
0047     % Denoise
0048     for tr = 1:size(bexp, 3)
0049         bexp(IA, :, tr) = denoise*bexp(IA, :, tr);
0050     end
0051     
0052     % Save denoised MEG data
0053     ix = findstr(old_file, p.trial_dirname);
0054     Nc =length(p.trial_dirname);
0055     new_file = [old_file(1:ix+Nc) prefix old_file(ix+Nc+1:end)];
0056     copyfile(old_file, new_file);
0057     vb_save(new_file, 'bexp');
0058     
0059     % Update fileinfo.filename
0060     fileinfo.filename{file} = new_file;
0061 end
0062 
0063 % Save fileinfo
0064 load(data_file, 'Measurement')
0065 save_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, [prefix '_all.info.mat']);
0066 save(save_file, 'fileinfo', 'Measurement')
0067 
0068