rir-solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[exercises](rir.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directly listening to a room impulse response (RIR) doesn't reveal much information (except probably for room acoustics experts).\n",
    "It's normally more helpful to use a bunch of dry recordings with different characteristics (speech, music, tonal, percussive, ...), convolve them with the given RIR and listen to the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import soundfile as sf\n",
    "import tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speech, fs = sf.read(\"data/xmas.wav\")\n",
    "rir, fs_rir = sf.read(\"data/rir_clap.wav\")\n",
    "assert fs == fs_rir\n",
    "speech_clap = signal.fftconvolve(speech, rir)\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",
    "sf.write(\"data/xmas_clap.wav\", speech_clap, fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<audio src=\"data/xmas.wav\" controls></audio>\n",
    "[data/xmas.wav](data/xmas.wav)\n",
    "\n",
    "<audio src=\"data/rir_clap.wav\" controls></audio>\n",
    "[data/rir_clap.wav](data/rir_clap.wav)\n",
    "\n",
    "<audio src=\"data/xmas_clap.wav\" controls></audio>\n",
    "[data/xmas_clap.wav](data/xmas_clap.wav)\n",
    "\n",
    "It doesn't sound exactly like the measured room, because the frequency response of the clapping device is not flat and its characteristics are part of the measured RIR and therefore also audible in the convolved signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(len(rir)) / fs\n",
    "plt.plot(t, rir)\n",
    "plt.xlabel(\"t / s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t, tools.db(rir))\n",
    "plt.xlabel(\"t / s\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_impulse_response(ir, fs=44100, db=True):\n",
    "    L = len(ir)\n",
    "    time = np.arange(L) / fs * 1000\n",
    "    \n",
    "    plt.figure()\n",
    "    if db:\n",
    "        plt.plot(time, tools.db(ir))\n",
    "        plt.ylabel('Level / dB')\n",
    "    else:\n",
    "        plt.plot(time, ir)\n",
    "        plt.ylabel('Amplitude')\n",
    "    plt.xlabel('t / ms')\n",
    "    plt.grid()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def energy_decay_curve(ir):\n",
    "    \"\"\"Normalized energy decay curve (EDC) of the impulse response.\n",
    "    \"\"\"\n",
    "    L = np.zeros_like(ir)\n",
    "    L[-1] = ir[-1]**2\n",
    "    for n in range(len(ir)-2, -1, -1):\n",
    "        L[n] = L[n+1] + ir[n]**2\n",
    "    return L / L[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the energy decay curve can be interpreted as a convolution of $h^2(t)$ and a rectangular window.\n",
    "Therefore, it can be computed more efficiently using the fast convolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fast_edc(ir):\n",
    "    L = len(ir)\n",
    "    window = np.ones_like(ir)\n",
    "    L = signal.fftconvolve(ir**2, window)[L-1:]\n",
    "    return L / L[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, this will be efficient only if the length of the impulse response is sufficiently long."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import io\n",
    "import soundfile as sf\n",
    "import tools\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"http://legacy.spa.aalto.fi/projects/poririrs/wavs/omni.zip\"\n",
    "filename = \"s1_r1_o.wav\"\n",
    "zf = zipfile.ZipFile(tools.HttpFile(url))\n",
    "pori, fs = sf.read(io.BytesIO(zf.read(filename)))\n",
    "print(pori.shape, fs)\n",
    "\n",
    "# you can also just download and unzip the file manually:\n",
    "#pori, fs = sf.read(filename)\n",
    "\n",
    "assert pori.shape[1] == 2  # stereo IR\n",
    "pori = pori.sum(axis=1)\n",
    "fs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sampling frequencies of the input signal and the impulse response have to match!\n",
    "\n",
    "It's very easy to convert between sampling frequencies with [SoX](http://sox.sourceforge.net/), e.g. like this:\n",
    "\n",
    "    sox xmas.wav -r 48000 xmas48k.wav"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speech48k, fs48k = sf.read(\"data/xmas48k.wav\")\n",
    "assert fs48k == 48000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speech_pori = signal.fftconvolve(speech48k, pori)\n",
    "speech_pori = tools.normalize(speech_pori, np.max(np.abs(speech)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sf.write(\"data/xmas_pori.wav\", speech_pori, fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<audio src=\"data/xmas.wav\" controls></audio>\n",
    "[data/xmas.wav](data/xmas.wav)\n",
    "\n",
    "<audio src=\"data/xmas_pori.wav\" controls></audio>\n",
    "[data/xmas_pori.wav](data/xmas_pori.wav)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(len(pori)) / fs\n",
    "plt.figure()\n",
    "plt.plot(t, pori)\n",
    "plt.xlabel(\"t / s\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: use custom plotting function\n",
    "plt.figure()\n",
    "plt.plot(t, tools.db(pori))\n",
    "plt.xlabel(\"t / s\")\n",
    "plt.ylabel(\"Level / dB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_impulse_response(pori, fs=fs, db=True)"
   ]
  },
  {
   "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": "mycomac",
   "language": "python",
   "name": "mycomac"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}