Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions doc/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@
"logo": {
"image_light": "_static/logo.png",
"image_dark": "_static/logo-dark.png",
"text": "meegkit",
"text": f"meegkit v{version}",
},
"show_toc_level": 1,
"external_links": [
Expand Down Expand Up @@ -126,7 +126,7 @@
"examples_dirs": "../examples", # path to your example scripts
"gallery_dirs": "auto_examples", # path to where to save gallery generated output
"filename_pattern": "/example_",
"ignore_pattern": "config.py",
"ignore_pattern": r"(config|run_all_notebooks)\.py",
"run_stale_examples": False,
}

Expand Down
34 changes: 32 additions & 2 deletions doc/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,31 @@ in mind that this is mostly development code, and as such is likely to change
without any notice. Also, while most of the methods have been fairly robustly
tested, bugs can (and should!) be expected.

The package is most useful for readers who want practical reference
implementations of denoising and component-analysis methods, together with
worked examples that show how to interpret the outputs.

The source code of the project is hosted on Github at the following address:
https://github.com/nbara/python-meegkit

To get started, follow the installation instructions `in the README <https://github.com/nbara/python-meegkit#installation>`_.
Quick start
-----------

Install the package with ``pip``:

.. code-block:: bash

pip install meegkit

Some ASR-related functionality requires optional dependencies. To install those
as well, use:

.. code-block:: bash

pip install 'meegkit[extra]'

For development, documentation building, or testing, see the fuller
installation guidance `in the README <https://github.com/nbara/python-meegkit#installation>`_.

Available modules
-----------------
Expand Down Expand Up @@ -49,7 +70,16 @@ Here is a list of the methods and techniques available in ``meegkit``:
Examples gallery
----------------

A number of example scripts and notebooks are available:
A number of example scripts and notebooks are available.

If you are new to the package, a good starting sequence is:

1. ``example_asr`` for a full artifact-removal workflow.
2. ``example_dss`` for a simple synthetic component-recovery example.
3. ``example_trca`` or ``example_ress`` for task-oriented spatial filtering.

Many examples are synthetic sanity checks with known ground truth, which makes
them useful for understanding what each method is expected to recover.


.. toctree::
Expand Down
17 changes: 17 additions & 0 deletions examples/README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -4,5 +4,22 @@ Examples gallery
Below is a list of example scripts showing basic usage for most methods in
MEEGkit.

The examples are intentionally mixed between two roles:

1. synthetic sanity checks, where the ground truth is known and the method can
be judged directly,
2. workflow-style demonstrations, where the goal is to show how to interpret
outputs on more realistic data.

Suggested starting points:

1. ``example_asr`` for artifact-subspace reconstruction.
2. ``example_dss`` for a compact repeated-signal example.
3. ``example_dss_line`` for line-noise removal.
4. ``example_trca`` for block-wise classification performance.

Examples such as ``example_trca`` and the notebook-executed gallery pages may
take longer to run than the smaller synthetic demos.

The examples can be browsed directly on `Github <https://github.com/nbara/python-meegkit/tree/master/examples>`_.

312 changes: 68 additions & 244 deletions examples/example_asr.ipynb

Large diffs are not rendered by default.

57 changes: 52 additions & 5 deletions examples/example_asr.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,26 @@
ASR example
===========

Denoise data using Artifact Subspace Reconstruction.
This example demonstrates a full ASR workflow on short EEG data:

1. Calibrate ASR on a mostly clean segment.
2. Apply ASR in 1-second windows.
3. Compare raw and cleaned traces and quantify amplitude reduction.

This is intended as a first-pass inspection example rather than a benchmark:
it shows how the calibration choice propagates to the cleaned output.

The most useful outputs are the calibration retention mask and the channel-wise
RMS attenuation summary.

Uses meegkit.ASR().

