To get instantaneous phases of specific epochs, we can bandpass filter a signal in a frequency band e.g. theta 4-8Hz and Hilbert transform. Can we avoid edge effects then, by simply first bandpass filtering and Hilbert transforming the whole continuous signal, and then epoch afterwards?
And to bandpass filter, can we use Eeglab’s eegfiltnew() instead of fir1()? The former seems more user friendly as you only need to enter your lower and upper bounds e.g. 4 to 8 Hz.
Hi Micky. Yes, you can filter the continuous data before epoching, but then you also need to think about cleaning, trial rejection, etc. And if you want to use a different frequency range, you’d need to start again from the beginning. So it would probably be better to process the broadband data and then filter; just make sure the epochs are long enough to absorb the edge effects.
As for the filter, I believe that eeglab’s eegfiltnew is actually a wrapper for fir1. But I could be wrong about that. Anyway, sure, go for the eeglab version.
I meant to apply the bandpass filter after all the cleaning, but now I’m confused.
During the cleaning, I already apply a High pass filter of 1 Hz and a Low pass of 45 Hz. Then I do the cleaning steps e.g re-referencing, ICA.
Now normally, after all this, I epoch. But my suggestion was thus to bandpass filter + Hilbert transform, and then epoch.
Is it bad to bandpass filter after already having filtered the data during cleaning?
From what I read in other synchrony papers, this is often done. And both Matlab and Python PLV functions tend to accept a cleaned data set and a range of frequencies for which you want PLV, so it seems to be the norm.
Also, eeglab uses firfilt(), which is indeed similar to fir1, so I’ll go ahead with that.
Oh I see. Yeah, that’s fine. Data cleaning is often done on epoched data, because a lot of noise occurs during experiment breaks, etc. But what you describe is also fine.
You can narrowband filter after bandpassing. No problem.
Great, thank you for your knowledge. I’ll try my method then, and if there is time, I’ll do it the usual way too to see how the ISCP values are influenced by those edge effects.
Update for future readers.
I tried 4 different orders to calculate ISPC of 5 seconds activity after 2 types of events (between 2 people’s brains). See the image at the bottom which shows the angles for 2 types of trial in 1 channel.
- Epoch (pop_epoch) > Filter (fir1) > angle(hilbert): this is the common way of doing it and introduces edge effects as can be seen in the image.
- Epoch > Filter (using eeglab’s eegfiltnew) > angle(hilbert)
- Filter (eeglab) > epoch > angle(hilbert)
- Filter (fir1) > angle(hilbert) > epoch: this way avoids edge effects completely by applying the filter and Hilbert transform on the continuous signal and ‘epochs’ the angles afterwards.
Most ISPC values were between 0.15 and 0.28.
From just looking at the absolute differences in 20 channels and 2 event types and both ISPC across trials and across times, I can say the following:
- ISPC values between method 1 and 2 differed a maximum of 0.06. This is a lot considering most ISPCs were around 0.20. I have no idea why this happens as the only difference is the fir1 filter vs eeglab’s filtnew.
- ISPC values between method 1 and 3 differed less, but still a noticeable amount of 0.5 at it’s max.
- ISPC values differed very little between method 1 and 4, which is surprising as this should have shown the largest differences as we are completely avoiding edge effects.
- Method 2 and 3 had the least differences.
The plots also show that method 1 and 4 look very similar and method 2 and 3 do. No idea what the pattern is here…
Obviously no conclusions should be made from this, as this is just from 1 session, but I hope one day someone will compare these orders on a larger scale statistically.
Hey, thanks for doing and reporting that, Mickey. Interesting way to test the different methods orders.
Did you compute ISPC across the entire time window? i.e., including the edge effects? That could explain the differences up to .06. What I recommend with edges is to make sure the epochs are wide enough that you can ignore the edges in the analyses. For example, computing ISPC between .5 and 4.5. Or possibly even a narrower window.
That said, I’m a bit surprised that the phase angle time series are visibly different in the middle of the segments. I suppose that’s the difference between the two filter kernel methods. Filtering is highly parameter-sensitive. I guess that inspecting the spectral gain functions of the two filters would show some differences. It’s always possible to match filters to be more similar.
The ISPCs were indeed across the entire window by differencing all the 20 angles (across trials) or 5000 angles (across time). After that, I just averaged all 20 or 5000 to get an average ISPC-trials and ISPC-time for each channel, and this was done for each event type (because I’m comparing 2 events/stimuli). Thus, in terms of ISPC, both methods end up using all the phase angles, just in different orders, so in the end I think it’s most important to look at the angles themselves.
As for the angles, the filter could indeed cause the difference. Specifically, for the fir1 filter, I set both the filter range (4 to 8 in the test above) and the order (4*EEG.srate/lower_frequency, so 4 cycles). While in EEGLAB’s pop_eegfiltnew, I only set the range, because the filter order is set automatically. The default is either 3.3 or 3*fix(srate/locutoff) (sources: github and sccn.uscd.edu). Thus the cause could be the filter order.
Aside from that, it’s still interesting how method 1 and 2 differed so little. Even thought they used the same filter, the plots do show that the edge effects (at least at the start) dissappeared. I guess then that the Hilbert transform might not introduce edge effects? Or somehow it does, but it’s edge effects go in the opposite way of the filter’s edge effects, ultimately cancelling them out. Again, this is just 1 specific case, so we’ll have to wait for a large scale analysis.
Finally, 1 tip in case anyone wants to easily run the Filter > Hilbert > Epoch order. Instead of manually coding your epoching on the continuous data after filtering and Hilberting, you can replace EEG.data with the angles you extracted as they will be the same shape of channels*timepoints. EEGLab does not care whether the values in EEG.data are amplitudes (in microVolt) or angles (in radians). You can then run the pop_epoch function on your EEG data as you normally would to get epochs. This saved me a lot of time and headache. For fun, you can even run the pop_eegplot to plot all the angles.
Hope this helps. Cheers.
Correction: from some EEGLAB list discussions and the github code, it seems the pop_eegfiltnew function does not use that simple equation afterall, but something a bit more complicated. While most of it goes over my head, I trust that EEGLAB’s function will set the filter order and other parameters better than I could using fir1 since it’s made for EEG data by the experts, thus I will be using that instead of my custom fir1(). I’ll also retest the 3 different orders using different frequencies and epoch lengths and report back later.
That all makes sense. The thing to keep in mind about edges is that they are not guaranteed to be fucked up; they are likely to contain distortions and so the best strategy is to ignore them, just in case. If the data naturally taper to small numerical values near the edges, then the edge effects will be very small.
As for the filters: Selecting an appropriate filter order is really tricky. There are guidelines that help select minimal orders, but there generally isn’t an “optimal” order – particularly for FIR filters, different orders provide different spectral resolutions and different computation times, but are not necessarily right or wrong (within some reasonable range). Anyway, yeah you can using eeglab’s filter function and sleep easy
Thanks again for the really nice discussion and for reporting your tests!
Small follow up question:
I need to normalize the average ISPC for both event types (each type has 20 trials) so I can compare them statistically using a parametric test.
What is better: a Fisher-Z transform or an ‘arcsine of the square root’ of the ISPC?
Do we perform the transform on the ISPC of each trial first and then average these values to get 1 average ISPC for each event type? Or can I apply the transform on the already calculated average ISPC of each event type (which is the average of the 20 trials of each type)?
I generally use Fisher-Z, but any other transform that gets the data closer to normal is fine.
As for the averaging, it’s often best to apply nonlinear transforms as late as possible. So I would first average the raw ISPC values and then apply the transform only during the statistics.
That makes sense, thank you!