Source code for mne_nirs.signal_enhancement._negative_correlation

# Authors: Robert Luke <mail@robertluke.net>
#
# License: BSD (3-clause)

import numpy as np
from mne import pick_types
from mne.io import BaseRaw
from mne.utils import _validate_type




def enhance_negative_correlation(raw):
    """
    Apply negative correlation enhancement algorithm.

    As described in :footcite:`cui2010functional`.
    This function will return a modified raw instance that has the
    enhancement applied to it.
    This function can only be run on data containing
    hbo and hbr types.

    Parameters
    ----------
    raw : instance of Raw
        Haemoglobin data.

    Returns
    -------
    raw : instance of Raw
        The modified raw instance.

    References
    ----------
    .. [1] Cui et al, Functional Near Infrared Spectroscopy (NIRS) signal
           improvement based on negative correlation between oxygenated and
           deoxygenated hemoglobin dynamics, NeuroImage
           https://doi.org/10.1016/j.neuroimage.2009.11.050
    """
    raw = raw.copy().load_data()
    _validate_type(raw, BaseRaw, "raw")

    hbo_channels = pick_types(raw.info, fnirs="hbo")
    hbr_channels = pick_types(raw.info, fnirs="hbr")

    if (not len(hbo_channels)) & (not len(hbr_channels)):
        raise RuntimeError(
            "enhance_negative_correlation should " "be run on haemoglobin data."
        )

    if len(hbo_channels) != len(hbr_channels):
        raise RuntimeError("Same number of hbo and hbr channels required.")

    for idx in range(len(hbo_channels)):
        if (
            raw.info["chs"][hbo_channels[idx]]["ch_name"][:-4]
            != raw.info["chs"][hbr_channels[idx]]["ch_name"][:-4]
        ):
            raise RuntimeError("Channels must alternate between HBO and HBR.")

    for idx in range(len(hbo_channels)):
        hbo = raw._data[hbo_channels[idx]]
        hbr = raw._data[hbr_channels[idx]]

        hbo = hbo - np.mean(hbo)
        hbr = hbr - np.mean(hbr)

        hbo_std = np.std(hbo)
        hbr_std = np.std(hbr)

        alpha = hbo_std / hbr_std

        raw._data[hbo_channels[idx]] = 0.5 * (hbo - alpha * hbr)
        raw._data[hbr_channels[idx]] = -(1 / alpha) * raw._data[hbo_channels[idx]]

    return raw