References
----------
.. [1] Mullen, T., Kothe, C., Chi, Y., Ojeda, A., Kerth, T., Makeig, S.,
Jung, T. P., & Cauwenberghs, G. (2015). Real-time neuroimaging and
cognitive monitoring using wearable dry EEG. IEEE Transactions on
Biomedical Engineering, 62(11), 2553-2567.
"""
import os

Expand All @@ -21,11 +38,15 @@
###############################################################################
# Calibration and processing
# -----------------------------------------------------------------------------
# We use the first 30 seconds as a calibration segment. In practice, this
# segment should be as artifact-free as possible because ASR thresholds are
# derived from it.

# Train on a clean portion of data
asr = ASR(method="euclid")
train_idx = np.arange(0 * sfreq, 30 * sfreq, dtype=int)
_, sample_mask = asr.fit(raw[:, train_idx])
selected_fraction = np.mean(sample_mask)

# Apply filter using sliding (non-overlapping) windows
X = sliding_window(raw, window=int(sfreq), step=int(sfreq))
Expand All @@ -36,16 +57,26 @@
raw = X.reshape(8, -1) # reshape to (n_chans, n_times)
clean = Y.reshape(8, -1)

# A simple quality metric: root-mean-square attenuation per channel.
rms_before = np.sqrt(np.mean(raw ** 2, axis=1))
rms_after = np.sqrt(np.mean(clean ** 2, axis=1))
rms_ratio = rms_after / np.maximum(rms_before, np.finfo(float).eps)

###############################################################################
# Plot the results
# -----------------------------------------------------------------------------
#
# Data was trained on a 40s window from 5s to 45s onwards (gray filled area).
# The algorithm then removes portions of this data with high amplitude
# artifacts before running the calibration (hatched area = good).
# The gray overlay marks the 30-second calibration region actually used by the
# code. The hatched overlay shows the subset of that region that ASR kept while
# estimating its clean-data statistics.
#
# What to look for:
# - After ASR, sharp bursts should be attenuated in many channels.
# - The RMS ratio (after/before) should generally be below 1.
# - Strong attenuation everywhere would suggest over-aggressive calibration.

times = np.arange(raw.shape[-1]) / sfreq
f, ax = plt.subplots(8, sharex=True, figsize=(8, 5))
f, ax = plt.subplots(8, sharex=True, figsize=(9, 6))
for i in range(8):
ax[i].fill_between(train_idx / sfreq, 0, 1, color="grey", alpha=.3,
transform=ax[i].get_xaxis_transform(),
Expand All @@ -59,8 +90,24 @@
ax[i].set_ylim([-50, 50])
ax[i].set_ylabel(f"ch{i}")
ax[i].set_yticks([])
ax[0].set_title("Raw and cleaned EEG traces")
ax[i].set_xlabel("Time (s)")
ax[0].legend(fontsize="small", bbox_to_anchor=(1.04, 1), borderaxespad=0)
plt.subplots_adjust(hspace=0, right=0.75)
plt.suptitle("Before/after ASR")

fig, axm = plt.subplots(1, 1, figsize=(7, 3))
axm.bar(np.arange(raw.shape[0]), rms_ratio)
axm.axhline(1.0, color="k", ls=":", lw=1)
axm.set_xlabel("Channel")
axm.set_ylabel("RMS ratio (after / before)")
axm.set_title("Channel-wise attenuation summary")
axm.set_xticks(np.arange(raw.shape[0]))
axm.grid(True, axis="y", ls=":", alpha=.4)
plt.tight_layout()

print(f"Median RMS ratio across channels: {np.median(rms_ratio):.3f}")
print(f"Fraction of calibration samples retained: {selected_fraction:.3f}")
print("Interpretation: if only a small fraction of the calibration window is")
print("retained, the chosen segment may not be clean enough for stable ASR.")
plt.show()
97 changes: 97 additions & 0 deletions examples/example_asr_calibration_sensitivity.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# ASR calibration sensitivity\n\nThis tutorial-style example illustrates how ASR results depend on the\ncalibration segment used during fitting.\n\nWe fit ASR on three candidate calibration windows and compare the resulting\nchannel-wise RMS attenuation. This helps users choose a calibration strategy in\nreal datasets.\n\nThe main question is whether the cleaning result is stable across plausible\ncalibration windows, or whether one window yields much stronger attenuation.\n\nUses `meegkit.asr.ASR()`.\n\n## References\n.. [1] Mullen, T., Kothe, C., Chi, Y., Ojeda, A., Kerth, T., Makeig, S.,\n Jung, T. P., & Cauwenberghs, G. (2015). Real-time neuroimaging and\n cognitive monitoring using wearable dry EEG. IEEE Transactions on\n Biomedical Engineering, 62(11), 2553-2567.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom meegkit.asr import ASR\nfrom meegkit.utils.matrix import sliding_window\n\nraw = np.load(os.path.join(\"..\", \"tests\", \"data\", \"eeg_raw.npy\"))\nsfreq = 250"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define candidate calibration windows\nIn practice, calibration quality strongly affects how aggressively ASR\nsuppresses artifacts. Here we compare three 20-second windows.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"windows_sec = {\n \"early (0-20s)\": (0, 20),\n \"middle (10-30s)\": (10, 30),\n \"late (20-40s)\": (20, 40),\n}\n\nX = sliding_window(raw, window=int(sfreq), step=int(sfreq))\nraw_win = X.reshape(raw.shape[0], -1)\nrms_before = np.sqrt(np.mean(raw_win ** 2, axis=1))\n\nratios = {}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit and apply ASR for each candidate window\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for label, (start_s, stop_s) in windows_sec.items():\n train_idx = np.arange(start_s * sfreq, stop_s * sfreq, dtype=int)\n\n asr = ASR(method=\"euclid\")\n asr.fit(raw[:, train_idx])\n\n Y = np.zeros_like(X)\n for i in range(X.shape[1]):\n Y[:, i, :] = asr.transform(X[:, i, :])\n\n clean = Y.reshape(raw.shape[0], -1)\n rms_after = np.sqrt(np.mean(clean ** 2, axis=1))\n ratios[label] = rms_after / np.maximum(rms_before, np.finfo(float).eps)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualize channel-wise attenuation by calibration choice\nWhat to look for:\n- Ratios below 1 indicate attenuation.\n- Large differences across windows suggest sensitivity to calibration period.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1, figsize=(8, 4))\nchannels = np.arange(raw.shape[0])\n\nfor label, ratio in ratios.items():\n ax.plot(channels, ratio, \"o-\", lw=1.5, label=label)\n\nax.axhline(1.0, color=\"k\", ls=\":\", lw=1)\nax.set_xticks(channels)\nax.set_xlabel(\"Channel\")\nax.set_ylabel(\"RMS ratio (after / before)\")\nax.set_title(\"ASR sensitivity to calibration window\")\nax.grid(True, axis=\"y\", ls=\":\", alpha=.4)\nax.legend(fontsize=\"small\")\nplt.tight_layout()\n\nfor label, ratio in ratios.items():\n print(f\"{label:>16}: median ratio = {np.median(ratio):.3f}\")\n\nbest_label = min(ratios, key=lambda key: np.median(ratios[key]))\nprint(f\"Most aggressive calibration in this example: {best_label}\")\nprint(\"Interpretation: large separation between curves means the ASR result\")\nprint(\"depends strongly on which segment is treated as clean calibration data.\")\n\nplt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.11"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
98 changes: 98 additions & 0 deletions examples/example_asr_calibration_sensitivity.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
"""
ASR calibration sensitivity
===========================

