apply_hilbert for Epochs class

[LegrandNico]

Wouldn't it be convenient to have an apply_hilbert method for Epochs class?

Artifact rejection modules like Autoreject directly output epoched data, which makes it difficult to have both clean data and the Hilbert transform at the end of the pipeline.

I was wondering if accessing the epochs as raw data to directly apply Hilbert and then epoch again with cropped interval could be a solution? But I don't know how it could affect the results...

[agramfort]

provided your Epochs are not too short yes it should work and could be a nice addition to MNE.

[LegrandNico]

Ok I will try this solution. If there is e.g + 200ms on each side would it be sufficient? Or maybe there is a more formal way to define the crop needed based on the Hilbert calculation?

[agramfort]

you want the phase or amplitude? if you need phase you need to filter the signals first so they are quite narrow band.

200ms depends on what frequency you are interested in. For high freqs I would say higher then 20Hz it should start to work.

HTH

[palday]

I'm happy to contribute this, but am currently travelling and won't get around to it for about a week.

In the meantime, here's some code I had been using to do arbitrary transformations on epochs with an example for the Hilbert transformation. It's a very naive approach that accesses private members, so consider yourself forewarned.

def transform_epochs(epochs, f, **kwargs):
    '''Perform an arbitrary transformation on the epochs data'''
    epochs._data = f(epochs._data, **kwargs)

    return epochs

from scipy.signal import hilbert

epochs = transform_epochs(epochs, hilbert)

[LegrandNico]

So I tried with this function:

def epochs_hilbert(epochs, frequencies, crop, baseline):
    """
    Hilbert transformation on epochs data

    Parameters
    ----------
    frequencies : tuple (Names, freq_min, freq_max)
        Name and values of the frequency band to explore.
    crop : tuple (tmin, tmax)
        New crop to apply after the Hilbert transform.
    baseline: int
        Baseline to apply after the Hilbert transform.

    Returns
    -------
    epochs : instance of Epochs
        The epochs object with transformed data.
    """

    band, fmin, fmax = frequencies

    epochs.filter(fmin, fmax, n_jobs=1,
                  l_trans_bandwidth=1,
                  h_trans_bandwidth=1,
                  fir_design='firwin')

    # Hilbert transformation
    epochs._data = hilbert(epochs._data)

    # Crop the new epochs
    epochs.crop(tmin = crop[0], tmax = crop[1])

    # Aplly baseline
    epochs.apply_baseline(baseline=baseline)

    # remove evoked response and get analytic signal (envelope)
    epochs.subtract_evoked()
    epochs = mne.EpochsArray(data=np.abs(epochs.get_data()),
                             info=epochs.info,
                             tmin=epochs.tmin)
    return epochs

When applying this on the MNE examples with:

# set epoching parameters
event_id, tmin, tmax    = 1, -3., 5.
baseline                = (-1.0, 0)
crop                    = (-1., 3.0)

epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=None,
                    reject=dict(grad=4000e-13, eog=350e-6), preload=True)

frequency_map = list()

for frequencies in freq_list:

    epochs_h = epochs.copy()

    epochs_h = epochs_hilbert(epochs_h, frequencies, crop=crop, baseline=baseline)

    frequency_map.append((frequencies, epochs_h.average()))

plot_hilbert(frequency_map, baseline)
plt.show()

I obtained similar but not identical results, especially for the Alpha band. Note that I first extract the epochs with large interval (+2 seconds on each side) and I return to the classic interval after that.

I am more concerned with the theta amplitude, do you think this method would bias the results?

[LegrandNico]

The plot_hilbert function I used is borrowed from the same example:

def plot_hilbert(frequency_map, baseline):
    """
    Plot results from Hilbert transformation

    """
    fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True)
    colors = plt.get_cmap('winter_r')(np.linspace(0, 1, 4))

    for ((freq_name, fmin, fmax), average), color, ax in zip(
            frequency_map, colors, axes.ravel()[::-1]):
        times = average.times * 1e3

        # Compute GFP
        gfp = np.sum(average.data ** 2, axis=0)
        gfp = mne.baseline.rescale(gfp, times, baseline=baseline)

        ax.plot(times, gfp, label=freq_name, color=color, linewidth=2.5)
        ax.axhline(0, linestyle='--', color='grey', linewidth=2)

        # Compute CI
        ci_low, ci_up = _bootstrap_ci(average.data, random_state=0,
                                      stat_fun=lambda x: np.sum(x ** 2, axis=0))

        ci_low = rescale(ci_low, average.times, baseline=baseline)
        ci_up  = rescale(ci_up, average.times, baseline=baseline)

        ax.fill_between(times, gfp + ci_up, gfp - ci_low, color=color, alpha=0.3)
        ax.grid(True)
        ax.set_ylabel('GFP')
        ax.annotate('%s (%d-%dHz)' % (freq_name, fmin, fmax),
                    xy=(0.95, 0.8),
                    horizontalalignment='right',
                    xycoords='axes fraction')
        ax.set_xlim(-1000, 3000)

    axes.ravel()[-1].set_xlabel('Time [ms]')

[jona-sassenhagen]

In principle, wouldn't autoreject not also work on Hilbert-transformed data?

[agramfort]

why would you want to do this?

[jona-sassenhagen]

It would solve the issue here ...

[LegrandNico]

I'm not sure that autoreject is designed to find artifacts on filtered data. Plus, it may be better to keep artifact correction first in the pipeline before ICA to have the opportunity save clean data, and after that run any analysis you want. Applying autoreject at the end on each Hilbert transform can be computationally really expensive.

[agramfort]

I'm not sure that autoreject is designed to find artifacts on filtered data. Plus, it may be better to keep artifact correction first in the pipeline before ICA to have the opportunity save clean data,

+1

[palday]

And I would also note that doing the Hilbert transform on Epochs, especially with envelope=False is really handy in terms of memory.

[LegrandNico]

So do you think the method I tried (access as raw data -> filter -> crop) is worth being added to MNE? It still quite basic but as far I see no other alternative...

[agramfort]

looking are your plot I see something very comparable to :

http://martinos.org/mne/stable/auto_examples/time_frequency/plot_time_frequency_global_field_power.html#sphx-glr-auto-examples-time-frequency-plot-time-frequency-global-field-power-py

do you expect to see more with the code you suggest?

[LegrandNico]

My idea here was just to test an apply_hilbert method for epochs, so I tried to reproduce the plot from the website with this one. Plots are similar but not exactly the same, due to the cropped interval. As a consequence, results are slightly different when apply_hilbert applies to epochs or raw data, which can be confusing. So, although this turnaround is possible, it is not perfect and I wonder if it should be added to MNE.

[agramfort]

if there is not added value to what we already have, I suggest to leave it out to limit future software maintenance.