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

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

.. _sphx_glr_auto_examples_spectrum_interpolation_plot_01_spectrum_interpolation.py:


Spectrum Interpolation: Power-Line Noise Removal.
=================================================

This example demonstrates :class:`~mne_denoise.spectrum_interpolation.SpectrumInterpolation`
on a synthetic recording contaminated by 60 Hz line noise and its harmonics.

Spectrum interpolation removes the line frequency by replacing the *amplitude*
of a thin band around it (and its harmonics) with the mean amplitude of the
neighbouring bins, while keeping the phase unchanged. Compared with a notch
filter, this edits only a narrow amplitude band and leaves the broadband
spectrum and phase intact.

Reference:
    Leske, S., & Dalal, S. S. (2019). Reducing power line noise in EEG and MEG
    data via spectrum interpolation. NeuroImage, 189, 763-776.

Authors: Sina Esmaeili (sina.esmaeili@umontreal.ca)
         Hamza Abdelhedi (hamza.abdelhedi@umontreal.ca)

.. GENERATED FROM PYTHON SOURCE LINES 23-25

Imports
-------

.. GENERATED FROM PYTHON SOURCE LINES 25-31

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy import signal

    from mne_denoise.spectrum_interpolation import SpectrumInterpolation








.. GENERATED FROM PYTHON SOURCE LINES 32-37

Synthetic data
--------------
We build an 8-channel, 10-second recording that contains a 10 Hz oscillation
we want to keep, strong 60 Hz line noise with 120 and 180 Hz harmonics, and
broadband background noise.

.. GENERATED FROM PYTHON SOURCE LINES 37-52

.. code-block:: Python


    sfreq = 1000.0
    duration = 10.0
    n_channels = 8
    n_times = int(sfreq * duration)
    t = np.arange(n_times) / sfreq

    rng = np.random.default_rng(42)
    neural = np.sin(2 * np.pi * 10 * t)

    data = rng.standard_normal((n_channels, n_times)) * 0.5
    data += neural[None, :]
    for harmonic, amplitude in [(60.0, 3.0), (120.0, 1.5), (180.0, 0.8)]:
        data += amplitude * np.sin(2 * np.pi * harmonic * t)[None, :]








.. GENERATED FROM PYTHON SOURCE LINES 53-59

Apply spectrum interpolation
----------------------------
The estimator follows the scikit-learn ``fit`` / ``transform`` API. We target
60 Hz and remove all of its harmonics below the Nyquist frequency. Spectrum
interpolation is best suited to continuous recordings or long segments;
inspect short epochs for edge effects.

.. GENERATED FROM PYTHON SOURCE LINES 59-65

.. code-block:: Python


    si = SpectrumInterpolation(sfreq=sfreq, line_freq=60.0)
    clean = si.fit_transform(data)

    print(f"Targeted frequencies (Hz): {si.freqs_}")





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

 .. code-block:: none

    Targeted frequencies (Hz): [ 60. 120. 180. 240. 300. 360. 420. 480.]




.. GENERATED FROM PYTHON SOURCE LINES 66-70

Compare power spectra
---------------------
The line frequency and its harmonics drop to the noise floor, while the 10 Hz
peak and the surrounding broadband spectrum are preserved.

.. GENERATED FROM PYTHON SOURCE LINES 70-85

.. code-block:: Python


    freqs, psd_before = signal.welch(data, fs=sfreq, nperseg=2048)
    _, psd_after = signal.welch(clean, fs=sfreq, nperseg=2048)

    fig, ax = plt.subplots(figsize=(9, 4.5))
    ax.semilogy(freqs, psd_before.mean(0), label="Before", color="0.6")
    ax.semilogy(freqs, psd_after.mean(0), label="After", color="C0")
    for harmonic in (60.0, 120.0, 180.0):
        ax.axvline(harmonic, color="C3", ls="--", lw=0.8, alpha=0.6)
    ax.set_xlim(0, 220)
    ax.set_xlabel("Frequency (Hz)")
    ax.set_ylabel("Power spectral density")
    ax.set_title("Spectrum interpolation removes line noise and harmonics")
    ax.legend()
    fig.tight_layout()



.. image-sg:: /auto_examples/spectrum_interpolation/images/sphx_glr_plot_01_spectrum_interpolation_001.png
   :alt: Spectrum interpolation removes line noise and harmonics
   :srcset: /auto_examples/spectrum_interpolation/images/sphx_glr_plot_01_spectrum_interpolation_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_examples_spectrum_interpolation_plot_01_spectrum_interpolation.py:

.. only:: html

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

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

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

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

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

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

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