This tutorial-style example illustrates how ASR results depend on the
calibration segment used during fitting.

We fit ASR on three candidate calibration windows and compare the resulting
channel-wise RMS attenuation. This helps users choose a calibration strategy in
real datasets.

The main question is whether the cleaning result is stable across plausible
calibration windows, or whether one window yields much stronger attenuation.

Uses `meegkit.asr.ASR()`.

References
----------
.. [1] Mullen, T., Kothe, C., Chi, Y., Ojeda, A., Kerth, T., Makeig, S.,
Jung, T. P., & Cauwenberghs, G. (2015). Real-time neuroimaging and
cognitive monitoring using wearable dry EEG. IEEE Transactions on
Biomedical Engineering, 62(11), 2553-2567.
"""
import os

import matplotlib.pyplot as plt
import numpy as np

from meegkit.asr import ASR
from meegkit.utils.matrix import sliding_window

raw = np.load(os.path.join("..", "tests", "data", "eeg_raw.npy"))
sfreq = 250

###############################################################################
# Define candidate calibration windows
# -----------------------------------------------------------------------------
# In practice, calibration quality strongly affects how aggressively ASR
# suppresses artifacts. Here we compare three 20-second windows.
windows_sec = {
"early (0-20s)": (0, 20),
"middle (10-30s)": (10, 30),
"late (20-40s)": (20, 40),
}

X = sliding_window(raw, window=int(sfreq), step=int(sfreq))
raw_win = X.reshape(raw.shape[0], -1)
rms_before = np.sqrt(np.mean(raw_win ** 2, axis=1))

ratios = {}

###############################################################################
# Fit and apply ASR for each candidate window
# -----------------------------------------------------------------------------
for label, (start_s, stop_s) in windows_sec.items():
train_idx = np.arange(start_s * sfreq, stop_s * sfreq, dtype=int)

asr = ASR(method="euclid")
asr.fit(raw[:, train_idx])

Y = np.zeros_like(X)
for i in range(X.shape[1]):
Y[:, i, :] = asr.transform(X[:, i, :])

clean = Y.reshape(raw.shape[0], -1)
rms_after = np.sqrt(np.mean(clean ** 2, axis=1))
ratios[label] = rms_after / np.maximum(rms_before, np.finfo(float).eps)

###############################################################################
# Visualize channel-wise attenuation by calibration choice
# -----------------------------------------------------------------------------
# What to look for:
# - Ratios below 1 indicate attenuation.
# - Large differences across windows suggest sensitivity to calibration period.
fig, ax = plt.subplots(1, 1, figsize=(8, 4))
channels = np.arange(raw.shape[0])

for label, ratio in ratios.items():
ax.plot(channels, ratio, "o-", lw=1.5, label=label)

ax.axhline(1.0, color="k", ls=":", lw=1)
ax.set_xticks(channels)
ax.set_xlabel("Channel")
ax.set_ylabel("RMS ratio (after / before)")
ax.set_title("ASR sensitivity to calibration window")
ax.grid(True, axis="y", ls=":", alpha=.4)
ax.legend(fontsize="small")
plt.tight_layout()

for label, ratio in ratios.items():
print(f"{label:>16}: median ratio = {np.median(ratio):.3f}")

best_label = min(ratios, key=lambda key: np.median(ratios[key]))
print(f"Most aggressive calibration in this example: {best_label}")
print("Interpretation: large separation between curves means the ASR result")
print("depends strongly on which segment is treated as clean calibration data.")

plt.show()
Loading