.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/asr/plot_07_riemannian_asr.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_asr_plot_07_riemannian_asr.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_asr_plot_07_riemannian_asr.py:


Riemannian ASR versus standard ASR.
===================================

Blum et al. (2019) replace ASR's sample-covariance calibration with a
**Riemannian geometric median**, making the clean-covariance estimate robust to
the occasional contaminated calibration window. ``mne-denoise`` exposes this as
``method="riemannian_windowed"`` --- it keeps the Riemannian robust calibration
but applies a standard per-window eigendecomposition at processing time, so
(unlike the MATLAB-parity ``method="riemannian"``) the ``cutoff`` knob still
works.

This example runs both backends on real EEG with eye blinks (the MNE *sample*
dataset) and compares frontal blink removal (cf. Blum 2019, Fig. 3). On this
relatively clean recording the two are comparable; rASR's documented edge ---
robustness when the calibration windows themselves are contaminated --- is by
design and is not stressed here.

References
----------
* Blum, S., Jacobsen, N. S. J., Bleichner, M. G., & Debener, S. (2019).
  A Riemannian Modification of Artifact Subspace Reconstruction for EEG
  Artifact Handling. Frontiers in Human Neuroscience, 13, 141.
  doi:10.3389/fnhum.2019.00141
* Chang, C.-Y., et al. (2020). Evaluation of Artifact Subspace
  Reconstruction... IEEE TBME, 67(4), 1114-1121. doi:10.1109/TBME.2019.2930186

.. GENERATED FROM PYTHON SOURCE LINES 30-32

Load real EEG with eye blinks
-----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 32-49

