Notebooks/Day 2/2 - Filtering and Spectra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "\n",
    "import mne\n",
    "import numpy as np\n",
    "\n",
    "from _helper_functions import plot_psd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral analysis of data - the `time_frequency` module\n",
    "\n",
    "Often in neuroscience, we are interested in determining the spectral composition of signals, represented as a power spectral density (PSD).\n",
    "\n",
    "Furthermore, we may want to isolate the spectral content of a signal in a particular range of frequencies, which involves filtering the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1 - Filtering data\n",
    "\n",
    "To examine activity at a limited range of frequencies, we perform spectral filtering.\n",
    "\n",
    "This can take the form of:\n",
    "- Lowpass filtering - retaining information content of a signal below a certain frequency.\n",
    "- Highpass filtering - retaining information content of a signal above a certain frequency.\n",
    "- Bandpass filtering - retaining information content of a signal within a certain frequency range.\n",
    "- Bandstop filtering - retaining information content of a signal outside a certain frequency range.\n",
    "\n",
    "<br>\n",
    "<img src=\"../../_static/figures/filter_types.png\" alt=\"Filter types\" width=\"30%\" height=\"30%\">\n",
    "\n",
    "Credit: [allaboutcircuits.com](https://www.allaboutcircuits.com/technical-articles/low-pass-filter-tutorial-basics-passive-RC-filter/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an understanding of how to perform such filtering in MNE, we start by simulating 5 seconds of data sampled at 200 Hz, consisting of sine waves at 5 Hz, 10 Hz, and 20 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation settings\n",
    "duration = 10  # seconds\n",
    "sfreq = 200  # sampling rate (Hz)\n",
    "\n",
    "# Timepoints of the simulated data\n",
    "times = np.linspace(start=0, stop=duration, num=sfreq * duration, endpoint=False)\n",
    "\n",
    "# Simulate data as sine waves of given frequencies\n",
    "chan_1 = np.sin(2 * np.pi * times * 5)  # 5 Hz signal\n",
    "chan_2 = np.sin(2 * np.pi * times * 10)  # 10 Hz signal\n",
    "chan_3 = np.sin(2 * np.pi * times * 20)  # 20 Hz signal\n",
    "\n",
    "# Combine channels into a single array\n",
    "data = np.array([chan_1, chan_2, chan_3])\n",
    "ch_names = [\"5Hz\", \"10Hz\", \"20Hz\"]  # channel names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To play around with the data, let us first store it in a `Raw` object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises - Spectral filtering**\n",
    "\n",
    "**Exercise:** Create an [`Info`](https://mne.tools/stable/generated/mne.Info.html) object for the 3 channels, setting the channel types to be EEG, and using the sampling frequency we specified above.\n",
    "\n",
    "Afterwards, use the `data` array and the `Info` object to create a [`RawArray`](https://mne.tools/stable/generated/mne.io.RawArray.html) object for the signals, called `raw`.\n",
    "\n",
    "*Hint:* use the [`create_info()`](https://mne.tools/stable/generated/mne.create_info.html) function to create the `Info` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now easily plot the data.\n",
    "\n",
    "Count the number of cycles in each channel per second. Do they match our expectations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The object containing the signals should be called `raw`\n",
    "raw.plot(duration=6, scalings=\"auto\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing power spectra - a brief introduction\n",
    "\n",
    "We can easily compute the power spectral density of channels in `Raw` objects.\n",
    "\n",
    "This is done by calling the [`compute_psd()`](https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.compute_psd) method.\n",
    "\n",
    "`compute_psd()` returns a [`mne.time_frequency.Spectrum`](https://mne.tools/stable/generated/mne.time_frequency.Spectrum.html#mne.time_frequency.Spectrum) object containing the power spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute PSD of the data\n",
    "spectrum = raw.compute_psd(fmax=30)\n",
    "spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the PSD using the [`plot()`](https://mne.tools/stable/generated/mne.time_frequency.Spectrum.html#mne.time_frequency.Spectrum.plot) method of the `Spectrum` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the PSD\n",
    "spectrum.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, there are distinct peaks in the power spectrum at 5, 10, and 20 Hz.\n",
    "\n",
    "We will explore PSD computation in more detail below after examining spectral filtering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spectral filtering\n",
    "\n",
    "Let us now look at how we can filter the data.\n",
    "\n",
    "##### Lowpass, highpass, bandpass, and bandstop filtering\n",
    "\n",
    "Lowpass, highpass, bandpass, and bandstop filtering is most easily done using the [`filter()`](https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.filter) method of `Raw` objects.\n",
    "\n",
    "Frequencies to filter using the `filter()` method are specified using the `l_freq` and `h_freq` parameters:\n",
    "- `l_freq` specifies the lowest frequency of information to retain (in Hz).\n",
    "- `h_freq` specifies the highest frequency of information to retain (in Hz).\n",
    "\n",
    "<br>\n",
    "\n",
    "In this way:\n",
    "- specifying only `l_freq` highpass filters the data.\n",
    "- specifying only `h_freq` lowpass filters the data.\n",
    "- specifying `l_freq` to be lower than `h_freq` bandpass filters the data.\n",
    "- specifying `l_freq` to be higher than `h_freq` bandstop filters the data.\n",
    "\n",
    "<br>\n",
    "<img src=\"../../_static/figures/filter_types_marked.png\" alt=\"Filter types with frequency parameters\" width=\"30%\" height=\"30%\">\n",
    "\n",
    "Adapted from: [allaboutcircuits.com](https://www.allaboutcircuits.com/technical-articles/low-pass-filter-tutorial-basics-passive-RC-filter/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example below shows how to highpass filter the data to remove the 5 Hz activity.\n",
    "\n",
    "**N.B.** Note that the `Raw` object is copied so that the original data is not modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy to preserve original data\n",
    "raw_copy = raw.copy()\n",
    "\n",
    "# Highpass filter at 8 Hz to exclude 5 Hz activity\n",
    "raw_copy.filter(l_freq=8, h_freq=None)\n",
    "\n",
    "# Compute the PSD of the new data and plot it\n",
    "raw_copy.compute_psd(fmax=30).plot()\n",
    "\n",
    "# View the filtered data\n",
    "raw_copy.plot(duration=6, scalings=\"auto\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the activity in the 5 Hz channel has been removed, along with the peak in the power spectrum at 5 Hz.\n",
    "\n",
    "On the other hand, the 10 and 20 Hz activity remains."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises - Lowpass, highpass, bandpass, and bandstop filtering**\n",
    "\n",
    "**Exercise:** Lowpass filter the data at 15 Hz to remove the 20 Hz activity, then plot the PSD and raw data to confirm that only the 20 Hz activity has been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Lowpass filter the data at 8 Hz to remove the 10 and 20 Hz activity, then plot the PSD and raw data to confirm that only the 5 Hz activity remains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Bandstop filter the data between 8 and 15 Hz to remove the 10 Hz activity, then plot the PSD and raw data to confirm that only the 10 Hz activity has been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Bandpass filter the data between 8 and 15 Hz to remove the 5 and 20 Hz activity, then plot the PSD and raw data to confirm that only the 10 Hz activity remains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the filters from the `filter()` method can operate over a large frequency range.\n",
    "\n",
    "Sometimes however, the attenuation of a very limited frequency range is desired, for example when removing line noise artefacts, see e.g.:\n",
    "- https://pressrelease.brainproducts.com/eeg-artifacts-handling-in-analyzer/#technical\n",
    "- https://labeling.ucsd.edu/tutorial/labels\n",
    "- https://mne.tools/stable/auto_tutorials/preprocessing/30_filtering_resampling.html#power-line-noise\n",
    "\n",
    "In these situations, a notch filter is often used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Notch filtering\n",
    "\n",
    "Notch filters have their own dedicated [`notch_filter()`](https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.notch_filter) method in the `Raw` object.\n",
    "\n",
    "Below, we use a notch filter to remove the 5 Hz activity alone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_copy = raw.copy()\n",
    "\n",
    "# Apply notch filter at 5 Hz\n",
    "raw_copy.notch_filter(freqs=5)\n",
    "\n",
    "# Plot the PSD and timeseries of the filtered data\n",
    "raw_copy.compute_psd(fmax=30).plot()\n",
    "raw_copy.plot(duration=6, scalings=\"auto\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises - Notch filtering**\n",
    "\n",
    "**Exercise:** Notch filter the data at 10 Hz, and visualise the PSD and raw data to confirm that only the 10 Hz activity has been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Notch filter the data at 20 Hz, and visualise the PSD and raw data to confirm that only the 20 Hz activity has been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Notch filter the data at 10 and 20 Hz in a single call to the `notch_filter()` method, and visualise the PSD and raw data to confirm that only the 5 Hz activity remains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, MNE provides a number of convenient tools for the spectral filtering of data.\n",
    "\n",
    "There are many options for specifying filter parameters to fine tune the filters for your needs which are discussed in more depth in the following tutorials:\n",
    "- Background on filtering: https://mne.tools/stable/auto_tutorials/preprocessing/25_background_filtering.html\n",
    "- Filtering and resampling: https://mne.tools/stable/auto_tutorials/preprocessing/30_filtering_resampling.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2 - Computing power spectral densities\n",
    "\n",
    "Up until now, we have computed PSDs using the `compute_psd()` method of `Raw` objects.\n",
    "\n",
    "Note that an equivalent method exists for `Epochs` objects: [`mne.Epochs.compute_psd()`](https://mne.tools/stable/generated/mne.Epochs.html#mne.Epochs.compute_psd).\n",
    "\n",
    "The `compute_psd()` methods of `Raw` and `Epochs` objects support PSD computations using the Welch and multitaper methods.\n",
    "\n",
    "There exist equivalent functions for computing PSDs using the Welch and multitaper methods from arrays of data:\n",
    "- [`mne.time_frequency.psd_array_welch()`](https://mne.tools/stable/generated/mne.time_frequency.psd_array_welch.html)\n",
    "- [`mne.time_frequency.psd_array_multitaper()`](https://mne.tools/stable/generated/mne.time_frequency.psd_array_multitaper.html)\n",
    "\n",
    "**N.B.** Performing the PSD computations on arrays requires the sampling frequency of the data (`sfreq` parameter) to be specified."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, using the sample data, we specify the multitaper method to use in `compute_psd()`.\n",
    "\n",
    "Using the `fmax` parameter, we only return the results until 50 Hz.\n",
    "\n",
    "We additionally take only the EEG channels and crop to the first 60 seconds to reduce computation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the sample data\n",
    "raw = mne.io.read_raw_fif(\n",
    "    os.path.join(mne.datasets.sample.data_path(), \"MEG\", \"sample\", \"sample_audvis_raw.fif\")\n",
    ")\n",
    "raw.pick(picks=\"eeg\", exclude=\"bads\")\n",
    "raw.crop(tmax=60)\n",
    "raw.load_data()\n",
    "\n",
    "# Compute PSD\n",
    "spectrum = raw.compute_psd(method=\"welch\", fmax=50, n_fft=2048)\n",
    "\n",
    "# Plot the PSD\n",
    "spectrum.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can extract the power values and the corresponding frequencies from the `Spectrum` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract PSD data\n",
    "psd = spectrum.get_data()\n",
    "\n",
    "# Extract frequencies in the PSD\n",
    "freqs = spectrum.freqs\n",
    "\n",
    "print(f\"PSD data has shape: {psd.shape}  # channels x frequencies\")\n",
    "print(f\"Frequencies has shape: {freqs.shape}  # frequencies\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a custom function `plot_psd()`, we can verify that these values match those plotted using the `plot()` method of the `Spectrum` object.\n",
    "\n",
    "For the `plot_psd()` function, we pass in the array of power spectral density values for a set of channels, alongisde the corresponding frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot PSD from arrays with custom function\n",
    "plot_psd(psd=psd, freqs=freqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing PSDs from standalone functions\n",
    "\n",
    "**Exercises - Computing PSDs**\n",
    "\n",
    "**Exercise:** Perform the equivalent computation using the `psd_array_welch()` function on the data array extracted from the `Raw` object,\n",
    "\n",
    "Remember to specify a maximum frequency of 50 Hz and an FFT length of 2,048.\n",
    "\n",
    "Use the custom `plot_psd()` function to visualise the results.\n",
    "\n",
    "Do the results match the output of `compute_psd()`?\n",
    "\n",
    "*Hint:* data can be extracted from `Raw` objects using the [`get_data()`](https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw.get_data) method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Again using `psd_array_welch()`, compute the PSDs for the frequency range from 5 Hz onwards (i.e. no 50 Hz limit), and visualise the results.\n",
    "\n",
    "*Hint:* use the `fmin` parameter to specify the starting frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Using `psd_array_welch()`, compute the PSDs for the frequency range 5 - 50 Hz, and visualise the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Indexing frequencies of results\n",
    "\n",
    "**Exercise:** Using `psd_array_welch()`, compute the PSDs for the entire frequency range, but only plot the results up to 50 Hz.\n",
    "\n",
    "*Hint:* Use the [`np.where()`](https://numpy.org/doc/stable/reference/generated/numpy.where.html) function to find where the appropriate frequency values are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Using `psd_array_welch()`, compute the PSDs for the entire frequency range, but only plot the results from 5 Hz onwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Using `psd_array_welch()`, compute the PSDs for the entire frequency range, but only plot the results from 5 - 50 Hz.\n",
    "\n",
    "*Hint:* Use the form `np.where((condition1) & (condition2))` when you want to index an array based on multiple conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CODE GOES HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Summary of PSD computation\n",
    "\n",
    "As you can see, the `compute_psd()` methods of `Raw` and `Epochs` objects are very convenient ways of computing PSDs, with equivalent standalone functions for computations on arrays.\n",
    "\n",
    "However, MNE also offers tools for more advanced time-frequency analyses based on epoched data. These include time-frequency representations ([TFRs](https://mne.tools/stable/documentation/glossary.html#term-tfr)) based on the multitaper, Morlet wavelet, or Stockwell transformation methods:\n",
    "- Multitaper:\n",
    "    - From `Epochs` objects: [`mne.time_frequency.tfr_multitaper()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_multitaper.html)\n",
    "    - From arrays: [`mne.time_frequency.tfr_array_multitaper()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_array_multitaper.html)\n",
    "- Morlet wavelet:\n",
    "    - From `Epochs` objects: [`mne.time_frequency.tfr_morlet()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_morlet.html)\n",
    "    - From arrays: [`mne.time_frequency.tfr_array_morlet()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_array_morlet.html)\n",
    "- Stockwell transformation:\n",
    "    - From `Epochs` objects: [`mne.time_frequency.tfr_stockwell()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_stockwell.html)\n",
    "    - From arrays: [`mne.time_frequency.tfr_array_stockwell()`](https://mne.tools/stable/generated/mne.time_frequency.tfr_array_stockwell.html)\n",
    "\n",
    "Time-frequency analyses are discussed in more detail here: https://mne.tools/stable/auto_tutorials/time-freq/20_sensors_time_frequency.html#time-frequency-analysis-power-and-inter-trial-coherence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3 - Spectral filtering to remove artefacts\n",
    "\n",
    "Spectral filtering is not only useful for isolating activity at some frequencies of interest, but it can also be used to remove artefacts from the data.\n",
    "\n",
    "The ability to remove technical artefacts was previously mentioned in the context of notch filtering line noise, but biological artefacts such as cardiac activity can also be identified using spectral filtering, and subsequently removed (see e.g. https://labeling.ucsd.edu/tutorial/labels)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cardiac artefacts can be clearly seen in the MEG channels of MNE's sample data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the sample data\n",
    "raw = mne.io.read_raw_fif(\n",
    "    os.path.join(mne.datasets.sample.data_path(), \"MEG\", \"sample\", \"sample_audvis_raw.fif\")\n",
    ")\n",
    "raw.del_proj()  # delete existing PCA projections\n",
    "\n",
    "# Pick some channels with strong artefacts and plot them\n",
    "artefact_picks = [152, 155, 158, 164, 167, 170, 272, 275, 278, 284, 287, 290]\n",
    "raw.plot(order=artefact_picks);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If our desire is to analyse neural data, not removing these non-neural artefacts could of course lead to erroneous conclusions.\n",
    "\n",
    "Thankfully, MNE has a convenient function for doing just that: [`mne.preprocessing.compute_proj_ecg()`](https://mne.tools/stable/generated/mne.preprocessing.compute_proj_ecg.html).\n",
    "\n",
    "`compute_proj_ecg()` involves:\n",
    "- Filtering the data within a given frequency range to isolate the cardiac activity.\n",
    "- Finding the peaks of cardiac activity.\n",
    "- Creating epochs around these peaks of activity.\n",
    "- Using these epochs to create [projection vectors](https://mne.tools/stable/documentation/glossary.html#term-projector) that can be used to minimise the cardiac artefacts in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find projections to minimise cardiac artefacts\n",
    "projs, _ = mne.preprocessing.compute_proj_ecg(raw=raw)\n",
    "\n",
    "# Apply projections to the data and plot the cleaned data\n",
    "raw.add_proj(projs=projs)\n",
    "raw.plot(order=artefact_picks);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An equivalent function exists for removing eye movement artefacts: [`mne.preprocessing.compute_proj_eog()`](https://mne.tools/stable/generated/mne.preprocessing.compute_proj_eog.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "Spectral filtering is an important part of many analyses in neuroscience, involving e.g. the extraction of activity at specific frequencies of interest and the removal of artefacts. The `filter()` and `notch_filter()` methods of `Raw` and `Epochs` objects provide convenient ways of performing such filtering, with equivalent standalone functions for working with arrays of data.\n",
    "\n",
    "Spectral activity can be represented as PSDs, which can be computed using the `compute_psd()` methods of `Raw` and `Epochs` objects, or the equivalent standalone functions for computations on arrays. More advanced spectral analyses are also offered in the form of TFR computations, with MNE's [Time-Frequency module](https://mne.tools/stable/api/time_frequency.html) also offering several other useful tools, such as for computing cross-spectral densities (CSDs)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional resources\n",
    "\n",
    "MNE tutorial on spectral analysis: https://mne.tools/stable/auto_tutorials/time-freq/20_sensors_time_frequency.html\n",
    "\n",
    "MNE tutorial on `Spectrum` and `EpochsSpectrum` classes: https://mne.tools/stable/auto_tutorials/time-freq/10_spectrum_class.html\n",
    "\n",
    "Video introducing the Fourier transform with some very nice visualisations:\n",
    "https://youtu.be/spUNpyF58BY?si=hUC2zG8dG6Zah8tP"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mne_course",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}