Home > functions > job > vb_job_filter_trial_data.m

vb_job_filter_trial_data

PURPOSE ^

band/lowpass filtering and bias removal for trial MEG/EEG data file

SYNOPSIS ^

function [data] = vb_job_filter_trial_data(root_dir, parm)

DESCRIPTION ^

  band/lowpass filtering and bias removal for trial MEG/EEG data file
    vb_job_filter_trial_data(parm)
    vb_job_filter_trial_data(root_dir, parm)
 --- Input
 root_dir       : root directory of data files
 parm.data_file : input data  file path/name
 parm.new_file  : output data file path/name
  
 --- preprocess parameter
 parm
  .bias_flg    : Bias correction flag  
               = 1: Bias correction for each channel & trial by mean
               = 2: Linear trend correction
               = [t1 t2]: Bias correction by time window [t1 t2]
  .fsamp       : Down sampling frequency [Hz]
                 lowpass or bandpass below Nyquist frequency is necessarily
  .common_flg  : Common reference flag (= 1/0 : ON/OFF)
  .trial       : trial index if trial extraction is required
  .trial_mean  = 1: trend mean is saved after filtering each trend

 --- If these fields are not given or empty,
        corresponding process are not applied
 --- Multiple filtering specification
 --- n-th filter specification
  .cutfreq{n}     : Cutoff frequency [Hz]
               = [f1]    for 'low','high'
               = [f1 f2] for 'band','stop'
  .filter_type{n} : 'low','high','band','stop'
     if this field is not given or empty
     filter_type = 'low'  for length(cutfreq) = 1
     filter_type = 'band' for length(cutfreq) = 2
  .online(n) : filter order of IIR bandpass filter
      = 0: FIR filter (eegfilt) [Default]
      > 0: Butterworth filter (filtfilt)
      < 0: Butterworth filter (online)

 2009-10-10 Masa-aki Sato
   2011-03-02 (Sako) added checking old data fields exist or not
   2011-06-01 (Sako) converted return values of vb_load_device to upper case

 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    [data] = vb_job_filter_trial_data(root_dir, parm)
