eegfilt

 eegfilt() - (high|low|band)pass filter EEG data using two-way least-squares 
              FIR filtering. Multiple data channels and epochs supported.
                    (Requires the Signal Processing Toolbox)
 Usage:
        >> [smoothdata] = eegfilt(data,srate,locutoff,hicutoff);
 or
   [smoothdata,filtwts] = eegfilt(data,srate,locutoff,hicutoff, ...
                                             epochframes,filtorder);
 Inputs:
    data       = (channels,frames*epochs) data to filter
    srate      = data sampling rate (Hz)
    locutoff   = low edge frequency in pass band (Hz)  {0 -> lowpass}
    hicutoff   = high edge frequency in pass band (Hz) {0 -> highpass}

    epochframes= frames per epoch (filter each epoch separately {def/0: all}
    filtorder  = length of the filter in points {default 3*fix(srate/locutoff)}

 Outputs:
    smoothdata = smoothed data
    filtwts    = filter coefficients [smoothdata <- filtfilt(filtwts,1,data)]

0001 % eegfilt() - (high|low|band)pass filter EEG data using two-way least-squares
0002 %              FIR filtering. Multiple data channels and epochs supported.
0003 %                    (Requires the Signal Processing Toolbox)
0004 % Usage:
0005 %        >> [smoothdata] = eegfilt(data,srate,locutoff,hicutoff);
0006 % or
0007 %   [smoothdata,filtwts] = eegfilt(data,srate,locutoff,hicutoff, ...
0008 %                                             epochframes,filtorder);
0009 % Inputs:
0010 %    data       = (channels,frames*epochs) data to filter
0011 %    srate      = data sampling rate (Hz)
0012 %    locutoff   = low edge frequency in pass band (Hz)  {0 -> lowpass}
0013 %    hicutoff   = high edge frequency in pass band (Hz) {0 -> highpass}
0014 %
0015 %    epochframes= frames per epoch (filter each epoch separately {def/0: all}
0016 %    filtorder  = length of the filter in points {default 3*fix(srate/locutoff)}
0017 %
0018 % Outputs:
0019 %    smoothdata = smoothed data
0020 %    filtwts    = filter coefficients [smoothdata <- filtfilt(filtwts,1,data)]
0021 
0022 % 4-22-97 Scott Makeig CNL / Salk Institute, La Jolla CA from bandpass.m
0023 % 5-08-97 fixed frequency bound computation -sm
0024 % 10-22-97 added MINFREQ tests -sm
0025 % 12-05-00 added error() calls -sm
0026 
0027 function [smoothdata,filtwts] = eegfilt(data,srate,locutoff,hicutoff,epochframes,filtorder)
0028 
0029 if nargin<4
0030     fprintf('');
0031     help eegfilt
0032     return
0033 end
0034 
0035 %if ~exist('firls')
0036 %   error('*** eegfilt() requires the signal processing toolbox. ***');
0037 %end
0038 
0039 [chans frames] = size(data);
0040 if chans > 1 & frames == 1,
0041     help eegfilt
0042     error('input data should be a row vector.');
0043 end
0044 nyq            = srate*0.5;  % Nyquist frequency
0045 MINFREQ = 0.1/nyq;
0046 
0047 minfac         = 3;    % this many (lo)cutoff-freq cycles in filter
0048 min_filtorder  = 15;   % minimum filter length
0049 trans          = 0.25; % fractional width of transition zones
0050 
0051 if locutoff>0 & hicutoff > 0 & locutoff > hicutoff,
0052     help eegfilt
0053     error('locutoff > hicutoff ???\n');
0054 end
0055 if locutoff < 0 | hicutoff < 0,
0056    help eegfilt
0057    error('locutoff | hicutoff < 0 ???\n');
0058 end
0059 
0060 if locutoff>nyq,
0061     help eegfilt
0062     error('locutoff cannot be > srate/2');
0063 end
0064 
0065 if hicutoff>=nyq
0066    hicutoff = 0; 
0067 end
0068 
0069 if nargin<6
0070    filtorder = 0;
0071 end
0072 
0073 if filtorder==0,
0074    if locutoff>0,
0075      filtorder = minfac*fix(srate/locutoff);
0076    elseif hicutoff>0,
0077      filtorder = minfac*fix(srate/hicutoff);
0078    end
0079      
0080    if filtorder < min_filtorder
0081         filtorder == min_filtorder;
0082     end
0083 end
0084 
0085 if nargin<5
0086        epochframes = 0;
0087 end
0088 if epochframes ==0,
0089     epochframes = frames;    % default
0090 end
0091 epochs = fix(frames/epochframes);
0092 if epochs*epochframes ~= frames,
0093     error('epochframes does not divide frames.\n');
0094 end
0095 
0096 if filtorder*3 > epochframes,   % Matlab filtfilt() restriction
0097     fprintf('filter order is %d. ',filtorder);
0098     error('epochframes must be 3 times the filtorder.');
0099 end
0100 
0101 if locutoff > 0 & hicutoff > 0,    % bandpass filter
0102 %  fprintf('eegfilt() - performing %d-point bandpass filtering.\n',filtorder);
0103   f=[MINFREQ (1-trans)*locutoff/nyq locutoff/nyq hicutoff/nyq (1+trans)*hicutoff/nyq 1]; 
0104   m=[0       0                      1            1            0                      0]; 
0105 elseif locutoff > 0                % highpass filter
0106  if locutoff/nyq < MINFREQ
0107     help eegfilt
0108 %    fprintf('eegfilt() - highpass cutoff freq must be > %g Hz\n\n',MINFREQ*nyq);
0109     return
0110  end
0111 % fprintf('eegfilter() - performing %d-point highpass filtering.\n',filtorder);
0112  f=[MINFREQ locutoff*(1-trans)/nyq locutoff/nyq 1]; 
0113  m=[   0             0                   1      1];
0114 elseif hicutoff > 0                %  lowpass filter
0115  if hicutoff/nyq < MINFREQ
0116     help eegfilt
0117 %    fprintf('eegfilt() - lowpass cutoff freq must be > %g Hz\n\n',MINFREQ*nyq);
0118     return
0119  end
0120 % fprintf('eegfilt() - performing %d-point lowpass filtering.\n',filtorder);
0121  f=[MINFREQ hicutoff/nyq hicutoff*(1+trans)/nyq 1]; 
0122  m=[     1           1              0                 0];
0123 else
0124  help eegfilt
0125  return
0126 end
0127 
0128 filtwts = firls(filtorder,f,m); % get FIR filter coefficients
0129 
0130 smoothdata = zeros(chans,frames);
0131 for e = 1:epochs                % filter each epoch, channel
0132     for c=1:chans
0133       smoothdata(c,(e-1)*epochframes+1:e*epochframes) ...
0134         = filtfilt(filtwts,1,data(c,(e-1)*epochframes+1:e*epochframes));
0135     end
0136 end
0137