[Eeglablist] Using the basic FIR filter to extract the EEG bands

[Andreas Widmann]

Hi Anastasios,

> I am using the basic FIR filtering function (from EEGLAB’s GUI) to extract the time-domain waveforms in each individual band (delta,theta,alpha,beta,gamma). Firstly, I band-pass filtered my non-epoched EEG data from 0.5 to 80Hz and then I filtered this data version 5 several times (0.5-4 , 4-8, 8-13, 13-30 , 30-80). I expected that the sum of these 5 time-domain versions would be equal to the initial waveform (0.5-80Hz), but this was not noticed. Do you have any reasons that this happened?

In principle your reasoning is correct that the sum of the bands should approximately equal the original signal (ignoring the small effects of passband and stopband ripple). There are, however, several issues with your implementation:

* Most importantly the transition bands of your filters do overlap with the passbands of the other filters. Therefore your signal in the transition bands will be actually amplified (by about a factor of 2) and not filtered. I suggest to first make yourself familiar with the concepts of cutoff frequency and passband and transition band (e.g. here: http://home.uni-leipzig.de/biocog/eprints/widmann_a2015jneuroscimeth250_34.pdf :) The problem is partly due to the EEGLAB basic FIR filter defining filters by passband edge rather than cutoff frequency (as in most other filter implementations; the implementation was done this way on explicit request of the EEGLAB developers).
* You use filters of different length/order. You have to manually specify a suitable and constant filter order (with basic FIR filter; note that in general also other approaches combining appropriate high-pass and low-pass filters are possible).
* You do filter the 0-0.5 and 80-80.5 bands twice.

Additionally, at signal edges or discontinuities (boundaries) there might be minor non-linearities which, however, cannot easily be resolved.

Find a demonstration of the issues in MATLAB code below. First, a test signal is defined, here an impulse (if you read the above paper you will learn why). In a first run the test signal is filtered with the filters you described, the frequency response is computed and displayed. Finally, the sums of the frequency responses and the sum of the band-filtered time domain signals are displayed. Notice how the signal is amplified in the overlapping transition and passbands.

In the second run I suggest a set of filters with (a) constant order and transition band (0.5 Hz) and (b) overlapping transition bands but not overlapping passbands and transition bands (i.e., higher pass-band edge of one filter plus transition band width equals the lower pass-band edge of the other filter) and (c) not filtering below 0.5 and above 80 Hz. The sum of the frequency responses does not exceed any of the frequency responses (ignoring expected ripple not exceeding 0.46%) and the difference between the original (0.5-80 Hz filtered) and the sum of the band-filtered signals does not exceed 10^-15.

Hope this helps! Best,
Andreas

EEG = eeg_emptyset;
EEG.srate = 500;
EEG.xmin = 0;
EEG.nbchan = 1;
EEG.trials = 1;
EEG.pnts = 3301;
EEG.xmax = 9.998;
EEG.data = zeros( EEG.nbchan, EEG.pnts );
EEG.data( 1, 1651 ) = 1;

% Firstly, I band-pass filtered my non-epoched EEG data from 0.5 to 80Hz
EEG = pop_eegfiltnew( EEG, 'locutoff',  0.5, 'hicutoff',  80, 'filtorder', 3300 );

b = {};
% then I filtered this data version 5 several times (0.5-4 , 4-8, 8-13, 13-30 , 30-80)
[ TMPEEG(1), com, b{ 1 } ] = pop_eegfiltnew( EEG, 'locutoff',  0.5, 'hicutoff',  4, 'filtorder', 3300 );
[ TMPEEG(2), com, b{ 2 } ] = pop_eegfiltnew( EEG, 'locutoff',  4,   'hicutoff',  8, 'filtorder',  826 );
[ TMPEEG(3), com, b{ 3 } ] = pop_eegfiltnew( EEG, 'locutoff',  8,   'hicutoff', 13, 'filtorder',  826 );
[ TMPEEG(4), com, b{ 4 } ] = pop_eegfiltnew( EEG, 'locutoff', 13,   'hicutoff', 30, 'filtorder',  508 );
[ TMPEEG(5), com, b{ 5 } ] = pop_eegfiltnew( EEG, 'locutoff', 30,   'hicutoff', 80, 'filtorder',  220 );

b_eff = zeros( 8192, 1 );
data = zeros( EEG.nbchan, EEG.pnts );
figure
subplot( 2, 2, 1 )
for iBand = 1:5
    [ h, f ] = freqz( b{ iBand }, 1, 8192, EEG.srate ); % Get frequency response
    plot( f, abs( h ) )
    hold all
    b_eff = b_eff + abs( h ); % Sum over bands in frequency domain
    data = data + TMPEEG( iBand ).data; % Sum over bands in time domain
end

plot( f, b_eff )
xlim( [ 0 100 ] ), ylim( [ 0 2 ] ), xlabel( 'Frequency (Hz)' ), ylabel( 'Magnitude' )
title( 'Frequency response' )
legend( '0.5-4', '4-8', '8-13', '13-30', '30-80', 'sum' )

% I expected that the sum of these 5 time-domain versions would be equal to the initial waveform (0.5-80Hz), but this was not noticed.
subplot( 2, 2, 2 )
plot( ( 0:EEG.pnts - 1 ) / EEG.srate, [ EEG.data' data' ] )
xlim( [ 2.5 4 ] ), xlabel( 'Time (s)' ), ylabel( 'Amplitude' )
title( 'Impulse response' )
legend( 'data', 'sum of band-filtered data' )

b = {};
% Order = 3300 at fs = 500 gives a transition band width of 0.5 Hz
[ TMPEEG(1), com, b{ 1 } ] = pop_eegfiltnew( EEG,                    'hicutoff',  3.75, 'filtorder', 3300 ); % Data are already 0.5 Hz high-pass filtered
[ TMPEEG(2), com, b{ 2 } ] = pop_eegfiltnew( EEG, 'locutoff',  4.25, 'hicutoff',  7.75, 'filtorder', 3300 );
[ TMPEEG(3), com, b{ 3 } ] = pop_eegfiltnew( EEG, 'locutoff',  8.25, 'hicutoff', 12.75, 'filtorder', 3300 );
[ TMPEEG(4), com, b{ 4 } ] = pop_eegfiltnew( EEG, 'locutoff', 13.25, 'hicutoff', 29.75, 'filtorder', 3300 );
[ TMPEEG(5), com, b{ 5 } ] = pop_eegfiltnew( EEG, 'locutoff', 30.25,                    'filtorder', 3300 ); % Data are already 80 Hz low-pass filtered

b_eff = zeros( 8192, 1 );
data = zeros( EEG.nbchan, EEG.pnts );
subplot( 2, 2, 3 )
for iBand = 1:5
    [ h, f ] = freqz( b{ iBand }, 1, 8192, EEG.srate ); % Get frequency response
    plot( f, abs( h ) )
    hold all
    b_eff = b_eff + abs( h ); % Sum over bands in frequency domain
    data = data + TMPEEG( iBand ).data; % Sum over bands in time domain
end

plot( f, b_eff )
xlim( [ 0 100 ] ), ylim( [ 0 2 ] ), xlabel( 'Frequency (Hz)' ), ylabel( 'Magnitude' )
title( 'Frequency response' )
legend( '4 lowpass', '4-8', '8-13', '13-30', '30 highpass', 'sum' )

subplot( 2, 2, 4 )
plot( ( 0:EEG.pnts - 1 ) / EEG.srate, [ EEG.data' data' ] )
xlim( [ 2.5 4 ] ), xlabel( 'Time (s)' ), ylabel( 'Amplitude' )
title( 'Impulse response' )
legend( 'data', 'sum of band-filtered data' )

max(abs(EEG.data-data))