update documentation for release 1.13.2

Author: github-actions[bot]

1.13.1/_downloads/0e1f06ac031754b1722bc7b181b7568d/10_peak_detection.ipynb

@@ -0,0 +1,281 @@
+"""
+Bandpower rolling window
+========================
+
+With a :class:`~mne_lsl.stream.StreamLSL`, we can compute the bandpower on a time
+rolling window. For this example, we will look at the alpha band power, between 8 and
+13 Hz.
+"""
+
+# sphinx_gallery_thumbnail_path = '_static/tutorials/bp-performance.png'
+
+import time
+import uuid
+
+import numpy as np
+from matplotlib import colormaps
+from matplotlib import pyplot as plt
+from mne.io import read_raw_fif
+from mne.time_frequency import psd_array_multitaper
+from numpy.typing import NDArray
+from scipy.integrate import simpson
+from scipy.signal import periodogram, welch
+
+from mne_lsl.datasets import sample
+from mne_lsl.player import PlayerLSL
+from mne_lsl.stream import StreamLSL
+
+# dataset used in the example
+raw = read_raw_fif(sample.data_path() / "sample-ant-raw.fif", preload=False)
+raw.crop(40, 60).load_data()
+raw
+
+# %%
+# Preprocessing
+# -------------
+#
+# In a real-time scenario, we would want to apply artifact rejection methods online to
+# estimate the bandpower on brain signals, not on artifacts. For this example, we will
+# only apply a bandpass filter to the data.
+#
+# Estimating the bandpower
+# ------------------------
+#
+# First, we will define the function estimating the bandpower on a time window. The
+# bandpower will be estimated by integrating the power spectral density (PSD) on the
+# frequency band of interest, using the composite Simpson's rule
+# (:func:`scipy.integrate.simpson`).
+
+
+def bandpower(
+    data: NDArray[np.float64],
+    fs: float,
+    method: str,
+    band: tuple[float, float],
+    relative: bool = True,
+    **kwargs,
+) -> NDArray[np.float64]:
+    """Compute the bandpower of the individual channels.
+
+    Parameters
+    ----------
+    data : array of shape (n_channels, n_samples)
+        Data on which the the bandpower is estimated.
+    fs : float
+        Sampling frequency in Hz.
+    method : 'periodogram' | 'welch' | 'multitaper'
+        Method used to estimate the power spectral density.
+    band : tuple of shape (2,)
+        Frequency band of interest in Hz as 2 floats, e.g. ``(8, 13)``. The
+        edges are included.
+    relative : bool
+        If True, the relative bandpower is returned instead of the absolute
+        bandpower.
+    **kwargs : dict
+        Additional keyword arguments are provided to the power spectral density
+        estimation function.
+        * 'periodogram': scipy.signal.periodogram
+        * 'welch'``: scipy.signal.welch
+        * 'multitaper': mne.time_frequency.psd_array_multitaper
+
+        The only provided arguments are the data array and the sampling
+        frequency.
+
+    Returns
+    -------
+    bandpower : array of shape (n_channels,)
+        The bandpower of each channel.
+    """
+    # compute the power spectral density
+    assert data.ndim == 2, (
+        "The provided data must be a 2D array of shape (n_channels, n_samples)."
+    )
+    if method == "periodogram":
+        freqs, psd = periodogram(data, fs, **kwargs)
+    elif method == "welch":
+        freqs, psd = welch(data, fs, **kwargs)
+    elif method == "multitaper":
+        psd, freqs = psd_array_multitaper(data, fs, verbose="ERROR", **kwargs)
+    else:
+        raise RuntimeError(f"The provided method '{method}' is not supported.")
+    # compute the bandpower
+    assert len(band) == 2, "The 'band' argument must be a 2-length tuple."
+    assert band[0] <= band[1], (
+        "The 'band' argument must be defined as (low, high) (in Hz)."
+    )
+    freq_res = freqs[1] - freqs[0]
+    idx_band = np.logical_and(freqs >= band[0], freqs <= band[1])
+    bandpower = simpson(psd[:, idx_band], dx=freq_res)
+    bandpower = bandpower / simpson(psd, dx=freq_res) if relative else bandpower
+    return bandpower
+
+
+# %%
+# Real-time estimation on a rolling window
+# ----------------------------------------
+#
+# Next, we can estimate the alpha band power on a rolling window of 4 seconds by running
+# an infinite loop that reads the data from the stream and computes the bandpower on the
+# last 4 seconds of data.
+#
+# .. note::
+#
+#     A chunk size of 200 samples is used to ensure stability in our documentation
+#     build, but in practice, a real-time application will likely publish new samples
+#     in smaller chunks and thus at a higher frequency. Due to the large chunk size,
+#     the acquisition delay of the connected stream is also increased to reduce the
+#     load on the CPU.
+
+source_id = uuid.uuid4().hex
+with PlayerLSL(raw, chunk_size=200, name="bandpower-example", source_id=source_id):
+    stream = StreamLSL(bufsize=4, name="bandpower-example", source_id=source_id)
+    stream.connect(acquisition_delay=0.1, processing_flags="all")
+    stream.pick("eeg").filter(1, 30)
+    stream.get_data()  # reset the number of new samples after the filter is applied
+
+    datapoints, times = [], []
+    while stream.n_new_samples < stream.n_buffer:
+        time.sleep(0.1)  # wait for the buffer to be entirely filled
+    while len(datapoints) != 30:
+        if stream.n_new_samples == 0:
+            continue  # wait for new samples
+        data, ts = stream.get_data()
+        bp = bandpower(data, stream.info["sfreq"], "periodogram", band=(8, 13))
+        datapoints.append(bp)
+        times.append(ts[-1])
+    stream.disconnect()
+
+# %%
+# Plot in function of time
+# ------------------------
+#
+# We can now plot the rolling-window bandpower in function of time, using the timestamps
+# of the last sample for each window on the X-axis. For simplicity, let's average all
+# channels together.
+
+f, ax = plt.subplots(1, 1, layout="constrained")
+ax.plot(times - times[0], [np.average(dp) * 100 for dp in datapoints])
+ax.set_xlabel("Time (s)")
+ax.set_ylabel("Relative α band power (%)")
+plt.show()
+
+# %%
+# Delay between 2 samples
+# -----------------------
+#
+# Let's also have a look at the delay between 2 data points, i.e. the overlap between
+# the windows.
+
+f, ax = plt.subplots(1, 1, layout="constrained")
+timedeltas = np.diff(times - times[0]) * 1000
+ax.hist(timedeltas, bins=15)
+ax.set_xlabel("Delay between 2 samples (ms)")
+plt.show()
+
+# %%
+# .. note::
+#
+#     Due to the low resources available on our CIs to build the documentation, some of
+#     those datapoints might have been computed with 2 acquisition window of delay
+#     instead of 1, yielding a delay between 2 samples of 2 acquisition windows instead
+#     of 1. In practice, with a large chunk size of 200 samples, we should get a delay
+#     between 2 computed time points to 200 samples, i.e. around 195.31 ms.
+#
+# Compare power spectral density estimation methods
+# -------------------------------------------------
+#
+# Let's compare both the bandpower estimation and the computation time of the different
+# methods to estimate the power spectral density.
+
+methods = ("periodogram", "welch", "multitaper")
+with PlayerLSL(raw, chunk_size=200, name="bandpower-example", source_id=source_id):
+    stream = StreamLSL(bufsize=4, name="bandpower-example", source_id=source_id)
+    stream.connect(acquisition_delay=0.1, processing_flags="all")
+    stream.pick("eeg").filter(1, 30)
+    stream.get_data()  # reset the number of new samples after the filter is applied
+
+    datapoints, times = {method: [] for method in methods}, []
+    while stream.n_new_samples < stream.n_buffer:
+        time.sleep(0.1)  # wait for the buffer to be entirely filled
+    while len(datapoints[methods[0]]) != 30:
+        if stream.n_new_samples == 0:
+            continue  # wait for new samples
+        data, ts = stream.get_data()
+        for method in methods:
+            bp = bandpower(data, stream.info["sfreq"], method, band=(8, 13))
+            datapoints[method].append(bp)
+        times.append(ts[-1])
+    stream.disconnect()
+
+f, ax = plt.subplots(1, 1, layout="constrained")
+for k, method in enumerate(methods):
+    ax.plot(
+        times - times[0],
+        [np.average(dp) * 100 for dp in datapoints[method]],
+        label=method,
+        color=colormaps["viridis"].colors[k * 60 + 20],
+    )
+ax.set_xlabel("Time (s)")
+ax.set_ylabel("Relative α band power (%)")
+ax.legend()
+plt.show()
+
+# %%
+# For the computation time and the overall loop execution speed, we need to run each
+# method on a separate loop.
+
+methods = ("periodogram", "welch", "multitaper")
+with PlayerLSL(raw, chunk_size=200, name="bandpower-example", source_id=source_id):
+    stream = StreamLSL(bufsize=4, name="bandpower-example", source_id=source_id)
+    stream.connect(acquisition_delay=0.1, processing_flags="all")
+    stream.pick("eeg").filter(1, 30)
+    stream.get_data()  # reset the number of new samples after the filter is applied
+
+    times = {method: [] for method in methods}
+    while stream.n_new_samples < stream.n_buffer:
+        time.sleep(0.1)  # wait for the buffer to be entirely filled
+    for k, method in enumerate(methods):
+        while len(times[methods[k]]) != 30:
+            if stream.n_new_samples == 0:
+                continue  # wait for new samples
+            data, ts = stream.get_data()
+            bp = bandpower(data, stream.info["sfreq"], method, band=(8, 13))
+            times[method].append(ts[-1])
+    stream.disconnect()
+
+timedeltas = {
+    method: np.diff(times[method] - times[method][0]) * 1000 for method in methods
+}
+timedeltas_average = {method: np.average(timedeltas[method]) for method in methods}
+
+for method in methods:
+    print(
+        f"Average delay between 2 samples for {method}: "
+        f"{timedeltas_average[method]:.2f} ms"
+    )
+
+# %%
+# .. note::
+#
+#     For this example, the average delay between 2 estimation of the bandpower is
+#     similar between all 3 methods because we are waiting for new samples which come in
+#     chunks of 200 samples, i.e. every 195.31 ms at the sampling frequency of 1024 Hz.
+#     The figure obtained for a chunk size of 1 sample and an acquisition delay of 1 ms
+#     is shown below.
+#
+#     .. code-block:: python
+#
+#         f, ax = plt.subplots(1, 1, layout="constrained")
+#         for k, method in enumerate(methods):
+#             ax.hist(
+#                 timedeltas[method],
+#                 bins=15,
+#                 label=method,
+#                 color=colormaps["viridis"].colors[k * 60 + 20],
+#             )
+#         ax.set_xlabel("Delay between 2 samples (ms)")
+#         ax.legend()
+#         plt.show()
+#
+# .. image:: ../../_static/tutorials/bp-performance.png
+#     :align: center

