vb_megfile_remove_environmental_noise

 Remove environmental noise by using temporal PCA of reference sensor signal
   
 --- Input
   megfile : .meg.mat (including MEG and REFERENCE data)
   new_meg_file : denoised meg file (.meg.mat) to be saved
   Twin_sec = [Nbck , Nfwd] : time shift length for backward and forward
   direction [sec]
     if Twin is scalar, Nbck = Nfwd = Twin is assumed
 --- Optional Input
   smin : >= 1 :  number of components for threshould PCA (default 40)
           < 1  : contribution rate for threshould PCA

 2008-1-25 Masa-aki Sato
 2015-4-20 Yusuke Fujiwara
 2015-6-12 Yusuke Takeda (Modified for pipeline of VBMEG)

0001 function vb_megfile_remove_environmental_noise(meg_file,new_meg_file,Twin_sec,smin)
0002 % Remove environmental noise by using temporal PCA of reference sensor signal
0003 %
0004 % --- Input
0005 %   megfile : .meg.mat (including MEG and REFERENCE data)
0006 %   new_meg_file : denoised meg file (.meg.mat) to be saved
0007 %   Twin_sec = [Nbck , Nfwd] : time shift length for backward and forward
0008 %   direction [sec]
0009 %     if Twin is scalar, Nbck = Nfwd = Twin is assumed
0010 % --- Optional Input
0011 %   smin : >= 1 :  number of components for threshould PCA (default 40)
0012 %           < 1  : contribution rate for threshould PCA
0013 %
0014 % 2008-1-25 Masa-aki Sato
0015 % 2015-4-20 Yusuke Fujiwara
0016 % 2015-6-12 Yusuke Takeda (Modified for pipeline of VBMEG)
0017 
0018 % Set parameters
0019 if length(Twin_sec)==1, Twin_sec = [Twin_sec, Twin_sec]; end;
0020 if length(Twin_sec)>2, error('Time shift parameter error'); end;
0021 if ~exist('smin','var'),  smin = 40; end;
0022 LINE_LEN = 15;
0023 
0024 % Prepare directory to save binary files
0025 new_bin_dir=new_meg_file(1:strfind(new_meg_file,'.meg.mat')-1);
0026 if exist(new_bin_dir, 'dir') ~= 7
0027     vb_mkdir(new_bin_dir);
0028 end
0029 
0030 % Load MEG channel information
0031 ch_labels = vb_megfile_get_channel_label_meg(meg_file);
0032 ch_len = size(ch_labels,1);
0033 meg_info  = vb_load_measurement_info(meg_file);
0034 precision = vb_meginfo_get_precision(meg_info);
0035 
0036 % Calculate temporal PCA of reference sensor data
0037 Twin=Twin_sec*meg_info.SampleFreq;
0038 loadspec_ref.ChannelType = 'REFERENCE';
0039 ref = vb_load_meg_data(meg_file,loadspec_ref);
0040 [zref,tx] = vb_temporal_pca(ref, Twin, smin);
0041 
0042 % Start denoising
0043 fprintf('Start denoising using reference sensor\n')
0044 for i_ch=1:ch_len
0045     fprintf('[%s]', ch_labels{i_ch});
0046     if ~rem(i_ch, LINE_LEN), fprintf('\n'); end
0047     
0048     load_spec.ChannelName = {ch_labels{i_ch}};
0049     load_spec.ChannelSwitch = true; % ChannelName is to read
0050     load_spec.ChannelType = 'MEG';
0051     load_spec.ActiveChannel = 0; % all channels
0052     
0053     % Denoise
0054     meg = vb_load_meg_data(meg_file, load_spec);
0055     denoised_meg=meg;
0056     denoised_meg(tx) = meg(tx)-(meg(tx)*zref')*zref;
0057     
0058     % Save denoised data
0059     new_file = sprintf('%s/%s.ch.meg.dat', new_bin_dir, ch_labels{i_ch});
0060     fid = fopen(new_file, 'wb');
0061     if fid == -1
0062         warning('(%s) cannot open file : %s\n', mfilename, new_file);
0063         continue;
0064     end
0065     fwrite(fid, denoised_meg, precision);
0066     fclose(fid);
0067 end
0068 
0069 % Copy and paste other files
0070 ch_labels_ext = vb_megfile_get_channel_label_extra(meg_file);
0071 ch_labels_ref = vb_megfile_get_channel_label_refmg(meg_file);
0072 ch_labels_extra = [ch_labels_ext; ch_labels_ref];
0073 ch_len_extra = size(ch_labels_extra, 1);
0074 
0075 for i_ch = 1:ch_len_extra
0076     load_spec.ChannelName = {ch_labels_extra{i_ch}};
0077     load_spec.ChannelSwitch = true; % ChannelName is to read
0078     load_spec.ChannelType = 'ALL';
0079     load_spec.ActiveChannel = 0; % all channels
0080     
0081     % Do nothing
0082     meg = vb_load_meg_data(meg_file, load_spec);
0083     
0084     % Save data
0085     new_file = sprintf('%s/%s.ch.meg.dat', new_bin_dir, ch_labels_extra{i_ch});
0086     fid = fopen(new_file, 'wb');
0087     if fid == -1
0088         warining('(%s) cannot open file : %s\n', mfilename, new_file);
0089         continue;
0090     end
0091     fwrite(fid, meg, precision);
0092     fclose(fid);
0093 end
0094 
0095 load(meg_file)
0096 a=strfind(new_bin_dir,'/');
0097 MEGinfo.saveman.data_dir = ['.' new_bin_dir(a(end):end)];
0098 bexp         = [];
0099 bexp_ext     = [];
0100 refmg        = [];
0101 fprintf('\n');
0102 fprintf('save: %s\n',new_meg_file);
0103 vb_fsave(new_meg_file, 'MEGinfo','bexp','bexp_ext','refmg', ...
0104   'CoordType','Measurement','PositionFile',...
0105   'Qpick','pick','ref_Qpick','ref_pick');