Analyzing_Rhythms_Lab_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: Analyzing Rhythms Part 1\n",
    "project:\n",
    "  type: website\n",
    "format:\n",
    "  html:\n",
    "    code-fold: false\n",
    "    code-tools: true\n",
    "jupyter: python 3\n",
    "number-sections: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data and look at it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load modules we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import loadmat\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.signal import spectrogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = loadmat('Data/EEG-1.mat')    # Load the EEG data\n",
    "# Data available here\n",
    "# https://github.com/Mark-Kramer/BU-MA665-MA666/tree/master/Data\n",
    "#\n",
    "EEG  = data['EEG'].reshape(-1) # Extract the EEG variable\n",
    "t    = data['t'][0]            # ... and the t variablea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t,EEG)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('EEG');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q:** What is the sampling interval (dt)?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = \"SOMETHING\"\n",
    "print(dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q.** What is the sampling frequency (f0)?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = \"SOMETHING\"\n",
    "print(f0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q.** What is the number of points in the data (N)?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = \"SOMETHING\"\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q.** What is the total time of the observation (T)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = \"SOMETHING\"\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Compute the spectrum \"by hand\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Q. What is the Nyquist frequency and frequency resolution?\n",
    "fNQ = \"SOMETHING\"\n",
    "df  = \"SOMETHING\"\n",
    "\n",
    "# Q. What is the frequency axis?\n",
    "fj  = np.arange(0,fNQ,df)\n",
    "\n",
    "# Then, compute the Fourier transform \"by hand\".\n",
    "x = EEG\n",
    "X = np.ndarray(np.size(fj), complex);\n",
    "for j in range( np.size(fj) ):\n",
    "    X[j] = \"SOMETHING\"\n",
    "\n",
    "# And the spectrum,\n",
    "Sxx = \"SOMETHING\"\n",
    "\n",
    "# Plot it,\n",
    "plot(fj, Sxx.real)\n",
    "plt.xlim([0, 100])                          # Select frequency range\n",
    "plt.xlabel('Frequency [Hz]')                # Label the axes\n",
    "plt.ylabel('Power [$\\mu V^2$/Hz]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the power spectrum using the FFT function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x   = EEG\n",
    "X   = np.fft.fft(x)                        # Compute Fourier transform of x\n",
    "Sxx = \"SOMETHING\"                          # Compute spectrum\n",
    "Sxx = Sxx[0:int(N / 2)].real               # Ignore negative frequencies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the frequency axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df  = \"SOMETHING\"                  # Determine frequency resolution\n",
    "fNQ = \"SOMETHING\"                  # Determine Nyquist frequency\n",
    "faxis = np.arange(0,fNQ,df)        # Construct frequency axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the spectrum versus frequency\n",
    "\n",
    "**Q.** What do you see?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(faxis, Sxx)\n",
    "plt.xlim([0, 100])                          # Select frequency range\n",
    "plt.xlabel('Frequency [Hz]')                # Label the axes\n",
    "plt.ylabel('Power [$\\mu V^2$/Hz]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the spectrum versus frequency on a **decibel scale**.\n",
    "\n",
    "**Q.** Now what do you see?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(faxis, [???])   # Plot the spectrum in decibels.\n",
    "plt.xlim([0, 100])                               # Select the frequency range.\n",
    "plt.ylim([-60, 0])                               # Select the decibel range.\n",
    "plt.xlabel('Frequency [Hz]')                     # Label the axes.\n",
    "plt.ylabel('Power [dB]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the spectrum versus frequency on a **logarithmic frequency axis**.\n",
    "\n",
    "**Q.** And now what do you see?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.[???]                                            # Log-log scale\n",
    "plt.xlim([df, 100])                                  # Select frequency range\n",
    "plt.ylim([-60, 0])                                   # ... and the decibel range.\n",
    "plt.xlabel('Frequency [Hz]')                         # Label the axes.\n",
    "plt.ylabel('Power [dB]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply Hanning taper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the Hanning taper to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_tapered  = ???\n",
    "plt.figure()\n",
    "plt.plot(t,x)\n",
    "plt.plot(t,x_tapered);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the Hanning taper and look at the spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf_tapered  = \"SOMETHING\"                        # Compute Fourier transform of x.\n",
    "Sxx_tapered = \"SOMETHING\"                        # Compute the spectrum,\n",
    "Sxx_tapered = np.real(Sxx_tapered[:int(N / 2)])  # ... and ignore negative frequencies.\n",
    "\n",
    "plt.figure()\n",
    "plt.semilogx(faxis,10*np.log10(Sxx))         # Plot spectrum of untapered signal.  \n",
    "plt.semilogx(faxis,10*np.log10(Sxx_tapered)) # Plot spectrum vs tapered signal.\n",
    "plt.xlim([faxis[1], 100])                    # Select frequency range,\n",
    "plt.ylim([-70, 20])                          # ... and the power range.\n",
    "plt.xlabel('Frequency [Hz]')                 # Label the axes\n",
    "plt.ylabel('Power [$\\mu V^2$/Hz]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectrogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the spectrogram.\n",
    "\n",
    "Fs = 1 / dt               # Define the sampling frequency,\n",
    "interval = int(Fs)        # ... the interval size,\n",
    "overlap = int(Fs * 0.95)  # ... and the overlap intervals\n",
    "\n",
    "                          # Compute the spectrogram\n",
    "f0, t0, Sxx0 = spectrogram(\n",
    "    EEG,                  # Provide the signal,\n",
    "    fs=Fs,                # ... the sampling frequency,\n",
    "    nperseg=interval,     # ... the length of a segment,\n",
    "    noverlap=overlap)     # ... the number of samples to overlap,\n",
    "plt.pcolormesh(t0, f0, 10 * np.log10(Sxx0),\n",
    "               cmap='jet')# Plot the result\n",
    "plt.colorbar()            # ... with a color bar,\n",
    "plt.ylim([0, 70])             # ... set the frequency range,\n",
    "plt.xlabel('Time [s]')       # ... and label the axes\n",
    "plt.ylabel('Frequency [Hz]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}