.. code-block:: Python

    import matplotlib.pyplot as plt
    import mne
    import numpy as np

    from mne_denoise.asr import ASR
    from mne_denoise.viz import plot_signal_overlay

    sample = mne.datasets.sample.data_path()
    raw = mne.io.read_raw_fif(
        sample / "MEG" / "sample" / "sample_audvis_raw.fif",
        preload=True,
        verbose="ERROR",
    )
    raw.pick(["eeg", "eog"]).crop(0, 60).resample(160, verbose="ERROR")
    raw.set_eeg_reference("average", verbose="ERROR")
    raw.filter(1.0, None, verbose="ERROR")  # ASR assumes high-pass-filtered data






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <script type="text/javascript">
        // must be `var` (not `const`) because this can get embedded multiple times on a page
    var toggleVisibility = (className) => {

        const elements = document.querySelectorAll(`.${className}`);

        elements.forEach(element => {
            if (element.classList.contains("mne-repr-section-header")) {
                return  // Don't collapse the section header row
            }
            element.classList.toggle("mne-repr-collapsed");
        });

        // trigger caret to rotate
        var sel = `.mne-repr-section-header.${className} > th.mne-repr-section-toggle > button`;
        const button = document.querySelector(sel);
        button.classList.toggle("collapsed");

        // adjust tooltip
        sel = `tr.mne-repr-section-header.${className}`;
        const secHeadRow = document.querySelector(sel);
        secHeadRow.classList.toggle("collapsed");
        secHeadRow.title = secHeadRow.title === "Hide section" ? "Show section" : "Hide section";
    }
    </script>

    <style type="text/css">
        /*
    Styles in this section apply both to the sphinx-built website docs and to notebooks
    rendered in an IDE or in Jupyter. In our web docs, styles here are complemented by
    doc/_static/styles.css and other CSS files (e.g. from the sphinx theme, sphinx-gallery,
    or bootstrap). In IDEs/Jupyter, those style files are unavailable, so only the rules in
    this file apply (plus whatever default styling the IDE applies).
    */
    .mne-repr-table {
        display: inline;  /* prevent using full container width */
    }
    .mne-repr-table tr.mne-repr-section-header > th {
        padding-top: 1rem;
        text-align: left;
        vertical-align: middle;
    }
    .mne-repr-section-toggle > button {
        all: unset;
        display: block;
        height: 1rem;
        width: 1rem;
    }
    .mne-repr-section-toggle > button > svg {
        height: 60%;
    }

    /* transition (rotation) effects on the collapser button */
    .mne-repr-section-toggle > button.collapsed > svg {
        transition: 0.1s ease-out;
        transform: rotate(-90deg);
    }
    .mne-repr-section-toggle > button:not(.collapsed) > svg {
        transition: 0.1s ease-out;
        transform: rotate(0deg);
    }

    /* hide collapsed table rows */
    .mne-repr-collapsed {
        display: none;
    }


    @layer {
        /*
        Selectors in a `@layer` will always be lower-precedence than selectors outside the
        layer. So even though e.g. `div.output_html` is present in the sphinx-rendered
        website docs, the styles here won't take effect there as long as some other rule
        somewhere in the page's CSS targets the same element.

        In IDEs or Jupyter notebooks, though, the CSS files from the sphinx theme,
        sphinx-gallery, and bootstrap are unavailable, so these styles will apply.

        Notes:

        - the selector `.accordion-body` is for MNE Reports
        - the selector `.output_html` is for VSCode's notebook interface
        - the selector `.jp-RenderedHTML` is for Jupyter notebook
        - variables starting with `--theme-` are VSCode-specific.
        - variables starting with `--jp-` are Jupyter styles, *some of which* are also
          available in VSCode. Here we try the `--theme-` variable first, then fall back to
          the `--jp-` ones.
        */
        .mne-repr-table {
            --mne-toggle-color: var(--theme-foreground, var(--jp-ui-font-color1));
            --mne-button-bg-color: var(--theme-button-background, var(--jp-info-color0, var(--jp-content-link-color)));
            --mne-button-fg-color: var(--theme-button-foreground, var(--jp-ui-inverse-font-color0, var(--jp-editor-background)));
            --mne-button-hover-bg-color: var(--theme-button-hover-background, var(--jp-info-color1));
            --mne-button-radius: var(--jp-border-radius, 0.25rem);
        }
        /* chevron position/alignment; in VSCode it looks ok without adjusting */
        .accordion-body .mne-repr-section-toggle > button,
        .jp-RenderedHTML .mne-repr-section-toggle > button {
            padding: 0 0 45% 25% !important;
        }
        /* chevron color; MNE Report doesn't have light/dark mode */
        div.output_html .mne-repr-section-toggle > button > svg > path,
        .jp-RenderedHTML .mne-repr-section-toggle > button > svg > path {
            fill: var(--mne-toggle-color);
        }
        .accordion-body .mne-ch-names-btn,
        div.output_html .mne-ch-names-btn,
        .jp-RenderedHTML .mne-ch-names-btn {
            -webkit-border-radius: var(--mne-button-radius);
            -moz-border-radius: var(--mne-button-radius);
            border-radius: var(--mne-button-radius);
            border: none;
            background-image: none;
            background-color: var(--mne-button-bg-color);
            color: var(--mne-button-fg-color);
            font-size: inherit;
            min-width: 1.5rem;
            padding: 0.25rem;
            text-align: center;
            text-decoration: none;
        }
        .accordion-body .mne-ch-names-btn:hover,
        div.output_html .mne.ch-names-btn:hover,
        .jp-RenderedHTML .mne-ch-names-btn:hover {
            background-color: var(--mne-button-hover-bg-color);
            text-decoration: underline;
        }
        .accordion-body .mne-ch-names-btn:focus-visible,
        div.output_html .mne-ch-names-btn:focus-visible,
        .jp-RenderedHTML .mne-ch-names-btn:focus-visible {
            outline: 0.1875rem solid var(--mne-button-bg-color) !important;
            outline-offset: 0.1875rem !important;
        }
    }
    </style>



    <table class="table mne-repr-table">
    







    <tr class="mne-repr-section-header general-26a0423c-6e32-49f6-a462-b9fffd7b2593"
         title="Hide section" 
        onclick="toggleVisibility('general-26a0423c-6e32-49f6-a462-b9fffd7b2593')">
        <th class="mne-repr-section-toggle">
            <button >
                <svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 512 512"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d="M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z"/></svg>
            </button>
        </th>
        <th colspan="2">
            <strong>General</strong>
        </th>
    </tr>


    <tr class="repr-element general-26a0423c-6e32-49f6-a462-b9fffd7b2593 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Filename(s)</td>
        <td>
        
            sample_audvis_raw.fif
        
        
        </td>
    </tr>

    <tr class="repr-element general-26a0423c-6e32-49f6-a462-b9fffd7b2593 ">
        <td class="mne-repr-section-toggle"></td>
        <td>MNE object type</td>
        <td>Raw</td>
    </tr>
    <tr class="repr-element general-26a0423c-6e32-49f6-a462-b9fffd7b2593 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Measurement date</td>
    
        <td>2002-12-03 at 19:01:10 UTC</td>
    
    </tr>
    <tr class="repr-element general-26a0423c-6e32-49f6-a462-b9fffd7b2593 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Participant</td>
    
        <td>Unknown</td>
    
    </tr>
    <tr class="repr-element general-26a0423c-6e32-49f6-a462-b9fffd7b2593 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Experimenter</td>
    
        <td>MEG</td>
    
    </tr>
    







    <tr class="mne-repr-section-header acquisition-b4eef280-1d75-440e-8a35-b98212880220"
         title="Hide section" 
        onclick="toggleVisibility('acquisition-b4eef280-1d75-440e-8a35-b98212880220')">
        <th class="mne-repr-section-toggle">
            <button >
                <svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 512 512"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d="M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z"/></svg>
            </button>
        </th>
        <th colspan="2">
            <strong>Acquisition</strong>
        </th>
    </tr>


    <tr class="repr-element acquisition-b4eef280-1d75-440e-8a35-b98212880220 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Duration</td>
        <td>00:01:00 (HH:MM:SS)</td>
    </tr>








    <tr class="repr-element acquisition-b4eef280-1d75-440e-8a35-b98212880220 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Sampling frequency</td>
        <td>160.00 Hz</td>
    </tr>


    <tr class="repr-element acquisition-b4eef280-1d75-440e-8a35-b98212880220 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Time points</td>
        <td>9,600</td>
    </tr>


    







    <tr class="mne-repr-section-header channels-b5ee876a-ecb5-42be-b5b9-2dbf317d2d13"
         title="Hide section" 
        onclick="toggleVisibility('channels-b5ee876a-ecb5-42be-b5b9-2dbf317d2d13')">
        <th class="mne-repr-section-toggle">
            <button >
                <svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 512 512"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d="M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z"/></svg>
            </button>
        </th>
        <th colspan="2">
            <strong>Channels</strong>
        </th>
    </tr>


    
    <tr class="repr-element channels-b5ee876a-ecb5-42be-b5b9-2dbf317d2d13 ">
        <td class="mne-repr-section-toggle"></td>
        <td>EEG</td>
        <td>
            <button class="mne-ch-names-btn sd-sphinx-override sd-btn sd-btn-info sd-text-wrap sd-shadow-sm" onclick="alert('Good EEG:\n\nEEG&nbsp;001, EEG&nbsp;002, EEG&nbsp;003, EEG&nbsp;004, EEG&nbsp;005, EEG&nbsp;006, EEG&nbsp;007, EEG&nbsp;008, EEG&nbsp;009, EEG&nbsp;010, EEG&nbsp;011, EEG&nbsp;012, EEG&nbsp;013, EEG&nbsp;014, EEG&nbsp;015, EEG&nbsp;016, EEG&nbsp;017, EEG&nbsp;018, EEG&nbsp;019, EEG&nbsp;020, EEG&nbsp;021, EEG&nbsp;022, EEG&nbsp;023, EEG&nbsp;024, EEG&nbsp;025, EEG&nbsp;026, EEG&nbsp;027, EEG&nbsp;028, EEG&nbsp;029, EEG&nbsp;030, EEG&nbsp;031, EEG&nbsp;032, EEG&nbsp;033, EEG&nbsp;034, EEG&nbsp;035, EEG&nbsp;036, EEG&nbsp;037, EEG&nbsp;038, EEG&nbsp;039, EEG&nbsp;040, EEG&nbsp;041, EEG&nbsp;042, EEG&nbsp;043, EEG&nbsp;044, EEG&nbsp;045, EEG&nbsp;046, EEG&nbsp;047, EEG&nbsp;048, EEG&nbsp;049, EEG&nbsp;050, EEG&nbsp;051, EEG&nbsp;052, EEG&nbsp;054, EEG&nbsp;055, EEG&nbsp;056, EEG&nbsp;057, EEG&nbsp;058, EEG&nbsp;059, EEG&nbsp;060')" title="(Click to open in popup)&#13;&#13;EEG&nbsp;001, EEG&nbsp;002, EEG&nbsp;003, EEG&nbsp;004, EEG&nbsp;005, EEG&nbsp;006, EEG&nbsp;007, EEG&nbsp;008, EEG&nbsp;009, EEG&nbsp;010, EEG&nbsp;011, EEG&nbsp;012, EEG&nbsp;013, EEG&nbsp;014, EEG&nbsp;015, EEG&nbsp;016, EEG&nbsp;017, EEG&nbsp;018, EEG&nbsp;019, EEG&nbsp;020, EEG&nbsp;021, EEG&nbsp;022, EEG&nbsp;023, EEG&nbsp;024, EEG&nbsp;025, EEG&nbsp;026, EEG&nbsp;027, EEG&nbsp;028, EEG&nbsp;029, EEG&nbsp;030, EEG&nbsp;031, EEG&nbsp;032, EEG&nbsp;033, EEG&nbsp;034, EEG&nbsp;035, EEG&nbsp;036, EEG&nbsp;037, EEG&nbsp;038, EEG&nbsp;039, EEG&nbsp;040, EEG&nbsp;041, EEG&nbsp;042, EEG&nbsp;043, EEG&nbsp;044, EEG&nbsp;045, EEG&nbsp;046, EEG&nbsp;047, EEG&nbsp;048, EEG&nbsp;049, EEG&nbsp;050, EEG&nbsp;051, EEG&nbsp;052, EEG&nbsp;054, EEG&nbsp;055, EEG&nbsp;056, EEG&nbsp;057, EEG&nbsp;058, EEG&nbsp;059, EEG&nbsp;060">
                59
            </button>

        
        
            and <button class="mne-ch-names-btn sd-sphinx-override sd-btn sd-btn-info sd-text-wrap sd-shadow-sm" onclick="alert('Bad EEG:\n\nEEG&nbsp;053')" title="(Click to open in popup)&#13;&#13;EEG&nbsp;053">
                1 bad
            </button>
        
        </td>
    </tr>

    
    <tr class="repr-element channels-b5ee876a-ecb5-42be-b5b9-2dbf317d2d13 ">
        <td class="mne-repr-section-toggle"></td>
        <td>EOG</td>
        <td>
            <button class="mne-ch-names-btn sd-sphinx-override sd-btn sd-btn-info sd-text-wrap sd-shadow-sm" onclick="alert('Good EOG:\n\nEOG&nbsp;061')" title="(Click to open in popup)&#13;&#13;EOG&nbsp;061">
                1
            </button>

        
        </td>
    </tr>


    <tr class="repr-element channels-b5ee876a-ecb5-42be-b5b9-2dbf317d2d13 ">
        <td class="mne-repr-section-toggle"></td>
        <td>Head & sensor digitization</td>
    
        <td>146 points</td>
    
    </tr>
    







    <tr class="mne-repr-section-header filters-040f3b59-9248-4597-85cd-b4fed666ff1e"
         title="Hide section" 
        onclick="toggleVisibility('filters-040f3b59-9248-4597-85cd-b4fed666ff1e')">
        <th class="mne-repr-section-toggle">
            <button >
                <svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 512 512"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d="M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z"/></svg>
            </button>
        </th>
        <th colspan="2">
            <strong>Filters</strong>
        </th>
    </tr>


    <tr class="repr-element filters-040f3b59-9248-4597-85cd-b4fed666ff1e ">
        <td class="mne-repr-section-toggle"></td>
        <td>Highpass</td>
        <td>1.00 Hz</td>
    </tr>


    <tr class="repr-element filters-040f3b59-9248-4597-85cd-b4fed666ff1e ">
        <td class="mne-repr-section-toggle"></td>
        <td>Lowpass</td>
        <td>80.00 Hz</td>
    </tr>


    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 50-53

