brainstorm_filt

(Toolbox/brainstorm_filt.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[smoothdata,filtwts] = brainstorm_filt(data,srate,locutoff,hicutoff,epochframes,filtorder)

Help Text

BRAINSTORM_FILT - Zero-phase filter dedicated to multi-channel MEG and EEG data
 function [smoothdata,filtwts] = brainstorm_filt(data,srate,locutoff,hicutoff,epochframes,filtorder)
Zero-phase filter dedicated to multi-channel MEG and EEG data
Necessitates the Matlab Signal Processing Toolbox

Cross-Reference Information

bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m

dataplot_cb                        C:\BrainStorm_2001\Toolbox\dataplot_cb.m
manage_corticalscouts              C:\BrainStorm_2001\Toolbox\manage_corticalscouts.m

Listing of function C:\BrainStorm_2001\Toolbox\brainstorm_filt.m

function [smoothdata,filtwts] = brainstorm_filt(data,srate,locutoff,hicutoff,epochframes,filtorder)
%BRAINSTORM_FILT - Zero-phase filter dedicated to multi-channel MEG and EEG data
% function [smoothdata,filtwts] = brainstorm_filt(data,srate,locutoff,hicutoff,epochframes,filtorder)
%Zero-phase filter dedicated to multi-channel MEG and EEG data
%Necessitates the Matlab Signal Processing Toolbox

%<autobegin> ---------------------- 09-Jul-2004 22:16:45 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%
% At Check-in: $Author: Mosher $  $Revision: 19 $  $Date: 7/09/04 8:42p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 09-Jul-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 09-Jul-2004 22:16:45 -----------------------


% Author: Sylvain Baillet, Ph.D. 
% ----------------------------- Script History ---------------------------------
% JCM 19-May-2004  Comments cleaning
% ----------------------------- Script History ---------------------------------



if nargin<4
    help brainstorm_filt
    return
end

if ~exist('firls')
    bst_message_window('\n This filter requires the Matlab Signal Processing toolbox.\ ***n\n');
    errordlg('This filter requires the Matlab Signal Processing toolbox')
    return
end

[chans frames] = size(data);
if chans > 1 & frames == 1,
    help brainstorm_filt
    errordlg('Input data should be a row vector.\n\n');
    return
end
nyq            = srate*0.5;  % Nyquist frequency
MINFREQ = 0.1/nyq;

minfac         = 3;    % this many (lo)cutoff-freq cycles in filter 
min_filtorder  = 15;   % minimum filter length
trans          = 0.25; % fractional width of transition zones

if locutoff>0 & hicutoff > 0 & locutoff > hicutoff,
    bst_message_window('\n locutoff cannot be larger than hicutoff\n');
    help brainstorm_filt
    return
end
if locutoff < 0 | hicutoff < 0,
    help brainstorm_filt
    return
end

if locutoff>nyq,
    bst_message_window(sprintf('\n Locutoff freq cannot be larger than srate/2.\n.',nyq));
    help brainstorm_filt
    return
end

if hicutoff>=nyq
    hicutoff = 0; 
end

if nargin<6
    filtorder = 0;
end

if filtorder==0,
    if locutoff>0,
        filtorder = minfac*fix(srate/locutoff);
    elseif hicutoff>0,
        filtorder = minfac*fix(srate/hicutoff);
    end
    
    if filtorder < min_filtorder
        filtorder == min_filtorder;
    end
end

if nargin<5
    epochframes = 0;
end
if epochframes ==0,
    epochframes = frames;    % default
end
epochs = fix(frames/epochframes);
if epochs*epochframes ~= frames,
    %fprintf('brainstorm_filt(): epochframes does not divide frames.\n');
    return
end

if filtorder*3 > epochframes,   % Matlab filtfilt() restriction
    errordlg('Number of time samples is too small to apply this filter.\n');
    return
end

if locutoff > 0 & hicutoff > 0,    % bandpass filter
    bst_message_window(sprintf('Filtering...%d-point bandpass filter.\n',filtorder));
    f=[MINFREQ (1-trans)*locutoff/nyq locutoff/nyq hicutoff/nyq (1+trans)*hicutoff/nyq 1]; 
    m=[0       0                      1            1            0                      0]; 
elseif locutoff > 0                % highpass filter
    if locutoff/nyq < MINFREQ
        help brainstorm_filt
        errodlg('Highpass cut-off freq must be > %g Hz\n\n',MINFREQ*nyq);
        return
    end
    bst_message_window(sprintf('Filtering...%d-point highpass filter.\n',filtorder));
    f=[MINFREQ locutoff*(1-trans)/nyq locutoff/nyq 1]; 
    m=[   0             0                   1      1];
elseif hicutoff > 0                %  lowpass filter
    if hicutoff/nyq < MINFREQ
        help brainstorm_filt
        errodlg('Lowpass cut-off freq must be > %g Hz\n\n',MINFREQ*nyq);
        return
    end
    bst_message_window(sprintf('Filtering...%d-point lowpass filter.\n',filtorder));
    f=[MINFREQ hicutoff/nyq hicutoff*(1+trans)/nyq 1]; 
    m=[     1           1              0                 0];
else
    errordlg(sprintf(...
        'High or low cut-off frequency might be larger that Nyquist largest admissible frequency value: %3.2f Hz.',nyq),...
        'Error in specifying cut-off frequencies')
    return
end

filtwts = firls(filtorder,f,m); % get FIR filter coefficients

smoothdata = zeros(chans,frames);
for e = 1:epochs                % filter each epoch, channel 
    for c=1:chans
        smoothdata(c,(e-1)*epochframes+1:e*epochframes) ...
            = filtfilt(filtwts,1,data(c,(e-1)*epochframes+1:e*epochframes));
    end
end