segmentation and computational order in ft_specest_irasa

[ruiliu-psy]

IRASA (Wen&Liu, 2016) is a method developed to separate the rhythmic and fractal components in MEG/EEG signals. To keep consistency of frequency resolution across all data segments, the current implementation of this method in FieldTrip (see example script and ft_specest_irasa.m by @StolkArjen) does the segmentation (to generate overlapped segments for each experimental trial) outside the IRASA procedure, and estimates the fractal component for each of these segments independently, then averages across all segments. However, mathematically, such processes reversed the computation order as in the script provided by the original authors; conceptually, the segments within and between experimental trials should not be traded equally (as in the Welch's method).

The following schemes explained the different procedures.

The computation order matters, also, because the computations of geometric mean (and median) and arithmetic mean do not commute. With a simulated dataset (pink noise + 10Hz oscillation), the following plot illustrates the difference of results by reversing the computational order. The blues lines refer to the mixed spectrum, and the red lines refer to the fractal component estimated by the original script, dash lines was generated with the same dataset and script except for reversing the order of computation.

With the same simulated dataset, the next plot shows the different results returned by FT (dashed lines) and the original script (solid lines).

So indeed FT offered a better frequency specificity and that's something we want to keep. Meanwhile, we would like to adjust the FT implementation so as to be closer to the original algorithm. The difficulty is how to combine ft_freqanalysis and the IRASA structure.

[StolkArjen]

Thanks for this clear report, Rui. You raise two points, 1) one related to taking the arithmetic mean when averaging across sub-segments, and 2) another related to the moment those sub-segments are being created. This second point is where the fieldtrip implementation primarily diverges from Wen & Liu's (the first point can be relatively easily addressed as it occurs outside ft_specest_irasa):

Wen & Liu both create and average across sub-segments inside the IRASA procedure (these sub-segments are thought to benefit the integrity of the IRASA output for reasons I won't go in here). In coordination with the original authors, the sub-segmentation step was moved outside ft_specest_irasa (i.e. using ft_redefinetrial in a prior step, cf. example script) for consistency with ft_specest_mtmfft. Put differently, if the data were sub-segmented inside ft_specest_irasa (typically 90% of the original window length), the frequency resolution of the fractal component computed by ft_specest_irasa would deviate from the frequency resolution of the original power spectrum computed by ft_specest_mtmfft. An effect of this altered frequency resolution (defined as fr = 1 / T, where T is window length) can be observed in your second plot: Wen & Liu's sub-segmentation ends up using smaller time windows than ft_specest_irasa, resulting in a coarser frequency resolution. To really see if the outputs differed, one would need to compare the two implementations under equal circumstances, including identical frequency resolutions.

Having said that, I agree it would be preferred to have the segmentation take place within rather than outside ft_specest_irasa, if only for reasons of parity with Wen & Liu. Perhaps @robertoostenveld can chime in here. I see two possibilities for achieving parity with Wen & Liu while maintaining consistency with ft_specest_mtmfft.

Option A (non-preferred): a pipeline that gives differently segmented input to ft_spectest_mtmfft (90% window length) and to ft_spectest_irasa (100%, internally sub-segmented into 90% chunks), such that the outputs have identical frequency resolution. This option is non-preferred because the consistency of the frequency resolutions would be hard to enforce on the input/output side.

Option B (preferred if valid): apply zero-padding to the sub-segments such that those sub-segments retain their original length and frequency resolution, and, thus, consistency with ft_specest_mtmfft. That way, the sub-segmentation can be moved inside ft_specest_irasa, as well as averaging over those sub-segments. However, this would only work if zero-padding (10% of the time window) didn't alter the spectral estimation.

[ruiliu-psy]

Thank Arjen for the thorough explanation!

Also thank Robert for his sharp input that we need to start with our research question. So before go into those details of the implementations, I would have to talk to my co-authors about their 'forest'. In case, we indeed need IRASA, there would be a bit more validation to do.

A bit more comments on the algorithm, for a reference, in case one is using it.

1. The discrepancy is not on when the sub-segments are created (both of the two implementations segmented the data in the very first step.), but on when the arithmetic means are computed, before or after geometric mean, because the computational order matters.
2. I agree with Arjen that it is essential to keep consistency of frequency resolution in fractal and original spectrum.
3. Not only for reasons of parity with Wen & Liu, I was also wondering why they did it like this. So it turns out that they tried on different datasets, and when the segments are noisy and/or short, they recommend to do the noise filtering (arithmetic mean) as early as possible.
4. It seems reversing the computational order would cause smaller estimates at higher frequency, so that the slope of fractal activity became larger (see the green solid and dashed lines, solid lines refer to simulated pink noise, dashed lines are estimates given by FT).
   

[schoffelen]

Is there any progress towards a satisfactory solution to this issue?

[ruiliu-psy]

Is there any progress towards a satisfactory solution to this issue?

Hi JM, sorry for taking so long. As Robert suggested, we spent some time on completing our hypothesis, leading to the motivation of adopting IRASA (or not). I now have a clearer view on this issue, will talk with my supervisor tomorrow and maybe also on the MEG hackathon next time. I'm also going to do the validations on possible interaction between the algorithm and strength of activity (seeing the increasing offsets between the green solid and dashed lines with frequencies, in the last plot), and on possible interaction between fractal and periodic components (the last plot on tested pure fractal activity).

[ruiliu-psy]

@schoffelen @StolkArjen

I've made a PR #1602 including a testing scripting that generates the following figures.
New features were summarized in the PR.

Red lines refers to the new implementation, addressed the offset issue at higher frequency (i.e. the solid and dashed red lines getting closer to each other, so better represented the 1/f slope). Might be clearer to see the offset in the earlier plot for multiple trials.