Find the blink-dominated EEG channel
------------------------------------
The channel most correlated with the EOG carries the strongest blinks.

.. GENERATED FROM PYTHON SOURCE LINES 53-60

.. code-block:: Python

    eeg_names = [raw.ch_names[i] for i in mne.pick_types(raw.info, eeg=True)]
    eeg = raw.get_data(picks="eeg")
    eog = raw.get_data(picks="eog")[0]
    corr = np.array([abs(np.corrcoef(ch, eog)[0, 1]) for ch in eeg])
    blink_ch = eeg_names[int(np.argmax(corr))]
    print(f"Most blink-correlated EEG channel: {blink_ch} (|r|={corr.max():.2f})")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Most blink-correlated EEG channel: EEG 002 (|r|=0.81)




.. GENERATED FROM PYTHON SOURCE LINES 61-63

Clean with standard and Riemannian-windowed ASR
-----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 63-69

.. code-block:: Python

    asr_std = ASR(cutoff=20.0, picks="eeg", method="standard", verbose=False)
    clean_std = asr_std.fit_transform(raw.copy())

    asr_rie = ASR(cutoff=20.0, picks="eeg", method="riemannian_windowed", verbose=False)
    clean_rie = asr_rie.fit_transform(raw.copy())








.. GENERATED FROM PYTHON SOURCE LINES 70-72