1.13.1/_downloads/15e9412ae22a63c550c15e151b3ffffc/20_bandpower.py

@@ -0,0 +1,107 @@
+"""
+Automatic vs Manual acquisition
+===============================
+
+.. include:: ./../../links.inc
+
+The :class:`~mne_lsl.stream.StreamLSL` object offers 2 mode of acquisition: automatic or
+manual. In automatic mode, the stream object acquires new chunks of samples at a
+regular interval. In manual mode, the user has to call the
+:meth:`~mne_lsl.stream.StreamLSL.acquire` to acquire new chunks of samples from the
+network. The automatic or manual acquisition is selected via the ``acquisition_delay``
+argument of :meth:`~mne_lsl.stream.StreamLSL.connect`:
+
+- a non-zero positive integer value will set the acquisition to automatic mode with the
+  specified delay in seconds.
+- ``0`` will set the acquisition to manual mode.
+
+Automatic acquisition
+---------------------
+
+When the stream is set to automatically acquire new samples at a regular interval, a
+background thread is created with :class:`concurrent.futures.ThreadPoolExecutor`. The
+background thread is periodically receives a job to acquire new samples from the
+network.
+
+.. important::
+
+    If the main thread is hogging all of the CPU resources, the delay between two
+    acquisition job might be longer than the specified delay. The background thread
+    will always do its best to acquire new samples at the specified delay, but it is not
+    able to do so if the CPU is busy.
+"""
+
+import uuid
+from time import sleep
+
+from matplotlib import pyplot as plt
+
+from mne_lsl.datasets import sample
+from mne_lsl.player import PlayerLSL
+from mne_lsl.stream import StreamLSL
+
+# create a mock LSL stream for this tutorial
+fname = sample.data_path() / "sample-ant-raw.fif"
+source_id = uuid.uuid4().hex
+player = PlayerLSL(fname, chunk_size=200, source_id=source_id).start()
+player.info
+
+# %%
+# .. note::
+#
+#     A ``chunk_size`` of 200 samples is used here to ensure stability and reliability
+#     while building the documentation on the CI. In practice, a ``chunk_size`` of 200
+#     samples is too large to represent a real-time application.
+
+stream = StreamLSL(bufsize=2, source_id=source_id).connect(acquisition_delay=0.1)
+sleep(2)  # wait for new samples
+print(f"New samples acquired: {stream.n_new_samples}")
+stream.disconnect()
+
+# %%
+# Manual acquisition
+# -------------------
+#
+# In manual acquisition mode, the user has to call the
+# :meth:`~mne_lsl.stream.StreamLSL.acquire` to get new samples from the network. In this
+# mode, all operation happens in the main thread and the user has full control over when
+# to acquire new samples.
+
+stream = StreamLSL(bufsize=2, source_id=source_id).connect(acquisition_delay=None)
+sleep(2)  # wait for new samples
+print(f"New samples acquired (before stream.acquire()): {stream.n_new_samples}")
+stream.acquire()
+print(f"New samples acquired (after stream.acquire()): {stream.n_new_samples}")
+
+# %%
+# However, it is also now up to the user to make sure he acquires new samples regularly
+# and does not miss part of the stream. The created :class:`~mne_lsl.lsl.StreamInlet`
+# has its buffer set to the same value as the :class:`~mne_lsl.stream.StreamLSL` object.
+
+stream.acquire()
+data1, ts1 = stream.get_data(picks="Cz")
+sleep(4)  # wait for 2 buffers
+stream.acquire()
+data2, ts2 = stream.get_data(picks="Cz")
+
+f, ax = plt.subplots(1, 1, layout="constrained")
+ax.plot(ts1 - ts1[0], data1.squeeze(), color="blue", label="acq 1")
+ax.plot(ts2 - ts1[0], data2.squeeze(), color="red", label="acq 2")
+ax.legend()
+ax.set_xlabel("Time (s)")
+ax.set_ylabel("EEG amplitude")
+plt.show()
+
+# %%
+# Free resources
+# --------------
+#
+# When you are done with a :class:`~mne_lsl.player.PlayerLSL` or
+# :class:`~mne_lsl.stream.StreamLSL`, don't forget to free the resources they both use
+# to continuously mock an LSL stream or receive new data from an LSL stream.
+
+stream.disconnect()
+
+# %%
+
+player.stop()