<div dir="ltr">
<div>Hi,</div><div><br></div><div>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). <br></div><div><br></div><div>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.</div><div><br></div><div>Good luck,</div><div><br></div><div>Matt G<br></div>
</div><br><div class="gmail_quote"><div dir="ltr">On Mon, Nov 19, 2018 at 2:10 PM Andreas Widmann <<a href="mailto:widmann@uni-leipzig.de">widmann@uni-leipzig.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Anastasios,<br>
<br>
> 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?<br>
<br>
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:<br>
<br>
* 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: <a href="http://home.uni-leipzig.de/biocog/eprints/widmann_a2015jneuroscimeth250_34.pdf" rel="noreferrer" target="_blank">http://home.uni-leipzig.de/biocog/eprints/widmann_a2015jneuroscimeth250_34.pdf</a> :) 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).<br>
* 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).<br>
* You do filter the 0-0.5 and 80-80.5 bands twice.<br>
<br>
Additionally, at signal edges or discontinuities (boundaries) there might be minor non-linearities which, however, cannot easily be resolved.<br>
<br>
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.<br>
<br>
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.<br>
<br>
Hope this helps! Best,<br>
Andreas<br>
<br>
EEG = eeg_emptyset;<br>
EEG.srate = 500;<br>
EEG.xmin = 0;<br>
EEG.nbchan = 1;<br>
EEG.trials = 1;<br>
EEG.pnts = 3301;<br>
EEG.xmax = 9.998;<br>
EEG.data = zeros( EEG.nbchan, EEG.pnts );<br>
EEG.data( 1, 1651 ) = 1;<br>
<br>
% Firstly, I band-pass filtered my non-epoched EEG data from 0.5 to 80Hz<br>
EEG = pop_eegfiltnew( EEG, 'locutoff', 0.5, 'hicutoff', 80, 'filtorder', 3300 );<br>
<br>
b = {};<br>
% then I filtered this data version 5 several times (0.5-4 , 4-8, 8-13, 13-30 , 30-80)<br>
[ TMPEEG(1), com, b{ 1 } ] = pop_eegfiltnew( EEG, 'locutoff', 0.5, 'hicutoff', 4, 'filtorder', 3300 );<br>
[ TMPEEG(2), com, b{ 2 } ] = pop_eegfiltnew( EEG, 'locutoff', 4, 'hicutoff', 8, 'filtorder', 826 );<br>
[ TMPEEG(3), com, b{ 3 } ] = pop_eegfiltnew( EEG, 'locutoff', 8, 'hicutoff', 13, 'filtorder', 826 );<br>
[ TMPEEG(4), com, b{ 4 } ] = pop_eegfiltnew( EEG, 'locutoff', 13, 'hicutoff', 30, 'filtorder', 508 );<br>
[ TMPEEG(5), com, b{ 5 } ] = pop_eegfiltnew( EEG, 'locutoff', 30, 'hicutoff', 80, 'filtorder', 220 );<br>
<br>
b_eff = zeros( 8192, 1 );<br>
data = zeros( EEG.nbchan, EEG.pnts );<br>
figure<br>
subplot( 2, 2, 1 )<br>
for iBand = 1:5<br>
[ h, f ] = freqz( b{ iBand }, 1, 8192, EEG.srate ); % Get frequency response<br>
plot( f, abs( h ) )<br>
hold all<br>
b_eff = b_eff + abs( h ); % Sum over bands in frequency domain<br>
data = data + TMPEEG( iBand ).data; % Sum over bands in time domain<br>
end<br>
<br>
plot( f, b_eff )<br>
xlim( [ 0 100 ] ), ylim( [ 0 2 ] ), xlabel( 'Frequency (Hz)' ), ylabel( 'Magnitude' )<br>
title( 'Frequency response' )<br>
legend( '0.5-4', '4-8', '8-13', '13-30', '30-80', 'sum' )<br>
<br>
% 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.<br>
subplot( 2, 2, 2 )<br>
plot( ( 0:EEG.pnts - 1 ) / EEG.srate, [ EEG.data' data' ] )<br>
xlim( [ 2.5 4 ] ), xlabel( 'Time (s)' ), ylabel( 'Amplitude' )<br>
title( 'Impulse response' )<br>
legend( 'data', 'sum of band-filtered data' )<br>
<br>
b = {};<br>
% Order = 3300 at fs = 500 gives a transition band width of 0.5 Hz<br>
[ TMPEEG(1), com, b{ 1 } ] = pop_eegfiltnew( EEG, 'hicutoff', 3.75, 'filtorder', 3300 ); % Data are already 0.5 Hz high-pass filtered<br>
[ TMPEEG(2), com, b{ 2 } ] = pop_eegfiltnew( EEG, 'locutoff', 4.25, 'hicutoff', 7.75, 'filtorder', 3300 );<br>
[ TMPEEG(3), com, b{ 3 } ] = pop_eegfiltnew( EEG, 'locutoff', 8.25, 'hicutoff', 12.75, 'filtorder', 3300 );<br>
[ TMPEEG(4), com, b{ 4 } ] = pop_eegfiltnew( EEG, 'locutoff', 13.25, 'hicutoff', 29.75, 'filtorder', 3300 );<br>
[ TMPEEG(5), com, b{ 5 } ] = pop_eegfiltnew( EEG, 'locutoff', 30.25, 'filtorder', 3300 ); % Data are already 80 Hz low-pass filtered<br>
<br>
b_eff = zeros( 8192, 1 );<br>
data = zeros( EEG.nbchan, EEG.pnts );<br>
subplot( 2, 2, 3 )<br>
for iBand = 1:5<br>
[ h, f ] = freqz( b{ iBand }, 1, 8192, EEG.srate ); % Get frequency response<br>
plot( f, abs( h ) )<br>
hold all<br>
b_eff = b_eff + abs( h ); % Sum over bands in frequency domain<br>
data = data + TMPEEG( iBand ).data; % Sum over bands in time domain<br>
end<br>
<br>
plot( f, b_eff )<br>
xlim( [ 0 100 ] ), ylim( [ 0 2 ] ), xlabel( 'Frequency (Hz)' ), ylabel( 'Magnitude' )<br>
title( 'Frequency response' )<br>
legend( '4 lowpass', '4-8', '8-13', '13-30', '30 highpass', 'sum' )<br>
<br>
subplot( 2, 2, 4 )<br>
plot( ( 0:EEG.pnts - 1 ) / EEG.srate, [ EEG.data' data' ] )<br>
xlim( [ 2.5 4 ] ), xlabel( 'Time (s)' ), ylabel( 'Amplitude' )<br>
title( 'Impulse response' )<br>
legend( 'data', 'sum of band-filtered data' )<br>
<br>
max(abs(EEG.data-data))<br>
<br>
_______________________________________________<br>
Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" rel="noreferrer" target="_blank">http://sccn.ucsd.edu/eeglab/eeglabmail.html</a><br>
To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu" target="_blank">eeglablist-unsubscribe@sccn.ucsd.edu</a><br>
For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu" target="_blank">eeglablist-request@sccn.ucsd.edu</a></blockquote></div>