linear_systems_II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Systems II\n",
    "\n",
    "[return to main page](index.ipynb)\n",
    "\n",
    "## Preparations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-11-26T13:10:06.061351Z",
     "iopub.status.busy": "2020-11-26T13:10:06.060911Z",
     "iopub.status.idle": "2020-11-26T13:10:06.491573Z",
     "shell.execute_reply": "2020-11-26T13:10:06.491116Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "> 0 Built-in Microphone, Core Audio (2 in, 0 out)\n",
       "< 1 Built-in Output, Core Audio (0 in, 2 out)\n",
       "  2 DisplayPort, Core Audio (0 in, 2 out)\n",
       "  3 Aggregate Device, Core Audio (0 in, 2 out)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import tools\n",
    "import numpy as np\n",
    "import soundfile as sf  # for reading a soundfile\n",
    "from scipy import signal\n",
    "\n",
    "# use the system OS sound device for convenience\n",
    "import sounddevice as sd  # for playback\n",
    "# sd.default.device = 0  # set default device -> not used, know what you're doing\n",
    "sd.query_devices()  # check if desired default device is set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-dimensional continuous Systems\n",
    "\n",
    "In the first part of this exercise, we had a quick review on linear systems in one dimension, namely the time $t$. Now, we add three additional dimension $\\mathbf x = [x,y,z]^T$ in order to cover the three-dimensional space. Generally, the input signal $q \\in \\mathbb{C}$ and the corresponding output signal $p \\in \\mathbb{C}$ of a system $\\mathcal H$ are related via:\n",
    "\n",
    "$$p(\\mathbf x, t) = p(x,y,z,t) = \\mathcal{H}\\{q(\\mathbf x,t)\\} = \\mathcal{H}\\{q(x,y,z,t)\\}\\,.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Time-Space-Invariant (LTSI) Systems\n",
    "\n",
    "As simple as it sounds, LTI-system are linear, time-invariant and **space-invariant**. We had already covered the first two properties in the [last exercise](linear_systems_I.ipynb#Linear-Time-Invariant-%28LTI%29-Systems)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Space-Invariance\n",
    "\n",
    "*Exercise*: Explain the term \"space-invariance\" in your own words.\n",
    "\n",
    "\n",
    "*Exercise*: What does this mean mathematically?\n",
    "\n",
    "$$\\mathcal{H}\\{q(\\mathbf x - \\mathbf x_0, t)\\} = \\,???\\, \\text{ for all } \\mathbf x_0 \\in \\mathbb{R}$$\n",
    "\n",
    "if $p(\\mathbf x, t) = \\mathcal{H} \\{ q(\\mathbf x,t) \\}$ is known.\n",
    "\n",
    "*Exercise*: Can you give an example, for which environment the space invariance does **not** hold?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Impulse Response\n",
    "\n",
    "#### ... of an LTSI System\n",
    "\n",
    "The impulse response $h(\\mathbf x, t)$ of an LTSI system characterises it completely. It is the system's response to a spatio-temporal Dirac impulse $\\delta(\\mathbf x)\\delta(t)$\n",
    "\n",
    "$$h(\\mathbf x, t) = h(x, y, z, t) = \\mathcal{H}\\{\\delta(\\mathbf x)\\delta(t)\\} = \\mathcal{H}\\{\\delta(x)\\delta(y)\\delta(z)\\delta(t)\\}\\,.$$\n",
    "\n",
    "We can describe the output signal $p(\\mathbf x, t)$ by a 4-dimensional convolution of the corresponding input signal $q(\\mathbf x, t)$ and the impulse response $h(\\mathbf x, t)$: \n",
    "\n",
    "$$ p(\\mathbf x, t) = \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} q(x_0, y_0, z_0, t_0)\\,h(x - x_0, y - y_0, z - z_0, t - t_0)\\,\\mathrm d x_0 \\mathrm d y_0 \\mathrm d z_0 \\mathrm d t_0$$\n",
    "\n",
    "Its short version reads\n",
    "\n",
    "$$ p(\\mathbf x, t) = q(\\mathbf x, t) *_{\\mathbf x, t} h(\\mathbf x, t) $$\n",
    "\n",
    "where $*$ is a common notation of the convolution.\n",
    "\n",
    "#### ... of an shift-variant LTI System\n",
    "\n",
    "For a shift-variant system, the impulse response is not independent from the shift in space:\n",
    "\n",
    "$$h(\\mathbf x, \\mathbf x_0, t)  = \\mathcal{H}\\{\\delta(\\mathbf x - \\mathbf x_0)\\delta(t)\\}$$\n",
    "\n",
    "We can no longer describe the resulting output signal as a convolution with respect to the three dimensions in space:\n",
    "\n",
    "$$ p(\\mathbf x, t) = \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} q(x_0, y_0, z_0, t_0)\\,h(x,x_0,y,y_0,z,z_0,t-t_0)\\,\\mathrm d x_0\\,\\mathrm d y_0\\,\\mathrm d z_0\\,\\mathrm d t_0$$\n",
    "\n",
    "The integral with respect to time is still a linear convolution, though.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Spatial Frequency Response\n",
    "\n",
    "For an LTSI-System, the spatial transfer characteristics in terms of its frequency response is defined as the three-dimensional Fourier transform with respect to each dimension in space:\n",
    "\n",
    "$$ \\tilde h(\\mathbf k, t) = \\tilde h(k_x, k_y, k_z, t) = \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} h(x,y,z,t) e^{+jk_x x} e^{+jk_y y} e^{+jk_z z} \\,\\mathrm d x\\,\\mathrm d y\\,\\mathrm d z$$\n",
    "\n",
    "with the wave vector $\\mathbf k = [k_x, k_y, k_z]^T$. The Fourier transform with respect to one dimension is equivalent to the temporal Fourier transform with e.g. $\\omega \\rightarrow -k_x$.\n",
    "\n",
    "### The Spatio-Temporal Frequency Response\n",
    "\n",
    "Applying the temporal Fourier transform to $\\tilde h(\\mathbf k, t)$ yields the **Spatio-Temporal Frequency Response** of an LTSI-System\n",
    "\n",
    "$$ \\tilde H(\\mathbf k, \\omega) = \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} h(x,y,z,t) e^{+jk_x x} e^{+jk_y y} e^{+jk_z z} e^{-j\\omega t} \\,\\mathrm d x\\,\\mathrm d y\\,\\mathrm d z\\,\\mathrm d t$$\n",
    "\n",
    "## Sound Propagation as a Linear System\n",
    "\n",
    "Remember the inhomogeneous wave equation?\n",
    "\n",
    "$$\\Delta p(\\mathbf x, t) - \\frac{1}{c^2} \\frac{\\partial^2}{\\partial t^2} p(\\mathbf x, t) = -q(\\mathbf x, t)$$\n",
    "\n",
    "Its solution $p(\\mathbf x, t)$ has to incorporate boundary conditions and the inhomogenity $q(\\mathbf x, t)$. It describes sound propagation and can be interpreted as an LTI-system with the source density $q(\\mathbf x, t)$ as the input signal and the sound pressure $p(\\mathbf x, t)$ as the output signal. Remember the Green's function? It corresponds to the special inhomogenity $q(\\mathbf x, t) = \\delta(\\mathbf x - \\mathbf x_0)\\delta(t-t_0)$ and fulfils the boundary conditions. The inhomogeneous wave equation for the Green's functions reads\n",
    "\n",
    "$$\\Delta g(\\mathbf x | \\mathbf x_0 , t) - \\frac{1}{c^2} \\frac{\\partial^2}{\\partial t^2} g(\\mathbf x | \\mathbf x_0 , t) = - \\delta(\\mathbf x - \\mathbf x_0)\\delta(t)\\,.$$\n",
    "\n",
    "*Exercise*: What is the impulse response $h(\\mathbf x,\\mathbf x_0,t)$ of the wave equation as a linear system?\n",
    "\n",
    "Under free-field conditions, i.e. no boundary conditions, the wave equation becomes an LTSI-System with free-field Green's function $g_0(\\mathbf x , t) = \\delta(t - |\\mathbf x|/c)$ as its spatio-temporal impulse responses. In rooms, the impulse response is often subsumed under the term **Room Impulse Response**.\n",
    "\n",
    "### Measuring a Room Impulse Response\n",
    "\n",
    "#### A Very Simplistic Procedure\n",
    "\n",
    "There are several possibilities to excite the room with a signal being close to a Dirac - bursting balloons, gunshots, electrical spark discharges, ...; for simplicity (and for safety reasons), one could use two wooden boards which are clapped together. The signal recorded by a microphone standing at the certain positions $\\mathbf x$ then yields the approximate room impulse response $h(\\mathbf x,\\mathbf x_0,t)$\n",
    "\n",
    "![Two Wooden Boards](data/wooden_boards.jpg)\n",
    "\n",
    "That's not ideal - ideally all frequencies should be excited equally - but it shall be sufficient for demonstration purposes. As we will get to know later, that there are *much better* ways to measure room impulse responses!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-11-26T13:10:06.496070Z",
     "iopub.status.busy": "2020-11-26T13:10:06.495576Z",
     "iopub.status.idle": "2020-11-26T13:10:13.101530Z",
     "shell.execute_reply": "2020-11-26T13:10:13.101930Z"
    }
   },
   "outputs": [],
   "source": [
    "# speech signal\n",
    "speech, fs = sf.read(\"data/xmas.wav\")\n",
    "# sloppy measured room impulse response\n",
    "# (measured with clap seen above, do not harm your ears!)\n",
    "rir, fs_rir = sf.read(\"data/rir_clap.wav\")\n",
    "\n",
    "assert fs == fs_rir  # compare sample rates of signal and rir\n",
    "speech_clap = signal.fftconvolve(speech, rir)  # convolve with respect to time\n",
    "# normalize to the same maximum value as the original speech signal:\n",
    "speech_clap = tools.normalize(speech_clap, np.max(np.abs(speech)))\n",
    "\n",
    "sd.play(speech, fs, blocking=True)\n",
    "sd.play(speech_clap, fs, blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### More elaborated procedures\n",
    "\n",
    "In slightly more modern *sweep method*, room is excited with a sine sweep, which is reproduced by means of a loudspeaker. The actual impulse response will be calculated from the excitation signal and the signal recorded by the microphone. Let's listen to some sweep signals (**watch the volume**!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-11-26T13:10:13.105844Z",
     "iopub.status.busy": "2020-11-26T13:10:13.105398Z",
     "iopub.status.idle": "2020-11-26T13:10:19.392534Z",
     "shell.execute_reply": "2020-11-26T13:10:19.392948Z"
    }
   },
   "outputs": [],
   "source": [
    "fstart = 65  # start with 65 Hertz, be nice to your onboard loudspeakers\n",
    "fstop = 15000  # stop a 15000 Hertz\n",
    "tlen = 3  # length of sweep signal\n",
    "t = np.arange(0, tlen*fs)/fs  # samples\n",
    "\n",
    "sweep_lin = signal.chirp(t, fstart, tlen, fstop, method='linear')\n",
    "sweep_log = signal.chirp(t, fstart, tlen, fstop, method='logarithmic')\n",
    "\n",
    "sd.play(sweep_lin, fs, blocking=True)\n",
    "sd.play(sweep_log, fs, blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fourier Analysis of a Sound Field\n",
    "\n",
    "### Spatial-Temporal Fourier Transform of a Plane Wave\n",
    "\n",
    "We already got to know the temporal spectrum of a [plane wave](physics_of_sound_I-solutions.ipynb#Plane-Wave) with a direction of propagation $\\mathbf n_{\\mathrm {pw}}$ which reads\n",
    "\n",
    "$$P_{pw}(\\mathbf x, \\omega) = \\mathrm{exp}\\left(-j\\frac{\\omega}{c} \\mathbf n_{\\mathrm {pw}} \\cdot \\mathbf x \\right) = \\mathrm{exp}\\left(-j\\frac{\\omega}{c} (n_{\\mathrm {pw},x} \\cdot x + n_{\\mathrm {pw},y} \\cdot y + n_{\\mathrm {pw},z} \\cdot z) \\right) \\,. $$\n",
    "\n",
    "*Exercise*: Calculate the three-dimensional, spatial Fourier Transform $\\tilde P_{pw}(\\mathbf k, \\omega)$ of the Plane Wave. Use the integral identity:\n",
    "\n",
    "$$ \\int_{-\\infty}^{\\infty} e^{+j a b}\\,\\mathrm d b\n",
    "    = \\int_{-\\infty}^{\\infty} \\mathrm{exp}\\left(+j a b\\right)\\,\\mathrm d b \n",
    "    = 2\\pi \\delta(a) \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourier Analysis using a Plane Wave Decomposition\n",
    "\n",
    "We now know that the value of $\\tilde P(\\mathbf k, \\omega)$ at $\\mathbf k = \\frac{\\omega}{c} \\mathbf n_{\\mathrm {pw}}$ corresponds to the amplitude of the plane wave propagation in the directions of $\\mathbf n_{\\mathrm {pw}}$. This can be straightforwardly used to compute the so-called plane wave composition:\n",
    "\n",
    "$$\n",
    "\\bar P(\\mathbf n_{\\mathrm {pw}}) = \\tilde P(\\mathbf k, \\omega)|_{\\mathbf k = \\frac{\\omega}{c} \\mathbf n_{\\mathrm {pw}}}\n",
    "$$\n",
    "\n",
    "After all this math, u can lay back and have a look at the demonstration, how to compute a plane wave decompositon out of two microphone signals.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions\n",
    "\n",
    "If you had problems solving some of the exercises, don't despair!\n",
    "Have a look at the [example solutions](linear_systems_II-solutions.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p xmlns:dct=\"http://purl.org/dc/terms/\">\n",
    "  <a rel=\"license\"\n",
    "     href=\"http://creativecommons.org/publicdomain/zero/1.0/\">\n",
    "    <img src=\"http://i.creativecommons.org/p/zero/1.0/88x31.png\" style=\"border-style: none;\" alt=\"CC0\" />\n",
    "  </a>\n",
    "  <br />\n",
    "  To the extent possible under law,\n",
    "  <span rel=\"dct:publisher\" resource=\"[_:publisher]\">the person who associated CC0</span>\n",
    "  with this work has waived all copyright and related or neighboring\n",
    "  rights to this work.\n",
    "</p>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mystiasp",
   "language": "python",
   "name": "mystiasp"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}