Before/after on the blink channel (Riemannian-windowed)
-------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 72-86

.. code-block:: Python

    plot_signal_overlay(
        raw,
        clean_rie,
        raw.times,
        pick=blink_ch,
        scale_after=False,
        before_label="raw",
        after_label="rASR-cleaned",
        x_label="Time (s)",
        y_label="Amplitude (V)",
        title=f"Riemannian-windowed ASR on {blink_ch}",
        show=False,
    )




.. image-sg:: /auto_examples/asr/images/sphx_glr_plot_07_riemannian_asr_001.png
   :alt: Riemannian-windowed ASR on EEG 002
   :srcset: /auto_examples/asr/images/sphx_glr_plot_07_riemannian_asr_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Figure size 2400x800 with 1 Axes>



.. GENERATED FROM PYTHON SOURCE LINES 87-93

Blink removal: decoupling the frontal EEG from the EOG (cf. Blum 2019, Fig. 3)
------------------------------------------------------------------------------
A blink-specific score: mean absolute correlation between the 8 most
blink-correlated EEG channels and the EOG, before vs after cleaning (lower is
better). Both backends collapse it from ~0.65 to <0.1 --- comparable here,
standard marginally more aggressive.

.. GENERATED FROM PYTHON SOURCE LINES 93-116

