dipole_fwhm.py

import numpy as np
import mne


def dipole_fwhm(dip, tmin=None, tmax=None, gof_thresh=0.):
    """Compute the time interval of full-width half-maximum of a dipole
    .. note:: This only considers the max(), not the min(), so if you're
              interested in large negative values, you could do
              ``dip.data = -dip.data``.
    Parameters
    ----------
    dip : instance of DipoleFixed
        The dipole to process.
    Returns
    -------
    t_start : float
        The first time point within 50% of the maximum response.
    t_stop : float
        The last time point within 50% of the maximum response.
    """
    assert isinstance(dip, mne.DipoleFixed)
    pick_data = mne.pick_types(dip.info, meg=False, dipole=True)
    pick_gof = mne.pick_types(dip.info, meg=False, gof=True)
    assert len(pick_data) == len(pick_gof) == 1
    data, gof = dip.data[[pick_data[0], pick_gof[0]]]
    idx = np.argmax(data)
    # At least 50%
    valid = (data >= 0.5 * data[idx])
    # And in the interval
    valid = np.logical_and(valid, mne.utils._time_mask(dip.times, tmin, tmax))
    # And sufficient GOF
    valid = np.logical_and(valid, gof >= gof_thresh)
    # make sure it's contiguous
    valid = np.where(valid)[0]
    breaks = np.where(np.diff(valid) != 1)[0]
    start = valid[breaks[valid[breaks] < idx] + 1]
    stop = valid[breaks[valid[breaks] >= idx]]
    start = start[-1] if len(start) > 0 else valid[0]
    stop = (stop[0] if len(stop) else valid[-1])
    return dip.times[[start, stop]]


gof_thresh = 0.3
tmin = 0.42
dip = mne.read_dipole(mne.datasets.testing.data_path() + '/dip/fixed_auto.fif')
tmin, tmax = dipole_fwhm(dip, tmin=tmin, gof_thresh=gof_thresh)

# do a little plot to verify it
fig = dip.plot()
for ai, ax in enumerate(fig.get_axes()):
    line = ax.get_lines()[0]
    xy = line.get_xydata()
    mask = np.logical_and(xy[:, 0] >= 1000 * tmin, xy[:, 0] <= 1000 * tmax)
    if ai == 0:
        y = 0.5 * np.max(xy[:, 1])
    else:
        y = gof_thresh
    ax.hlines(y, xy[0, 0], xy[-1, 0], color='0.7', linestyle=':')
    ax.vlines(1000 * tmin, xy[:, 1].min(), xy[:, 1].max(), color='0.7',
              linestyle=':')
    ax.fill_between(xy[mask][:, 0], 0, xy[mask][:, 1], color='r')