0002 %  band/lowpass filtering and bias removal for trial MEG/EEG data file
0003 %    vb_job_filter_trial_data(parm)
0004 %    vb_job_filter_trial_data(root_dir, parm)
0005 % --- Input
0006 % root_dir       : root directory of data files
0007 % parm.data_file : input data  file path/name
0008 % parm.new_file  : output data file path/name
0009 %
0010 % --- preprocess parameter
0011 % parm
0012 %  .bias_flg    : Bias correction flag
0013 %               = 1: Bias correction for each channel & trial by mean
0014 %               = 2: Linear trend correction
0015 %               = [t1 t2]: Bias correction by time window [t1 t2]
0016 %  .fsamp       : Down sampling frequency [Hz]
0017 %                 lowpass or bandpass below Nyquist frequency is necessarily
0018 %  .common_flg  : Common reference flag (= 1/0 : ON/OFF)
0019 %  .trial       : trial index if trial extraction is required
0020 %  .trial_mean  = 1: trend mean is saved after filtering each trend
0021 %
0022 % --- If these fields are not given or empty,
0023 %        corresponding process are not applied
0024 % --- Multiple filtering specification
0025 % --- n-th filter specification
0026 %  .cutfreq{n}     : Cutoff frequency [Hz]
0027 %               = [f1]    for 'low','high'
0028 %               = [f1 f2] for 'band','stop'
0029 %  .filter_type{n} : 'low','high','band','stop'
0030 %     if this field is not given or empty
0031 %     filter_type = 'low'  for length(cutfreq) = 1
0032 %     filter_type = 'band' for length(cutfreq) = 2
0033 %  .online(n) : filter order of IIR bandpass filter
0034 %      = 0: FIR filter (eegfilt) [Default]
0035 %      > 0: Butterworth filter (filtfilt)
0036 %      < 0: Butterworth filter (online)
0037 %
0038 % 2009-10-10 Masa-aki Sato
0039 %   2011-03-02 (Sako) added checking old data fields exist or not
0040 %   2011-06-01 (Sako) converted return values of vb_load_device to upper case
0041 %
0042 % Copyright (C) 2011, ATR All Rights Reserved.
0043 % License : New BSD License(see VBMEG_LICENSE.txt)
0044 
0045 if nargin==1,
0046     parm = root_dir;
0047     datafile = [parm.data_file];
0048     filename = [parm.new_file];
0049 else
0050     root_dir = vb_rm_trailing_slash(root_dir);
0051     datafile = [root_dir '/' parm.data_file];
0052     filename = [root_dir '/' parm.new_file];
0053 end
0054 
0055 % load trial data
0056 [data, channel_info] = vb_load_meg_data(datafile);
0057 
0058 % channel index of data
0059 ch_id = channel_info.ID; 
0060 
0061 % Sample Frequency of Original data
0062 info  = vb_load_meg_info(datafile);
0063 parm.freq  = info.SampleFreq;
0064 
0065 % trial extraction
0066 if isfield(parm,'trial') && ~isempty(parm.trial)
0067     if any(parm.trial < 1 | parm.trial > size(data,3)),
0068         error('trial index is not correct')
0069     end
0070     data = data(:,:,parm.trial);
0071 end
0072 
0073 %tic
0074 if isfield(parm,'cutfreq') 
0075     if iscell(parm.cutfreq)
0076         Nfilter = length(parm.cutfreq);
0077         % Bias remove
0078         Parm = parm;
0079         Parm.cutfreq     = [];
0080         Parm.filter_type = [];
0081         Parm.online      = [];
0082         Parm.fsamp       = [];
0083         Parm.common_flg  = []; 
0084         % Parm.bias_flg    = [];
0085         
0086         [data] = vb_filter_trial_data(data,Parm);
0087         
0088         Parm.bias_flg    = []; 
0089         % Multiple filtering operation
0090         for n=1:Nfilter
0091             Parm.cutfreq     = parm.cutfreq{n};
0092             Parm.filter_type = parm.filter_type{n} ;
0093             Parm.online      = parm.online(n) ;
0094             [data] = vb_filter_trial_data(data,Parm);
0095         end
0096         
0097         % common ref & down sampling
0098         Parm = parm;
0099         Parm.cutfreq     = [];
0100         Parm.filter_type = [];
0101         Parm.online      = [];
0102         %Parm.fsamp       = [];
0103         %Parm.common_flg  = [];
0104         Parm.bias_flg    = []; 
0105         
0106         [data] = vb_filter_trial_data(data,Parm);
0107     else
0108         [data] = vb_filter_trial_data(data,parm);
0109     end
0110 else
0111     [data] = vb_filter_trial_data(data,parm);
0112 end
0113 %vb_ptime(toc)
0114 
0115 if isfield(parm,'trial_mean') && ~isempty(parm.trial_mean) ...
0116     && parm.trial_mean==1,
0117     data = mean(data ,3);
0118 end
0119 
0120 [Nch,T,Ntr] = size(data);
0121 
0122 % load file info and data as structure
0123 data_info = load(datafile);
0124 
0125 measurement = vb_load_device(datafile);
0126 Measurement = upper(measurement);
0127 
0128 switch    Measurement
0129 case    'MEG'
0130     if ~isfield(data_info, 'bexp') || isempty(data_info.bexp)
0131         data_info.bexp = data;
0132     else
0133         if T~=info.Nsample
0134             data_info.bexp(:,T+1:end,:) = [];
0135         end
0136         if Ntr~=info.Nrepeat
0137             data_info.bexp(:,:,Ntr+1:end) = [];
0138         end
0139         %data_info.bexp = zeros(info.Nchannel, T, Ntr);
0140         data_info.bexp(ch_id,:,:) = data;
0141     end
0142     
0143     % Trial
0144     trial = [];
0145     if isfield(parm,'trial') && ~isempty(parm.trial)
0146         trial = parm.trial;
0147     end
0148     if isfield(parm,'trial_mean') && ~isempty(parm.trial_mean) ...
0149         && parm.trial_mean==1,
0150         trial = 1;
0151     end
0152     if ~isempty(trial)
0153         data_info.MEGinfo.Nrepeat = length(trial);
0154         if isfield(data_info.MEGinfo ,'Trial')
0155             data_info.MEGinfo.Trial = data_info.MEGinfo.Trial(trial);
0156         end
0157         if isfield(data_info.MEGinfo ,'ActiveTrial')
0158             data_info.MEGinfo.ActiveTrial = ones(1,length(trial));
0159         end
0160     end
0161     
0162     % SampleFrequency
0163     if T ~= info.Nsample
0164         data_info.MEGinfo.Nsample = T;
0165         data_info.MEGinfo.SampleFreq = parm.fsamp;
0166         data_info.MEGinfo.Pretrigger = ...
0167                 round(info.Pretrigger*parm.fsamp/info.SampleFreq);
0168     end
0169 
0170 case    'EEG'
0171     if ~isfield(data_info, 'eeg_data') || isempty(data_info.eeg_data)
0172         data_info.eeg_data = data;
0173     else
0174         if T~=info.Nsample
0175             data_info.eeg_data(:,T+1:end,:) = [];
0176         end
0177         if Ntr~=info.Nrepeat
0178             data_info.eeg_data(:,:,Ntr+1:end) = [];
0179         end
0180         %data_info.data = zeros(info.Nchannel, T, Ntr);
0181         data_info.eeg_data(ch_id,:,:) = data;
0182     end
0183     
0184     % Trial
0185     trial = [];
0186     if isfield(parm,'trial') && ~isempty(parm.trial)
0187         trial = parm.trial;
0188     end
0189     if isfield(parm,'trial_mean') && ~isempty(parm.trial_mean) ...
0190         && parm.trial_mean==1,
0191         trial = 1;
0192     end
0193     if ~isempty(trial)
0194         data_info.EEGinfo.Nrepeat = length(trial);
0195         if isfield(data_info.EEGinfo ,'Trial')
0196             data_info.EEGinfo.Trial = data_info.EEGinfo.Trial(trial);
0197         end
0198         if isfield(data_info.EEGinfo ,'ActiveTrial')
0199             data_info.EEGinfo.ActiveTrial = ones(1,length(trial));
0200         end
0201     end
0202     
0203     % SampleFrequency
0204     if T ~= info.Nsample
0205         data_info.EEGinfo.Nsample = T;
0206         data_info.EEGinfo.SampleFrequency = parm.fsamp;
0207         data_info.EEGinfo.Pretrigger = ...
0208                 round(info.Pretrigger*parm.fsamp/info.SampleFreq);
0209     end
0210 
0211 end
0212 
0213 %tic
0214 vb_save_struct(filename, data_info);
0215 %vb_ptime(toc)
0216 
0217 return
0218 
0219 
0220 %    SampleFreq: 500
0221 %    Pretrigger: 0
0222 %       Nsample: 251
0223 
0224 %       Nrepeat: 30
0225 %       Trial: [30x1 struct]
0226 % ActiveTrial: [1x30 double]
0227 
0228 %     Nchannel: 400
0229 %      MEGch_id: [400x1 double]
0230 %    MEGch_name: {400x1 cell}
0231 % ActiveChannel: [400x1 double]
0232 %   ChannelInfo: [1x1 struct]
0233 
0234 % load sensor position
0235 % [pick, Qpick, sensor_weight] = vb_load_sensor(datafile);
0236 % MEG
0237 %    data_info.pick  = pick;
0238 %    data_info.Qpick = Qpick;
0239 %    data_info.MEGinfo.sensor_weight = sensor_weight;
0240 % EEG
0241 %    data_info.EEGinfo.Coord = pick;
0242 
0243 clear data_info
0244 
0245 [data2, channel_info2] = vb_load_meg_data(datafile);
0246 
0247 err = max(abs(data(:)-data2(:)))/mean(data(:))
0248 err = sum(abs(ch_id - channel_info2.ID))
0249 return
0250 
0251 [Nch,T,N]=size(data);
0252 data = reshape(data, [Nch T*N]);
0253 plot(data')

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005