common_average_eeg

 Take common average

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

0001 function common_average_eeg(p)
0002 % Take common average
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 paremeters
0010 prefix = 'c';% Prefix of resultant EEG data
0011 
0012 % Load fileinfo
0013 data_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, 'rri_all.info.mat');
0014 load(data_file,'fileinfo')
0015 Nfile=length(fileinfo.filename);
0016 
0017 % Load indices of channels to take common average
0018 [~, ch_info, time_info] = vb_load_meg_data(data_file);
0019 Nch = length(ch_info.Active);
0020 load(fileinfo.filename{1}, 'EEGinfo')
0021 ix_in_eeg_data = zeros(1, Nch);
0022 for ii = 1:Nch
0023     ix_in_eeg_data(ii) = find(EEGinfo.ChannelInfo.ID == ch_info.ID(ii));
0024 end
0025 
0026 for file=1:Nfile
0027     
0028     % Load EEG data
0029     old_file=fileinfo.filename{file};
0030     load(old_file)
0031     
0032     % Take common reference
0033     fparm.common_flg=1;
0034     fparm.freq=EEGinfo.SampleFrequency;
0035     for tr=1:size(eeg_data,3)
0036         eeg_data(ix_in_eeg_data,:,tr)=vb_filter_raw_data(eeg_data(ix_in_eeg_data,:,tr), fparm);
0037     end
0038     
0039     % Save denoised EEG data
0040     ix = strfind(old_file, p.trial_dirname);
0041     Nc =length(p.trial_dirname);
0042     new_file = [old_file(1:ix+Nc) prefix old_file(ix+Nc+1:end)];
0043     copyfile(old_file, new_file);
0044     vb_save(new_file, 'eeg_data');
0045     
0046     % Update fileinfo
0047     fileinfo.filename{file} = new_file;
0048 end
0049 
0050 % Save new fileinfo file
0051 load(data_file, 'Measurement')
0052 save_file = fullfile(p.proj_root, p.eeg_dirname, p.trial_dirname, [prefix 'rri_all.info.mat']);
0053 save(save_file, 'fileinfo', 'Measurement')
0054 
0055 % Chack if average across channels is 0
0056 if 0
0057     a = vb_load_meg_data(save_file);
0058     ma = squeeze(mean(a,1));
0059     figure
0060     plot(time_info.time, ma)
0061     xlim([time_info.time(1) time_info.time(end)])
0062     xlabel('Time [s]')
0063 end
0064