.. code-block:: Python

    top = np.argsort(corr)[-8:]


    def mean_eog_coupling(data_eeg):
        return float(np.mean([abs(np.corrcoef(ch, eog)[0, 1]) for ch in data_eeg[top]]))


    coupling = {
        "raw": mean_eog_coupling(eeg),
        "standard": mean_eog_coupling(clean_std.get_data(picks="eeg")),
        "riemannian_windowed": mean_eog_coupling(clean_rie.get_data(picks="eeg")),
    }
    for name, value in coupling.items():
        print(f"  {name:20s} mean |corr with EOG| = {value:.3f}")

    fig, ax = plt.subplots(figsize=(5.5, 4))
    ax.bar(list(coupling), list(coupling.values()), color=["C3", "0.6", "C0"])
    ax.set_ylabel("mean |corr with EOG| (8 frontal channels)")
    ax.set_title("Blink coupling before / after ASR (lower = blinks removed)")
    ax.tick_params(axis="x", labelrotation=15)
    fig.tight_layout()

    plt.show()



.. image-sg:: /auto_examples/asr/images/sphx_glr_plot_07_riemannian_asr_002.png
   :alt: Blink coupling before / after ASR (lower = blinks removed)
   :srcset: /auto_examples/asr/images/sphx_glr_plot_07_riemannian_asr_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

      raw                  mean |corr with EOG| = 0.654
      standard             mean |corr with EOG| = 0.467
      riemannian_windowed  mean |corr with EOG| = 0.210





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 8.555 seconds)


.. _sphx_glr_download_auto_examples_asr_plot_07_riemannian_asr.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_07_riemannian_asr.ipynb <plot_07_riemannian_asr.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_07_riemannian_asr.py <plot_07_riemannian_asr.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_07_riemannian_asr.zip <plot_07_riemannian_asr.zip>`