3-Linear-Time-Invariant-Systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Week 5: Linear Time-Invariant Systems\n",
    "\n",
    "Last week, we discussed the inhomogeneous Poisson model and how it can represent a point process that has a probability distribution that only depends on an underlying **rate** or **intensity**.\n",
    "\n",
    "This week, we'll discuss how to model the rate as a function of an external stimulus.\n",
    "\n",
    "After the lecture, work on the cells marked **Q** with help from your instructor. Save the completed notebook as an HTML file and submit on Canvas, one submission per group.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Group Members**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Readings\n",
    "\n",
    "Before coming to class, you should have finished reading Chapter 1 in Dayan and Abbott.\n",
    "\n",
    "For a more in-depth treatment, consult the following papers:\n",
    "\n",
    "- Aljadeff et al. Analysis of Neuronal Spike Trains, Deconstructed. Neuron. [doi:10.1016/j.neuron.2016.05.039](https://doi.org/10.1016/j.neuron.2016.05.039)\n",
    "- Touryan and Dan. Analysis of sensory coding with complex stimuli. Curr Opin Neurobiol [doi:10.1016/S0959-4388(00)00232-4](https://doi.org/10.1016/S0959-4388(00)00232-4)\n",
    "- Schwartz et al. Spike-triggered neural characterization. J Vis [doi:10.1167/6.4.13](https://doi.org/10.1167/6.4.13)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import sys\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "from comp_neurosci_uva import signal, data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Neurons are integrators\n",
    "\n",
    "Biological neurons may receive synaptic inputs from thousands of other cells.\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_purkinje.png\" alt=\"Purkinje cell\" style=\"width: 300px;\"/>\n",
    "\n",
    "Synapses can be excitatory or inhibitory, and can vary in strength.\n",
    "\n",
    "Roughly speaking, the neuron sums up these inputs and decides whether to spike. Its firing rate is therefore a function of the firing rates of its inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The linear neuron model\n",
    "\n",
    "We can conceptualize this process as a **weighted sum** of the inputs.\n",
    "\n",
    "![linear neuron](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_linear_neuron.png)\n",
    "\n",
    "In equation form, the diagram above looks like:\n",
    "\n",
    "$$y = \\sum_{i=0}^{N-1} w_i x_i$$\n",
    "\n",
    "For now, we are not going to think about spikes, but about rates. The symbols in this equation mean:\n",
    "\n",
    "- $x_i$ : the rate of presynaptic neuron $i$\n",
    "- $w_i$ : the synaptic weight of neuron $i$\n",
    "- $y$ : the output rate\n",
    "\n",
    "Note that this weighted sum is equivalent to the [dot product](https://en.wikipedia.org/wiki/Dot_product) of the vector of weights, $\\mathbf{w}$, with the vector of input rates, $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear sensory models\n",
    "\n",
    "It's quite difficult to know or control the firing rates of all the neurons that might be providing inputs to our target neuron. However, if we're recording from a living animal, we can control the stimulus. So for sensory neurons, without any loss of generality, we can think of the inputs to our model as dimensions of a stimulus. For a visual stimulus, this would be the individual pixels of an image:\n",
    "\n",
    "We can also think of the inputs to our model neuron as corresponding to dimensions of a stimulus.\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/mnist-input.png\" alt=\"visual perceptron\" style=\"width: 500px;\"/>\n",
    "\n",
    "That is, the value of $x_i$ would correspond to the intensity of the image at pixel $i$, and the whole image would be represented as an $N$-dimensional vector $\\mathbf{x} = (x_1, x_2, \\ldots, x_N)$.\n",
    "\n",
    "As before, the output rate of the target neuron, $y$ is simply a weighted sum over the pixels of the image.\n",
    "\n",
    "However, biological neurons are firmly embedded in the physical world, and one consequence is that they are responding to events that happened in the past. That is, they have a **memory**.\n",
    "\n",
    "How do we incorporate time into the model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear time-invariant (LTI) systems\n",
    "\n",
    "A simple approach is to allow some of the inputs to be delayed or **lagged**.\n",
    "\n",
    "Let's consider the case of a univariate input $x(t)$, with time discretized over some interval $\\Delta t$.\n",
    "\n",
    "If each input corresponds to an increasing delay, then:\n",
    "\n",
    "$$y(t) = \\sum_{i=0}^\\infty x(t - i\\cdot\\Delta t) w_i$$\n",
    "\n",
    "At time $t$, $w_0$ is the weight of the input at $t$, $w_1$ is the weight of the input at $t - \\Delta t$, and so forth.\n",
    "\n",
    "This should remind you of the basic formula for a linear neuron model, but instead of $i$ indexing different dimensions of the stimulus (or different synapses), it indexes past moments in time.\n",
    "\n",
    "If the weights $w_i$ remain constant, this model is an example of a **linear**, **time-invariant** system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Properties of LTI systems\n",
    "\n",
    "Consider a dynamical system with input x(t) and output y(t).\n",
    "\n",
    "![dynamical system](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_dynamical_system.png)\n",
    "\n",
    "The system is linear if it obeys the laws of superposition and scaling. That is, if\n",
    "\n",
    "\\begin{align}\n",
    "a(t) & \\rightarrow b(t) \\\\\n",
    "c(t) & \\rightarrow d(t)\n",
    "\\end{align}\n",
    "\n",
    "Then the following must be true:\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha a(t) & \\rightarrow \\alpha b(t) \\\\\n",
    "\\alpha a(t) + \\beta c(t) & \\rightarrow \\alpha b(t) + \\beta d(t)\n",
    "\\end{align}\n",
    "\n",
    "Furthemore, for the system to be time-invariant, shifting the input in time will lead to the same output, shifted in time by an equal amount:\n",
    "\n",
    "\\begin{align}\n",
    "a(t + \\Delta) & \\rightarrow b(t + \\Delta)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Impulse Response Functions\n",
    "\n",
    "How can we characterize the transformation that occurs between the input and output of an LTI system?\n",
    "\n",
    "It turns out, we only need to know one thing: the system's **impulse response function**. This is the response you would obtain to a small pulse of unit amplitude, $\\delta(t) \\rightarrow h(t)$.\n",
    "\n",
    "For example, here's the impulse response function of a tuning fork. When you hit it, it starts oscillating at 100 Hz and then decays (rather quickly, for illustration)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0.0, 1.0, 1000)\n",
    "h = np.sin(100 * t) * np.exp(- 5 * t)\n",
    "plt.plot(t, h)\n",
    "plt.xlabel(\"Time (s)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Once we know $h(t)$, we can compute the output $y(t)$ to **any** input $x(t)$. Why? Because any signal can be represented as a sum of time-shifted unit impulses of varying amplitude:\n",
    "\n",
    "$$x(t) \\equiv \\sum_i \\delta(t - \\tau_i) x(\\tau_i)$$\n",
    "\n",
    "Recall that $\\delta(t)$ is equal to $\\infty$ at $t=0$, is zero everywhere else, and that the area under $\\delta(t)$ is equal to 1. $\\delta(t - \\tau_i)$ simply shifts the impulse to $\\tau_i$, and $x(\\tau_i)$ scales it by the value of $x(t)$ at $t = \\tau_i$. \n",
    "\n",
    "Although this equation seems trivial, it's the basis of **discretization** and an important operation called **convolution**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NB: the python keyword `lambda` allows us to define simple functions inline\n",
    "f = lambda t: np.sin(2 * np.pi * t)\n",
    "plt.plot(t, f(t))\n",
    "tau = np.linspace(0.0, 1.0, 10)\n",
    "plt.vlines(tau, 0, f(tau))\n",
    "plt.plot(tau, f(tau), 'ko')\n",
    "plt.xlabel(\"Time (s)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convolution\n",
    "\n",
    "So how do we predict the response of our system from $h(t)$?\n",
    "\n",
    "Because the system is time-invariant, we know the response to an impulse shifted by $\\tau_i$ is just $h(t)$ shifted by the same amount:\n",
    "\n",
    "\\begin{align}\n",
    "\\delta(t - \\tau_i) & \\rightarrow h(t - \\tau_i)\n",
    "\\end{align}\n",
    "\n",
    "If this impulse is scaled by $x(\\tau_i)$, then the output is scaled equally:\n",
    "\n",
    "\\begin{align}\n",
    "x(\\tau_i) \\delta(t - \\tau_i) & \\rightarrow x(\\tau_i) h(t - \\tau_i)\n",
    "\\end{align}\n",
    "\n",
    "Finally, because of superposition, if we add together many scaled, time-shifted impulses, the output is simply the sum of the scaled, time-shifted responses:\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_i x(\\tau_i) \\delta(t - \\tau_i) & \\rightarrow \\sum_i x(\\tau_i) h(t - \\tau_i)\n",
    "\\end{align}\n",
    "\n",
    "So,\n",
    "\n",
    "$$y(t) \\equiv \\sum_i h(t - \\tau_i) x(\\tau_i)$$\n",
    "\n",
    "Here's an illustration for the sine function in the previous cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(9, 6))\n",
    "tot = np.zeros(t.size + h.size)\n",
    "axes[0].vlines(tau, 0, f(tau))\n",
    "axes[0].plot(tau, f(tau), 'ko')\n",
    "axes[0].set_title(r\"$x(\\tau_i)$\")\n",
    "# this loop is implementing the summation in the convolution equation\n",
    "for tau_i in tau:\n",
    "    hf = h * f(tau_i)\n",
    "    axes[1].plot(t + tau_i, hf)\n",
    "    idx = t.searchsorted(tau_i)\n",
    "    tot[idx:idx+hf.size] += hf\n",
    "axes[2].plot(np.linspace(0, 2.0, tot.size), tot)\n",
    "axes[1].set_title(r\"$h(t - \\tau_i)x(\\tau_i)$\")\n",
    "axes[2].set_title(r\"$\\sum h(t - \\tau_i)x(\\tau_i)$\")\n",
    "axes[2].set_xlabel(\"Time (s)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how constructive and destructive interference produces a rather complex pattern."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "**Q** Why does the signal decay to zero?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convolution (II)\n",
    "\n",
    "As the spacing between samples becomes infinitesimally small, the convolution sum becomes an integral:\n",
    "\n",
    "$$y(t) = \\int_{-\\infty}^\\infty h(t - \\tau) x(\\tau) d\\tau$$\n",
    "\n",
    "The shorthand for convolution is $*$ : $y(t) = (h * t)(t)$\n",
    "\n",
    "Convolution is commutative:\n",
    "\n",
    "\\begin{align}\n",
    "(h * t)(t) & = (t * h)(t) \\\\\n",
    "\\sum_i h(t - \\tau_i) x(\\tau_i) & = \\sum_i h(\\tau_i) x(t - \\tau_i)\n",
    "\\end{align}\n",
    "\n",
    "Convolution can be done in any domain.\n",
    "\n",
    "One way of interpreting the convolution sum is that it tells us that the output is computed by taking a *weighted sum of the present and past input values*. We can see this by writing out the sum:\n",
    "\n",
    "$$\\sum_i h(\\tau_i) x(t - \\tau_i) = h(0)x(t) + h(1)x(t - 1) + \\cdots$$\n",
    "\n",
    "The system is **causal** if $h(\\tau)$ is only greater than zero for $\\tau \\geq 0$.\n",
    "\n",
    "In most physical systems, the impulse response decays away with time, so there is a point where we can consider $h(\\tau)$ to be essentially zero and truncate the function.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Here are some visual illustrations of convolution from [Wikipedia](https://en.wikipedia.org/wiki/Convolution). Essentially you are taking one of the functions, flipping it in time, sliding it past the other function, and computing the area where the two functions overlap. Convolution is also called **filtering**.\n",
    "\n",
    "![convolution_animation](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_convolution_box.gif)\n",
    "\n",
    "![convolution_diagram](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_convolution_static.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LTI neuron models\n",
    "\n",
    "We now have our first model of how sensory neurons respond to stimuli.\n",
    "\n",
    "In essense, we are representing the neuron as a linear filter that computes a weighted sum of the stimulus as it varies in time.\n",
    "\n",
    "Although very simple, LTI models (and linear filters) can generate surprisingly complex behavior. The exercises in this notebook will help you explore some of this complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model systems\n",
    "\n",
    "We will investigate the properties of two different LTI models. Their impulse response functions are:\n",
    "\n",
    "\\begin{align}\n",
    "h_1(t) & = \n",
    "\\frac{t}{\\tau_1^2} e^{(-t/\\tau_1)} \\\\\n",
    "h_2(t) & = \n",
    "\\frac{t}{\\tau_1^2} e^{(-t/\\tau_1)} - \\frac{t}{\\tau_2^2} e^{(-t/\\tau_2)}\n",
    "\\end{align}\n",
    "\n",
    "Both filters are causal, so $h_1(t) = h_2(t) = 0$ for all $t < 0$.\n",
    "\n",
    "Functions with the general form of $t \\exp (-t)$ are called **alpha** functions. To get you started, I've defined a Python *function* that will generate alpha kernels for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alpha(tau, duration, dt):\n",
    "    \"\"\"An alpha function kernel with time constant tau, scaled to \n",
    "    \n",
    "    tau: the time constant of the kernel (in units of duration/dt)\n",
    "    duration: the duration of the support for the kernel\n",
    "    dt: the sampling interval of the kernel\n",
    "    \n",
    "    Returns a tuple (h(t), h(t))\n",
    "    \"\"\"\n",
    "    t = np.arange(0, duration, dt)\n",
    "    k = t / tau**2 * np.exp(-t / tau)\n",
    "    return (k, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Like mathematical functions, Python functions take one or more **arguments**. The result of applying the function to those arguments is the **return value**. Functions can only return a single value, but we can easily get around this by returning a **tuple**. In Python, a tuple is a kind of list that can't be modified. You can unpack the tuple into separate variables when you call a function by using **deconstruction**, as illustrated in the code cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1, t = alpha(50, 1000, 1)\n",
    "plt.plot(t, h1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.** Let $\\tau_1 = 50$ ms and $\\tau_2 = 100$ ms. Plot $h_1(t)$ and $h_2(t)$ for $0 < t < 1000$ ms, with a time step of 1 ms. Use the cell above as a model for your code. Note that $h_2(t)$ is the *sum* of $h_1(t)$ and another alpha function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**2.** Consider three input signals:\n",
    "\n",
    "\\begin{align}\n",
    "s_1(t) & = \\sin(2 \\pi \\omega_1 t) \\\\\n",
    "s_2(t) & = \\sin(2 \\pi \\omega_2 t) \\\\\n",
    "s_3(t) & = \\mathrm{sign}\\; [s_1(t)]\n",
    "\\end{align}\n",
    "\n",
    "Let $\\omega_1 = 0.3$ Hz and $\\omega_2 = 3$ Hz. For $s_3$, `sign` means that the value is 1.0 if $s_1(t) > 0$ and -1.0 if $s_1(t) \\leq 0$.\n",
    "\n",
    "Generate and plot 10 s of data for each signal, using a time step of 1 ms. Keep your time units consistent!\n",
    "\n",
    "Hint: Use `plt.subplots` to generate a nice grid of plots (see above for examples)\n",
    "\n",
    "Another hint: Don't reinvent the wheel! See if there's a numpy function that will compute `sign` for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**3A.** Convolve $s_1$ with the $h_1$ impulse response functions and plot the results. Do this by writing your own loop. See the model above for an example of how to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**3B.** Now use `np.convolve` to compute the convolution of $s_1$ and $h_1$. Do you get the same result as when you did it manually? If not, you might need to play with the `mode` argument to `np.convolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**3C.** Now do the convolution for each combination of signal ($s_1$, $s_2$, and $s_3$) and kernel ($h_1$ and $h_2$) and plot the results.\n",
    "\n",
    "Hint: Use `plt.subplots` to generate a grid of 6 panels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3D.** What do you notice about the response amplitudes? What differences do you see between the outputs of the two filters?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**4.** Now let's consider a white noise input:\n",
    "\n",
    "$$s_4(t) \\sim N(0,1)$$\n",
    "\n",
    "$N$ means that each sample is drawn from a normal distribution with mean 0 and standard deviation 1.0. (Hint: use `np.random.randn`)\n",
    "\n",
    "Compute and plot $(h_1 * s_4)$ and $(h_2 * s_4)$, with $s_4(t)$ evaluated over a 10 s interval with resolution of 1 ms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**5.** It's a bit hard to compare the results of the convolution in the time domain, so let's see what the spectrum looks like.\n",
    "\n",
    "We'll discuss spectral analysis in some detail later in the course, but for now I've provided the code you need. Just change the variable names to match what you used in the previous question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "# s4 should be your white noise signal\n",
    "freq, S4 = signal.welch(s4, nperseg=10000, fs=1000)\n",
    "plt.plot(freq, S4, 'k:', label=r\"$s_4$\")\n",
    "# the variable h1s4 should contain the convolution of s_4 with h_1\n",
    "freq, H1S4 = signal.welch(h1s4, nperseg=10000, fs=1000)\n",
    "plt.plot(freq, H1S4, label=r\"$h_1 * s_4$\")\n",
    "# h2s4 = h_2 * s_4\n",
    "freq, H2S4 = signal.welch(h2s4, nperseg=10000, fs=1000)\n",
    "plt.plot(freq, H2S4, label=r\"$h_2 * s_4$\")\n",
    "plt.xlim(0, 20)\n",
    "plt.legend()\n",
    "plt.xlabel(\"Frequency (Hz)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5A** What do you notice about the spectra? Do they seem noisy? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5B** Use a 1000-s input signal to to get a better estimate of the spectrum.\n",
    "\n",
    "You can copy the code in the cell above to calculate and plot the spectra, but you'll need to write the code for generating the longer signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5C.** Now that we have a nice plot, describe the following:\n",
    "\n",
    "- the shape of the spectrum for the input signal\n",
    "- the shape of the spectra for the two convolutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation\n",
    "\n",
    "Is it possible to estimate the impulse response function from the output of an LTI system when then input is *not* an impulse?\n",
    "\n",
    "Yes! The opposite operation of convolution is called **correlation** or **cross-correlation**\n",
    "\n",
    "$$(a \\star b)(t) = \\sum_i a(t + \\tau_i) b(t)$$\n",
    "\n",
    "Notice how similar the definition is to that of convolution. However, there is a critical sign difference: in convolution the \"sliding\" function is inverted in time, in correlation it is not.\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l6_correlation_static.png\" alt=\"correlation_diagram\" style=\"width: 300px;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6.** Compute and plot the correlation between $s_4(t)$ and $(s_4 * h_1)(t)$, and between $s_4(t)$ and $(s_4 * h_2)(t)$.\n",
    "\n",
    "Use `np.correlate` with the argument `mode=\"same\"`. Note that although *we* know $h_1$ and $h_2$ are causal, the correlation function does not. The output of `np.correlate` therefore contains both the causal ($t \\geq 0$) and acausal ($t < 0$) components.\n",
    "\n",
    "Try using 10 seconds of data first and then 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the outputs of the correlation compare to the original $h_1$ and $h_2$ kernels?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens if you reverse the order of the arguments to `correlate`?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.11.4"
  },
  "rise": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}