Bias correction, Filtering, Down sampling and Common reference
[data] = vb_filter_data(data,parm)
--- Input
data : input signal [Nchannel x Tsample x Ntrial]
--- preprocess parameter
parm
.freq : Sampling frequency [Hz]
--- optional fields
.bias_flg : Bias correction flag
= 0/1/2 : OFF/Bias/Linear) : Bias correction by all data
= [t1 t2]: Bias correction by time window [t1 t2] (msec)
time is specified by [msec] from the begining of data
.cutfreq : Cutoff frequency [Hz] for filtering
= [f1] for 'low','high'
= [f1 f2] for 'band','stop'
.fsamp : Down sampling frequency [Hz]
lowpass or bandpass below Nyquist frequency is necessarily
.common_flg : Common reference flag (= 1/0 : ON/OFF)
--- If these fields are not given or empty,
corresponding process are not applied
--- Type of filter specification
.filter_type : 'low','high','band','stop'
if this field is not given or empty, filtering is not done
.impulse_response :
= 1: Finite impluse filter (eegfilt)
= 2: Butterworth filter
--- following fields are Butterworth filter only
.order : filter order of Butterworth filter
.online = 1: online filter with time delay
= 0: 'filtfilt' is used for filter with no time delay
if this is not given, 'filtfilt' is used
2012-9-19 Masa-aki Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [data, fs] = vb_filter_data(data,parm) 0002 % Bias correction, Filtering, Down sampling and Common reference 0003 % [data] = vb_filter_data(data,parm) 0004 % --- Input 0005 % data : input signal [Nchannel x Tsample x Ntrial] 0006 % --- preprocess parameter 0007 % parm 0008 % .freq : Sampling frequency [Hz] 0009 % --- optional fields 0010 % .bias_flg : Bias correction flag 0011 % = 0/1/2 : OFF/Bias/Linear) : Bias correction by all data 0012 % = [t1 t2]: Bias correction by time window [t1 t2] (msec) 0013 % time is specified by [msec] from the begining of data 0014 % .cutfreq : Cutoff frequency [Hz] for filtering 0015 % = [f1] for 'low','high' 0016 % = [f1 f2] for 'band','stop' 0017 % .fsamp : Down sampling frequency [Hz] 0018 % lowpass or bandpass below Nyquist frequency is necessarily 0019 % .common_flg : Common reference flag (= 1/0 : ON/OFF) 0020 % --- If these fields are not given or empty, 0021 % corresponding process are not applied 0022 % --- Type of filter specification 0023 % .filter_type : 'low','high','band','stop' 0024 % if this field is not given or empty, filtering is not done 0025 % .impulse_response : 0026 % = 1: Finite impluse filter (eegfilt) 0027 % = 2: Butterworth filter 0028 % --- following fields are Butterworth filter only 0029 % .order : filter order of Butterworth filter 0030 % .online = 1: online filter with time delay 0031 % = 0: 'filtfilt' is used for filter with no time delay 0032 % if this is not given, 'filtfilt' is used 0033 % 0034 % 2012-9-19 Masa-aki Sato 0035 % 0036 % Copyright (C) 2011, ATR All Rights Reserved. 0037 % License : New BSD License(see VBMEG_LICENSE.txt) 0038 0039 [Nch ,T, Ntr] = size(data); 0040 0041 fs = parm.freq; 0042 0043 % Common reference 0044 if isfield(parm,'common_flg') && ~isempty(parm.common_flg) ... 0045 && parm.common_flg==1 0046 0047 bias = mean(data,1); 0048 data = data - repmat(bias, [Nch 1 1]); 0049 end 0050 0051 % Bias correction 0052 if isfield(parm,'bias_flg') && ~isempty(parm.bias_flg) 0053 if length(parm.bias_flg) == 2 , 0054 % Bias correction by time window [t1 t2] 0055 twin = round(parm.bias_flg * fs * 0.001); 0056 twin = min(max(twin , 1) , T); 0057 bias = mean(data(:,twin(1):twin(2),:) , 2); 0058 data = data - repmat(bias, [1 T 1]); 0059 elseif parm.bias_flg > 0 0060 % Bias correction by all data 0061 for m=1:Ntr 0062 data(:,:,m) = vb_channel_bias_remove(data(:,:,m),parm.bias_flg); 0063 end 0064 end 0065 end 0066 0067 % Filtering 0068 if isfield(parm,'cutfreq') && ~isempty(parm.cutfreq) && ... 0069 isfield(parm,'filter_type') && ~isempty(parm.filter_type), 0070 0071 % filter_type 0072 filter_type = parm.filter_type; 0073 % cutoff frequency 0074 fcut = parm.cutfreq; 0075 Nfreq = length(fcut); 0076 0077 if Nfreq > 2, error('Cutoff frequency is wrong'); end; 0078 0079 if parm.impulse_response==1 0080 % FIR Lowpass filter 0081 switch filter_type 0082 case 'low' 0083 data = eegfilt(data, fs, 0, fcut ); 0084 case 'high' 0085 data = eegfilt(data, fs, fcut, 0 ); 0086 case 'band' 0087 data = eegfilt(data, fs, fcut(1), fcut(2) ); 0088 case 'stop' 0089 Norder = fix(3*(fs/(fcut(2)-fcut(1)))); 0090 [B, A] = fir1(Norder, fcut/(fs/2), 'stop'); 0091 for m=1:Ntr 0092 for n=1:Nch 0093 data(n,:,m) = filtfilt(B,A, data(n,:,m) ); 0094 end 0095 end 0096 end 0097 else 0098 % Butterworth filter 0099 Norder = parm.order; % filter order 0100 switch filter_type 0101 case 'low' 0102 [B,A] = butter(Norder, fcut(1)/(fs/2) ,'low'); 0103 case 'high' 0104 [B,A] = butter(Norder, fcut(1)/(fs/2) ,'high'); 0105 case 'band' 0106 [B,A] = butter(Norder, fcut(1:2)/(fs/2) ); 0107 case 'stop' 0108 [B,A] = butter(Norder, fcut(1:2)/(fs/2) , 'stop'); 0109 end 0110 0111 if isfield(parm,'online') && parm.online==1 0112 for m=1:Ntr 0113 for n=1:Nch 0114 data(n,:,m) = filter(B,A, data(n,:,m) ); 0115 end 0116 end 0117 else 0118 for m=1:Ntr 0119 for n=1:Nch 0120 data(n,:,m) = filtfilt(B,A, data(n,:,m) ); 0121 end 0122 end 0123 end 0124 end 0125 end 0126 0127 % Down sampling 0128 if isfield(parm,'fsamp') && ~isempty(parm.fsamp) && fs ~=parm.fsamp 0129 fsnew = parm.fsamp; 0130 %fprintf('Down sampling to %6.1f [Hz]\n',fsnew) 0131 0132 data = vb_convert_freq(data, fs, fsnew); 0133 fs = fsnew; 0134 end