Assignments/A2-FourierTransform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Name\n",
    "- PID\n",
    "- COGS118C - Assignment 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook has [45 + 5 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": {},
   "source": [
    "# Lab 2 - Fourier Analysis & Coding Your Own DFT\n",
    "In this lab, we will learn about Fourier analysis the hands-on way: through building your own Discrete Fourier transform (DFT). You will employ complex sinusoids and dot products from Lab 0, and apply your DFT to a real signal to measure neural oscillations via the power spectrum. There is a bonus section at the end to explain aliasing in the frequency domain.\n",
    "\n",
    "Key concepts:\n",
    "\n",
    "- Cosine and Sine waves\n",
    "- Dot product, revisited\n",
    "- Complex exponential, magnitude and phase\n",
    "- Wave number and frequency\n",
    "- Power spectrum\n",
    "- Aliasing\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": [
    "# Discrete Fourier Transform\n",
    "It's useful to keep the formula for DFT in mind, as we step through each component for building it in code.\n",
    "\n",
    "![](img/DFT_formula.png)\n",
    "from [Wikipedia][WP_DFT].\n",
    "\n",
    "[WP_DFT]:https://en.wikipedia.org/wiki/Discrete_Fourier_transform\n",
    "\n",
    "To unpack that scary-looking equation, notice first that the summation expression is performing a dot product between two discrete-time signals: our signal $x(n)$ and the complex exponential in polar form. \n",
    "\n",
    "The second line simply rewrites the complex exponential using Euler's formula, and because of the linearity of the dot product, we can treat the sine and cosine separately.\n",
    "\n",
    "k is the wave number, which is analogous to frequency.\n",
    "\n",
    "** At its core, the DFT turns an array of numbers, our signal in time, into another array of numbers, its \"frequency domain representation\".**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# make the imports\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the LFP signal\n",
    "Today, we'll be working with LFP recorded in the rat hippocampus. This dataset comes from an openly accessible neuroscience database. For more information on this particular dataset, see [here](https://crcns.org/data-sets/hc/hc-2/about-hc-2).\n",
    "\n",
    "As was with last assignment, we will load the .mat file into python, and unpack the dictionary (I've done all this for you). In this file, LFP values are actually stored as integers (representing ADC levels), so we'll just label the y-axis as voltage without specifying a particular unit (like mV). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = io.loadmat('./data/LFP.mat', squeeze_me=True)\n",
    "print(data.keys())\n",
    "\n",
    "# unpack the variables\n",
    "fs = data['fs'] # sampling rate \n",
    "print('Sampling rate = %iHz'%fs)\n",
    "\n",
    "lfp = data['lfp'][0,:]/1000 # this file contains two channels, we'll only work with the first one\n",
    "lfp_short = lfp[:int(2*fs)] # make a variable that has only the first two seconds of the LFP\n",
    "t_short = np.arange(0, len(lfp_short)/fs, 1/fs) # create the corresponding time vector\n",
    "\n",
    "plt.figure(figsize=(15,3))\n",
    "plt.plot(t_short,lfp_short, 'k', alpha=0.8)\n",
    "plt.xlim([0,2])\n",
    "plt.xlabel('Time (s)');plt.ylabel('Voltage');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Oscillations\n",
    "You will notice that the LFP signal has a very prominent rhythmic component, with peaks near, for example, 0.25, 0.37, and 0.5 seconds. Therefore, it has a **periodicity of about 0.125 seconds**, or a **frequency of 8Hz**. This is a very famous brain oscillation, called the theta oscillation. \n",
    "\n",
    "Fourier, or frequency domain analysis, is particularly useful for characterizing this type of rhythmic brain signal. The method has a long history from electrical and systems engineering, both for describing signals and processes that generate signals (linear time invariant systems). But first, we'll construct the primary ingredient of the Fourier Transform: **sine waves**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [10] Q1: Sine and Cosine Waves\n",
    "Before we start analyzing the neural data, we have to go over the concept of sine waves. As we alluded to above, the (Discrete) Fourier Transform is a series of dot products, between sines/cosines of different frequencies and our signal. This might be a little confusing, if you are mainly familiar with computing the cosine of an angle, which gives you a single number back. For example, the cosine of 90 degrees (or $\\pi/2$) is 0. \n",
    "\n",
    "Here, we want a sine *wave*, i.e., to compute cosine and sine over time to produce a vector, and you can think of that as computing the cosine of different \"angles\" at the same time points as your LFP signal. But what are these angles?\n",
    "\n",
    "It's helpful to picture a clock in your head. A clock's seconds hand rotates around the circle, which has 360 degrees, or $2\\pi$ radians. Every second, it moves clockwise by a little bit. 1/60th of $2\\pi$, to be exact. Therefore, at every time point, you can compute the cosine of the **angle between the seconds hand and 12 o'clock** (due North). \n",
    "\n",
    "**(Note that, mathematically, the reference for 0 degree is actually 3 o'clock, or due East. But we'll use 12 o'clock for this question because it's a little more intuitive.)**\n",
    "\n",
    "---\n",
    "[1] 1.1: Let `t=0` be the start of the minute (pointing at 12 o'clock). What is the angle (in radians) between the seconds hand and the reference (12 o'clock) at `t=15` seconds? Save that in the variable `angle_t15`. Hint: `np.pi` has the value for $\\pi$.\n",
    "\n",
    "[1] 1.2: What is the cosine and sine of that angle? Answer this without coding.\n",
    "\n",
    "[1] 1.3: Complete the code to compute and print the quantities above and make sure this is what you expected.\n",
    "\n",
    "[1] 1.4: What is the angle at `t=30` seconds? What is the cosine and sine of that angle?\n",
    "\n",
    "[1] 1.5: After how many seconds will the angle be $2\\pi$? In other words, in how many seconds will the seconds hand have completed a full revolution around the clock? \n",
    "\n",
    "[1] 1.6 What is the frequency of rotation of the seconds hand on a clock (in cycles per second, Hz)? Store that in the variable `clock_freq`.\n",
    "\n",
    "[4] 1.7: The above questions gave you two points where time is mapped to angle. In code below, compute the corresponding angle (**in radians**) for all time points in the time vector `t_clock`, and store it in the variable `angle_clock`. Note, you should use the variables `t_clock`, `clock_freq`, and the constant $\\pi$. Then, using this angle array, compute the cosine and sine of the entire vector and store them in `cos_clock` and `sin_clock`. Finally, plot both the cosine and sine values over time. Remember to label your lines and axes. Hint: both cosine and sine should go through just a single cycle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 1.2:**\n",
    "\n",
    "**Response for 1.4:**\n",
    "\n",
    "**Response for 1.5:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(np.pi) # print the value for pi\n",
    "\n",
    "# compute the cosine and sine at t=15 seconds\n",
    "angle_t15 = # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "# print the cosine and sine of that angle\n",
    "print(# _FILL_IN_YOUR_CODE_HERE)\n",
    "\n",
    "clock_freq = # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "t_clock= np.arange(0,60,1/fs)\n",
    "angle_clock = # _FILL_IN_YOUR_CODE_HERE\n",
    "cos_clock = # _FILL_IN_YOUR_CODE_HERE\n",
    "sin_clock = # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "plt.figure(figsize=(15,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] Q2: Cosine and Sine Dot Product\n",
    "You should be familiar with how to construct a cosine and sine wave of a particular frequency at this point. To start analyzing our brain data, we'll use a cosine/sine of a different frequency: 8Hz. \n",
    "\n",
    "[1] 2.1: Make an 8Hz cosine and sine wave (store in variables `cos` and `sin`) using the time vector `t_short` and the frequency variable `freq`. This should be analogous to what you did above.\n",
    "\n",
    "[1] 2.2: Plot, in the same graph, the 2 seconds of LFP signal (already done), as well as the cosine and sine you have just generated. Remember to label the traces and axes. \n",
    "\n",
    "[1] 2.3: How many cycles should an 8Hz oscillation go through in 2 seconds? Count the peaks in your plot to make sure the answer is what you expect.\n",
    "\n",
    "[2] 2.4: Compute the dot product between the 2 second chunk of LFP and the cosine vector, and store in the variable `a_coef`; repeat for the LFP and the sine vector, store in variable `b_coef`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 2.3:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "freq = 8 # Hz\n",
    "cos = # _FILL_IN_YOUR_CODE_HERE\n",
    "sin = # _FILL_IN_YOUR_CODE_HERE\n",
    "plt.figure(figsize=(15,3))\n",
    "plt.plot(t_short, lfp_short, 'k', label='LFP', alpha=0.5)\n",
    "plt.xlim([0,2]); plt.legend()\n",
    "\n",
    "a_coef = # _FILL_IN_YOUR_CODE_HERE\n",
    "b_coef = # _FILL_IN_YOUR_CODE_HERE\n",
    "print(a_coef, b_coef)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] Q2 - continued\n",
    "We'll be generating sine waves and performing the dot product very often from this point on, so we better make them into functions. Complete the functions below. Note that you essentially already wrote all the code in the cell above, we just want to wrap them inside functions for convenience.\n",
    "\n",
    "[1] 2.5: For the function `make_cos_sin`: given a time vector `t`, and a frequency of interest `freq`, return the corresponding cosine and sine waves.\n",
    "\n",
    "[2] 2.6: For the function `dot_cos_sin`: given a signal, the associated time vector, and the frequency of interest, first make cosine and sine waves by calling the function `make_cos_sin`, then compute and return the dot product between `signal` and the cosine and sine wave.\n",
    "\n",
    "[1] 2.7: Call your function `dot_cos_sin` to repeat the above analysis (for 8Hz waves). Confirm that you get the same answer as above. Print the value.\n",
    "\n",
    "[1] 2.8: Perform the same computation, but for 30Hz cosine/sine waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_cos_sin(t, freq):    \n",
    "    # return the corresponding cosine and sine waves for time vector t, at frequency freq\n",
    "    # _FILL_IN_YOUR_CODE_HERE    \n",
    "    return cos, sin\n",
    "\n",
    "def dot_cos_sin(signal, t, freq):\n",
    "    cos, sin = # _FILL_IN_YOUR_CODE_HERE\n",
    "    a_coef = # _FILL_IN_YOUR_CODE_HERE\n",
    "    b_coef = # _FILL_IN_YOUR_CODE_HERE\n",
    "    return a_coef, b_coef\n",
    "\n",
    "# compute the dot products for 8Hz cosine and sine waves\n",
    "freq=8\n",
    "# _FILL_IN_YOUR_CODE_HERE\n",
    "print()\n",
    "\n",
    "# compute the dot products for 30Hz cosine and sine waves\n",
    "freq=30\n",
    "# _FILL_IN_YOUR_CODE_HERE\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3 + 2 Bonus] Q3: Complex Exponential, Power (Amplitude), and Phase\n",
    "Why did we just do that? One way to conceptualize the DFT is that it's a \"sinusoidal sieve\": by passing a signal through the DFT, we split up the signal into different frequencies, and measure the **independent** contribution of each of those component frequencies. This is what a prism does, separating white light into lights of different frequencies. What we did in Q2 was to construct a small part of the sieve, creating the \"hole\" at 8Hz and 30Hz. Passing the signal through (via dot product) \"filters\" out everything except the component of the signal at 8Hz. It's like wearing red goggles: it blocks/absorbs all other colors except red. \n",
    "\n",
    "---\n",
    "We performed two dot products with the cosine and sine waves, separately computing `a_coef` and `b_coef`. But if you refer to the first line of the DFT equation, you see that it's actually performing just a single dot product: the signal dotted with (projected onto) the complex exponential. Using Euler's formula, we can separate the complex exponential into a real (cosine) and an imaginary (sine) component. The computed coefficients, `a_coef` and `b_coef`, are the real and imaginary component of that dot product result, which form a 2D (complex) vector.\n",
    "\n",
    "![](./img/complex_exp.png)\n",
    "Adapted from [Better Explained](https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/).\n",
    "\n",
    "What is this complex number, `(a_coef + i*b*coef)`? This complex number gives us the signals power (or amplitude) and phase delay at a particular frequency, such as 8Hz. Power is how much the 8Hz component contributes to (or is in) the signal, while phase is the angular shift we have to apply to the cosine (real) component to find maximum alignment with the signal, or the \"starting point\". To get power and phase, we have to convert that complex number from rectangular to polar coordinate, which we did in Lab 0. \n",
    "\n",
    "**Note**: power is squared amplitude (or squared magnitude) by convention. The FT returns amplitude and phase, so that's what we will be doing.\n",
    "\n",
    "---\n",
    "[BONUS: 2] 3.1: Construct the complex exponential using `np.exp` (instead of the components separately as above), and perform the dot product directly with the signal `lfp_short`, confirm that you get the same complex number as `(a_coef + i*b*coef)`, which I've printed for you.\n",
    "\n",
    "[2] 3.2: Complete the function `compute_amp_phase` for computing the amplitude (magnitude) and phase (angle) from the coefficients. You may use the numerical formula for magnitude and angle in A0, or use `np.abs()` and `np.angle()` on the complex number itself.\n",
    "\n",
    "[1] 3.3: Using the above function, find the magnitude and phase of the 8Hz component of the signal. Store them in the variables `amp` and `phase`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# BONUS: computing the dot product from the complex exponential\n",
    "complex_exp = # _FILL_IN_YOUR_CODE_HERE\n",
    "# print the dot product\n",
    "print() # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "# computing the dot product separately\n",
    "freq = 8\n",
    "a_coef, b_coef = dot_cos_sin(lfp_short, t_short, freq)\n",
    "print(a_coef-1j*b_coef) # printing the coefficients as a complex number, note the minus in front of sin\n",
    "\n",
    "def compute_amp_phase(a_coef, b_coef):\n",
    "    # _FILL_IN_YOUR_CODE_HERE    \n",
    "    return amp, phase\n",
    "\n",
    "amp, phase = # _FILL_IN_YOUR_CODE_HERE\n",
    "print(amp,phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Intermission\n",
    "Well done! You have built from scratch the atomic element of the Fourier Transform in python. To recap, we've gone over the concept of a time-varying sinusoid, created a cosine and sine wave (components of the complex exponential) that have a specific frequency (8Hz), projected our signal onto these \"sinusoidal bases\" to get the Fourier coefficients, and computed the magnitude and phase of the signal at 8Hz. To construct the full DFT, we just have to perform this sequence of operations:\n",
    "- create complex exponential/sinusoidal bases at a particular frequency\n",
    "- perform dot product to get complex coefficients\n",
    "- compute magnitude and phase\n",
    "\n",
    "and iterate over all the frequencies of interest. \n",
    "\n",
    "### The last questions are: at *how many* frequencies, and *which* frequencies?\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [8] Q4: Wave Number & Frequencies of DFT\n",
    "In the questions above, we've cherrypicked a particular frequency (8Hz), and used that sieve to query for the contribution of the 8Hz component in our signal. To build the full DFT, you have to complete the sieve. A sieve merely separates components of a material, it does not make anything disappear. The same is true for the DFT. The real ingenuity behind the Fourier Transform is that a signal - *any signal* - can be decomposed into contributions from sinusoidal components (or bases), if you pick the right frequencies. In other words, from the signal's Fourier transform, you can reconstruct the original signal with 100% accuracy, using the *inverse* Fourier Transform.\n",
    "\n",
    "One idea is to simply do this for every single possible frequency, 8Hz, 8.1Hz, 8.01Hz, etc. But we quickly realize that there are an infinite number of frequencies along the real number line, and that's probably an overkill anyway. If our discrete time signal is only 1000-points long, using a million sinusoids to represent that isn't very clever.\n",
    "\n",
    "The intuitive answer (and the right answer) is that **we need as many frequencies as we have samples in the signal itself**. In the LFP signal above (`lfp_short`), there are 2500 points, so we need 2500 frequencies of sine waves to deconstruct that signal. \n",
    "\n",
    "---\n",
    "The unintuitive part is that the DFT operation doesn't know (or care) about frequency in terms of real time, i.e., Hz. You can compute the DFT of any sequence of numbers, and it will return the coefficients from the dot product. **What it cares about is the \"wave number\", `k`, of the sinusoids, which always goes from 0 to N-1 (see eq.1).** The sampling rate of your signal, then, determines how fast in real clock time those wave numbers correspond to in units of Hz. So, in reality, the 8Hz above is actually computed from a wave number.\n",
    "\n",
    "![](./img/DFT_formula.png)\n",
    "\n",
    "To think about the wave numbers (somewhat) intuitively, think of a clock again, but a weird one. You have N of these clocks (N is the length of your sigal). All of these clocks have the same faces, which is the full circle equally divided into N equal partitions. Instead of always moving forward 1/60th of a circle, each clock moves forward by a different amount at each \"time step\". The 0th clock (k=0) moves 0/N steps per time step, so it always stays still (this is an especially weird one). The 1st clock (k=1) moves clockwise by 1/N steps on each time step. Hence, in N time steps, it will have completed just a single revolution around the clock. \n",
    "\n",
    "** In general, the kth clock will complete k revolutions around the clock after N steps. Altogether, these clocks are all the possible ones that complete an integer number of revolutions around the clock after N steps, from 0 times to N-1 times.**\n",
    "\n",
    "To compute their physical frequency in Hz, we need to know how much each time step is, in terms of real time. This is determined by our signal we are trying to analyze, and is simply the variable `dt` that we've been working with, or `1/fs`.\n",
    "\n",
    "---\n",
    "[1] 4.1: For `N=2000` and `fs=1000`, how long (in seconds) does the first sine wave (`k=1`) take to complete a single cycle/ revolution around the circle? (The weird clock might be helpful here: how many steps does it take to complete one cycle, and how much time passes with each step?) Provide your reasoning in words or equation.\n",
    "\n",
    "[1] 4.2: What is its frequency, in Hz (cycles per second)? Hint: this is just the inverse of the above answer\n",
    "\n",
    "[1] 4.3: What is the frequency of the 16th wave (`k=16`) in Hz? Follow the same process as the last two questions.\n",
    "\n",
    "[1] 4.4: What is the frequency of the 16th wave (`k=16`) in Hz, when `N=4000` and `fs=500`?\n",
    "\n",
    "[3] 4.5: From those examples, provide the general formula for wave frequency in Hz, given its wave number `k`, total number of samples `N`, and sampling frequency `fs`. Hint: it's a simple algebraic expression using only these 3 terms. Complete the function below for computing this in code, and use that to confirm your answer for the last 3 questions. Print the results.\n",
    "\n",
    "[1] 4.6: In Q2, we used the 8Hz sine wave, with `fs=1250` and `N=2500`. What is its corresponding wave number, `k`? (k starts from 0.)\n",
    "\n",
    "---\n",
    "(`np.fft.fftfreq` is numpy's function for computing the frequency vector (in Hz) given `N` and `fs` or `dt`. I'm printing it beside your function calls to help you make sure you have the right answer.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 4.1:**\n",
    "\n",
    "**Response for 4.2:**\n",
    "\n",
    "**Response for 4.3:**\n",
    "\n",
    "**Response for 4.4:**\n",
    "\n",
    "**Response for 4.5:**\n",
    "\n",
    "**Response for 4.6:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def freq_from_k(k, N, fs):\n",
    "    # k can be a number or an array, the same expression should work\n",
    "\n",
    "    # _FILL_IN_YOUR_CODE_HERE\n",
    "    return \n",
    "\n",
    "k, N, fs = \n",
    "print(freq_from_k(k, N, fs), np.fft.fftfreq(N,1/fs)[k])\n",
    "\n",
    "k, N, fs = \n",
    "print(freq_from_k(k, N, fs), np.fft.fftfreq(N,1/fs)[k])\n",
    "\n",
    "k, N, fs = \n",
    "print(freq_from_k(k, N, fs), np.fft.fftfreq(N,1/fs)[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [+2 Bonus] Negative Frequencies & Aliasing - A Caveat\n",
    "In the plot below, I plotted the computed wave frequency of all the wave numbers using your function (blue), as well as using numpy's function (orange). From k=0 to 1250, the two lines should overlap (if your function is doing the right thing). After 1250Hz, however, your function should keep returning increasing wave frequencies for increasing wave number, while numpy returns **negative** frequencies.\n",
    "\n",
    "In short, in the strange clock example, a negative frequency would be when the hand rotates **counterclockwise**. For wave number `k=N-1`, for example, after a single time step, you can picture that as the hand jumping clockwise **almost** all the way around the full circle, to the last tick mark. OR, you can picture it as the hand jumping a single step counterclockwise, to -1/N. Because everything we work with is in discrete time, these two scenarios are indistinguishable, and the boundary at which this phenomenon occurs is at `k=N/2`, which corresponds to `fs/2`. This is known as the **Nyquist frequency**, and it's where sinusoidal rotations start to appear as if they are rotating backwards. See [here](https://www.youtube.com/watch?v=VNftf5qLpiA) for a visual demonstration. This is also what gives rise to **aliasing**, which we will explore in detail in the last section.\n",
    "\n",
    "[BONUS: 2] Uncomment the code below. Write a modified version of your `freq_from_k` function inside `freq_from_k_nyq` such that it returns negative wave frequencies for frequencies greater than the Nyquist frequency. Call the function and store the result in `your_freqs_nyq`, and plot it in the same plot. It should overlap completely with the numpy result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# def freq_from_k_nyq(k, N, fs):\n",
    "# _FILL_IN_YOUR_CODE_HERE\n",
    "#     return\n",
    "\n",
    "N, fs = 2500, 1250\n",
    "k_s = np.arange(0,N)\n",
    "your_freqs = freq_from_k(k_s, N, fs)\n",
    "np_freqs = np.fft.fftfreq(N, 1/fs)\n",
    "# your_freqs_nyq = # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(k_s,your_freqs, label='Your Wave Frequencies')\n",
    "plt.plot(k_s,np_freqs, label='numpy Wave Frequencies')\n",
    "# plt.plot() # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('Wave Number (k)')\n",
    "plt.ylabel('Wave Frequency (Hz)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# [6] Q5: DFT - Bringing Everything Together\n",
    "\n",
    "Now, you have all the components for writing your own Discrete Fourier Transform! To recap: \n",
    "- you can compute the corresponding cosine and sine wave frequencies given the wave numbers (0 to N-1)\n",
    "- you have a function that returns the wave frequencies (`freq_from_k`, or just call `np.fft.fftfreq`)\n",
    "- you have a function that creates the cosine and sine waves, given a frequency, and perform the dot product to get the coefficients\n",
    "\n",
    "### Your job is now to combine these parts to make the function `my_dft`.\n",
    "\n",
    "---\n",
    "[1] 5.1: Compute the wave frequencies using either your function `freq_from_k`, or `np.fft.fftfreq`, for the given `N` and `fs`. Store that in the variable `freqs`. This should be a vector of all frequencies of interest, in Hz.\n",
    "\n",
    "[2] 5.2: Complete the function `my_dft`, which computes the Fourier coefficients for `signal`, given the time vector `t` and frequency vector `freqs`. Hint: you will probably need some kind of loop here, or list comprehension.\n",
    "\n",
    "[1] 5.3: Call your function `my_dft` to compute the Fourier Transform of our snippet of LFP data, `lfp_short`, and store it in the variable `F_lfp`.\n",
    "\n",
    "[2] 5.4: Complete the 2 lines of plotting code that plots the real and imaginary components of your Fourier Transform as the y-values, while the x-values is the frequency axis (`freqs`). They should overlap completely with the numpy result. Depending on how you implemented `my_dft`, you may return the coefficients as a 2-dimensional array, or a 1-dimensional array of complex numbers, so index your array accordingly. Also, depending which function you called in part 2 to compute your frequency vector, this plot might look different, for the reason explained in the last section. Label your axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fs = 1250\n",
    "N = len(lfp_short)\n",
    "\n",
    "freqs = # _FILL_IN_YOUR_CODE_HERE\n",
    "\n",
    "def my_dft(signal, t, freqs):\n",
    "    # _FILL_IN_YOUR_CODE_HERE\n",
    "    return \n",
    "\n",
    "F_lfp = # _FILL_IN_YOUR_CODE_HERE\n",
    "F_lfp_np = np.fft.fft(lfp_short) # calling np.fft to make sure you have the right answers\n",
    "\n",
    "# plotting code\n",
    "plt.figure(figsize=(15,8))\n",
    "plt.subplot(2,1,1)\n",
    "# _FILL_IN_YOUR_CODE_HERE\n",
    "plt.plot(_FILL_IN_YOUR_CODE_HERE, 'o-', alpha=0.5, label='My DFT')\n",
    "plt.plot(freqs, F_lfp_np.real, '.-', alpha=0.5, label='numpy FFT')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(_FILL_IN_YOUR_CODE_HERE, 'o-', alpha=0.5) # you might need a negative in front of the imag component\n",
    "plt.plot(freqs, F_lfp_np.imag, '.-', alpha=0.5)\n",
    "\n",
    "# print the first 10 coefficients to make sure they're the same\n",
    "# the signs will be flipped for the imaginary component\n",
    "print(F_lfp[:5,:])\n",
    "print(F_lfp_np[:5])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [3 + 1 Bonus] Q6: Power Spectrum\n",
    "The above two plots are the real and imaginary coefficients spectra. While looking at them is interesting, we are often interested in another spectrum, the **power spectrum** of a signal, which is the squared magnitude of the complex vector (defined by coefficients `a_coef` and `b_coef`) at every frequency. This representation tells us which components are the strongest (has the highest power) in our data.\n",
    "\n",
    "[1] 6.1: Compute the array of power values from your Fourier coefficients `F_lfp`, and store it in the variable `my_power`.\n",
    "\n",
    "[1] 6.2: Plot the power spectrum (the x-axis is, again, frequency `freqs`), and zoom into the first 100Hz. Remember to label your axes. It should overlap completely with the numpy-computed power.\n",
    "\n",
    "[1] 6.3: By looking at the power spectrum of the LFP data, what frequency contributes the most (has the most power)? Is this consistent with time series plot in Q2? Justify your response.\n",
    "\n",
    "[BONUS: 1] 6.4: Confirm that your Fourier transform follows [Parseval's theorem](https://en.wikipedia.org/wiki/Parseval%27s_theorem) by computing and printing the values from both sides of this equation:\n",
    "![](./img/parseval.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 6.3:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "my_power = # _FILL_IN_YOUR_CODE_HERE\n",
    "np_power = np.abs(F_lfp_np)**2\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(# _FILL_IN_YOUR_CODE_HERE)\n",
    "plt.plot(freqs, np_power, '.', alpha=0.5, label='numpy Power')\n",
    "plt.legend()\n",
    "plt.xlim([0,100])\n",
    "\n",
    "# Parseval's theorem\n",
    "# _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [5] Q6 - continued\n",
    "[3] 6.5: Repeat the above analysis, but for the full LFP signal (`lfp`, I've re-defined for you). Plot the power spectrum (and remember to label your axes). You may use `np.fft.fft()` instead for a speed boost.\n",
    "\n",
    "[2] 6.6: Comment on similarities and differences between this plot and the above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Response for 6.6:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fs = data['fs'] # sampling rate\n",
    "lfp = data['lfp'][0,:]/1000\n",
    "\n",
    "# _FILL_IN_YOUR_CODE_HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Congratulations!\n",
    "You have just performed your first frequency (Fourier) analysis on real neural recordings. Computing and looking at the power spectrum is one of the most common ways of characterizing neural signals, especially for identifying strong rhythmic components in the data (neural oscillations). Almost every lab on campus that works with EEG, ECoG, or LFPs will at some point compute a power spectrum to look for rhythmic signals, including Profs. Lara Rangel, Doug Nitz, Brad Voytek, and many more. Neural oscillations is a huge field of neuroscience research, and there's plenty of work to be done here!"
   ]
  },
  {
   "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!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Alising & Nyquist Sampling Theorem in Frequency Domain\n",
    "Some bonus material: we now have the machinary to examine what exactly happens in the frequency domain when a signal is undersampled, thus causing aliasing. Let's start with a finely sampled sinusoid, with an oscillation frequency of 14Hz, sampled with a sampling rate of 1000Hz. \n",
    "\n",
    "We can then create a downsampled version of our signal by taking every 50th sample, which gives us an effective sampling rate of 1000/50 = 20Hz. By the Nyquist sampling theorem, we know that for a signal of frequency B, we have to sample at least 2B Hz. In other words, if the sampling rate is fs, the highest frequency signal we can resolve is fs/2, which is called the Nyquist frequency. \n",
    "\n",
    "Now, our downsampled signal is sampled at 20Hz, which means we can resolve a signal of at most 10Hz. Remember the weird clocks? Well, if a 14Hz signal sampled at 20Hz is like rotating around the circle clockwise 14 times per second, it will produce the same sequence of clock positions as if it's rotating around the circle counterclockwise 6 times per second, or -6Hz. The math is $14-20 = -6$\n",
    "\n",
    "Explore different values for the variables `cos_freq` and `ds_factor` and re-run the code to see which frequencies will not be aliased, i.e., the blue and orange traces in the power spectrum plot overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "t = np.arange(0,5,1/fs)\n",
    "cos_freq = 14\n",
    "signal = np.cos(2*np.pi*cos_freq*t) # this is our sinusoidal signal\n",
    "\n",
    "# downsampling\n",
    "ds_factor = 50\n",
    "fs_ds = fs/ds_factor\n",
    "print(fs_ds)\n",
    "t_ds = t[::ds_factor]\n",
    "signal_ds = signal[::ds_factor]\n",
    "\n",
    "plt.figure(figsize=(15,3))\n",
    "plt.plot(t, signal)\n",
    "plt.plot(t_ds, signal_ds, 'o-')\n",
    "plt.xlabel('Time (s)');plt.ylabel('Voltage');\n",
    "plt.xlim([0,1])\n",
    "\n",
    "# compute power spectrum of original signal\n",
    "f_signal = np.fft.fftfreq(len(signal), 1/fs)\n",
    "power_signal = np.abs(np.fft.fft(signal))\n",
    "print(f_signal[np.where(power_signal==max(power_signal))])\n",
    "\n",
    "# compute power spectrum of downsampled signal\n",
    "f_signal_ds = np.fft.fftfreq(len(signal_ds), 1/(fs_ds))\n",
    "power_signal_ds = np.abs(np.fft.fft(signal_ds))*(len(signal)/len(signal_ds))\n",
    "print(f_signal_ds[np.where(power_signal_ds==max(power_signal_ds))])\n",
    "\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(f_signal, power_signal, '-o', label='Original Signal')\n",
    "plt.plot(f_signal_ds, power_signal_ds, '-o', label='Downsampled Signal')\n",
    "plt.xlabel('Frequency (Hz)');plt.ylabel('Power');\n",
    "plt.xlim([-40,40]); plt.legend()"
   ]
  }
 ],
 "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
}