'Subtract an evoked response from each epoch.
Can be used to exclude the evoked response when analyzing induced
activity, see e.g. [1].
References
[1] David et al. "Mechanisms of evoked and induced responses in
MEG/EEG", NeuroImage, vol. 31, no. 4, pp. 1580-1591, July 2006.
Parameters
evoked : instance of Evoked | None
The evoked response to subtract. If None, the evoked response
is computed from Epochs itself.
Returns
self : instance of Epochs
The modified instance (instance is also modified inplace).'

def subtract_evoked(self, evoked=None):

    logger.info('Subtracting    Evoked    from    Epochs')
    if (evoked is None):
        picks = _pick_data_channels(self.info, exclude=[])
        evoked = self.average(picks)
    picks = pick_channels(evoked.ch_names, include=self.ch_names)
    if (len(picks) < len(self.ch_names)):
        sel_ch = [evoked.ch_names[ii] for ii in picks]
        diff_ch = list(set(self.ch_names).difference(sel_ch))
        diff_idx = [self.ch_names.index(ch) for ch in diff_ch]
        diff_types = [channel_type(self.info, idx) for idx in diff_idx]
        bad_idx = [diff_types.index(t) for t in diff_types if (t in _DATA_CH_TYPES_SPLIT)]
        if (len(bad_idx) > 0):
            bad_str = ',    '.join([diff_ch[ii] for ii in bad_idx])
            raise ValueError(('The    following    data    channels    are    missing    in    the    evoked    response:    %s' % bad_str))
        logger.info(('                The    following    channels    are    not    included    in    the    subtraction:    %s' % ',    '.join(diff_ch)))
    if ((len(self.times) != len(evoked.times)) or (np.max(np.abs((self.times - evoked.times))) >= 1e-07)):
        raise ValueError('Epochs    and    Evoked    object    do    not    contain    the    same    time    points.')
    if ((not self.proj) and evoked.proj):
        warn('Evoked    has    SSP    applied    while    Epochs    has    not.')
    if (self.proj and (not evoked.proj)):
        evoked = evoked.copy().apply_proj()
    ep_picks = [self.ch_names.index(evoked.ch_names[ii]) for ii in picks]
    if self.preload:
        self._data[:, ep_picks, :] -= evoked.data[picks][None, :, :]
    else:
        if (self._offset is None):
            self._offset = np.zeros((len(self.ch_names), len(self.times)), dtype=np.float)
        self._offset[ep_picks] -= evoked.data[picks]
    logger.info('[done]')
    return self
