butterfilt() - Filters EEG data with low pass, high pass, band-pass,
or band-stop 3rd order Butterworth filter. Independent
component (IC) activations and features are re-computed after
filtering and baselining if the EEG variable contains ICs.
Usage:
>> EEG=butterfilt_linear(EEG,flt,baseline,n_pad);
Required Inputs:
EEG - EEGLAB "EEG" struct variable
filt - [low_boundary high_boundary] a two element vector indicating the frequency
cut-offs for a 3rd order Butterworth filter that will be applied to each
trial of data. If low_boundary=0, the filter is a low pass filter. If
high_boundary=(sampling rate)/2, then the filter is a high pass filter. If both
boundaries are between 0 and (sampling rate)/2, then the filter is a band-pass filter.
If both boundaries are between 0 and -(sampling rate)/2, then the filter is a band-
stop filter (with boundaries equal to the absolute values of the low_boundary and
high_boundary). Note, using this option requires the signal processing toolbox
function butter.m. You should probably use the 'baseline' option as well since the
mean prestimulus baseline may no longer be 0 after the filter is applied
Optional Inputs:
baseline - [low_boundary high_boundary] a time window (in msec) whose mean amplitude in
each trial will be removed from each trial (e.g., [-100 0]) after filtering.
{default: no further baselining of data}
n_pad - [integer] number of time points to pad each epoch of
data before filtering. For example, if n_pad is 512 and
you have 256 time points of data per epoch, each epoch of
data will be grown to 1280 time points (i.e., 512*2+256) by
adding 512 padding points before and after the 256 time points
of data. The values of the padding points form a line from the
value of the first and last data point (i.e., if you wrapped
the padded waveform around a circle the pre-padding and
post-padding points would form a line). This method of
padding is the way the Kutaslab function filter.c works.
If n_pad is not specified, butterfilt.m will follow the Kutaslab
convention of n_pad=2*(the number of time points per
epoch). If you wish to avoid zero padding, set n_pad to 0.
Example:
1) A 15 Hz low pass filter with no zero padding
>> EEG=butterfilt(EEG,[0 15],[-100 0],0);
2) A .2 Hz high pass filter with default, Kutaslab-like, zero padding.
Note, the sampling rate in this example is 250 Hz.
>> EEG=butterfilt(EEG,[0.2 125],[-100 0]);
3) A .2 to 15 Hz band pass filter with no zero padding
>> EEG=butterfilt(EEG,[0.2 15],[-100 0],0);
Author:
David Groppe
Kutaslab, 12/20090001 function EEG=butterfilt(EEG,flt,baseline,n_pad) 0002 % butterfilt() - Filters EEG data with low pass, high pass, band-pass, 0003 % or band-stop 3rd order Butterworth filter. Independent 0004 % component (IC) activations and features are re-computed after 0005 % filtering and baselining if the EEG variable contains ICs. 0006 % 0007 % Usage: 0008 % >> EEG=butterfilt_linear(EEG,flt,baseline,n_pad); 0009 % 0010 % Required Inputs: 0011 % EEG - EEGLAB "EEG" struct variable 0012 % filt - [low_boundary high_boundary] a two element vector indicating the frequency 0013 % cut-offs for a 3rd order Butterworth filter that will be applied to each 0014 % trial of data. If low_boundary=0, the filter is a low pass filter. If 0015 % high_boundary=(sampling rate)/2, then the filter is a high pass filter. If both 0016 % boundaries are between 0 and (sampling rate)/2, then the filter is a band-pass filter. 0017 % If both boundaries are between 0 and -(sampling rate)/2, then the filter is a band- 0018 % stop filter (with boundaries equal to the absolute values of the low_boundary and 0019 % high_boundary). Note, using this option requires the signal processing toolbox 0020 % function butter.m. You should probably use the 'baseline' option as well since the 0021 % mean prestimulus baseline may no longer be 0 after the filter is applied 0022 % 0023 % Optional Inputs: 0024 % baseline - [low_boundary high_boundary] a time window (in msec) whose mean amplitude in 0025 % each trial will be removed from each trial (e.g., [-100 0]) after filtering. 0026 % {default: no further baselining of data} 0027 % n_pad - [integer] number of time points to pad each epoch of 0028 % data before filtering. For example, if n_pad is 512 and 0029 % you have 256 time points of data per epoch, each epoch of 0030 % data will be grown to 1280 time points (i.e., 512*2+256) by 0031 % adding 512 padding points before and after the 256 time points 0032 % of data. The values of the padding points form a line from the 0033 % value of the first and last data point (i.e., if you wrapped 0034 % the padded waveform around a circle the pre-padding and 0035 % post-padding points would form a line). This method of 0036 % padding is the way the Kutaslab function filter.c works. 0037 % If n_pad is not specified, butterfilt.m will follow the Kutaslab 0038 % convention of n_pad=2*(the number of time points per 0039 % epoch). If you wish to avoid zero padding, set n_pad to 0. 0040 % 0041 % 0042 % Example: 0043 % 1) A 15 Hz low pass filter with no zero padding 0044 % >> EEG=butterfilt(EEG,[0 15],[-100 0],0); 0045 % 0046 % 2) A .2 Hz high pass filter with default, Kutaslab-like, zero padding. 0047 % Note, the sampling rate in this example is 250 Hz. 0048 % >> EEG=butterfilt(EEG,[0.2 125],[-100 0]); 0049 % 0050 % 3) A .2 to 15 Hz band pass filter with no zero padding 0051 % >> EEG=butterfilt(EEG,[0.2 15],[-100 0],0); 0052 % 0053 % Author: 0054 % David Groppe 0055 % Kutaslab, 12/2009 0056 % 0057 0058 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0059 % 0060 % 1/10/2010 - Function now checks for ICA mixing matrix before attempting to 0061 %recompute IC features 0062 % 0063 % 11/5/2010 - n_pad option added 0064 0065 %%%%%%%%%%%%%%%% FUTURE ADDITIONS %%%%%%%%%%%%%%%%% 0066 %-add verblevel? 0067 0068 if nargin<2, 0069 error('butterfilt.m requires at least two inputs.'); 0070 end 0071 if nargin<3, 0072 baseline=[]; 0073 end 0074 if nargin<4 || isempty(n_pad), 0075 n_pad=2*EEG.pnts; 0076 end 0077 0078 if length(flt)~=2, 0079 error('''filt'' parameter argument should be a two element vector.'); 0080 elseif max(flt)>(EEG.srate/2), 0081 error('''filt'' parameters need to be less than or equal to (sampling rate)/2 (i.e., %f).',EEG.srate/2); 0082 elseif (flt(2)==(EEG.srate/2)) && (flt(1)==0), 0083 error('If second element of ''filt'' parameter is EEG.srate/2, then the first element must be greater than 0.'); 0084 elseif abs(flt(2))<=abs(flt(1)), 0085 error('Second element of ''filt'' parameters must be greater than first in absolute value.'); 0086 elseif (flt(1)<0) || (flt(2)<0), 0087 if (flt(1)>=0) || (flt(2)>=0), 0088 error('BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.'); 0089 end 0090 if min(flt)<=(-EEG.srate/2), 0091 error('''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',EEG.srate/2); 0092 end 0093 end 0094 0095 fprintf('\nFiltering data with 3rd order Butterworth filter: '); 0096 if (flt(1)==0), 0097 %lowpass filter the data 0098 [B A]=butter(3,flt(2)*2/EEG.srate,'low'); 0099 fprintf('lowpass at %.2f Hz\n',flt(2)); 0100 elseif (flt(2)==(EEG.srate/2)), 0101 %highpass filter the data 0102 [B A]=butter(3,flt(1)*2/EEG.srate,'high'); 0103 fprintf('highpass at %.2f Hz\n',flt(1)); 0104 elseif (flt(1)<0) 0105 %bandstop filter the data 0106 flt=-flt; 0107 [B A]=butter(3,flt*2/EEG.srate,'stop'); 0108 fprintf('stopband from %.2f to %.2f Hz\n',flt(1),flt(2)); 0109 else 0110 %bandpass filter the data 0111 [B A]=butter(3,flt*2/EEG.srate); 0112 fprintf('bandpass from %.2f to %.2f Hz\n',flt(1),flt(2)); 0113 end 0114 %preallocate memory 0115 total_pad=n_pad*2; 0116 total_tpts=total_pad+EEG.pnts; 0117 padded=zeros(1,total_tpts); 0118 0119 0120 fprintf('Adding %d time points of zeros before and after each epoch of data.\n',n_pad); 0121 fprintf('Padded values are linear interpolations between the first and last value of each epoch of data.\n'); 0122 s=size(EEG.data); 0123 if n_pad, 0124 for chan=1:s(1), 0125 fprintf('Now filtering channel: %s\n',EEG.chanlocs(chan).labels); 0126 for trial=1:s(3), 0127 padded(n_pad+1:n_pad+EEG.pnts)=squeeze(EEG.data(chan,:,trial)); %put real data in the middle 0128 inc=(EEG.data(chan,EEG.pnts,trial)-EEG.data(chan,1,trial))/(total_pad+1); 0129 %add padding before the real data 0130 for a=n_pad:-1:1, 0131 padded(a)=padded(a+1)+inc; 0132 end 0133 %add padding after the real data 0134 for a=n_pad+EEG.pnts+1:total_tpts, 0135 padded(a)=padded(a-1)-inc; 0136 end 0137 %Note, those for loops are actually faster than using the 0138 %MATLAB linspace function. Go figure. 0139 filtered=filtfilt(B,A,padded); 0140 EEG.data(chan,:,trial)=filtered(n_pad+1:n_pad+EEG.pnts); 0141 end 0142 end 0143 else 0144 %no padding 0145 for chan=1:s(1), 0146 for trial=1:s(3), 0147 EEG.data(chan,:,trial)=filtfilt(B,A,squeeze(EEG.data(chan,:,trial))); 0148 end 0149 end 0150 end 0151 0152 %Find out if EEG variable has independent components 0153 if ~isempty(EEG.icawinv) && ~isempty(EEG.icasphere), 0154 ics_present=1; 0155 else 0156 ics_present=0; 0157 end 0158 0159 0160 %% Mean Baseline Each Trial (if requested) %% 0161 if ~isempty(baseline), 0162 %check argument values for errors 0163 if baseline(2)<baseline(1), 0164 error('First element of ''baseline'' argument needs to be less than or equal to second argument.'); 0165 elseif baseline(2)<EEG.times(1), 0166 error('Second element of ''baseline'' argument needs to be greater than or equal to epoch start time %.1f.',EEG.times(1)); 0167 elseif baseline(1)>EEG.times(end), 0168 error('First element of ''baseline'' argument needs to be less than or equal to epoch end time %.1f.',EEG.times(end)); 0169 end 0170 0171 %convert msec into time points 0172 if baseline(1)<EEG.times(1), 0173 strt_pt=1; 0174 else 0175 strt_pt=find_crspnd_pt(baseline(1),EEG.times,1:length(EEG.times)); 0176 strt_pt=ceil(strt_pt); 0177 end 0178 if baseline(end)>EEG.times(end), 0179 end_pt=length(EEG.times); 0180 else 0181 end_pt=find_crspnd_pt(baseline(2),EEG.times,1:length(EEG.times)); 0182 end_pt=floor(end_pt); 0183 end 0184 fprintf('\nRemoving pre-stimulus mean baseline from %.1f to %.1f msec.\n\n',EEG.times(strt_pt),EEG.times(end_pt)); 0185 EEG.data=rmbase(EEG.data,s(2),strt_pt:end_pt); 0186 if ics_present, 0187 %re-compute ICA activations 0188 EEG.icaact=EEG.icaweights*EEG.icasphere*EEG.data; 0189 EEG.icaact=reshape(EEG.icaact,s(1),s(2),s(3)); 0190 end 0191 EEG.data=reshape(EEG.data,s(1),s(2),s(3)); 0192 else 0193 if ics_present, 0194 %re-compute ICA activations 0195 EEG.icaact=EEG.icaweights*EEG.icasphere*reshape(EEG.data,s(1),s(2)*s(3)); 0196 EEG.icaact=reshape(EEG.icaact,s(1),s(2),s(3)); 0197 end 0198 end 0199 EEG.saved='no'; 0200 0201 if ~isempty(EEG.icawinv), 0202 % normalize ICs to unit topography length and recompute IC features 0203 EEG=cmpt_ic_art_ftrs(EEG); 0204 end 0205 0206 fprintf('EEG has been filtered. You must save the EEG variable to preserve changes.\n'); 0207 0208 if ischar(EEG.history), 0209 %if EEG.history is a string, turn it into a cell array 0210 temp=EEG.history; 0211 EEG.history=[]; 0212 EEG.history{1}=temp; 0213 end 0214 if isempty(baseline), 0215 EEG.history{length(EEG.history)+1}=sprintf('EEG=butterfilt(EEG,[%f %f],[],%d);',flt(1),flt(2),n_pad); 0216 else 0217 EEG.history{length(EEG.history)+1}=sprintf('EEG=butterfilt(EEG,[%f %f],[%d %d],%d);',flt(1),flt(2),baseline(1),baseline(2),n_pad); 0218 end 0219 0220 0221 %% %%%%%%%%%%%%%%%%%%%%% function find_crspnd_pt() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0222 % 0223 function y_pt=find_crspnd_pt(targ,vals,outtrials) 0224 %function id=find_crspnd_pt(targ,vals,outtrials) 0225 % 0226 % Inputs: 0227 % targ - desired value of sorting variable 0228 % vals - a vector of observed sorting variables (possibly smoothed) 0229 % outtrials - a vector of y-axis values corresponding to trials in the 0230 % ERPimage (this will just be 1:n_trials if there's no 0231 % smoothing) 0232 % 0233 % Output: 0234 % y_pt - y-axis value (in units of trials) corresponding to "targ". 0235 % If "targ" matches more than one y-axis pt, the median point is 0236 % returned. If "targ" falls between two points, y_pt is linearly 0237 % interpolated. 0238 % 0239 % Note: targ and vals should be in the same units (e.g., milliseconds) 0240 0241 %find closest point above 0242 abv=find(vals>=targ); 0243 if isempty(abv), 0244 %point lies outside of vals range, can't interpolate 0245 y_pt=[]; 0246 return 0247 end 0248 abv=abv(1); 0249 0250 %find closest point below 0251 blw=find(vals<=targ); 0252 if isempty(blw), 0253 %point lies outside of vals range, can't interpolate 0254 y_pt=[]; 0255 return 0256 end 0257 blw=blw(end); 0258 0259 if (vals(abv)==vals(blw)), 0260 %exact match 0261 ids=find(vals==targ); 0262 y_pt=median(outtrials(ids)); 0263 else 0264 %interpolate point 0265 0266 %lst squares inear regression 0267 B=regress([outtrials(abv) outtrials(blw)]',[ones(2,1) [vals(abv) vals(blw)]']); 0268 0269 %predict outtrial point from target value 0270 y_pt=[1 targ]*B; 0271 0272 end 0273 0274 0275