Assignments/A1-SamplingTimeSeries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Name\n",
    "- PID\n",
    "- COGS118C - Assignment 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook has [30 + 3 bonus] points in total \n",
    "##### The number of points for each question is denoted by []. Make sure you've answered all the questions and that the point total add up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---\n",
    "# Lab 1 - Time Series, Sampling, and Epoched Analysis (ERPs)\n",
    "In this lab, we will cover the first stages of signal processing: sampling data. This includes digitization and sampling theorem. We will generate and plot some signals. \n",
    "\n",
    "Then, we'll perform our first kind of neural signal analysis: event-related potentials.\n",
    "\n",
    "Key concepts:\n",
    "- visualizing time-series\n",
    "- digitization/quantization\n",
    "- sampling\n",
    "- (more) indexing arrays\n",
    "- epoching\n",
    "- event-related potentials (ERPs): noise and averaging\n",
    "\n",
    "**Answers for questions requiring written responses can be entered in the cell immediately below the question, so that when you write your response, it doesn't screw up the formatting of the question.** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analog signals\n",
    "Real world signals are continuous in time and amplitude (up to quantum-level limits, anyway). These are referred to as **\"analog\"** signals (Google it). Soundwaves that we produce when we speak or when we play a violin, for example, are analog signals. \n",
    "\n",
    "Equivalently, there are \"analog devices\" that produce, receive, and/or operate on analog signals. These often involve \"analog\" circuits.\n",
    "\n",
    "# [1] Q1:\n",
    "[1] 1.1: Give 3 examples of analog devices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 1.1:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Digital signals\n",
    "People used to analyze signals using analog circuits. This is pretty hardcore, and requires extensive hands-on knowledge about circuitry. Once you want to analyze the signal on a \"digital\" computer, however, you have to \"digitize\" the signal. This requires an **\"analog-to-digital converter\"** or ADC for short. \n",
    "\n",
    "---\n",
    "A tangent (without delving too much into how a computer works): all modern computers operate with binary transistors, which use a combination of \"bits\" to represent all other types of information. In the analog world, there are an infinite number of number between 0 and 1, so there is a limit to how accurately we can represent small decimals (or really big numbers). Python uses [floating point](https://0.30000000000000004.com/). Everything you see on your screen, at the lowest level, is converted into a numerical **binary** representation, even strings (see [ASCII](https://www.cs.cmu.edu/~pattis/15-1XX/common/handouts/ascii.html) table, for example).\n",
    "\n",
    "---\n",
    "Anyway, to digitize an analog signal, you have to discretely sample, both in value (voltage, brightness, etc) and in time. The former is usually called **digitization or quantization**, while **sampling** usually refers to the latter. It's like drawing a grid over your continuous signals and interpolating its values only at where the grid crosses.\n",
    "![sampling](img/WvD_fig1_6.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Let's get into it\n",
    "Without further ado: let's load up some EEG signals and explore. But first, make the necessary python module imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import io # this submodule let's us load the signal we want\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scipy loads .mat file into a dictionary\n",
    "# the details are not crucial, we just have to unpack them into python variables\n",
    "EEG_data = io.loadmat('data/EEG_exp.mat', squeeze_me = True)\n",
    "\n",
    "# print all the variables that exist in the dictionary\n",
    "print(EEG_data.keys())\n",
    "\n",
    "# this contains the EEG data\n",
    "EEG = EEG_data['EEG']\n",
    "# this contains the sampling rate, in Hz (or samples/second)\n",
    "fs = EEG_data['fs']\n",
    "\n",
    "# let's plot the signal\n",
    "plt.figure(figsize=(15,3))\n",
    "plt.plot(EEG)\n",
    "# ALWAYS label your plot axes in this course (and ever)\n",
    "plt.xlabel('Sample Number')\n",
    "plt.ylabel('Voltage (uV)')\n",
    "\n",
    "# now let's zoom in to see more detail\n",
    "plt.figure(figsize=(15,3))\n",
    "plt.plot(EEG, '.') # '.' means plot the data points as individual dots without linking them\n",
    "plt.xlim([0,1000]) # this limits the x-axis shown\n",
    "plt.xlabel('Sample Number')\n",
    "plt.ylabel('Voltage (uV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3] Q2: Digitization\n",
    "As you can see above, the signal we loaded is already a digitally sampled time series (a little over 70,000 samples), represented by discrete points in the second plot. To study the effect of quantization, let's simulate what would happen if we further quantized the signal, with a (prehistoric) 4-bit ADC.\n",
    "\n",
    "[1] 2.1: How many possible values can a 4-bit ADC represent? Remember, this means that the ADC has 4 binary 'bits' that it can use, thus giving you a total of how many levels? Compute this number in code and store that value in the variable `num_levels` below.\n",
    "\n",
    "[1] 2.2: Let's say our ADC has a total range between -32uV to 32uV. What is the voltage resolution of our ADC then? In other words, what is the finest voltage difference our ADC can distinguish between two samples? Compute this number in code and store that value in the variable `delta_v` below.\n",
    "\n",
    "[1] 2.3: Run the next two cells, they should produce a graph where the orange trace looks very quantized (kind of square). This is not good, because then we cannot distinguish small fluctuations in our signals, which, as we will see later in the course, are very important. **Re-run** the next two cells, but experiment with different values for `num_bits`. Just based on visual inspection of the plot, what is the minimum number of bits that you would want your ADC to have in this case, assuming the blue trace is a faithful representation of your signal? There's no one right answer, but justify your response."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 2.3:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bits = 4\n",
    "min_v, max_v = -32,32\n",
    "\n",
    "num_levels = # _FILL_IN_YOUR_CODE_HERE \n",
    "delta_v = # _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the quantization vector, these are the new possible values that your signal can take\n",
    "ADC_levels = np.arange(min_v,max_v,delta_v)+delta_v/2\n",
    "\n",
    "# quantize the EEG signal with our crappy ADC with the function np.digitize\n",
    "# note that we have to scale the redigitized signal to its original units\n",
    "EEG_quant = np.digitize(EEG,bins=ADC_levels)*delta_v+min_v\n",
    "\n",
    "plt.figure(figsize=(15,4))\n",
    "plt.plot(EEG, label='Original EEG')\n",
    "plt.plot(EEG_quant, label='Quantized EEG', alpha=0.8)\n",
    "plt.xlim([0,1000]); plt.ylim([-15, 15]);\n",
    "plt.legend()\n",
    "plt.xlabel('Sample Number')\n",
    "plt.ylabel('Voltage (uV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Sample Number vs. Time\n",
    "Notice that in all the plots above, the x-axis is \"sample number\", which simply correponds to the position each value is in the array `EEG`. We want to create a corresponding time vector, which marks at what clock time each value is sampled at. \n",
    "\n",
    "Sometimes your data will include a time vector. But for the sake of this exercise, you are asked to create the time vector based on the information/variables you have."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [6] Q3: Sampling in Time\n",
    "[1] 3.1: Given the sampling rate, what is the sampling **period**? In other words, how much time elapses between each consecutive sample? Compute this number as a function of `fs` and store it in the variable `dt` below.\n",
    "\n",
    "[1] 3.2: How long in total is this signal, in absolute time? Compute and store this in the variable `T_exp` below.\n",
    "\n",
    "[1] 3.3: Construct the corresponding time vector for the EEG data, assuming that the first sample came at t=0 and evenly spaced samples at `dt`. Store that in the variable `t_EEG` below. Hint: check out the function `np.arange()`.\n",
    "\n",
    "[2] 3.4: Re-plot the signal as a line chart, but with the x-axis as time (using the time vector you created above), and zoom into the first 1 second of the data. **Take note to label your plots carefully, with units!**\n",
    "\n",
    "[1] 3.5: To simulate **downsampling** in time, plot every **10th** value of the EEG data by indexing the array (check Google/StackExchange for how to do this). Remember, this applies both to the time vector and your EEG data. **Make sure to label your data and display the legend as Q2 above.**\n",
    "\n",
    "[BONUS: 1] 3.6: Sometimes it's useful to downsample your signal in time to conserve memory. As we did above, by taking every 10th value in our data, we essentially reduce the data size 10-fold. However, this is **NOT** the entirely right way to downsample your data. What issue do we introduce when we simply do that? (Hint: the answer can be as short as one word, and Google is your friend here.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 3.6:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dt = # _FILL_IN_YOUR_CODE_HERE \n",
    "T_exp = # _FILL_IN_YOUR_CODE_HERE \n",
    "t_EEG = # _FILL_IN_YOUR_CODE_HERE \n",
    "\n",
    "# Plotting the signal and its downsampled version\n",
    "plt.figure(figsize=(15,3))\n",
    "\n",
    "plt.plot(t_EEG, EEG, label='EEG')\n",
    "plt.plot(_FILL_IN_YOUR_CODE_HERE, _FILL_IN_YOUR_CODE_HERE, '.-', label=_FILL_IN_YOUR_CODE_HERE)\n",
    "plt.xlim([0,1]); plt.ylim([-15, 15]);\n",
    "plt.legend()\n",
    "# _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Event-Related Analysis\n",
    "The above data actually comes from an event-style EEG experiment. The participant is shown visual stimuli at regular intervals, aimed to trigger a reliable brain response for each type of stimuli (cat vs. dog pics, for example). This is a very common type of study design in neuroscience (and psychology). \n",
    "\n",
    "In this case, we will need to know when a stimulus was presented, and what type of stimulus it was. This information is stored in the variable `trial_info`, where the **first column has the stimulus onset time (in seconds), and the second column has the type of stimulus shown (1,2, or 3).** These are often extra streams of data sent through the \"trigger channel\" by the stimulus-presenting computer directly to the recording equipment, in order to synchronize with the EEG data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "trial_info = EEG_data['trial_info']\n",
    "\n",
    "# print the first 10 events\n",
    "print(trial_info[:10,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Process for Analyzing Event-Related Data\n",
    "These types of experiments follow a pretty standard analysis process. \n",
    "1. Import and pre-process your data (already done; we'll skip the pre-processing for now)\n",
    "2. Given the stimulus presentation timestamps (first column of `trial_info` above), find the corresponding indices in your EEG data by matching to the `t_EEG` time vector.\n",
    "3. Cut out an **epoch** (window of data) around the stimulus presentation time, which usually includes:\n",
    "    - pre-stimulus baseline (~0.5 seconds before stimulus presentation)\n",
    "    - stimulus presentation (t = 0)\n",
    "    - stimulus-driven response (or event-related response, 0-1 second after stimulus presentation)\n",
    "4. Baseline subtraction: subtract each epoch by its mean pre-stimulus value to account for any slow drifts over time.\n",
    "5. Group epochs based on stimulus type, and average epochs of the same type.\n",
    "6. Plot the average response (s)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Q4: Step 2 - Find Matching Timestamps in EEG Data\n",
    "Given the event times in `trial_info`, which we will assume to be the stimulus onset time for this experiment, we have to find the corresponding timestamp in the EEG data. Note that the timestamps may not always match exactly, as they could have different sampling rates. In those cases, you will have to settle for finding the **closest** timestamps. Currently, however, life was made easy for us by virtue of the fact that the EEG data (and timestamps) and the stimulus event timestamps are synchronously sampled at 1000Hz.\n",
    "\n",
    "In this case, we can directly convert the event timestamp into an integer index, since we know the sampling frequency and starting time. \n",
    "\n",
    "[1] 4.1: If the EEG timestamp starts at `t=0`, which is indexed by `i=0`, and is sampled at `fs=1000`, at which index will the EEG timestamp be equal to **3.050 seconds**? Compute and store this in the variable `trial_index` below. Note that to index an array, the number has to be an integer, which I've converted for you. (You will notice that the value is *a LITTLE* off. That's a precision issue and We can ignore that for now.)\n",
    "\n",
    "[3] 4.2: Following this logic, write a function that will find the corresponding index in the EEG data/timestamp for every event timestamp. Return that as an array of integers (`my_arr.astype(int)` will convert an array to all integers). You may use a for loop, list comprehension, or a simple (one-line) array calculation for this. Confirm that the timestamps match what you expect by printing the first 10 events (I've done this for you)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "trial_index = #_FILL_IN_YOUR_CODE_HERE \n",
    "print(t_EEG[np.array(trial_index).astype(int)]) # access the value at the corresponding index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_EEG_indices(event_timestamps, fs):\n",
    "    # _FILL_IN_YOUR_CODE_HERE\n",
    "    return \n",
    "\n",
    "# call your function to compute the corresponding indices\n",
    "EEG_indices = compute_EEG_indices()\n",
    "\n",
    "# print your solution and the actual event times to compare, they should be identical\n",
    "print(t_EEG[EEG_indices[:10]])\n",
    "print(trial_info[:10,0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [6] Q5: Step 3 - Grabbing Epochs\n",
    "Now that we have the corresponding indices in the EEG data, we know exactly where the **onset** of each stimulus is. The next thing we have to do is to grab a chunk of data surrounding the onset time, which we define to be `t=0` for every trial. That means you will want to grab a little bit of data before and after that time. \n",
    "\n",
    "[3] 5.1: Write a function that will, given an array of `data`, the sampling rate `fs`, and an `index`, grab a window of data surrounding that index, defined by `len_pre` and `len_post` in **seconds**. Note that `len_pre` should be negative to reflect that it's before the stimulus onset time. I've started this function for you below. Again, there are multiple ways to accomplish this, but the simplest solution can accomplish this in a single line.\n",
    "\n",
    "[1] 5.2: Use this function to grab an epoch for the **10th trial** (remember that's stored in `EEG_indices` already), with a pre-stimulus window of 0.5 seconds and a post-stimulus window of 1 second.\n",
    "\n",
    "[1] 5.3: Create a time vector `t_epoch` that corresponds to the timestamps for that epoch, relative to the stimulus onset time as zero. In other words, this time vector should start at `len_pre` and end at `len_post`, and has the same sampling frequency.\n",
    "\n",
    "[1] 5.4: Plot the epoch of data you grabbed. Note that the x-axis should be time. **Label your axes!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def grab_epoch(data, index, fs, len_pre, len_post):\n",
    "    # _FILL_IN_YOUR_CODE_HERE \n",
    "    return \n",
    "\n",
    "\n",
    "# _FILL_IN_YOUR_CODE_HERE \n",
    "len_pre = -0.5 #second\n",
    "len_post = 1 #second\n",
    "epoch = grab_epoch(_FILL_IN_YOUR_CODE_HERE)\n",
    "print(epoch[:5])\n",
    "\n",
    "t_epoch = # _FILL_IN_YOUR_CODE_HERE \n",
    "\n",
    "# plotting\n",
    "plt.figure(figsize=(6,4))\n",
    "# _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [4] Q6: Step 4 - Grab All & Baseline Correct (Bonus)\n",
    "[2] 6.1: If you grab an epoch for every trial and store that in a 2D numpy matrix, what should the dimensions of that matrix be, i.e., how many rows and how many columns? What do those numbers correspond to? Hint: you should organize your data such that there are more columns than rows in this particular case.\n",
    "\n",
    "[2] 6.2: Write a function that grabs **all** epochs (every trial) and store that in a 2D numpy matrix. There are a few ways to do this, but they will likely all use `grab_epoch()` somehow. Confirm that it has the same shape that you expect from above. Hint: you can append your epochs indefinitely to a python list using `list.append()`, and use `np.array()` to automatically convert that into a 2D matrix.\n",
    "\n",
    "[BONUS: 2] 6.3: Baseline all your epochs by subtracting the pre-stimulus epoch mean (-0.5 to 0 seconds) of each epoch from itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 6.1:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_all_epochs(data, indices, fs, len_pre, len_post):\n",
    "    # _FILL_IN_YOUR_CODE_HERE \n",
    "    # get all epochs\n",
    "\n",
    "    # baselining (if you want, it can also be a separate function)\n",
    "\n",
    "    return all_epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "epoched_EEG = get_all_epochs(EEG, EEG_indices, fs, len_pre, len_post)\n",
    "print(epoched_EEG.shape)\n",
    "\n",
    "# plot all the epochs and average\n",
    "plt.plot(t_epoch, epoched_EEG.T, '-k', alpha=0.01)\n",
    "plt.plot(t_epoch, np.mean(epoched_EEG,axis=0), label='Average Response')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Voltage (uV)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [6] Q7: Step 5 & 6 - Group Based on Trial Type\n",
    "In the plot above, I simply averaged over all the epochs to produce the average response (blue). However, as you will recall, there are several different types of trials (second column in `trial_info`). We should group epochs of the same trial type, and average over those. \n",
    "\n",
    "[5] 7.1: You have full flexibility for this part, with the only requirement being to produce a plot with 3 average responses corresponding to the 3 different trial types. Remember to label your plot axes and include a legend for which trace corresponds to which stimulus type. You will be evaluated on 3 things: whether you have successfully separated the epochs into their respective groupings, how well your code is commented to explain what you're doing, and whether you plot is correct and labeled. \n",
    "\n",
    "Since I have not given you a template for making a function, it may be useful to plan out what you want to do beforehand by writing pseudo code (i.e., plain English). Decide what strategy you will take (loops vs. list comprehension vs. others), and whether you want to separate the averaging and the plotting. You already know all the concepts required to tackle this problem (indexing, averaging, plotting), the challenge is putting them together. \n",
    "\n",
    "[1] 7.2: Briefly describe your results, e.g., what's similar and what's different between the conditions? Which stimulus produced the largest response.\n",
    "\n",
    "---\n",
    "Your plot should look something like:\n",
    "\n",
    "![example_erp](img/example_erp.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 7.2:** ANSWER 14 & -14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# That's All!\n",
    "There! You just performed your first neural data analysis. This type of stimulus-locked experiment design and analysis is very common in neuroscience, especially human and animal electrophysiology. Here at UCSD, Profs. Marta Kutas, Seana Coulson, and Sarah Creel all deploy these types of human EEG experiments to probe various aspects of the neural correlates of cognition. You will also see how it's applied in the paper we will discuss this week, and other commercial applications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# End Survery\n",
    "Please take a few minutes to fill out the following as it will help us to improve the following assignments & lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Content:\n",
    "What was one thing you learned from this lab & associated lectures?\n",
    "\n",
    "ANSWER: \n",
    "\n",
    "---\n",
    "What was one thing that you still found confusing after the lab, and need clarification?\n",
    "\n",
    "ANSWER:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Style:\n",
    "What was one thing you enjoyed about the formatting of this assignment (e.g., clarity, structure, guidance, etc.)?\n",
    "\n",
    "Answer:\n",
    "\n",
    "---\n",
    "What was one thing that you thought could use improvements on?\n",
    "\n",
    "Answer:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thank you!"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}