Source code for borsar.csd

import numpy as np


[docs]def current_source_density(inst, G, H, smoothing=1.0e-5, head_radius=1.):
    '''
    Compute current source density for given mne object instance.

    Note that this function works in-place.

    Parameters
    ----------
    inst : mne object instance
        Raw or Epochs data.
    G : numpy matrix
        G matrix.
    H : numpy matrix
        H matrix.
    smoothing : float
        CSD smoothing. Defaults to 1.0e-5.
    head_radius : float
        Radius of the head sphere.

    Returns
    -------
    inst : mne object instance
        The data with CSD reference.
    '''
    import mne
    from mne.utils import _get_inst_data

    data = _get_inst_data(inst)
    got_epochs = isinstance(inst, mne.Epochs)
    if got_epochs:
        n_epochs, n_channels, n_times = data.shape
    else:
        n_channels, n_times = data.shape

    # ensure the data is average referenced
    inst.set_eeg_reference('average', projection=False)

    # apply current source density
    data = _current_source_density(data, G, H, smoothing, head_radius)

    if not isinstance(inst, mne.Evoked):
        inst._data = data
    else:
        inst.data = data

    return inst


def _current_source_density(data, G, H, smoothing=1.0e-5, head_radius=1.):
    '''Python implementation of CSD.m from Matlab CSD toolbox.'''
    # FIXME add checks for G and H sizes

    n_channels, n_times = data.shape

    # average reference
    data -= data.mean(axis=0, keepdims=True)

    # check if bads exist and interpolate if necessary (or throw an error)

    # add lambda smoothing to the diagonal
    G = G + np.diag(np.ones(G.shape[0]) * smoothing)
    G_inv = np.linalg.inv(G)
    G_inv_row_sum = G_inv.sum(axis=1)
    G_inv_sum = G_inv_row_sum.sum()

    for idx in range(n_times):
        Cp = np.dot(G_inv, data[:, idx])
        c0 = Cp.sum() / G_inv_sum
        C = Cp - np.dot(c0, G_inv_row_sum.T)

        # for ch_idx in range(n_channels):
        #     data[ch_idx, idx] = (C * H[ch_idx].T).sum() / head_radius
        data[:, idx] = np.dot(C, H) / head_radius

    return data