.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_star_dss.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_example_star_dss.py>`
        to download the full example code

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

.. _sphx_glr_auto_examples_example_star_dss.py:


Example demonstrating STAR + DSS
================================

This example shows how one can effectively combine STAR and DSS to recover
signal components which would not have been discoverable with either these
two techniques alone, due to the presence of strong artifacts.

This example replicates figure 1 in [1]_.

References
----------
.. [1] de Cheveigné A (2016) Sparse Time Artifact Removal, Journal of
   Neuroscience Methods, 262, 14-20, doi:10.1016/j.jneumeth.2016.01.005

.. GENERATED FROM PYTHON SOURCE LINES 17-28

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.optimize import leastsq

    from meegkit import dss, star
    from meegkit.utils import demean, normcol, tscov

    # import config  # noqa

    rng = np.random.default_rng(9)








.. GENERATED FROM PYTHON SOURCE LINES 29-33

Create simulated data
-----------------------------------------------------------------------------
Simulated data consist of N channels, 1 sinusoidal target, N-3 noise sources,
with temporally local artifacts on each channel.

.. GENERATED FROM PYTHON SOURCE LINES 33-70

.. code-block:: Python


    # Source
    n_chans, n_samples = 10, 1000
    f = 2
    target = np.sin(np.arange(n_samples) / n_samples * 2 * np.pi * f)
    target = target[:, np.newaxis]
    noise = rng.standard_normal((n_samples, n_chans - 3))

    # Create artifact signal
    SNR = np.sqrt(1)
    x0 = normcol(np.dot(noise, rng.standard_normal((noise.shape[1], n_chans)))) + \
        SNR * target * rng.standard_normal((1, n_chans))
    x0 = demean(x0)
    artifact = np.zeros(x0.shape)
    for k in np.arange(n_chans):
        artifact[k * 100 + np.arange(20), k] = 1
    x = x0 + 10 * artifact


    def _sine_fit(x):
        """Fit a sinusoidal trend."""
        guess_mean = np.mean(x)
        guess_std = np.std(x)
        guess_phase = 0
        t = np.linspace(0, 4 * np.pi, x.shape[0])

        # Optimization function, in this case, we want to minimize the difference
        # between the actual data and our "guessed" parameters
        def func(y):
            return np.mean(x - (y[0] * np.sin(t + y[1]) + y[2])[:, None], 1)

        est_std, est_phase, est_mean = leastsq(
            func, [guess_std, guess_phase, guess_mean])[0]
        data_fit = est_std * np.sin(t + est_phase) + est_mean
        return np.tile(data_fit, (x.shape[1], 1)).T









.. GENERATED FROM PYTHON SOURCE LINES 71-73

1) Apply STAR
-----------------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 73-75

.. code-block:: Python

    y, w, _ = star.star(x, 2)





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

 .. code-block:: none

    proportion artifact free: 0.70
    proportion artifact free: 0.70
    proportion artifact free: 0.70
    depth: 1
    fixed channels: 10
    fixed samples: 299
    ratio: 1.01
    power ratio: 0.40




.. GENERATED FROM PYTHON SOURCE LINES 76-78

2) Apply DSS on raw data
-----------------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 78-83

.. code-block:: Python

    c0, _ = tscov(x)
    c1, _ = tscov(x - _sine_fit(x))
    [todss, _, pwr0, pwr1] = dss.dss0(c0, c1)
    z1 = normcol(np.dot(x, todss))








.. GENERATED FROM PYTHON SOURCE LINES 84-87

3) Apply DSS on STAR-ed data
-----------------------------------------------------------------------------
Here the bias function is the original signal minus the sinusoidal trend.

.. GENERATED FROM PYTHON SOURCE LINES 87-92

.. code-block:: Python

    c0, _ = tscov(y)
    c1, _ = tscov(y - _sine_fit(y))
    [todss, _, pwr0, pwr1] = dss.dss0(c0, c1)
    z2 = normcol(np.dot(y, todss))








.. GENERATED FROM PYTHON SOURCE LINES 93-95

Plots
-----------------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 95-112

.. code-block:: Python

    f, (ax0, ax1, ax2, ax3) = plt.subplots(4, 1, figsize=(7, 9))
    ax0.plot(target, lw=.5)
    ax0.set_title("Target")

    ax1.plot(x, lw=.5)
    ax1.set_title(f"Signal + Artifacts (SNR = {SNR})")

    ax2.plot(z1[:, 0], lw=.5, label="Best DSS component")
    ax2.set_title("DSS")
    ax2.legend(loc="lower right")

    ax3.plot(z2[:, 0], lw=.5, label="Best DSS component")
    ax3.set_title("STAR + DSS")
    ax3.legend(loc="lower right")

    f.set_tight_layout(True)
    plt.show()



.. image-sg:: /auto_examples/images/sphx_glr_example_star_dss_001.png
   :alt: Target, Signal + Artifacts (SNR = 1.0), DSS, STAR + DSS
   :srcset: /auto_examples/images/sphx_glr_example_star_dss_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_examples_example_star_dss.py:

.. only:: html

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

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

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

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

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


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_