{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ThinkDSP\n",
    "\n",
    "This notebook contains solutions to exercises in Chapter 9: Differentiation and Integration\n",
    "\n",
    "Copyright 2015 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get thinkdsp.py\n",
    "\n",
    "import os\n",
    "\n",
    "if not os.path.exists('thinkdsp.py'):\n",
    "    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/thinkdsp.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from thinkdsp import decorate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "The goal of this exercise is to explore the effect of `diff` and `differentiate` on a signal.  Create a triangle wave and plot it.  Apply the `diff` operator and plot the result.  Compute the spectrum of the triangle wave, apply `differentiate`, and plot the result.  Convert the spectrum back to a wave and plot it.  Are there differences between the effect of `diff` and `differentiate` for this wave?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Solution:* Here's the triangle wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkdsp import TriangleSignal\n",
    "\n",
    "in_wave = TriangleSignal(freq=50).make_wave(duration=0.1, framerate=44100)\n",
    "in_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diff of a triangle wave is a square wave, which explains why the harmonics in a square wave drop off like $1/f$, compared to the triangle wave, which drops off like $1/f^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = in_wave.diff()\n",
    "out_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we take the spectral derivative, we get \"ringing\" around the discontinuities: https://en.wikipedia.org/wiki/Ringing_(signal)\n",
    "\n",
    "Mathematically speaking, the problem is that the derivative of the triangle wave is undefined at the points of the triangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave2 = in_wave.make_spectrum().differentiate().make_wave()\n",
    "out_wave2.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "The goal of this exercise is to explore the effect of `cumsum` and `integrate` on a signal.  Create a square wave and plot it.  Apply the `cumsum` operator and plot the result.  Compute the spectrum of the square wave, apply `integrate`, and plot the result.  Convert the spectrum back to a wave and plot it.  Are there differences between the effect of `cumsum` and `integrate` for this wave?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Solution:* Here's the square wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkdsp import SquareSignal\n",
    "\n",
    "in_wave = SquareSignal(freq=50).make_wave(duration=0.1, framerate=44100)\n",
    "in_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cumulative sum of a square wave is a triangle wave.  After the previous exercise, that should come as no surprise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = in_wave.cumsum()\n",
    "out_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectral integral is also a triangle wave, although the amplitude is very different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum = in_wave.make_spectrum().integrate()\n",
    "spectrum.hs[0] = 0\n",
    "out_wave2 = spectrum.make_wave()\n",
    "out_wave2.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we unbias and normalize the two waves, they are visually similar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave.unbias()\n",
    "out_wave.normalize()\n",
    "out_wave2.normalize()\n",
    "out_wave.plot()\n",
    "out_wave2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And they are numerically similar, but with only about 3 digits of precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave.max_diff(out_wave2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "The goal of this exercise is the explore the effect of integrating twice.  Create a sawtooth wave, compute its spectrum, then apply `integrate` twice.  Plot the resulting wave and its spectrum.  What is the mathematical form of the wave?  Why does it resemble a sinusoid? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the sawtooth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkdsp import SawtoothSignal\n",
    "\n",
    "in_wave = SawtoothSignal(freq=50).make_wave(duration=0.1, framerate=44100)\n",
    "in_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first cumulative sum of a sawtooth is a parabola:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = in_wave.cumsum()\n",
    "out_wave.unbias()\n",
    "out_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second cumulative sum is a cubic curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = out_wave.cumsum()\n",
    "out_wave.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrating twice also yields a cubic curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum = in_wave.make_spectrum().integrate().integrate()\n",
    "spectrum.hs[0] = 0\n",
    "out_wave2 = spectrum.make_wave()\n",
    "out_wave2.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the result looks more and more like a sinusoid.  The reason is that integration acts like a low pass filter.  At this point we have filtered out almost everything except the fundamental, as shown in the spectrum below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave2.make_spectrum().plot(high=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4 \n",
    "\n",
    "The goal of this exercise is to explore the effect of the 2nd difference and 2nd derivative.  Create a `CubicSignal`, which is defined in `thinkdsp`.  Compute the second difference by applying `diff` twice.  What does the result look like.  Compute the second derivative by applying `differentiate` twice.  Does the result look the same?\n",
    "\n",
    "Plot the filters that corresponds to the 2nd difference and the 2nd derivative and compare them.  Hint: In order to get the filters on the same scale, use a wave with framerate 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Solution:* Here's the cubic signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkdsp import CubicSignal\n",
    "\n",
    "in_wave = CubicSignal(freq=0.0005).make_wave(duration=10000, framerate=1)\n",
    "in_wave.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first difference is a parabola and the second difference is a sawtooth wave (no surprises so far):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = in_wave.diff()\n",
    "out_wave.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_wave = out_wave.diff()\n",
    "out_wave.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we differentiate twice, we get a sawtooth with some ringing.  Again, the problem is that the deriviative of the parabolic signal is undefined at the points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum = in_wave.make_spectrum().differentiate().differentiate()\n",
    "out_wave2 = spectrum.make_wave()\n",
    "out_wave2.plot()\n",
    "decorate(xlabel='Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The window of the second difference is -1, 2, -1.  By computing the DFT of the window, we can find the corresponding filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkdsp import zero_pad\n",
    "from thinkdsp import Wave\n",
    "\n",
    "diff_window = np.array([-1.0, 2.0, -1.0])\n",
    "padded = zero_pad(diff_window, len(in_wave))\n",
    "diff_wave = Wave(padded, framerate=in_wave.framerate)\n",
    "diff_filter = diff_wave.make_spectrum()\n",
    "diff_filter.plot(label='2nd diff')\n",
    "\n",
    "decorate(xlabel='Frequency (Hz)',\n",
    "                 ylabel='Amplitude ratio')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for the second derivative, we can find the corresponding filter by computing the filter of the first derivative and squaring it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "PI2 = np.pi * 2\n",
    "\n",
    "deriv_filter = in_wave.make_spectrum()\n",
    "deriv_filter.hs = (PI2 * 1j * deriv_filter.fs)**2\n",
    "deriv_filter.plot(label='2nd deriv')\n",
    "\n",
    "decorate(xlabel='Frequency (Hz)',\n",
    "                 ylabel='Amplitude ratio')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the two filters look like on the same scale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_filter.plot(label='2nd diff')\n",
    "deriv_filter.plot(label='2nd deriv')\n",
    "\n",
    "decorate(xlabel='Frequency (Hz)',\n",
    "                 ylabel='Amplitude ratio')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both are high pass filters that amplify the highest frequency components.  The 2nd derivative is parabolic, so it amplifies the highest frequencies the most.  The 2nd difference is a good approximation of the 2nd derivative only at the lowest frequencies, then it deviates substantially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}