{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Extracting and visualizing subject head movement\n\nContinuous head movement can be encoded during MEG recordings by use of\nHPI coils that continuously emit sinusoidal signals. These signals can then be\nextracted from the recording and used to estimate head position as a function\nof time. Here we show an example of how to do this, and how to visualize\nthe result.\n\n## HPI frequencies\n\nFirst let's load a short bit of raw data where the subject intentionally moved\ntheir head during the recording. Its power spectral density shows five peaks\n(most clearly visible in the gradiometers) corresponding to the HPI coil\nfrequencies, plus other peaks related to power line interference (60 Hz and\nharmonics).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Eric Larson <larson.eric.d@gmail.com>\n#          Richard H\u00f6chenberger <richard.hoechenberger@gmail.com>\n#          Daniel McCloy <dan@mccloy.info>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from os import path as op\n\nimport mne\n\ndata_path = op.join(mne.datasets.testing.data_path(verbose=True), \"SSS\")\nfname_raw = op.join(data_path, \"test_move_anon_raw.fif\")\nraw = mne.io.read_raw_fif(fname_raw, allow_maxshield=\"yes\").load_data()\nraw.compute_psd().plot(picks=\"data\", exclude=\"bads\", amplitude=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can use `mne.chpi.get_chpi_info` to retrieve the coil frequencies,\nthe index of the channel indicating when which coil was switched on, and the\nrespective \"event codes\" associated with each coil's activity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "chpi_freqs, ch_idx, chpi_codes = mne.chpi.get_chpi_info(info=raw.info)\nprint(f\"cHPI coil frequencies extracted from raw: {chpi_freqs} Hz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Estimating continuous head position\n\nFirst, let's extract the HPI coil amplitudes as a function of time:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "chpi_amplitudes = mne.chpi.compute_chpi_amplitudes(raw)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Second, let's compute time-varying HPI coil locations from these:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "chpi_locs = mne.chpi.compute_chpi_locs(raw.info, chpi_amplitudes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, compute head positions from the coil locations:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "head_pos = mne.chpi.compute_head_pos(raw.info, chpi_locs, verbose=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that these can then be written to disk or read from disk with\n:func:`mne.chpi.write_head_pos` and :func:`mne.chpi.read_head_pos`,\nrespectively.\n\n## Visualizing continuous head position\n\nWe can plot as traces, which is especially useful for long recordings:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mne.viz.plot_head_positions(head_pos, mode=\"traces\", totals=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Or we can visualize them as a continuous field (with the vectors pointing\nin the head-upward direction):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mne.viz.plot_head_positions(head_pos, mode=\"field\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These head positions can then be used with\n:func:`mne.preprocessing.maxwell_filter` to compensate for movement,\nor with :func:`mne.preprocessing.annotate_movement` to mark segments as\nbad that deviate too much from the average head position.\n\n\n## Computing SNR of the HPI signal\n\nIt is also possible to compute the SNR of the continuous HPI measurements.\nThis can be a useful proxy for head position along the vertical dimension,\ni.e., it can indicate the distance between the HPI coils and the MEG sensors.\nUsing `~mne.chpi.compute_chpi_snr`, the HPI power and SNR are computed\nseparately for each MEG sensor type and each HPI coil (frequency), along with\nthe residual power for each sensor type. The results can then be visualized\nwith `~mne.viz.plot_chpi_snr`. Here we'll just show a few seconds, for speed:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw.crop(tmin=5, tmax=10)\nsnr_dict = mne.chpi.compute_chpi_snr(raw)\nfig = mne.viz.plot_chpi_snr(snr_dict)"
      ]
    }
  ],
  "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.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}