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

[Matt Gerhold]

 Hi,

There are many ways to skin a cat. I think Andreas' solution fits the
questions. Personally, I would use a time-frequency representation and then
either average (band-average) or sum (total band power) the bands of
interest. The time-frequency distribution will converse the total power
from the time-domain but distributed across frequency (see Parceval's
theorem).

Depending on how much data you have and whether it's single-trial or
multiple realisations, the spectral power estimate will have some
estimation bias and variability. If you have multiple realisations to
average, the degree to which they are statistically independent will
decrease the variability of the estimate.

Good luck,

Matt G

On Mon, Nov 19, 2018 at 2:10 PM Andreas Widmann <widmann at uni-leipzig.de>
wrote:

> 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))
>
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20181119/80d46587/attachment.html>