filtering and down sampling for trial data
[data] = vb_filter_trial_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)
= [t1 t2]: Bias correction by time window [t1 t2]
.cutfreq : Cutoff frequency [Hz]
= [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
filter_type = 'low' for length(cutfreq) = 1
filter_type = 'band' for length(cutfreq) = 2
.online : filter order of IIR bandpass filter
= 0: FIR filter (eegfilt)
> 0: Butterworth filter (filtfilt)
< 0: Butterworth filter (online)
2009-8-10 Masa-aki Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [data] = vb_filter_trial_data(data,parm) 0002 % filtering and down sampling for trial data 0003 % [data] = vb_filter_trial_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) 0012 % = [t1 t2]: Bias correction by time window [t1 t2] 0013 % .cutfreq : Cutoff frequency [Hz] 0014 % = [f1] for 'low','high' 0015 % = [f1 f2] for 'band','stop' 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 % --- If these fields are not given or empty, 0020 % corresponding process are not applied 0021 % --- Type of filter specification 0022 % .filter_type : 'low','high','band','stop' 0023 % if this field is not given or empty 0024 % filter_type = 'low' for length(cutfreq) = 1 0025 % filter_type = 'band' for length(cutfreq) = 2 0026 % .online : filter order of IIR bandpass filter 0027 % = 0: FIR filter (eegfilt) 0028 % > 0: Butterworth filter (filtfilt) 0029 % < 0: Butterworth filter (online) 0030 % 2009-8-10 Masa-aki Sato 0031 % 0032 % Copyright (C) 2011, ATR All Rights Reserved. 0033 % License : New BSD License(see VBMEG_LICENSE.txt) 0034 0035 [Nch ,T, Ntr] = size(data); 0036 0037 fsamp = parm.freq; 0038 0039 % Bias correction 0040 if isfield(parm,'bias_flg') && ~isempty(parm.bias_flg) 0041 if length(parm.bias_flg) == 2 , 0042 % Bias correction by time window [t1 t2] 0043 %fprintf('Bias correction by first %d sample\n',parm.bias_flg) 0044 bias = mean(data(:,parm.bias_flg(1):parm.bias_flg(2),:) , 2); 0045 data = data - repmat(bias, [1 T 1]); 0046 elseif parm.bias_flg > 0 0047 % Bias correction by all data 0048 %fprintf('Bias correction mode = %d\n',parm.bias_flg) 0049 for m=1:Ntr 0050 data(:,:,m) = vb_channel_bias_remove(data(:,:,m),parm.bias_flg); 0051 end 0052 end 0053 end 0054 0055 % filtering 0056 if isfield(parm,'cutfreq') && ~isempty(parm.cutfreq), 0057 % cutoff frequency 0058 fcut = parm.cutfreq; 0059 Nfreq = length(fcut); 0060 0061 if Nfreq > 2, error('Cutoff frequency is wrong'); end; 0062 0063 if isfield(parm,'filter_type') && ~isempty(parm.filter_type), 0064 filter_type = parm.filter_type; 0065 else 0066 if Nfreq==1, 0067 filter_type = 'low'; 0068 elseif Nfreq==2 0069 filter_type = 'band'; 0070 end 0071 end 0072 0073 % Lowpass filter 0074 if isfield(parm,'online') 0075 % fprintf('Online Lowpass filter %6.1f [Hz]\n', fcut) 0076 % Online Butterworth filter 0077 Norder = parm.online; % filter order 0078 if Norder > 0 0079 %fprintf('butter filtfilt lowpass filter %6.1f [Hz]\n',fcut) 0080 [B,A] = butter(Norder, fcut/(fsamp/2) ,filter_type); 0081 0082 for m=1:Ntr 0083 for n=1:Nch 0084 data(n,:,m) = filtfilt(B,A, data(n,:,m) ); 0085 end 0086 end 0087 elseif Norder < 0 0088 %fprintf('Online butter lowpass filter %6.1f [Hz]\n',fcut) 0089 [A, B, C, D, Z0] = ... 0090 vb_online_butter_init(fcut,fsamp,abs(Norder),Nch,filter_type); 0091 if Nch > 10 0092 for m=1:Ntr 0093 Z = Z0; 0094 [data(:,:,m) ,Z] = ... 0095 vb_online_filter_loop(A, B, C, D, data(:,:,m), Z); 0096 end 0097 else 0098 for m=1:Ntr 0099 Z = Z0; 0100 [data(:,:,m) ,Z] = ... 0101 vb_online_filter_loop_mex(A, B, C, D, data(:,:,m), Z); 0102 end 0103 end 0104 end 0105 end 0106 if ~isfield(parm,'online') || parm.online==0 0107 % FIR Lowpass filter 0108 % fprintf('FIR Lowpass filter %6.1f [Hz]\n', fcut) 0109 switch filter_type 0110 case 'low' 0111 data = eegfilt(data, fsamp, 0, fcut ); 0112 case 'high' 0113 data = eegfilt(data, fsamp, fcut, 0 ); 0114 case 'band' 0115 data = eegfilt(data, fsamp, fcut(1), fcut(2) ); 0116 case 'stop' 0117 Norder = fix(3*(fsamp/(fcut(2)-fcut(1)))); 0118 [B, A] = fir1(Norder, fcut/(fsamp/2), 'stop'); 0119 for m=1:Ntr 0120 for n=1:Nch 0121 data(n,:,m) = filtfilt(B,A, data(n,:,m) ); 0122 end 0123 end 0124 end 0125 end 0126 end 0127 0128 % Down sampling 0129 if isfield(parm,'fsamp') && ~isempty(parm.fsamp) && fsamp ~=parm.fsamp 0130 fsnew = parm.fsamp; 0131 %fprintf('Down sampling to %6.1f [Hz]\n',fsnew) 0132 0133 data = vb_convert_freq(data, fsamp, fsnew); 0134 fsamp = fsnew; 0135 end 0136 0137 % Common reference 0138 if isfield(parm,'common_flg') && ~isempty(parm.common_flg) && parm.common_flg==1 0139 % fprintf('Common reference\n') 0140 bias = mean(data,1); 0141 data = data - repmat(bias, [Nch 1 1]); 0142 end