{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## ThinkDSP\n", "\n", "This notebook contains solutions to exercises in Chapter 6: Discrete Cosine Transform\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": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "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", "In this chapter I claim that `analyze1` takes time proportional\n", "to $n^3$ and `analyze2` takes time proportional to $n^2$. To\n", "see if that's true, run them on a range of input sizes and time\n", "them. In IPython, you can use the magic command `%timeit`.\n", "\n", "If you plot run time versus input size on a log-log scale, you\n", "should get a straight line with slope 3 for `analyze1` and\n", "slope 2 for `analyze2`. You also might want to test `dct_iv`\n", "and `scipy.fftpack.dct`.\n", "\n", "I'll start with a noise signal and an array of power-of-two sizes" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "(16384,)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from thinkdsp import UncorrelatedGaussianNoise\n", "\n", "signal = UncorrelatedGaussianNoise()\n", "noise = signal.make_wave(duration=1.0, framerate=16384)\n", "noise.ys.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function takes an array of results from a timing experiment, plots the results, and fits a straight line." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import linregress\n", "\n", "loglog = dict(xscale='log', yscale='log')\n", "\n", "def plot_bests(ns, bests): \n", " plt.plot(ns, bests)\n", " decorate(**loglog)\n", " \n", " x = np.log(ns)\n", " y = np.log(bests)\n", " t = linregress(x,y)\n", " slope = t[0]\n", "\n", " return slope" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "PI2 = np.pi * 2\n", "\n", "def analyze1(ys, fs, ts):\n", " \"\"\"Analyze a mixture of cosines and return amplitudes.\n", "\n", " Works for the general case where M is not orthogonal.\n", "\n", " ys: wave array\n", " fs: frequencies in Hz\n", " ts: times where the signal was evaluated \n", "\n", " returns: vector of amplitudes\n", " \"\"\"\n", " args = np.outer(ts, fs)\n", " M = np.cos(PI2 * args)\n", " amps = np.linalg.solve(M, ys)\n", " return amps" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "def run_speed_test(ns, func):\n", " results = []\n", " for N in ns:\n", " print(N)\n", " ts = (0.5 + np.arange(N)) / N\n", " freqs = (0.5 + np.arange(N)) / 2\n", " ys = noise.ys[:N]\n", " result = %timeit -r1 -o func(ys, freqs, ts)\n", " results.append(result)\n", " \n", " bests = [result.best for result in results]\n", " return bests" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the results for `analyze1`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([ 64, 128, 256, 512, 1024, 2048, 4096])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ns = 2 ** np.arange(6, 13)\n", "ns" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "64\n", "270 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1000 loops each)\n", "128\n", "834 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1000 loops each)\n", "256\n", "3.63 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)\n", "512\n", "18 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)\n", "1024\n", "75.7 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)\n", "2048\n", "343 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n", "4096\n", "1.98 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" ] }, { "data": { "text/plain": [ "2.1523010349686165" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3RUdf7/8eeHEAgtoYQOIUBC70aaoAgWFCN2RVfXBqJi+/50RWXXtrbV77oiWHBtLIqsNImiqAj2QhESCAFCIIUAAUJCepvP7w/Y72FZcBMyk3tn5vU4x3PM9c6dl2TOvLh33vczxlqLiIiI29RzOoCIiMiJqKBERMSVVFAiIuJKKigREXElFZSIiLhSfacD/JbIyEgbHR3tdAwREfGhdevWHbDWtj5+uysLyhgTD8THxMSwdu1ap+OIiIgPGWPST7TdlZf4rLUJ1topERERTkcRERGHuLKgREREVFAiIuJKKigREXElFZSIiLiSCkpERFxJBSUiIq6kghIREVdSQYmISI15PJaf0w769DlcWVDGmHhjzJz8/Hyno4iIyHHyisu5de5arp7zE5t2++592pUFpZUkRETcKTErjwkzv+Pb7ft5cmJf+nYI99lzuXItPhERcRdrLfN+zuDJhGRaN2vIh1NHMqhzc58+pwpKRER+U1FZJY8sSWLphmzG9GzNi1cNokWTBj5/XhWUiIicVGpOAbfPW8+O/YXcf14P7hgTQ716pk6eWwUlIiIntGxjNtMXJdIoNIR/3DKMM2Ii6/T5VVAiIvJvyiqreOqTLcz9MZ24Li2Yde0Q2kWE1XkOFZSIiPyfrEPF3PneejZm5TN5dFf+ML4XoSHODHyroEREBIBVKTncu2ADHo/ltd+dxvh+7RzNo4ISEQlyVR7Li19sY9aqVHq3D+fV64YQHdnE6VgqKBGRYHagsIy75//KDzsOcnVcZx6f2Jew0BCnYwEqKBGRoLVmVy7T3l9PXnEFf7liAFfFdXY60r9RQYmIBBlrLX//difPfpZC5xaNePuOofTx4ZJFp0oFJSISRA6XVvDAhxtZsXkf4/u24y9XDiA8LNTpWCdUZwVljGkCvAKUA6utte/V1XOLiAhszs7njvfWs/tQCTMm9OaWUV0xpm5WhTgVtRpuN8a8ZYzJMcZsOm77eGPMVmNMqjFm+tHNlwELrbWTgYtr87wiIlIzC9ZkcOkrP1BaUcUHU4Zz6+huri4nqP0Z1DvALGDuvzYYY0KA2cC5QBawxhizDOgEJB3draqWzysiItVQUl7Fnz7axIfrshgVE8nfrhlEZNOGTseqlloVlLX2G2NM9HGbhwKp1to0AGPMB8BEjpRVJ2ADLv0eKhGRQLLzQBG3z1vH1n0F3D0ulnvGxRJSRwu9eoMvPoPqCGQe83MWMAyYCcwyxkwAEk72YGPMFGAKQFRUlA/iiYgEvk+T9vDAwkTqhxjevvF0xvRs43SkGvNFQZ2onq21tgi46b892Fo7B5gDEBcXZ72cTUQkoFVUeXj20xTe/G4nAzs355XrhtCxeSOnY50SXxRUFnDs3V6dgOyaHMAYEw/Ex8TEeDOXiEhA25tfyp3vr2dd+iFuHBnNwxf2pkF9//1ExRfJ1wCxxpiuxpgGwDXAspocwFqbYK2dEhER4YN4IiKB57vtB5gw81tS9hzm5UmDeezivn5dTlD7MfP5wI9AT2NMljHmFmttJTANWAFsAf5prd1c+6giInI8j8cyc+V2rn/rZ1o1bcBH00YRP7CD07G8orZTfJNOsn05sPxUj6tLfCIi/11uUTn3LtjAN9v2c8mgDjx9WX8aNwicBYJcef6nS3wiIr/t14xDXDTzW37acZCnLu3Hi1cPCqhyAq3FJyLiV6y1vPvDLp5avoW24WEsun0k/TsF5l/mVVAiIn6isKyS6YsS+ThxD+N6teGvVw0iorE7F3r1BlcWlD6DEhH5d9v2FTB13jp2HSjiwfG9uO3MbtTzo1UhToU+gxIRcbklv2Yxcdb3HC6p5L1bh3P7mO4BX07g0jMoERGB0ooqnvg4mfd/zmBo15bMmjSYNuFhTseqM64sKF3iE5Fgl5lbzO3vrWPT7sNMPas795/Xg/ohrrzo5TOu/L/VJT4RCWZfJu9jwsxvST9YzBs3xDH9gl5BV07g0jMoEZFgVFnl4YXPt/Ha1zvo1zGcV649jahWjZ2O5RgVlIiIC+QUlHLX+7/y885cJg2N4tH4PoSFhjgdy1EqKBERh/2UdpC75v9KQWkFf71qIJcN6eR0JFdwZUFpSEJEgoHHY3n9mzSeX5FCdKsmzLtlGD3bNXM6lmu48lM3DUmISKDLL65gyj/W8dxnKVzQvz3L7hqlcjqOK8+gREQCWVJWPne8v469+aU8Ft+H34+MxpjAv/G2plRQIiJ1xFrL+79k8PiyZCKbNmDBbSMYEtXC6ViupYISEakDxeWVzFiyicW/7mZ0bCQvXTOYlk0aOB3L1VRQIiI+tmN/IbfPW8f2nELuO6cH08bGEBIEa+nVlisLSlN8IhIoPk7M5sGFiTQMDWHuzUMZHdva6Uh+Q1N8IiI+UF7p4bFlm5n2/q/0bNeMT+4epXKqIVeeQYmI+LPdeSXc+d56NmTmccuorky/oBehQbiWXm2poEREvMRay4drs3jyk2SshVevG8IF/ds7HctvqaBERLxgd14JDy1O4ptt+xnatSXPXzGALq2aOB3Lr6mgRERq4V/3Nj2zPAWPtTwxsS+/G9YlKL7x1tdUUCIipygzt5gHFyXyw46DjOzeiucuH0DnlsH79Rje5sqC0pi5iLiZx2P5x0/pPPdZCvWM4elL+zNpaGctV+Rlriwoa20CkBAXFzfZ6SwiIsfadaCIPyxK5JeduZzZozXPXNafjs0bOR0rILmyoERE3KbKY3n7+5288PlWQkPq8ZcrBnDlaZ101uRDKigRkf8iNaeQPyzcyPqMPMb1asNTl/anXUSY07ECngpKROQkKqs8/P27nfz1i200Cg3hxasHcsmgjjprqiMqKBGRE9i2r4AHPtzIxqx8zu/blicv6UebZjprqksqKBGRY1RUeXj96x3MXJlK07D6vDxpMBcNaK+zJgeooEREjkrOPswDCzeyOfswEwa054mL+9KqaUOnYwUtFZSIBL3ySg+zV6Uye1UqzRuH8trvhjC+n9bQc5oKSkSCWlJWPg8s3EjK3gIuGdSBR+P70kLfdOsKriworSQhIr5WVlnFS19u5/Vv0ohs2oC/3xDHOX3aOh1LjuHKgtJKEiLiSxsy83jgw41szynkytM6MeOiPkQ0CnU6lhzHlQUlIuILpRVVvPjFNt74No224WG8c9PpjOnZxulYchIqKBEJCuvSc3ngw0TSDhQxaWhnHrqwN+FhOmtyMxWUiAS04vJKXlixjbd/2EmHiEbMu2UYo2IjnY4l1aCCEpGA9VPaQR5clEj6wWJuGNGFP4zvRdOGetvzF/pNiUjAKSqr5LnPUpj7YzpdWjXmgynDGd6tldOxpIZUUCISUL5PPcCDixLZnVfCzWd05f7ze9C4gd7q/JF+ayISEApKK3h6eQrzf8mgW2QTPrxtBHHRLZ2OJbWgghIRv7d6aw4PLU5i3+FSppzZjf85twdhoSFOx5JaUkGJiN/KL67gyU+SWbgui9g2TXnl9pEMjmrhdCzxEhWUiPilL5P38fCSJA4WlXPn2d25e1wsDevrrCmQqKBExK8cKirn8YTNLN2QTa92zXjz96fTv1OE07HEB1RQIuI3Ptu0lxlLN5FXXM4942K58+wYGtSv53Qs8ZE6KyhjTDfgESDCWntFXT2viPi/g4Vl/GnZZj5J3EPfDuHMvXkofTqEOx1LfKxaf/UwxrxljMkxxmw6bvt4Y8xWY0yqMWb6bx3DWptmrb2lNmFFJLhYa0nYmM25L37D55v3cv95PVh65xkqpyBR3TOod4BZwNx/bTDGhACzgXOBLGCNMWYZEAI8c9zjb7bW5tQ6rYgEjZyCUv64dBMrNu9jYKcInr9yOD3aNnM6ltShahWUtfYbY0z0cZuHAqnW2jQAY8wHwERr7TPARacayBgzBZgCEBUVdaqHERE/Za3low3ZPJawmeLyKqZf0ItbR3Wlfog+awo2tfmNdwQyj/k56+i2EzLGtDLGvAYMNsY8dLL9rLVzrLVx1tq41q1b1yKeiPibfYdLmTx3Lfcu2EC3yCYsv3s0U8/qrnIKUrUZkjAn2GZPtrO19iAwtRbPJyIBylrLh+uyePLjZCqqPMyY0JubzuhKSL0Tvc1IsKhNQWUBnY/5uROQXbs4Rxhj4oH4mJgYbxxORFwsO6+EhxYn8fW2/QyNbslfrhhAdGQTp2OJC9SmoNYAscaYrsBu4BrgWm+EstYmAAlxcXGTvXE8EXEfay3zf8nk6eVb8FjL4xf35frhXainsyY5qloFZYyZD4wBIo0xWcCj1to3jTHTgBUcmdx7y1q72WdJRSRgZOYWM31xIt+nHmRk91Y8d/kAOrds7HQscZnqTvFNOsn25cByryZCl/hEApXHY3nv53Se+TSFesbw1KX9uHZoFMborEn+kyuXOtIlPpHAs3VvAQ8tTmR9Rh6jYyN59vIBdGzeyOlY4mKuLCgRCRylFVXMXLmdOd+kEd4olL9eNZBLB3fUWZP8V64sKF3iEwkM327fzyNLNpGRW8wVp3Xi4Qt707JJA6djiZ9wZUHpEp+IfztQWMaTHyfz0YZsukU2Yf7k4Yzo3srpWOJnXFlQIuKfPB7LP9dm8synKZSUV3HPuFjuOLu7vkhQTokKSkS8IjWngIcXb+KXXbkM69qSpy7tT0ybpk7HEj/myoLSZ1Ai/qO0oopXVqXy6tc7aNygPn+5fABXxnXSEITUmisLSp9BifiHH1IP8MjSTew8UMRlgzvy8ITeRDZt6HQsCRCuLCgRcbfconL+/Ekyi9fvpkurxsy7ZRijYiOdjiUBRgUlItVmrWXhuiyeXr6FgtJKpp0dw7SxMYSFaghCvE8FJSLVsmN/IY8sSeKntFziurTg6cv66xtuxadcWVAakhBxj7LKKl5bncbsVamEhdbjmcv6c3VcZ606Lj7nyoLSkISIO/ycdpCHlySxY38R8QM78MeLetOmWZjTsSRIuLKgRMRZecXlPLM8hQVrM+ncshHv3HQ6Y3q2cTqWBBkVlIj8H2stSzfs5s8fbyGvpIKpZ3XnnnGxNGqgIQipeyooEQFg14EiZizdxHepBxjUuTnzLutP7/bhTseSIObKgtKQhEjdKa/08Ma3acxcuZ0GIfV4cmJfrh3WhRANQYjDXFlQGpIQqRtrd+Xy8JIktu0r5ML+7Xg0vi9twzUEIe7gyoISEd/KL67g2c9SmP9LBh2bN+LN38cxrndbp2OJ/BsVlEgQsdbyceIeHk9IJreojMmju3LvOT1o0lBvBeI+elWKBInM3GJmLN3E19v2M6BTBO/cdDr9OkY4HUvkpFRQIgGuosrDm9/t5G9fbiPEGB6N78MNI6I1BCGup4ISCWC/ZhziocVJpOwt4Lw+bXl8Yl/aRzRyOpZItaigRALQ4dIKXlixlX/8lE7bZmG8fv1pnN+3ndOxRGrElQWl+6BETo21ls827eWxhM3kFJTx+xHR3H9+T5pqCEL8kCtftboPSqTmdueV8Kelm1iZkkPfDuHMuT6OgZ2bOx1L5JS5sqBEpPoqqzy888Mu/vrFNqyFGRN6c+PIaOqH1HM6mkitqKBE/FhiVh4PLU5ic/ZhxvZqwxMT+9KpRWOnY4l4hQpKxA8VllXyv59v5d0fdhHZtCGvXjeE8f3aYYxGxyVwqKBE/Mznm/fy6LLN7D1cyu+GdeGB8T0JDwt1OpaI16mgRPzEnvwSHv1oM58n76NXu2bMvm4IQ6JaOB1LxGdUUCIuV+WxzP1xFy+s2EqVtUy/oBe3jOpKqIYgJMCpoERcbNPufB5ekkRiVj5n9mjNnyf2I6qVhiAkOKigRFyouLySF7/Yxlvf76JF4wbMnDSY+AHtNQQhQcWVBaWVJCSYfZWyjz8u3czuvBImDY1i+vheRDTWEIQEH1cWlFaSkGCUU1DKY8s2szxpL7FtmvLh1BGcHt3S6VgijnFlQYkEm8827eWhxYkUlVfxwPk9mTy6Gw3qawhCgpsKSsRBhWWVPL5sMx+uy6J/xwhevHoQMW2aOh1LxBVUUCIOWbsrl/v+uYHdh0qYdnYMd4+L1VmTyDFUUCJ1rKLKw0tfbueV1al0bNGIf942gjh91iTyH1RQInVox/5C7luwgcSsfK48rRN/iu9DMy1TJHJCKiiROmCtZd5P6Ty1fAuNQkN47XdDGN+vvdOxRFxNBSXiYzkFpfxhYSKrt+7nzB6teeGKAbQJD3M6lojrqaBEfGjF5r08tDiJorJKHr+4LzeM6KLVIESqSQUl4gOFZZU8mZDMgrWZ9OsYzt+uHkRMm2ZOxxLxKyooES9bl36I+xZsIOtQMXeM6c695/TQ+LjIKVBBiXhJRZWHl1duZ9aqVDo0b8SC27RUkUhtqKBEvCDt6Pj4xqx8Lh/Siccu1vi4SG3VWUEZYy4BJgBtgNnW2s/r6rlFfMVay3s/Z/DnT5IJCw3h1euGcEF/jY+LeEO1LowbY94yxuQYYzYdt328MWarMSbVGDP9t45hrV1qrZ0M3AhcfcqJRVxif0EZt7y7lhlLN3F6dEtW3HumyknEi6p7BvUOMAuY+68NxpgQYDZwLpAFrDHGLANCgGeOe/zN1tqco/8+4+jjRPzWF8n7mL4okcKySh6L78MNI6KpV0/j4yLeVK2CstZ+Y4yJPm7zUCDVWpsGYIz5AJhorX0GuOj4Y5gjN388C3xqrV1/sucyxkwBpgBERUVVJ55InSkqq+TJj5P5YE0mfdqH88E1g4htq/FxEV+ozWdQHYHMY37OAob9xv53AecAEcaYGGvtayfayVo7B5gDEBcXZ2uRT8Sr1mccGR/PyC3m9jHduU/j4yI+VZuCOtH1jJMWirV2JjCzFs8n4oiKKg8vf5XK7FWptAsP44PJwxnWrZXTsUQCXm0KKgvofMzPnYDs2sU5whgTD8THxMR443Aip+zY8fHLhnTksYv7Eq7xcZE6UZvrE2uAWGNMV2NMA+AaYJk3QllrE6y1UyIiIrxxOJEas9by/s8ZTJj5HbsOFjP72iH89apBKieROlStMyhjzHxgDBBpjMkCHrXWvmmMmQas4Mjk3lvW2s0+SypSR/YXlDF9USIrU3IYHRvJ81cMpF2EVh8XqWvVneKbdJLty4HlXk2ELvGJc75M3seDixIpKKvkTxf14caRGh8XcYorlzqy1iYACXFxcZOdziLBoaiskj9/soX5v2TQp304868ZRA+Nj4s4ypUFJVKXfj06Pp6eW8zUs7pz37mxNKwf4nQskaDnyoLSJT6pC5VVHmatSuXlr46Mj8+fPJzhGh8XcQ1XFpQu8Ymv7TxQxH0LNrAhM49LB3fk8YkaHxdxG1cWlIivWGv5YE0mTyQkExpieHnSYOIHdnA6loicgApKgsaBwjKmL0riyy37OCOmFS9cOZD2EY2cjiUiJ+HKgtJnUOJtK7ccGR8/XFrJjAm9ufmMrhofF3E5V650qZUkxFuKyyt5eEkSt7y7lsimDUmYNopbR3dTOYn4AVeeQYl4w4bMPO5bsIFdB4u47cxu/M95PTQ+LuJHVFAScCqrPMxetYOZX22nbbOGvH/rcEZ01/i4iL9RQUlAST9YxL0LNvBrRh4TB3XgiYn9iGik8XERf+TKgtKQhNSUtZYFazJ54uNk6tczzJw0mIs1Pi7i11xZULpRV2riYGEZ0xcn8UXyPkZ2PzI+3qG5xsdF/J0rC0qkulal5PDAwkQOl1RofFwkwKigxC+VlFfx1PJk5v2UQa92zZh361B6tQt3OpaIeJEKSvzOxqPj42kHipg8uiv/77yehIVqfFwk0LiyoDQkISdSWeXh1dU7eGnldlo3a8j7tw5jZEyk07FExEdcWVAakpDjpR88svr4+ow8Lh7YgScn9iOiscbHRQKZKwtK5F+qPJZ3f9jFC59vJaSe4aVrBjFxUEenY4lIHVBBiWtt3VvAg4sS2ZCZx5ierXnq0v501Pi4SNBQQYnrlFVWMXvVDl5dnUqzsFBeumYQFw/sgDEaHxcJJioocZV16bk8uCiJ1JxCLhnUgT/F96VlkwZOxxIRB6igxBUKyyp5/rMU5v6UTvvwMN6+6XTO7tnG6Vgi4iBXFpTGzIPLqpQcHlmSxJ7Dpfx+RDT3n9+Tpg1d+dIUkTrkyncBjZkHh4OFZTzxcTIfbcgmpk1TFk4dyWldWjgdS0RcwpUFJYHNWsvSDbt5IiGZwrJK7hkXyx1nd9eXCYrIv1FBSZ3KOlTMjKWbWL11P4OjmvPc5QPo0baZ07FExIVUUFInqjyWuT/u4vkVWwF4NL4PN4yIJkQrj4vISaigxOe27Ttyw+2vGXmc1aM1T13aj04tGjsdS0RcTgUlPlNWWcUrq3bwyupUmjasz4tXD+SSQR11w62IVIsKSnxiXfohpi9KZHtOIRMHdeBPF/WhVdOGTscSET+ighKvKiqr5PkVW3n3x120Dw/jrRvjGNurrdOxRMQPqaDEa1ZvzeGRJZvIzi/hhuFdeGB8L91wKyKnTO8eUmu5ReU8kbCZpRuy6d66CQunjuC0Li2djiUifs6VBaWljvyDtZZlG7N5PCGZwyUV3D02hjvHxuiGWxHxClcWlJY6cr/deSXMWJLEqq37Gdi5Oc9d3p9e7cKdjiUiAcSVBSXu5fFY/vFTOn/5LAWPhT9e1IcbR+qGWxHxPhWUVFtqTgEPLkpiXfohRsdG8vSl/encUjfciohvqKDkvyqv9PDq6h3MXpVK44Yh/O+VA7lsiG64FRHfUkHJb/o14xDTFyWxdV8B8QM78Gh8HyJ1w62I1AEVlJxQUVklL3y+lXd+2EW78DDe/H0c43rrhlsRqTsqKPkPX2/bz8OLk9idV8L1w7vwh/E9aRYW6nQsEQkyKij5P4eKynny42QW/7qbbq2b8OHUEZwerRtuRcQZKijBWktC4h4eX7aZ/JIK7hobw51nxxAWqhtuRcQ5Kqggl51Xwh+XbmJlSg4DO0Uw79Zh9G6vG25FxHkqqCDl8Vje+zmd5z7bSqXHw4wJvbnpjK664VZEXEMFFYRScwqZviiRtemHGBVz5IbbqFa64VZE3EUFFUTKKz28/vUOXv4qlUYNQnj+igFccVon3XArIq5UZwVljOkN3ANEAiutta/W1XMLbMjMY/qiRFL2FjBhQHsei+9L62a64VZE3KtaBWWMeQu4CMix1vY7Zvt44CUgBPi7tfbZkx3DWrsFmGqMqQe8UavUUm3F5ZX87+fbePv7nbRu1pA3bojj3D664VZE3K+6Z1DvALOAuf/aYIwJAWYD5wJZwBpjzDKOlNUzxz3+ZmttjjHmYmD60WOJj327fT8PL0kiM7eE64ZF8eAFvQjXDbci4ieqVVDW2m+MMdHHbR4KpFpr0wCMMR8AE621z3DkbOtEx1kGLDPGfAK8f6J9jDFTgCkAUVFR1Yknx8krLufJj7ewaH0W3SKbsGDKcIZ1a+V0LBGRGqnNZ1Adgcxjfs4Chp1sZ2PMGOAyoCGw/GT7WWvnAHMA4uLibC3yBR1rLZ8k7eGxZZvJK67gzrO7c9fYWN1wKyJ+qTYFdaLRr5MWirV2NbC6Fs8nv2FP/pEbbr/ckkP/jhHMvXkYfTrohlsR8V+1KagsoPMxP3cCsmsX5whjTDwQHxMT443DBbQqj2X+Lxk8+2kKlR4Pj1zYm5vOiKZ+SD2no4mI1EptCmoNEGuM6QrsBq4BrvVGKGttApAQFxc32RvHC0RVHsvHidm8/FUqqTmFnBHTiqcv7U+XVk2cjiYi4hXVHTOfD4wBIo0xWcCj1to3jTHTgBUcmdx7y1q72WdJBYDKKg8fbchm9qpU0g4U0bNtM2ZdO5gJ/dvrhlsRCSjVneKbdJLty/mNgYdTpUt8/6miysOS9buZtSqVjNxierVrxqvXDeH8vu2op/XzRCQAGWvdOygXFxdn165d63QMR5VXeli0PovZq1LJOlRCv47h3D02lnN6t1UxiUhAMMass9bGHb9da/G5VFllFf9cm8Wrq1LJzi9lYKcInpjYl7N7ttGlPBEJCq4sqGC+xFdaUcWCNZm8unoHew+XMjiqOU9f1p+zerRWMYlIUNElPpcoKa/i/V8yeP3rHeQUlHF6dAvuGdeDM2JaqZhEJKDpEp9LFZdX8t5PGbz+TRoHCssY3q0lf7tmECO6qZhEJLipoBxSWFbJP35M541v08gtKmdUTCR3jR2sNfNERI5yZUEF8mdQBaUVzD1aTHnFFZzZozX3jIvhtC4tnY4mIuIqriyoQFxJIr+kgne+38Wb36VxuLSSsb3acNfYGAZHtXA6moiIK7myoAJJXnE5b32/i7e/30lBaSXn9mnL3WNj6d8pwuloIiKupoLykdyict78Lo13f0insKyS8X3bcde4GPp2UDGJiFSHKwvKnz+DOlhYxhvf7mTuj7soqajiwv7tuWtsDL3a6asvRERqwpUF5Y+fQe0vKGPONzuY91MGpZVVxA/owLSxMfRo28zpaCIifsmVBeVP9h0u5fWv03jv53QqqjxcMqgjd5wdQ0ybpk5HExHxayqoU7Qnv4TXVu9g/ppMqjyWSwd35M6zY+gaqe9jEhHxBhVUDe3OK+GVVal8uDYLj7VccVon7hgTQ1Srxk5HExEJKCqoasrMLeaV1aksXJcFwJVxnbljTHc6tVAxiYj4gisLyk1TfOkHi5i9KpXF63dTzxgmDY1i6lnd6dC8kdPRREQCmisLyg1TfGn7C5m1KpWPNmRTv57h+hFduO3M7rSLCHMqkohIUHFlQTkpNaeAWV+lsmxjNg3q1+OmkdFMOasbbZqpmERE6pIK6qitewt4+avtfJK0h0ahIUw+sxuTR3cjsmlDp6OJiASloC+oLXsO8/JX21metJcmDUK4/azu3Dq6Gy2bNHA6mohIUAvagtq0O5+ZK7fzefI+mjWsz91jY7h5VFeaN1YxiaoqXwAAAAOWSURBVIi4QdAV1MbMPGau3M7KlBzCw+pz3zk9uPGMaCIahTodTUREjuHKgvLFmPn6jEPMXLmd1Vv307xxKPef14MbRkYTHqZiEhFxI2OtdTrDScXFxdm1a9fW6hhrd+Xy0srtfLv9AC2bNGDy6G5cP6ILTRu6sptFRIKOMWadtTbu+O0B/S595/vr+SRxD5FNG/Dwhb24blgXmqiYRET8QkC/Ww/v1orBnZtz3bAuNGoQ4nQcERGpgYAuqOuHd3E6goiInKJ6TgcQERE5ERWUiIi4kgpKRERcSQUlIiKupIISERFXcmVBGWPijTFz8vPznY4iIiIOcWVBWWsTrLVTIiIinI4iIiIOcWVBiYiIqKBERMSVXL1YrDFmP5Beg4dEADX54MoX+1dnn0jgQA2e1x/V9M/WV3ydw1vHr81xTuWxNXlMdffVa/+IYHjte/PYEUBza23r//gv1tqA+QeY4/T+1dxnrdN/Vm77XfhrDm8dvzbHOZXH1uQx1d1Xr33vvibcnMObx/6tYwXaJb4EF+xf02MGKrf8Ofg6h7eOX5vjnMpja/KY6u7rlt+509zy5+DLHN489kmP5epLfIHKGLPWnuC7T0QCnV77UhOBdgblL+Y4HUDEIXrtS7XpDEpERFxJZ1AiIuJKKigREXElFZSIiLiSCkpERFxJBeUCxphLjDFvGGM+Msac53QekbpgjOltjHnNGLPQGHO703nEfVRQPmKMecsYk2OM2XTc9vHGmK3GmFRjzHQAa+1Sa+1k4EbgagfiinhFDV/3W6y1U4GrAN0bJf9BBeU77wDjj91gjAkBZgMXAH2AScaYPsfsMuPofxfxV+9Qg9e9MeZi4DtgZd3GFH+ggvIRa+03QO5xm4cCqdbaNGttOfABMNEc8RzwqbV2fV1nFfGWmrzuj+6/zFo7EriubpOKP6jvdIAg0xHIPObnLGAYcBdwDhBhjImx1r7mRDgRHznh694YMwa4DGgILHcgl7icCqpumRNss9bamcDMug4jUkdO9rpfDayu2yjiT3SJr25lAZ2P+bkTkO1QFpG6ote9nBIVVN1aA8QaY7oaYxoA1wDLHM4k4mt63cspUUH5iDFmPvAj0NMYk2WMucVaWwlMA1YAW4B/Wms3O5lTxJv0uhdv0mrmIiLiSjqDEhERV1JBiYiIK6mgRETElVRQIiLiSiooERFxJRWUiIi4kgpKRERcSQUlIiKu9P8BqZtL8oVjmyoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "bests = run_speed_test(ns, analyze1)\n", "plot_bests(ns, bests)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimated slope is close to 2, not 3, as expected. One possibility is that the performance of `np.linalg.solve` is nearly quadratic in this range of array sizes.\n", "\n", "Here are the results for `analyze2`:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def analyze2(ys, fs, ts):\n", " \"\"\"Analyze a mixture of cosines and return amplitudes.\n", "\n", " Assumes that fs and ts are chosen so that M is orthogonal.\n", "\n", " ys: wave array\n", " fs: frequencies in Hz\n", " ts: times where the signal was evaluated \n", "\n", " returns: vector of amplitudes\n", " \"\"\"\n", " args = np.outer(ts, fs)\n", " M = np.cos(PI2 * args)\n", " amps = np.dot(M, ys) / 2\n", " return amps" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "64\n", "171 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 10000 loops each)\n", "128\n", "670 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1000 loops each)\n", "256\n", "2.58 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)\n", "512\n", "10.3 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)\n", "1024\n", "45.3 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)\n", "2048\n", "165 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)\n", "4096\n", "625 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n" ] }, { "data": { "text/plain": [ "1.9837162383836813" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3RVVeL28e8mkBBa6EQIIUAIIQREiDQdsYCCgiBiwzqOgjOjM7ZRQEDFAuqoPx3bYHfGDqGEIjZsiAo4kEaAQCihhZoEQurd7x9h5kUETLnJOffe57OWa3kP5577rHDM47lnn72NtRYRERG3qeN0ABERkRNRQYmIiCupoERExJVUUCIi4koqKBERcaW6Tgc4lZYtW9qoqCinY4iISA1atWrVXmttq+O3u7qgoqKiWLlypdMxRESkBhljtpxouyu/4jPGjDDGzMzNzXU6ioiIOMSVBWWtTbLWjgsLC3M6ioiIOMSVBSUiIuLKgtJXfCIi4sqC0ld8IiLiyoISERFRQYmIiCu5sqB0D0pERFxZULoHJSLibh6P5YdN+2r0M1xZUCIi4l4bdudz5T+Xc/XMH0jdXnPfdLl6qiMREXGPwpIyXlqayctfb6RhSF2eGtOT7m2b1NjnubKgjDEjgBHR0dFORxEREWD5xn08MCeFTXsPc9kZ7Zh8STdaNAqp0c90ZUFZa5OApISEhFudziIiEsgOFhTz+KK1fLQym/bNQ3nn5r6cE/OricdrhCsLSkREnGWtZf6aHTyyIJ0DBSXcNqgzf72gC6HBQbWWQQUlIiK/sG1/AQ/MTeWb9Xs4PSKMd27uR1wN3ms6GRWUiIgAUFrm4Y1lWTzz2XqCjOHBEXHcMCCKoDrGkTyuLCgNkhARqV3J2QeZMDuF9J15DO7Wmmkj42nbNNTRTK4sKA2SEBGpHYeLSnn60/W89X0WLRuF8PK1vRkaH44xzlw1HcuVBSUiIjXvi7W7mTovje0Hj3Bd/0juGxpLk/r1nI71PyooEZEAk5NXyMNJ6SxM2UmX1o2YddsAEqKaOx3rV1RQIiIBwuOxfLBiG9MXr6Wo1MO9F8Yw7pzOBNd156x3KigRkQCwYXc+ExNTWLnlAP07Nefxy3rQqVUjp2OdkisLSqP4RES840Tz543pE+GKQRC/xZUFpVF8IiLV98OmfUxKLJ8/b1SvtkweHkfLGp4/z5tcWVAiIlJ1BwuKmb4ogw9XbqN981Devrkvg2pp/jxvUkGJiPiJ4+fPGz+oE3deEFOr8+d5kwpKRMQPbNtfwOS5qXx9dP68t2/uS/e2vr0quQpKRMSH/Xf+vGc/24AxOD5/njepoEREfFRy9kEmJqaQtiOPC2JbM21UPO0cnj/Pm1RQIiI+5vj58166tjfDXDJ/nje5sqD0HJSIyIl9mbGbKXPL58+7tl/5/Hlhoe6ZP8+bXFlQeg5KROSXcvKPzp+X7O7587zJlQUlIiLlfjF/XomHe4bEMH6Qe+fP8yYVlIiIS2XmlM+ft2LzAfp1bM7jo3vQ2eXz53mTCkpExGUKS8p46auNvPxVJg2C6/Lk5T25IsE35s/zJhWUiIiL/LBpH5PmpLBpz2FG9mrLFB+bP8+bVFAiIi5w7Px5Ec18d/48b1JBiYg4yFpLUvJOpiWllc+fd04n/jq4Cw2C9etZPwEREYccO39eTz+ZP8+bVFAiIrWstMzDm8s288xn6zEGpg6P48aB/jF/njepoEREalFKdi4TEpP9dv48b3JlQWmqIxHxN4eLSnnms/W8uSyLFn48f543ubKgNNWRiPiTY+fPG9svkvv9eP48b3JlQYmI+INj58+Lbt2Ij28bwJl+Pn+eN6mgRES87L/z581YvJbCEg93D4lh/KBOhNT1zaXXnaKCEhHxonW78pk8N3Dnz/MmFZSIiBccLirluS828Pp3WTSuH7jz53mTCkpEpBqstSxJ2820pDR25BZyZUIEE4Z1o3nDYKej+TwVlIhIFW3bX8CD89P4MiOH2PDGPH/NGX6/iGBtUkGJiFRScamHV7/dxD++3EAdY3jg4m7cdFYU9YL8fxHB2qSCEhGphO837mXK3FQ27jnM0O7hTB0RR1vNBFEjVFAiIhWwJ7+IxxamM3f1Dto3D+XNm87kvNjWTsfyayooEZFTKPNY3vtxC08uWUdhSRm3nxfNn8+LJjRYzzTVNBWUiMhJpGTnMnluCmuycxnYuQWPjIrXM021SAUlInKcvMISnl6yjn/9sIXmDUN47upeXHp6Wz3TVMtUUCIiR1lrmb9mB48sWMu+w0Vc378D91zYVRO7OqTWCsoY0wl4AAiz1o6prc8VEamIjXsOMXVeKssy99GjXRhv3JRAz4imTscKaBUqKGPMG8BwIMdaG3/M9qHAc0AQ8Jq1dsbJjmGt3QT8wRgzq3qRRUS8p7CkjBeXZvLPrzcRUq8Oj4zszth+HbS6rQtU9ArqLeAF4J3/bjDGBAEvAkOAbGCFMWY+5WU1/bj332ytzal2WhERL1q6LocH56WxdX8Bo3q1ZdIl3WjduL7TseSoChWUtfYbY0zUcZv7AplHr4wwxnwAjLTWTqf8aqtKjDHjgHEAkZGRVT2MiMhJ7cw9wrSkdBan7qJTq4a8d0s/Bka3dDqWHKc696DaAduOeZ0N9DvZzsaYFsBjwBnGmIlHi+xXrLUzgZkACQkJthr5RER+obTMw1vfb+aZz9ZT5rHce2EMt56jdZrcqjoFdaIvaE9aKNbafcBt1fg8EZEqW7VlPw/MSSVjVz7ndW3Fw5fGE9migdOx5BSqU1DZQPtjXkcAO6oXp5wxZgQwIjo62huHE5EAduBwMTMWZ/Dhym2cFlafV67rzUXdw/VMkw+oTkGtALoYYzoC24GrgbHeCGWtTQKSEhISbvXG8UQk8Hg8llmrspm+eC15haXc+ruO3Dk4hoYhevzTV1R0mPn7wLlAS2NMNvCgtfZ1Y8ztwBLKR+69Ya1Nq7GkIiIVlLErj8lzUlm55QAJHZrx6GXxxIY3cTqWVFJFR/Fdc5Lti4BFXk0kIlJFxy673uTosutj+kRQR880+SRXXuvqHpSIVEb5suu7eDgpnZ25hVyV0J77h8Vq2XUf58qC0j0oEamobfsLmDovlaXr9hAb3ph/aNl1v+HKghIR+S1FpWW8+s0m/vFlJkF1DJMv6caNA7Xsuj9xZUHpKz4ROZXvM/cyeV4qm/YcZlh8+bLrp4Vp2XV/48qC0ld8InIiOfmFPL5wrZZdDxCuLCgRkWOVeSzv/riFp44uu37H+eXLrtevpymK/JkKSkRcLTn7IJPnppKsZdcDjisLSvegRCT3SAlPf1q+7HoLLbsekFxZULoHJRK4jl92/Yb+Hbhby64HJFcWlIgEpo17DjFlbirfb9xHz4gw3rzpTHpEhDkdSxyighIRx2nZdTkRFZSIOGppRg5T56eybf8RLbsuv+DKgtIgCRH/t+Ng+bLrn6Rp2XU5MWOte1dVT0hIsCtXrnQ6hoh4UUmZh7eWbebZz8uXXf/LBV245Xcdtex6ADPGrLLWJhy/3ZVXUCLin9J25HLPR2vI2JXP+bGtefjS7rRvrmXX5cRUUCJS48o8lle/3cTTn66jaYNgXrmuDxd1b6NnmuSUVFAiUqOyDxRwz0dr+DFrP0O7h/P46B5ap0kqRAUlIjXCWsu81TuYMjcVj7U8OaYnV/SJ0FWTVJgrC0qj+ER8W25BCQ/MTWFB8k76dGjGs1f2IrKF7jVJ5biyoDTVkYjv+j5zL/d8vIY9+UXce2EMtw3qTF0tIihV4MqCEhHfU1hSxt+XrOO177Lo1KohiX8aSM+Ipk7HEh+mghKRasvYlcedH6wmY1c+1/fvwKSLuxEarOeapHpUUCJSZR6P5fXvsnhqyTqahNbTCrfiVSooEamSHQePcM9Ha1i+aR9D4towY3QPWjQKcTqW+BEVlIhU2vw1O5g8J4VSj+WJy3twZUJ7DR8Xr1NBiUiF5R4pYeq8VOat3sEZkU159speRLVs6HQs8VOuLCg9ByXiPss37uOej1azO7+IuwbH8OfzNHxcapYrzy5rbZK1dlxYmFbSFHFaUWkZ0xetZexrPxBctw6zbhvAXwd3UTlJjXPlFZSIuMO6Xfnc+eFq1u7M45q+kUy+pBsNQ/RrQ2qHzjQR+RWPx/Lm95t54pMMGofU5bUbEhgc18bpWBJgVFAi8gu7cgu59+M1fJe5lwtiWzPj8p60aqzh41L7VFAi8j8Lk3cyaU4KxaUeHrssnrF9IzV8XByjghIR8gpLeGh+Gok/b+f0iDCevaoXnVo1cjqWBDgVlEiA+ylrP3d9uJqduUf4ywVduOP8aOpphJ64gApKJEAVl3p49vP1vPL1RiKbN+Dj2wbSp0Mzp2OJ/I8KSiQAZebk89cPVpO2I4+rz2zPlOFxGj4urqMzUiSAWGt5Z/kWHl+0loYhdfnn9X24qHu407FETsiVBaWpjkS8LyevkL/NSubr9Xs4t2srnhzTk9aN6zsdS+SkXFlQWvJdxLs+Sd3JxMQUjpSU8cjI7lzXv4OGj4vrubKgRMQ7DhWV8vD8ND5elU2PduXDx6Nba/i4+AYVlIifWrVlP3d9uIbsAwXcfl40f7mgC8F1NXxcfIcKSsTPlJR5eP6LDby4NJN2zUL5aPwAEqKaOx1LpNJUUCJ+ZOOeQ9z14WqSs3MZ0yeCB0fE0bh+PadjiVSJCkrED1hr+fePW3lsYTr16wXx8rW9GdbjNKdjiVSLCkrEx+3JL+L+2cl8mZHD77q05O9XnE6bJho+Lr5PBSXiwz5L382E2ckcKirloRFx3DAgijp1NHxc/IMKSsQHHS4q5dGF6bz/0zbiTmvCB1f3okubxk7HEvEqFZSIj/l56wHu/nA1W/YX8MdzO3PX4BgNHxe/pIIS8RGlZR7+8WUmLyzNJLxJfT64tT/9OrVwOpZIjVFBifiArL2HuevD1azedpDRZ7TjoZHdaaLh4+LnVFAiLmat5YMV25iWlE5w3Tq8MPYMhvds63QskVqhghJxqb2HipgwO4XP1+7m7Ojy4ePhYRo+LoFDBSXiQl9m7Oa+WcnkFZYyZXgcvx+o4eMSeGqtoIwxo4BLgNbAi9baT2vrs0V8RUFxKY8tXMu7P24lNrwx797Sn67hGj4ugalCY1ONMW8YY3KMManHbR9qjFlnjMk0xkw41TGstXOttbcCNwFXVTmxiJ9ave0gw5//jvd+2sq4czox7/azVE4S0Cp6BfUW8ALwzn83GGOCgBeBIUA2sMIYMx8IAqYf9/6brbU5R/998tH3iQhQWFLGs5+t59VvNxHepD7v3tKPgZ1bOh1LxHEVKihr7TfGmKjjNvcFMq21mwCMMR8AI62104Hhxx/DlC/fOQNYbK39+WSfZYwZB4wDiIyMrEg8EZ+1YvN+7puVTNbew4ztF8nEYbGafVzkqOrcg2oHbDvmdTbQ7xT73wEMBsKMMdHW2ldOtJO1diYwEyAhIcFWI5+IaxUUl/LkJ+t4e/lm2jUN5d1b+nFWtK6aRI5VnYI60ZCikxaKtfZ54PlqfJ6IX/g+cy/3Jyazbf8RbhoYxd8u6krDEA2oFTledf6ryAbaH/M6AthRvTjljDEjgBHR0dHeOJyIK+QXljB9cQbv/biVji0b8tH4AfTtqJVuRU6mOgW1AuhijOkIbAeuBsZ6I5S1NglISkhIuNUbxxNx2lfrcpiYmMLuvELGndOJu4fEUL9ekNOxRFytQgVljHkfOBdoaYzJBh601r5ujLkdWEL5yL03rLVpNZZUxAflFpTwyMJ0Zq3KpkvrRrz0x4GcEdnM6VgiPqGio/iuOcn2RcAiryZCX/GJf/gsfTcPzElh3+Fibj8vmjsuiCakrq6aRCrKlXdm9RWf+LL9h4t5aH4a89fsoNtpTXjjpjOJbxfmdCwRn+PKghLxRdZaFqXsYuq8VPIKS7h7SAy3DeqsxQRFqkgFJeIFOfmFTJ2bxidpu+gZEca7Y/oRG97E6VgiPs2VBaV7UOIrrLXMXb2dh5PSKSguY8KwWG45uyN1g3TVJFJdriwo3YMSX7Art5AH5qTwRUYOvSOb8uSY04lu3cjpWCJ+w5UFJeJm1lo+WrmNRxespcTjYcrwOG4aGEWQ1msS8SoVlEglZB8oYGJiCt9u2Eu/js15ckxPOrRo6HQsEb/kyoLSPShxG4/H8u6PW5ixOAOAR0bFc23fSK1yK1KDXFlQugclbrJ572Hun53Mj1n7+V2Xlkwf3YOIZg2cjiXi91xZUCJuUOaxvLksi79/uo56QXV48vKeXJEQQfnSZiJS01RQIieQmXOI+2at4eetB7kgtjWPXdaD8LD6TscSCSgqKJFjlJZ5ePXbLJ79fD0NgoP4v6t6MbJXW101iTjAlQWlQRLihIxdefzt42RStucyLD6caSPjadU4xOlYIgHLlQWlQRJSm4pLPbz81UZeWLqBJvXr8dK1vbm4x2lOxxIJeK4sKJHakro9l3s/XkPGrnxG9mrLgyO607xhsNOxRAQVlASootIynv9iA698vYkWDYN59YYEhsS1cTqWiBxDBSUB5z9bD/C3Wclk5hziij4RTL4kjrAG9ZyOJSLHUUFJwDhSXMYzn63j9e+yCG9Sn7dv7sugmFZOxxKRk3BlQWkUn3jbT1n7uW/WGjbvK+DafpFMGBZL4/q6ahJxM1cWlEbxibccLirlyU8yeHv5Fto3D+W9W/oxMLql07FEpAJcWVAi3rAscy/3z05m+8Ej3DQwivuGdqVBsE55EV+h/1rF7+QVljB9UQbv/7SVTi0b8tH4AZwZ1dzpWCJSSSoo8StL1+UwKTGF3XmFjD+nE3cNiaF+vSCnY4lIFaigxC8cLChm2oJ0En/eTkybRrx83Vn0at/U6VgiUg0qKPF5S9J2MXluKvsPF3PH+dHcfn40IXV11STi61RQ4rP2HSrioaR0ktbsIO60Jrx505nEtwtzOpaIeIkrC0rPQcmpWGtZkLyTB+enkV9Ywj1DYrjt3M7UC6rjdDQR8SJXFpSeg5KTyckvZMrcVJak7eb0iDCeHNOfruGNnY4lIjXAlQUlcjxrLYk/b2fagnSOlJQxcVgsfzi7I3V11STit1RQ4nrbDx5h8pwUlq7bQ0KHZjwxpiedWzVyOpaI1DAVlLhWmcfyzvLNPLVkHQBTh8dx48Aogupo+XWRQKCCElfK2JXHhNkprN52kEExrXh0VDztmzdwOpaI1CIVlLhKYUkZL3yZyStfb6RJaD2eu7oXl57eFmN01SQSaFRQ4ho/btrHxMQUNu09zOje7Zh8SZyWXxcJYCoocVzukRJmLC6f3LV981Deubkv52ghQZGAp4ISR32SupOp89LYe6iIW3/XkbuGxGhJDBEBVFDikN15hUydV/7AbdxpTXj9xjPpEaFpikTk/3NlQWmqI//l8Vje+2krTyzOoLjMw4SjD9xqmiIROZ4rC0pTHfmnzJxDTExMZsXmAwzs3ILHL+tBVMuGTscSEZdyZUGJfyku9fDK1xt54ctMQoODeGpMT8b0idDQcRE5JRWU1KhVWw4wMTGZ9bsPMbznaTw4ojutGoc4HUtEfIAKSmrEoaJSnvokg3d+2EJ4k/q8fmMCF3Rr43QsEfEhKijxui/W7mby3FR25RVy44Ao7r2oK41CdKqJSOXot4Z4zZ78Ih5OSmNB8k5i2jTixWsH0juymdOxRMRHqaCk2qy1fLwqm8cWruVIcRn3DIlh/KDOBNfV0HERqToVlFTL5r2HmTQnhe837qNvVHMeH92D6NZaq0lEqk8FJVVSUubhtW+z+L/P1xMcVIfHLovnmjMjqaO1mkTES1RQUmkp2bncPzuZ9J15XNS9DdNGxtOmSX2nY4mIn1FBSYUVFJfyzKfreWNZFi0bhfDKdb0ZGn+a07FExE+poKRCvlm/h0lzUsg+cISx/SK5f2gsYaH1nI4lIn5MBSWntP9wMY8uSCfxP9vp1KohH47rT79OLZyOJSIBQAUlJ2StZd7qHUxbkE7ekRLuOD+aP58XTf16QU5HE5EAUWsFZYzpBvwVaAl8Ya19ubY+Wypn2/4CJs9N5ev1e+jVvikzLu9BbHgTp2OJSICpUEEZY94AhgM51tr4Y7YPBZ4DgoDXrLUzTnYMa+1a4DZjTB3g1WqllhpR5rG8uSyLpz9dTx0DD42I4/oBUQRp6LiIOKCiV1BvAS8A7/x3gzEmCHgRGAJkAyuMMfMpL6vpx73/ZmttjjHmUmDC0WOJi6TvyGNiYjJrsnM5P7Y1j4yKp13TUKdjiUgAq1BBWWu/McZEHbe5L5Bprd0EYIz5ABhprZ1O+dXWiY4zH5hvjFkIvFfV0OI9hSVlPP/FBmZ+s4mw0Hr845ozGN7zNK3VJCKOq849qHbAtmNeZwP9TrazMeZcYDQQAiw6xX7jgHEAkZGR1Ygnv+X7jXuZlJjC5n0FjOkTwQMXd6NZw2CnY4mIANUrqBP9L7Y92c7W2q+Ar37roNbamcBMgISEhJMeT6out6CExxet5cOV24hs3oB//6EfZ3dp6XQsEZFfqE5BZQPtj3kdAeyoXhypSdZaFqXs4sH5aRwoKGb8oE7ceUEMocEaOi4i7lOdgloBdDHGdAS2A1cDY70RyhgzAhgRHR3tjcMJsDP3CFPmpvL52hzi2zXhrd+fSXy7MKdjiYicVIUW7DHGvA8sB7oaY7KNMX+w1pYCtwNLgLXAR9baNG+EstYmWWvHhYXpF2h1eTyWfy3fzJBnvuG7zL08cHE35v7pLJWTiLheRUfxXXOS7Ys4xYAHcdaG3flMSExh1ZYD/K5LSx4b1YPIFg2cjiUiUiGunOpIX/FVT1FpGS8t3chLX2XSMKQuT19xOqN7t9PQcRHxKa4sKGttEpCUkJBwq9NZfM3KzfuZkJhCZs4hRvZqy5ThcbRsFOJ0LBGRSnNlQUnl5ReW8MQnGfz7h620axrKm78/k/O6tnY6lohIlbmyoPQVX+V8mraLqfPS2J1fyM1ndeSeC2NoGOLKv1oRkQpz5W8xfcVXMXvyi3hofhoLU3YSG96YV67vQ6/2TZ2OJSLiFa4sKDk1ay2JP29n2oJ0jhSXce+FMYwf1Jl6QRV6akBExCeooHxM9oECJs1J5Zv1e+jToRlPXN6T6NaNnI4lIuJ1Kigf4fFY/vXDFp74JAOAhy/tzvX9O1BHazWJiJ9yZUFpkMQvZeYcYsLsZFZuOcA5Ma14/LJ4IprpgVsR8W+uLCgNkihXUuZh5jebeO7zDYQGB+mBWxEJKK4sKIGU7Fzum53M2p15XNLjNB66tDutGuuBWxEJHCoolyksKePZz9fz2rdZNG8YzCvX9WFofLjTsUREap0KykV+2LSPiYkpZO09zFUJ7Zl0cTfCGtRzOpaIiCNcWVCBNkgiv7CEGYszePfHrbRvHsq7t/TjrGitcCsigc2VBRVIgyS+zNjNA3NS2Z1XyC1nd+TuC2NoEOzKvxYRkVql34QO2XeoiGkL0pm3egcxbRrx0rUDOSOymdOxRERcQwVVy6y1zF+zg4eT0skvLOHOwV3407nRBNfVNEUiIsdSQdWinblHmDwnlS8ycji9fVOevLwnXcMbOx1LRMSVVFC1wOOxvL9iK9MXZVDq8TD5km78/qyOBGmaIhGRk3JlQfnTKL6svYeZMDuZH7P2M7BzC2aM7klkC01TJCLyW1xZUP4wiq+0zMPr32XxzGfrCa5bhycu78GVCe01TZGISAW5sqB8XfqOPO6fnUzK9lyGxLXh0VHxtGlS3+lYIiI+RQXlRUWlZbzwZSYvf7WRpg3q8eLY3lzcI1xXTSIiVaCC8pJVW/Zz36xkNu45zOje7ZhySRzNGgY7HUtExGepoKrpcFEpTy1Zx9vLN9M2LJS3fn8m53Zt7XQsERGfp4Kqhq/X72FSYgo7co9wQ/8O/G1oLI1C9CMVEfEGV/42dfsw84MFxTyyYC2zf86mU6uGfDx+AAlRzZ2OJSLiV1xZUG4dZm6tZXHqLqbOS+VAQQm3nxfN7edHU79ekNPRRET8jisLyo1y8gqZMi+VJWm7iW/XhLdv7kv3tmFOxxIR8VsqqN9greXjldk8sjCd4lIPE4bFcsvZHakbpMldRURqkgrqFLbuK2DinGSWZe6jb8fmzBjdg06tGjkdS0QkIKigTqDMY3lzWRZPf7qeoDqGR0fFM7ZvJHU0uauISK1RQR1n3a587p+dzOptBzk/tjWPjoqnbdNQp2OJiAQcFdRRxaUeXvoqkxeXZtIopC7PXd2LS09vq2mKREQcooICVm87yP2zklm3O59LT2/LgyPiaNEoxOlYIiIBLaALqqC4lGc+Xc8by7Jo3bg+r92QwOC4Nk7HEhERArigvs/cy4TEFLbuL2Bsv0gmDIulSf16TscSEZGjAq6gco+UMH3RWj5YsY2oFg14/9b+DOjcwulYIiJyHFcWVE3Nxfdp2i4mz01l76Eixp/TiTsHxxAarGmKRETcyJUF5e25+PbkF/FQUhoLk3cSG96Y125MoGdEU28cWkREaogrC8pbrLXM+c92pi1Ip6CojHsvjGH8oM7U0zRFIiKu57cFZa3lz+/9zKKUXfSObMqTY3oS3bqx07FERKSC/LagjDGcFd2SvlHNuX5AFEGapkhExKf4bUEBXNuvg9MRRESkinQzRkREXEkFJSIirqSCEhERV1JBiYiIK6mgRETElVRQIiLiSiooERFxJRWUiIi4kgpKRERcyVhrnc5wUsaYPcCWSrwlDMh1eP+K7NMS2FuJz/VFlf3Z1pSazuGt41fnOFV5b2XeU9F9de6XC4Rz35vHDgOaWmtb/epPrLV+8w8w0+n9K7jPSqd/Vm77u/DVHN46fnWOU5X3VuY9Fd1X5753zwk35/DmsU91LH/7ii/JBftX9pj+yi0/h5rO4a3jV+c4VXlvZd5T0X3d8nfuNLf8HGoyhzePfdJjuforPn9ljFlprU1wOodIbdO5L5Xhb1dQvmKm0wFEHKJzXypMV1AiIuJKuoISERFXUkGJiIgrqaBERMSVVFAiIuJKKigXMMaMMsa8aoyZZ4y50Ok8IrXBGNPNGPOKMcAdcpAAAAGDSURBVGaWMeaPTucR91FB1RBjzBvGmBxjTOpx24caY9YZYzKNMRMArLVzrbW3AjcBVzkQV8QrKnner7XW3gZcCejZKPkVFVTNeQsYeuwGY0wQ8CIwDIgDrjHGxB2zy+Sjfy7iq96iEue9MeZS4Dvgi9qNKb5ABVVDrLXfAPuP29wXyLTWbrLWFgMfACNNuSeAxdban2s7q4i3VOa8P7r/fGvtQODa2k0qvqCu0wECTDtg2zGvs4F+wB3AYCDMGBNtrX3FiXAiNeSE570x5lxgNBACLHIgl7icCqp2mRNss9ba54HnazuMSC052Xn/FfBV7UYRX6Kv+GpXNtD+mNcRwA6HsojUFp33UiUqqNq1AuhijOlojAkGrgbmO5xJpKbpvJcqUUHVEGPM+8ByoKsxJtsY8wdrbSlwO7AEWAt8ZK1NczKniDfpvBdv0mzmIiLiSrqCEhERV1JBiYiIK6mgRETElVRQIiLiSiooERFxJRWUiIi4kgpKRERcSQUlIiKu9P8AFFzg6SK/VCgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "bests2 = run_speed_test(ns, analyze2)\n", "plot_bests(ns, bests2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results for `analyze2` fall in a straight line with the estimated slope close to 2, as expected.\n", "\n", "Here are the results for the `scipy.fftpack.dct`" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import scipy.fftpack\n", "\n", "def scipy_dct(ys, freqs, ts):\n", " return scipy.fftpack.dct(ys, type=3)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "64\n", "7.02 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 100000 loops each)\n", "128\n", "8.33 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 100000 loops each)\n", "256\n", "8.3 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 100000 loops each)\n", "512\n", "9.42 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 100000 loops each)\n", "1024\n", "13.3 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 100000 loops each)\n", "2048\n", "24.4 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 10000 loops each)\n", "4096\n", "41.9 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 10000 loops each)\n" ] }, { "data": { "text/plain": [ "0.4112565850438262" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXgV9b3H8fc3G2FN2MKWIEvY10AAxVJwxY26K4hWrUXctba23lpre2trtVdbrV5Fr16vqARQVFTcalnEupBA2AlEEBK2hC0hhJDl/O4fCW0aA+SEJHOWz+t58jw5c+bMfBLHfPidmfMbc84hIiISaCK8DiAiIlIbFZSIiAQkFZSIiAQkFZSIiAQkFZSIiASkKK8DNLYOHTq4Hj16eB1DRESOISMjY49zrmPN5SFfUD169CA9Pd3rGCIicgxmtrW25XqLT0REApIKSkREApIKSkREApIKSkREApIKSkREApIKSkREApIKSkREApIKSkREApIKSkRE6mVv0RGydh1stO2roERExG9lFT5ue205U174kkNHyhtlHyooERHx28PvreOrLft48KIBtGzWOLPmqaBERMQvaV9v4/++2Mq0cT25NCWx0fYTlAVlZhPM7DMze87MJnidR0QkXGRs3ceD76xhXJ8O/OK8/o26rzoXlJlFmtkKM3uvvjszs5fMLM/M1tTy3HlmlmVm2WZ2/wk25YAiIBbIrW8eERGpu50Fh5k+cznd4pvz9JQRREU27hjHnzcO7wbWA21qPmFmCcBh59zBasuSnXPZNVZ9GXgaeKXG6yOBZ4BzqCycZWY2H4gEHqmxjR8BnznnFptZJ+AJYKofP4eIiPippKyC6TMzOFxazuvTxhDXIrrR91mn+jOzROBC4H+Oscp44B0zi61afxrwVM2VnHNLgH21vH40kO2c2+ycKwXSgIudc6udcxfV+MpzzvmqXrcfaHaMzJPM7PmCgoK6/IgiInIMzjl+OW81q3IL+PPVw+nbqXWT7Leu47O/AD8HfLU96ZybC3wIpJnZVCpHOVf5kaMbkFPtcW7VslqZ2WVmNgOYSeWIrLZM7zrnbo6Li/MjhoiI1PTi0i3MW7Gde8/py7mDOjfZfk/4Fp+ZXQTkOecyjndBgnPuMTNLA54FejvnivzIYbVt8jj7mgfM82P7IiJSD59tyucPC9Zz3qDO3HFGcpPuuy4jqNOBH5jZt1S+9Xammb1acyUzGwcMBt4CHvIzRy6QVO1xIrDDz22IiEgD2rr3EHe8voI+Ca15/KphRETUNpZoPCcsKOfcfzjnEp1zPYDJwN+dc9dWX8fMUoAXgIuBG4F2ZvawHzmWAX3MrKeZxVTtZ74frxcRkQZUdKScaa+kYwYv/DC10T6MezwNdY1gC+BK59w3VRcwXA9srbmSmc0CvgD6mVmumd0E4JwrB+4APqLySsE5zrm1DZRNRET84PM57p2dyTf5h3jmmhF0b9/Ckxx+VaJzbhGwqJbln9d4XEbliKrmelOOs+0FwAJ/8oiISMN78tNNfLxuN7++aCCnJ3fwLEdQziQhIiKN48M1O3ny001cMTKRG0/v4WkWFZSIiACwYVch985ZybCkeB6+ZDBmTXtRRE0qKBERYf+hUqa9kk6rZlE8f91IYqMjvY7k3zkoEREJPeUVPu6YtZzdBUdIm34qndrEeh0JUEGJiIS9Rz7YwOfZe3nsiqGM6N7W6zj/pLf4RETC2BsZuby4dAs3jO3BValJJ35BE1JBiYiEqcycA/zyrdWM7d2eBy4c4HWc71BBiYiEobzCEqbPTCehdTOevmYE0Y18b6f60DkoEZEwc6S8gumvZlB4uJx5t42lXcsYryPVSgUlIhJGnHM8+PYaVmw7wLNTRzCgy3fuQRswAm9MJyIijeaVL7YyJz2XO89M5vwhXbyOc1wqKBGRMPGPb/bwn++t4+wBCfzk7L5exzkhFZSISBjI2VfM7a8tp2eHlvz56uFNfm+n+lBBiYiEuOLSyns7VfgcL/wwldax0V5HqhNdJCEiEsKcc/xs7ko27j7I/944mp4dWnodqc40ghIRCWHPLMxmwepd3H9+f8b37eh1HL+ooEREQtTf1u3m8U82csnwrkwb18vrOH5TQYmIhKDsvIPcMzuTQV3b8MfLh3p+b6f6UEGJiISYgsNlTHslg9joCJ6/LjUg7u1UH7pIQkQkhFT4HHfNWkHu/mJen3YqXeObex2p3lRQIiIh5LGPNrB4Yz5/uHQIo3q08zrOSdFbfCIiIeKdzO3MWLyZa0/tzjVjunsd56SpoEREQsDq3AJ+/sYqRvdsx68vGuR1nAahghIRCXL5B49w88x02reM4b+njiAmKjT+tOsclIhIECst93HbaxnsLy7ljVvG0qFVM68jNRgVlIhIEPvtu2tZ9u1+npqSwuBucV7HaVBBOQ40swlm9pmZPWdmE7zOIyLihVe/3MprX23jlvG9+cGwrl7HaXAnLCgzizWzr81spZmtNbPf1ndnZvaSmeWZ2ZpanjvPzLLMLNvM7j/BphxQBMQCufXNIyISrL7eso/fzF/LhH4duW9iP6/jNIq6vMV3BDjTOVdkZtHAUjP7wDn35dEVzCwBOOycO1htWbJzLrvGtl4GngZeqb7QzCKBZ4BzqCycZWY2H4gEHqmxjR8BnznnFptZJ+AJYGodfg4RkZCw/cBhbn01g+7tWvDk5BQig+DeTvVxwoJyzh0drQBEV325GquNB241swuccyVmNg24FLigxraWmFmPWnYzGsh2zm0GMLM04GLn3CPARceJtx+o9YygmU0CJiUnJx/n5SIiweVwaQXTZ6ZTWu7j+R+mEtc8OO7tVB91OgdlZpFmlgnkAZ84576q/rxzbi7wIZBmZlOpHOVc5UeObkBOtce5VcuOlecyM5sBzKRyRPYdzrl3nXM3x8WF1klDEQlfzjl+8eYq1u4o5C+Th5Oc0MrrSI2qTlfxOecqgOFmFg+8ZWaDnXNraqzzWNXI51mgt3OuqLZtHUNt49Oao7Tq+5oHzPNj+yIiQW/Gks3MX7mD+yb246wBnbyO0+j8uorPOXcAWAScV/M5MxsHDAbeAh7yM0cukFTtcSKww89tiIiErEVZeTz64QYuHNqF2yb09jpOk6jLVXwdq0ZOmFlz4GxgQ411UoAXgIuBG4F2ZvawHzmWAX3MrKeZxQCTgfl+vF5EJGRtzi/izlkr6N+5DX+6Ijjv7VQfdRlBdQEWmtkqKovkE+fcezXWaQFc6Zz7xjnnA64HttbckJnNAr4A+plZrpndBOCcKwfuAD4C1gNznHNr6/tDiYiEioMlZUx7JZ3oyAiev24kLWLCZ36FulzFtwpIOcE6n9d4XEbliKrmelOOs40FwIIT5RERCRc+n+OetEy+3VvMqzeNIaldC68jNamgnElCRCQcPPHJRj7dkMdDkwZyWu/2XsdpciooEZEA9P6qnTy9MJvJo5K47tRTvI7jCRWUiEiAWbejkJ/NXcmI7vH89uJBYXNRRE0qKBGRALLvUCnTXkknrnk0z107kmZRkV5H8kz4XA4iIhLgyioq7+2UX3SEudNPI6FNrNeRPKURlIhIgPj9++v5cvM+/njZEIYlxXsdx3MqKBGRADB72TZe/se3/Ph7PblsRKLXcQKCCkpExGMZW/fzq7fXMK5PB+4/v7/XcQKGCkpExEO7Ckq45dUMusQ1569TUoiK1J/lo3SRhIiIR0rKKu/tVHyknNd+PIb4FjFeRwooKigREQ845/jlW6tZmVvAjOtG0rdTa68jBRyNJUVEPPDS598yb/l27jm7DxMHdfY6TkBSQYmINLGlm/bw+/fXMXFQJ+46s4/XcQKWCkpEpAlt3XuI219fTp+E1jx+1XAiIsJzGqO6UEGJiDSRoiPlTHslHYDnfziSVs10GcDxqKBERJqAz+f46ZxMsvOKeOaaEZzSvqXXkQKeCkpEpAk89fdNfLR2Nw9cOJDv9engdZygoIISEWlkH63dxV/+tonLRyTyo9N7eB0naKigREQa0abdB7l3dibDEuP4/aWDw/beTvWhghIRaSQlZRXc/vpymsdEMuO6VGKjw/feTvWhS0hERBrJw++vY+PuIv7vR6PpHBfe93aqD42gREQawYdrdvHql9uYNq4n4/t29DpOUFJBiYg0sB0HDvOLN1cxpFsc903U7TPqSwUlItKAKnyOe2ZnUlbh46kpKcRE6c9sfekclIhIA3pmYTZfb9nH41cOo2cHfRj3ZKjaRUQaSPq3+3jy001cPLwrl43o5nWcoKeCEhFpAAWHy7g7LZOu8bE8fIk+79QQ9BafiMhJcs7xy3mr2V1YwtxbTqN1bLTXkUKCRlAiIidpTnoO76/eyb3n9iWle1uv44QMFZSIyEnIzjvIb+av4/Tk9tzy/d5exwkpKigRkXoqKavgzlmZNI+J5AndfLDB6RyUiEg9/fGDDazfWciL16fSqY2mMmpoGkGJiNTDp+t38/I/vuWGsT04a0Anr+OEJBWUiIifdheWcN8bqxjQpQ33n6+pjBqLCkpExA8+n+PeOZkcLq3gr1NSdAuNRqSCEhHxw4wlm/k8ey8PTRpIckIrr+OENBWUiEgdrdi2n8c/zuLCIV24elSS13FCngpKRKQOCkvKuCttBZ3axPKHy4ZoKqMmoMvMRUROwDnHg2+vYfv+w8yZfhpxzTWVUVPQCEpE5ATmLd/OO5k7uOfsvqT2aOd1nLChghIROY7N+UU8+M4aRvdsx+1nJHsdJ6wEZUGZ2QQz+8zMnjOzCV7nEZHQVFru4+60TGKiInhy8nAiNZVRkzphQZlZkpktNLP1ZrbWzO6u787M7CUzyzOzNbU8d56ZZZlZtpndf4JNOaAIiAVy65tHROR4/uvjLFZvL+DRy4fSJa6513HCTl0ukigHfuqcW25mrYEMM/vEObfu6ApmlgAcds4drLYs2TmXXWNbLwNPA69UX2hmkcAzwDlUFs4yM5sPRAKP1NjGj4DPnHOLzawT8AQwtQ4/h4hInS3emM/zSzZz7andmTios9dxwtIJR1DOuZ3OueVV3x8E1gM172U8HnjHzGIBzGwa8FQt21oC7KtlN6OBbOfcZudcKZAGXOycW+2cu6jGV55zzlf1uv1As9pym9kkM3u+oKDgRD+iiMi/yT94hJ/OyaRvp1b86sKBXscJW36dgzKzHkAK8FX15c65ucCHQJqZTaVylHOVH5vuBuRUe5zLd0uweo7LzGwGMJPKEdl3OOfedc7dHBcX50cMEQl3Pp/jp3NXcrCknL9OGaGpjDxU589BmVkr4E3gHudcYc3nnXOPmVka8CzQ2zlX5EeO2s48umOt7JybB8zzY/siInXy0udbWLIxn99dMph+nVt7HSes1WkEZWbRVJbTa1XlUNs644DBwFvAQ37myAWqzxuSCOzwcxsiIidldW4Bj364gXMHduLaMd29jhP26nIVnwEvAuudc08cY50U4AXgYuBGoJ2ZPexHjmVAHzPraWYxwGRgvh+vFxE5KYeOlHNX2grat2zGo5cP1VRGAaAuI6jTgeuAM80ss+rrghrrtACudM59U3UBw/XA1pobMrNZwBdAPzPLNbObAJxz5cAdwEdUXoQxxzm3tt4/lYiInx6av5Zv9x7iL5OH07ZljNdxhDqcg3LOLaX2c0TV1/m8xuMyKkdUNdebcpxtLAAWnCiPiEhDeydzO29k5HLXmcmc2qu913GkSlDOJCEi0lC27S3mgbfWMPKUttx1Vh+v40g1KigRCVtlFT7uTFuBGTw5eThRkfqTGEh0uw0RCVt//mQjK3MO8Mw1I0hs28LrOFKD/rkgImHp8+w9PLv4GyaPSuLCoV28jiO1UEGJSNjZW3SEn8zOpFeHlvx6kqYyClR6i09Ewopzjp+/sYoDxWX8742jaBGjP4OBSiMoEQkr//ePb/l0Qx7/cUF/BnXVXJ2BTAUlImFj3Y5C/rBgA2f2T+CGsT28jiMnoIISkbBQXFrOnbOWE98imj9doamMgoHefBWRsPC799axec8hXr1pDO1b1XobOQkwGkGJSMh7f9VOZn2dwy3je3N6cgev40gdqaBEJKTl7i/m/nmrGJYUz73n9PU6jvhBBSUiIau8wsc9aZk4B09NHk60pjIKKjoHJSIh66m/Z5O+dT9PTh7OKe1beh1H/KR/TohISPpq816e/vsmLhvRjYuHd/M6jtSDCkpEQs6B4lLumZ1J93Yt+M+LB3sdR+pJb/GJSEhxzvGLN1exp+gI8249nVbN9GcuWGkEJSIh5bWvtvHR2t38fGJ/hiRqKqNgpoISkZCRtesgv3tvHeP6dOCm7/X0Oo6cJBWUiISEkrIK7pq1gtaxUTx+1TAiIjSVUbDTm7MiEhJ+//56snYf5OUbR5HQOtbrONIANIISkaD30dpdzPxyKz/+Xk8m9EvwOo40EBWUiAS1nQWH+cWbqxjcrQ33ndfP6zjSgFRQIhK0KnyOe9IyKS338dTkFJpFRXodSRqQzkGJSND674XZfLVlH3+6Yii9OrbyOo40MI2gRCQoZWzdx18+3cQPhnXlipGJXseRRqCCEpGgU3C4jLtmZdI1PpaHLx2su+OGKL3FJyJBxTnHL+etZndhCXNvOY02sdFeR5JGohGUiASVOek5vL96Jz85py8p3dt6HUcakQpKRIJGdl4Rv5m/jrG923PL+N5ex5FGpoISkaBQUlbBnbNWEBsdwZ+vHk6kpjIKeToHJSJB4dEPN7B+ZyH/88NUOrXRVEbhQCMoEQl4f9+wm//9/FtuGNuDswd28jqONBEVlIgEtLzCEn42dxX9O7fm/vP7ex1HmpAKSkQCls/n+MmcTIpLy3n6mhRiozWVUThRQYlIwJqxZDOfZ+/loUmDSE5o7XUcaWIqKBEJSJk5B3j84ywuGNKZyaOSvI4jHlBBiUjAOVhSxl2zVtCpTSyPXDpUUxmFKV1mLiIB58G315C7v5g5008jroWmMgpXGkGJSECZtzyXtzN3cPdZfUnt0c7rOOIhFZSIBIwtew7x4NtrGN2zHXecmex1HPGYCkpEAkJpuY+7Zq0gKjKCv2gqI0HnoEQkQPzXx1ms3l7Ac9eOoGt8c6/jSABQQYmIp5xzvLh0C88v2czUMd05b3AXryNJgFBBiYhnikvLuf/N1cxfuYNzB3biwYsGeh1JAogKSkQ8sXXvIabPzCBr90Hum9iPW8f3JkLnnaQaFZSINLmFWXncPWsFZsbLN45mfN+OXkeSAKSCEpEm4/M5nl6YzZ//tpEBndsw47qRJLVr4XUsCVAqKBFpEoUlZdw7eyV/W7+bS1O68YdLh9A8RrOTy7GpoESk0W3afZDpMzPYtq+Y30wayPVje2h+PTkhFZSINKoFq3fys7kraRETxWs/HsOYXu29jiRBQgUlIo2ivMLHnz7OYsbizaR0j+fZqSPpHBfrdSwJIiooEWlw+w6VctesFSzN3sPUMd359aSBNIvS+SbxjwpKRBrUmu0FTJ+ZQX7RER67fChX6WaDUk8qKBFpMG9k5PLAW6tp3zKGudNPY1hSvNeRJIipoETkpJWW+3j4/XW88sVWxvZuz1+npNC+VTOvY0mQU0GJyEnJKyzhtteWk751Pzd/vxc/n9iPqEjdyUdOngpKROot/dt93PracopKyvnrlBQmDevqdSQJISooEfGbc46ZX27lP99dR2Lb5rx60xj6dW7tdSwJMSooEfFLSVkFD7y1hjeX53JW/wSeuHo4cc2jvY4lIUgFJSJ1lru/mFtezWDN9kLuObsPd53ZR7fIkEajghKROlm6aQ93zlpOuc/x4vWpnDWgk9eRJMSpoETkuJxzzFiymcc+3EByQitmXJdKzw4tvY4lYUAFJSLHVHSknPvmruSDNbu4cGgXHrt8KC2b6c+GNA0daSJSq835RUyfmcE3+UU8cMEAfjyup26RIU1KBSUi3/HJut3cOzuT6KgIXr1pDGOTO3gdScKQCkpE/qnC53jybxt56u/ZDE2M49lrR9ItvrnXsSRMqaBEBICC4jLunr2CRVn5XDkykd9dMpjYaN0iQ7yjghIR1u8sZPrMDHYWHOb3lw7mmtHddb5JPKeCEglz72Ru5/43V9OmeRSzp5/GiO5tvY4kAqigRMJWWYWPP36wgReXbmF0j3Y8PTWFhNa6JbsEDhWUSBjaU3SE219bzldb9nHD2B48cOEAonWLDAkwKiiRMJOZc4BbX81gf3Epf756GJemJHodSaRWKiiRMJL29TZ+/c5aEto0481bxzKoa5zXkUSOSQUlEgaOlFfwm/lrmfV1DuP6dOCvU1KIbxHjdSyR41JBiYS4nQWHueXV5azMOcDtZ/Tm3nP6EalbZEgQCKqCMrMJwO+AtUCac26Rp4FEAtyXm/dyx+vLKSnz8dy1IzlvcGevI4nUWZNdtmNmL5lZnpmtqbH8PDPLMrNsM7v/BJtxQBEQC+Q2VlaRYOec48WlW5j6P18R1zyat28/XeUkQacpR1AvA08DrxxdYGaRwDPAOVQWzjIzmw9EAo/UeP2PgM+cc4vNrBPwBDC1CXKLBJXi0nL+Y95q3sncwcRBnfivK4fROla3ZJfg02QF5ZxbYmY9aiweDWQ75zYDmFkacLFz7hHgouNsbj/QrDFyigSzrXsPMX1mBlm7D3LfxH7cNqG3piySoOX1OahuQE61x7nAmGOtbGaXAROBeCpHY8da72bgZoDu3bs3SFCRQLcwK4+7Z63AzHj5xtGM79vR60giJ8Xrgqrtn3buWCs75+YB8060Uefc88DzAKmpqcfcnkgo8PkczyzM5om/bWRA5zbMuG4kSe1aeB1L5KR5XVC5QFK1x4nADo+yiASdwpIyfjpnJZ+s282lKd34w6VDaB6jW2RIaPC6oJYBfcysJ7AdmAxc420kkeCwafdBps/MYNu+Yn4zaSDXj+2h800SUpqsoMxsFjAB6GBmucBDzrkXzewO4CMqr9x7yTm3tqkyiQSrD1bv5GdzV9I8JorXfjyGMb3aex1JpME15VV8U46xfAGwoKlyiASzCp/jTx9l8dzib0jpHs+zU0fSOU63yJDQ5PVbfCJSB/sOlbJkYz6vf72Nr7fsY+qY7vx60kCaRel8k4QuFZRIAPL5HGt3FLIwK4+FWXlk5hzAOejQKobHLh/KVaOSTrwRkSCnghIJEAWHy1i6aQ8Ls/JYlJXPnqIjmMGwxHjuOasvZ/TvyOCucURoolcJEyooEY8458jafZCFG/JZmJVHxtb9VPgccc2jGd+3I2f078j3+3SkfStNmiLhSQUl0oQOHSnn8+w9LMzKZ1FWHjsLSgAY2KUNt4zvxRn9EhieFE+Ubr8uooISaUzOOTbvOcTCDZVv2329ZR+lFT5aNYvie8kduOfsjozvm6Ar8URqoYISaWAlZRV8sXkvizbksTArn237igHok9CKG07vwYR+HUk9pR0xURoliRxPyBaUmU0CJiUnJ3sdRcJAzr5iFmVVFtI/vtlDSZmP2OgITu/dgWnf78WEvh01P56In8y50J5LNTU11aWnp3sdQ0JMabmP9G/3VV0Gnk92XhEAp7RvwRn9EjijfwJjerYjNlqfUxI5ETPLcM6l1lwesiMokYa2q6CkapSUx9JNezhUWkFMZARjerXjmtHdOaN/Aj07tPQ6pkjIUEGJHEN5hY8VOQdYWHUuaf3OQgC6xsVySUo3zuiXwGm929Oymf43EmkM+j9LpJo9RUdYnFX5uaQlG/MpLCknMsJIPaUt95/fnzP7J9AnoZVmDRdpAiooCWs+n2PV9oKqy8DzWLW9AOegY+tmTBzUmTP6J/C9Ph1oExvtdVSRsKOCkrBzoLiUJZv2sGhDHos35rP3UClmkJIUz71n9+WM/gkM7NJGUwqJeEwFJSHPOcf6nQer5rirnFLI56Bti6NTCiUwrk9H2rWM8TqqiFSjggpT63cWsmXPISLMiDCIjLDK7yOqHpthZlXL+bfvK19T7fHR1x59LsKItH9/LtIMi+Bf31ffZ9X2G9LBkjI+z977z6vudhceAWBItzjuOCOZCf0TGJYYT6RGSSIBSwUVRgpLypifuYPZy3JYvb3A6zj/xqpKsbLg/lWCR0vuX4VZ7bljFB7AN/lFlFU4WjeL4vt9OzKhX0fG9+tIQmtNKSQSLFRQIc45R/rW/aR9ncP7q3dQUuajf+fW/PYHgxjdsx3Ogc+5qq/KO7Y656jwVT4++lzlcv75fW3PVS7/17Z8zuHz/Wu7te3n3/Z59HW+o9utnu27z9X2mqPPHX3rbuQpbYnWxKsiQUkFFaL2FB1h3vJc0pblsDn/EK2aRXHZiEQmj0piSLc4XSYtIgFPBRVCKnyOzzblM3tZDp+s2025z5F6SltuvaI3Fw7tQosY/ecWkeChv1ghIHd/MXPSc3kjPYcdBSW0axnDjaf34OpRSSQntPY6nohIvaigglRpuY9P1u0mbdk2lmbvAWBcn4786qKBnD2gk27lICJBL2QLKlRvt5Gdd5DZy3J4c/l29h0qpWtcLHed2YcrUxNJbKvbOYhI6AjZgnLOvQu8m5qaOs3rLCeruLSc91btZM6yHNK37icqwjhnYCeuHpXEuD4d9VkeEQlJIVtQwc45x6rcAtKW5fDuyh0UHSmnV8eW/PKC/lw2IpEOrZp5HVFEpFGpoALMgeJS3l6xnbRlOWzYdZDY6AguHNKVyaOTSD2lrS4PF5GwoYIKAD6f48ste5m9LIcP1uyitNzHkG5xPHzJYH4wvKtm0haRsKSC8lBeYQlzM3KZk57D1r3FtI6NYvKoJK5KTWJwtziv44mIeEoF1cTKK3wsysonbVkOC7PyqPA5xvRsxz1n9+H8wV2IjY70OqKISEBQQTWRrXsPMSc9h7npueQdPEKHVs2YNq4XV6Um0qtjK6/jiYgEHBVUIyopq+CjtbuYvSyHf3yzlwiDCf0SuHpUEmf2T9AkpiIix6GCagQbdhWS9nUOb63YTsHhMhLbNuen5/TlitREusQ19zqeiEhQUEE1kKIj5by7cgdpy3JYmXOAmMgIzh3UicmjujO2d3vdPlxExE8qqJPgnGP5tgPMXraN969HB5AAAAOBSURBVFbtpLi0gr6dWvHgRQO5NKWbbiEuInISVFD1sO9QKfOW5zJ7WQ6b8opoERPJpKFduXp0EilJ8fowrYhIA1BB1ZHP5/j8mz2kLcvh47W7KKtwDE+K54+XDeGiYV1p1Uy/ShGRhqS/qiews+Awc9MrP0ybu/8w8S2iufbUU7h6VBL9O7fxOp6ISMhSQR3Hox9uYMbib/A5OD25PT8/rz/nDuykD9OKiDSBkC2ohrgf1LDEeG6bkMxVqUl0b697LYmINCVzznmdoVGlpqa69PR0r2OIiMgxmFmGcy615nJNZSAiIgFJBSUiIgFJBSUiIgFJBSUiIgFJBSUiIgFJBSUiIgFJBSUiIgFJBSUiIgFJBSUiIgEp5GeSMLN8YGsdV48DCvzYvD/r12XdE63TAdhTx/0FG39/98G0/4badn23E+jHNYTuse31cd2YGRryuI53znX8zjPOOX1VfQHPN9b6dVn3ROsA6V7/jgLldx9M+2+obdd3O4F+XFetE5LHttfHdWNmaIrjWm/x/bt3G3H9uqzr7/5Didc/e2Puv6G2Xd/t6Lj2TiD87I2VodGP65B/iy+UmFm6q2VCRZFgp2NbaqMRVHB53usAIo1Ex7Z8h0ZQIiISkDSCEhGRgKSCEhGRgKSCEhGRgKSCEhGRgKSCCmJmdomZvWBm75jZuV7nEWkIZjbAzJ4zszfM7Fav84h3VFABxsxeMrM8M1tTY/l5ZpZlZtlmdj+Ac+5t59w04Abgag/iitSJn8f1eufcLcBVgD4bFcZUUIHnZeC86gvMLBJ4BjgfGAhMMbOB1Vb5VdXzIoHqZfw4rs3sB8BS4NOmjSmBRAUVYJxzS4B9NRaPBrKdc5udc6VAGnCxVXoU+MA5t7yps4rUlT/HddX6851zY4GpTZtUAkmU1wGkTroBOdUe5wJjgDuBs4E4M0t2zj3nRTiReqr1uDazCcBlQDNggQe5JECooIKD1bLMOeeeAp5q6jAiDeRYx/UiYFHTRpFApLf4gkMukFTtcSKww6MsIg1Fx7UclwoqOCwD+phZTzOLASYD8z3OJHKydFzLcamgAoyZzQK+APqZWa6Z3eScKwfuAD4C1gNznHNrvcwp4g8d11Ifms1cREQCkkZQIiISkFRQIiISkFRQIiISkFRQIiISkFRQIiISkFRQIiISkFRQIiISkFRQIiISkP4f8gMgzYDaeiwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "bests3 = run_speed_test(ns, scipy_dct)\n", "plot_bests(ns, bests3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This implementation of dct is even faster. The line is curved, which means either we haven't seen the asymptotic behavior yet, or the asymptotic behavior is not a simple exponent of $n$. In fact, as we'll see soon, the run time is proportional to $n \\log n$.\n", "\n", "The following figure shows all three curves on the same axes." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUVfr48c9JJ4UUUighjRI6EUIvgg1UEBBXxI4iurv63b7rftdd/brruuu6uz/bqnRBxN6ioK7SQ2+hSAspJCQkgfQySWbm/P64IYaQQAgzmcnkeb9eeTlz59xzz40DD+fcc56jtNYIIYQQzsbN0Q0QQgghmiIBSgghhFOSACWEEMIpSYASQgjhlCRACSGEcEoejm6APYSGhuqYmBhHN0MIIUQL7Nmz56zWOqzxcZcMUDExMezevdvRzRBCCNECSqnMpo7LEJ8QQginJAFKCCGEU3KpAKWUmq6UWlhSUuLopgghhLhKLhWgtNZJWusFgYGBjm6KEEKIq+RSAUoIIYTrkAAlhBDCKUmAEkII4ZQkQAkhhGg1i9V+WzZJgGqHli9fzuOPP26z+o4ePcqYMWPw9vbmxRdftFm9QgjXlXmugv9ZvY9ff5Bit2u4VCYJpdR0YHrv3r0d3ZR2JSQkhJdffplPP/3U0U0RQji5/DITr3yXyuqdp/BwVzw8PhatNUopm1/LpXpQ7WGa+cyZMxk+fDgDBw5k4cKFAPj7+/OHP/yBoUOHMnr0aPLy8gBISkpi1KhRXHPNNdxwww31x88rKysjNjaW2tpaAEpLS4mJiSEzM5OEhIT6H3d3dzIzMykoKGD27NmMGDGCESNGkJycDEB4eDgjRozA09OzDX8TQoj2pNRUyz++Psq1L2xg9c5T3DWyJ5t+M5nfTOlnl+AELtaDaqn/SzrM9zmlNq1zQPfOPD194GXLLV26lJCQEKqqqhgxYgSzZ8+moqKC0aNH89xzz/Hb3/6WRYsW8dRTTzF+/Hi2b9+OUorFixfzwgsv8M9//rO+roCAACZNmsSXX37JzJkzeffdd5k9ezbR0dHs378fgNdee42NGzcSHR3N3XffzS9+8QvGjx/PqVOnmDJlCkeOHLHp70EI4VpMtRZWbMvgPxtOUlxZy/Sh3fnVjX2JCfWz+7U7ZIBypJdffplPPvkEgKysLE6cOIGXlxfTpk0DYPjw4fz3v/8FIDs7mzlz5pCbm0tNTQ2xsbEX1Td//nxeeOEFZs6cybJly1i0aFH9Z8nJySxevJjNmzcD8O233/L999/Xf15aWkpZWRkBAQF2u18hRPtktlj5cE82/+/bE5wpNXFt3zB+MyWeQT3aboSqQwaolvR07GHDhg18++23bNu2DV9fXyZNmoTJZMLT07O+i+zu7o7ZbAbgiSee4Je//CW33XYbGzZs4JlnnrmoznHjxpGRkcHGjRuxWCwMGjQIgNzcXB5++GE+//xz/P39AbBarWzbto1OnTq1zQ0LIdodrTVfHTrDP745RlpBBQk9g/j3nATG9OrS5m1xqWdQzq6kpITg4GB8fX05evQo27dvv2z5Hj16APDWW281W+7+++9n7ty5zJs3D4Da2lruvPNO/v73v9O3b9/6cjfddBOvvvpq/fvzw4BCCAGQnHqWma8l8+NVe3FTijfvG84nPxnrkOAEEqDa1NSpUzGbzQwZMoQ//vGPjB49+pLln3nmGX70ox8xYcIEQkNDmy13zz33UFRUxNy5cwHYunUru3bt4umnn66fKJGTk8PLL7/M7t27GTJkCAMGDOCNN94A4MyZM0RGRvKvf/2Lv/zlL0RGRlJaattndEII53Uwu4T7luzgnsU7KCir5h93DOHrn09kysCudpsA0RJKa/stsnKUxMRE3ZE2LPzwww/57LPPWLlypaObIoRoR9IKyvnnN8f58mAuwb6e/HRyb+4dHY2Pp3ubtkMptUdrndj4uNM/g1JK+QH/AWqADVrrVQ5uklN54oknWLt2LWvWrHF0U4QQ7cSZEhMvfXeC93dn4e3hxv9c34dHJsQS4ONcS00cEqCUUkuBaUC+1npQg+NTgZcAd2Cx1vpvwO3Ah1rrJKXUe4AEqAZeeeUVRzdBCNFOFFfW8PrGkyxPzsCqNfeNjuank3sTFuDt6KY1yVE9qOXAq8CK8weUUu7Aa8CNQDawSyn1ORAJHKwrZmnbZgohRPtXVWNh2dZ03thwkrJqM7MSevCLG/vSM8TX0U27JIcEKK31JqVUTKPDI4FUrXUagFLqXWAGRrCKBPZziUkdSqkFwAKAqKgo2zdaCCHamVqLlfd2ZfHydyfIL6vm+n7h/HpKPP27dXZ001rEmZ5B9QCyGrzPBkYBLwOvKqVuBZKaO1lrvRBYCMYkCTu2UwghnJrVqvnyYC7//OYYGecqSYwO5rV7hjEiJsTRTbsizhSgmprLqLXWFcC8FlUgyWKFEB2Y1ppNJ87ywldHOZxTSnxEAEseSOS6fuEOnS7eWs60Diob6NngfSSQcyUVtIdksbZg6+02Vq1axZAhQxgyZAhjx44lJcV+6fOFEPax71QRcxdt54GlOympquXfc4ay5mcTuL5/RLsMTuBcPahdQB+lVCxwGrgLuPtKKpAeVOvExsayceNGgoODWbt2LQsWLGDHjh2ObpYQogVS88v4x9fH+PpwHl38vHhm+gDmjorC26Nt1zLZg0N6UEqp1cA2IF4pla2UelhrbQYeB74GjgDva60PX0m97aEH5YzbbYwdO5bg4GAARo8eTXZ2dlv9OoQQrXS6uIrffJDCTf/eRHLqOX55Y182/nYyD46LdYngBI6bxTe3meNrgFavOG1xD2rtk3Dm4KXLXKmug+Hmv122mLNvt7FkyRJuvvlm2/5uhBA2U1hRw3/Wp7JieyZomDculp9O7k2In5ejm2ZzzjTEd9W01klAUmJi4iOObktznHm7jfXr17NkyRK2bNlin5sXQrRaRbWZJVvSWbgpjcoaM7OHRfLzG/vSI8hBuxPUVsH3n4HZBMMftMslXCpAtVgLejr24MzbbRw4cID58+ezdu1aunRxTOZiIcTFasxWVu88xSvrTnC2vIabBkTwmynx9Ilw0D5uZ0/A7mWQ8g5UFUHP0RKgWsLZJ0nYe7uNP/7xj8Dlt9v4zW9+AxjbbSQkJHDq1Cluv/12Vq5ceUF5IYTjWK2az1JO86//HiersIpRsSEsvL8fw6KC274x5ho4mmQEpozN4OYB/aZB4kMQO9Ful3WmaeZXzdknSTjrdhvPPvss586d4yc/+QkJCQkkJl6UVFgI0Ua01qw7msctL2/mF++l0NnHk7ceGsm7C0a3fXAqTIf/Pg3/HgAfPgTFmXD9n+AX38Odb0HctWDHKeyy3YYLkO02hHANuzMK+ftXR9mVUURMF19+dVM8tw7uhptbG65jstTC8a9g91I4uQ6UO8TfDMPnQa/rwM32/Zp2u93GlXD2IT57kO02hGj/jp4p5cWvj/HtkXzCArz5y8xBzBnRE0/3NhzkKs6CvSuMn/Iz0LkHTPpfGHYfdO7edu1oQHpQQgjhIFmFlfz7v8f5ZP9p/L09+PGkXjw4NgZfrzbqO1gtkPqt0Vs68Q1oDX1uNHpLfW4C97ZpR4foQQkhRHtwtryaV9elsmpHJm5KsWBiHD++thdBvm20lqnsDOxdCXvfgpIs8AuH8b+EYfdDcHTbtKEFJEAJIUQbKTPVsmhzOos3p1FttnJnYk9+dn0fugb62P/iViukbzB6S0fXgLZA3CS46S/Q71Zwd67ddMHFAlRHfAYlhHB+ploLq3ac4rX1qRRW1HDr4G788qa+9Arzt//FK87Cvrdhz3IoSodOITDmp8bapS697H/9qyDPoIQQwk5qzFY+3Xeal747weniKib0CeU3U+IZEhlk3wtrDZnJRm/p+8/BWgvR44xnSwNuAw/n2uJdnkE5gZdffpnXX3+dhIQECgoKOHv2LL///e/Jzc1lwYIF+PrabvvljIwMpk2bxqFDh1p8zjPPPIO/vz+//vWvm61z69at3H33FSWZF6LDKTPVsnrnKZZuyeBMqYkhkYG8cMcQxvVufj2jTVQWQsq7sGcZnD0OPoEwYr7RWwrvZ99r24EEqDb0n//8h7Vr15KXl8fvfve7+oSuMTEx3HvvvTYNUPaQkZHBO++8IwFKiGbklZpYmpzOO9tPUVZtZkxcF/42ezDX9g2z355MWkP2LqO3dPgTIzde5AiY+ToMmAlezv33yqVIgGojjz32GGlpadxyyy0cP34cf39/EhISmDdvHjk5OUyePJnQ0FDWr1+Pv78/jz76KOvXryc4OJh3332XsLAwFi1axMKFC6mpqaF3796sXLkSX19f8vLy6usHeP311+ne/Yd1C2lpacyePZuFCxcyYsSIC9r13HPPsWLFCnr27ElYWBjDhw8HIDU1lccee4yCggLc3d354IMPePLJJzly5AgJCQk88MAD/OIXv2i7X6AQTux4XhkLN6Xx2f7TWKyaWwZ3Y8HEOPsO5ZlK4MD7Rvqh/MPgFQAJ90DiPGN3BRfgUs+gGkySeOTEiRPNlvv7zr9ztPCoTa/dL6Qfvxv5u0uWiYmJYffu3Rw6dIgXX3yRL7744oLj59MZKaV4++23ueeee3j22WfJz8/n1Vdf5dy5c/WJXJ966ikiIiJ44oknmDNnDmPGjOHnP/85FouF8vJyioqKmDZtGh999BF33XUXy5YtIyEh4YL27NmzhwcffJAdO3ZgNpsZNmwYjz32GL/+9a8ZNWoUTz75JLNmzcJkMmG1Wtm5c+cF7RaiI9NaszO9kDc3pbHuaD4+nm7MSezJw+PjiOpix17L6b3GEN7BD6G2EroNNXLiDboDvNtg0oUddIhnUO1hu42WcHNzY86cOQDce++93H777QAcOnSIp556iuLiYsrLy5kyZQoA69atY8WKFYCRDT0wMJCioiIKCgqYMWMGH330EQMHDrzoOps3b2bWrFn1Q4u33XYbYGyEePr0aWbNmgWAj08bTIEVop2wWDXfHD7DG5vSSMkqJsTPi1/c0Jf7xkTbb0+m6nI49KHRW8rdD56+MGi2EZh6DLPPNZ2ASwWolrpcT8fZnB+7fvDBB/n0008ZOnQoy5cvZ8OGDZc8LzAwkJ49e5KcnNxkgGpYd0Ou1KsWwlZMtRY+2JPNks1pZJyrJLqLL3+ZOYg7hkfi42mnHWzPHDJ6SynvQU0ZhA+AW16EIXcaEyBcnEtlM2+vAgICKCsrq39vtVr58MMPAXjnnXcYP348YPRsunXrRm1tLatWraovf/311/P6668DYLFYKC0tBcDLy4tPP/2UFStW8M4771x03YkTJ/LJJ59QVVVFWVkZSUlJAHTu3JnIyEg+/fRTAKqrq6msrLyonUJ0BEUVNbz07QnG/W0df/z0EIG+Xrx+zzDW/WoS946Otn1wqq2C/e/A4hvhjXFGxod+t8JD38CPt8LIRzpEcIIO2oNyNgsWLODmm2+mW7durF+/Hj8/Pw4fPszw4cMJDAzkvffeA+DPf/4zo0aNIjo6msGDB9cHi5deeokFCxawZMkS3N3def311+nWrRsAfn5+fPHFF9x44434+fkxYsQI5s+fz5o1axg2bBhz5swhISGB6OhoJkyYUN+mlStX8uijj/KnP/0JT09PPvjgA4YMGYKHhwdDhw7lwQcflEkSwqVlFVayeHMa7+3OwlRr5fp+4SyYGMfI2BD7zMgrOG70lva/A6Zi6NIbpvwVhs4F3xDbX68dcKlJEue194W6/v7+lJeXO7oZQnRIB7NLeHPTSdYczMXdTTEzoQcLJsbZZwdbczUcqdsIMHMLuHlC/+nGs6WY8Xbda8mZdIhJEpLqSAjRGlprNh4v4M2NaWxLO0eAtwePTIxj3thY++TJO3fSSD20fxVUnoOgaLjhGUi4F/zDbH+9dkp6UEKIDqvWYiUpJYeFm9I4eqaMrp19eGh8DHNHRhHgY+PkqZZaOLbG6C2lrf9hI8DEhyBusl02AmwvOkQPSgghWqLMVMu7O7NYmpxObomJ+IgA/vmjoUwf2h0vDxsHiuJTsOct2LcSyvOgcyRM/gNccx907mbba7kYCVBCiA4jr9TEsuQMVu3IpMxkZnRcCH+9fTCTbJ2KyFwDx9cagenkOuNYn5uM3lKfG8HNTtPSXYwEKCGEy0vNN1IRfbLPSEV08+BuPGqPVERnTxhbpqeshooCY9v0a38L19wLQVG2vVYHIAFKCOGStNbsyihi4aaTfHvESEU0d2QU822diqi2Cr7/zAhMmck/PFsa9gD0vl56S1dBApQQwqVYrJr/fn+GNzamsb8uFdHPb+jD/WNibJuK6MxBYwjvwPtQXQLBsXD900bC1oAI212nA5MAJYRwCaZaCx/uyWZxXSqiqBBf/jxjIHcM70knLxv1YkylcOgj2PsW5OwDd29jA8Bh90P0+A49E88eJEAJIdq1oooaVm7P5K2tGZyrqGFoZCD/uWcYUwZ2xd3NBhMfzu+3tPctOPSxkUE8fABM/buRE6+DZnloC04foJRSccAfgECt9R2Obo8QwjlkFVayZEs67+3KoqrWwnV1qYhG2SoV0fndafeugIIj4OkHg+8wni31GN5hsjw4kl0DlFJqKTANyNdaD2pwfCrwEuAOLNZa/625OrTWacDDSqkP7dlWIUT7cOh0CW9uSmPNwVzcFMyoS0XU1xapiKxWyNhs9JaOJIGlxghG01+GQbeDtx3SHYlm2bsHtRx4FVhx/oBSyh14DbgRyAZ2KaU+xwhWzzc6/yGtdb6d2yiEcHJaazadOMvCTSdJTjVSEc0fH8u8cTZKRVR2xkg7tHclFKUb2cKHzzOeLXUddPnzhV3YNUBprTcppWIaHR4JpNb1jFBKvQvM0Fo/j9HbahWl1AJgAUBUlKw3EMIV1FqsfHEghzc3GqmIIjp78/ub+zF3VBSdrzYVkcUMqd8aQ3jHvwJtgZgJMPl/jYStnp1scxOi1RzxDKoHkNXgfTYwqrnCSqkuwHPANUqp39cFsotorRcCC8HIxWe75goh2lp5tZl3d55i6ZZ0ckpM9I3w58UfDeU2W6QiKso00g7tWwVlOeAXDmOfMHpLXXrZ5gaETTgiQDX1ZLHZgKK1Pgc81qKKJZu5EO1afqmJZVszWLU9k1KTmVGxITw3azCT4q8yFZG5Go5+afSW0jYYx3rfALe8AH2ngruNE8MKm3BEgMoGejZ4Hwnk2KJirXUSkJSYmPiILeoTQrSN1PxyFtWlIjJbrdw8qBsLJsYxtOdVpiIqOPZD6qHKcxDYEyY9aaQeCoy0TeOF3TgiQO0C+iilYoHTwF3A3baoWHpQQrQfWms2nzjLsuR01h8rwMfTjTkjejJ/QizRXfxaX3FNJXz/qZHlIWs7uHlA/C3G9PBekyX1UDti1/2glFKrgUlAKJAHPK21XqKUugX4fxgz95ZqrZ+z5XVlPyghnFdVjYWP92WzPDmDE/nlhPp7c9/oaO4dHUUXf+/WV5ybYgSlgx9AdamxZfqw+40t0/3DbXcDwuYcsh+U1npuM8fXAGtsfT3pQQnhvHKKq1ixLZN3d52iuLKWQT068687h3LrkG54e7SyV2MqgYMfGuuWclPAwwcGzDB6S9FjZTFtOyc76goh7EZrzd5TxSxNTuerQ2fQWjN1UFfmjYslMTq4dRMftIasHcazpcOfGKmHIgbD8AeMTA+dgm1/I8KuZEddIUSbqTFbWXsol6XJGaRkFdPZx1hYe9+YaCKDW7nVRcXZH1IPnT0GXv5GLrxhD0D3a6S35IJcKkDJEJ8QjnWuvJrVO0+xYlsm+WXVxIX58ecZA7l9WCR+3q3468ZqhfQNRlA68gVYayFyJNz2KgycBd7+Nr8H4TxkiE8IcdWOnill2ZYMPtl/mhqzlYl9w3hoXAwT+4Th1pqM4qU5xkLafSug+JQxbDd0LlxzH0QMsP0NCIeSIT4hhE1ZrJp1R/NZlpzO1pPn8PF040fDI5k3Lobe4a1Iqmoxw4lvjAkPJ74BbYXYicYmgP2mgacNcu6JdsWlApQM8Qlhf2WmWj7Ync1b2zLIPFdJ90Afnry5H3eN6EmQbyt2rC05DbuXGD2m8jPgHwHjfg7D7oOQOJu3X7QfMsQnhGiRzHMVLN+awQe7symvNjM8OpiHxsUyZWAEHu6tyI+XmwJbX4XDHxu9pT43GeuW+kwBd5f6t7O4DBniE0JcMa0129LOsXRLBt8dzcNdKaYN6ca8cbGtS0NktULqf2HrK8a+S17+MHIBjHoMgqNtfwOiXXOpACVDfELYhqnWwuf7c1ianM7RM2WE+Hnx+OTe3Ds6mojOrXgWVGuCA+/BtteMKeIB3eHGZ40p4p2uMt+ecFkyxCeEqJdXauLt7Zms2nGKwooa+nUN4KFxsdyW0B0fz1Zke6g4C7uWwK5FUFEAXQfDmCeMKeIerXheJVySDPEJIZqVklXMsuR0vjiQi0VrbugfwbxxMYyJ69K6bA9nTxi9pZTVYDYZz5fGPG7MypMFtaKFJEAJ0UGZLVa+OnyGZckZ7Mkswt/bg/vHxPDA2OjWZRPXGjKTjYkPx9eCuzcMvQvG/BTC4m1/A8LlSYASooMprqxh9c4sVm7LIKfERHQXX56ePoA7hkcS0Jpt1C218P1nxsSH3P3g2wWufRJGzAf/MJu3X3QcLhWgZJKEEM07kVfGsq0ZfLw3G1OtlbG9uvDsjEFM7heOe2uyPZhKjBRE29+A0mzo0gem/T+j1+TZyfY3IDocmSQhhAuzWjUbTxSwdEs6m0+cxcvDjVkJPXhwXAz9u3VuXaXFWbDjDWPvpZoyiB4PY58wnjO5tWI9lOjwZJKEEB1IRbWZj/dms2xrBmkFFYQHePPrm/oyd+RVbAp4ei9sexUOf2q8HzgLxj5uZBIXwg4kQAnhQrKLKlmxLZPVO09RZjIzNDKQl+5K4OZB3fDyaEXvxmqF418ZgSkzGbw7w+gfGwtrg3ra/gaEaEAClBDtnNaaXRlFLEtO5+vDZ1BKMXVQVx4aF8uwqKDWTROvqTSmiG//D5xLhcCecNNzRioin1YODQpxhSRACdFOVZstfJGSy7Kt6Rw6XUpgJ08WTOzF/WOi6R7UykkK5fmwcxHsWgxVhcbw3ewlMGCm5McTbU6+cUK0MwVl1azakcnb209xtrya3uH+PDdrELdfE0knr1ZkewDIP2oM4x14Hyw1EH+zsbA2eqwsrBUO41IBSqaZC1d26HQJy5IzSErJocZiZXJ8GPPGxTKhT2jrhvG0hvSNxsLa1P+Chw9ccw+M/imEyp8h4XgyzVwIJ7cj7RyvrEtlS+pZfL3cuWN4JA+MjaFXWCu3OzfXGFtcbHsVzhwEvzAjo3jiw+DXxbaNF6IFZJq5EO2I1potqWd55btUdmYUEurvzZM392PuyCgCO7Ui2wNAVRHsWQ473oSyXAjrB7e9AoPvlN1qhVOSACWEE9Ha2Eb9lXWp7M8qpmtnH56ZPoC7Rka1Lps4QFEGbH8d9q6E2gqIvdYITL1vkOdLwqlJgBLCCVitmq8Pn+GVdal8n1tKZHAn/jprMLOH98Dbo5WBKWsXbHsFjiSBcoNBdxiJW7sNsW3jhbATCVBCOJDFqvniQA6vrkvlRH45saF+/OOOIcy8pgeerdlG3WqBo18az5eydoBPIIz9Hxj1KHTubvsbEMKOJEAJ4QC1Fiuf7DvN6xtOkn62gr4R/rx0VwLThnRvXeLWmgrYt8pYWFuUDkHRMPXvcM294N3KyRRCOJgEKCHaULXZwod7snl9w0myi6oY2L0zb9w7jJsGdMWtNYGp7Iwx6WH3UjAVQ+QIuOEZ6D8d3Fo5NCiEk2gXAUopNRO4FQgHXtNaf+PgJglxRapqLLy76xRvbkzjTKmJhJ5BPDtjIJPjw1u3hinvsLF+6eAHYDVD/2nGVupRo2zfeCEcxO4BSim1FJgG5GutBzU4PhV4CXAHFmut/9ZcHVrrT4FPlVLBwIuABCjRLlRUm3l7eyaLNqdxtryGkbEhvPijoYzr3Yqt1LWGk98ZgSltPXj6QuI8I3lrSJx9bkAIB7psgFJKjQHuBSYA3YAq4BDwJfC21rrkMlUsB14FVjSo0x14DbgRyAZ2KaU+xwhWzzc6/yGtdX7d66fqzhPCqZWaankrOYMlyekUV9YyoU8oj0/uzai4ViyEramEQx8aGwPmHwb/rnD9n2D4PPANsX3jhXASlwxQSqm1QA7wGfAckA/4AH2BycBnSql/aa0/b64OrfUmpVRMo8MjgVStdVrddd4FZmitn8fobTVuhwL+BqzVWu9t2a0J0faKKmpYmpzO8q0ZlJnMXN8vnJ9e15thUcFXXllhGuxaAvveNp4vRQyCma/DoNng0co9nYRoRy7Xg7pPa3220bFyYG/dzz+VUqGtuG4PIKvB+2zgUoPnTwA3AIFKqd5a6zcaF1BKLQAWAERFRbWiSUK0XkFZNYs3p7FyeyaVNRZuHtSVn07uzaAegVdWkdUKJ9fBzoVw4htjokP/6UYqoqgxsrBWdCiXDFDng5NSyg+o0lpblVJ9gX4YvZnaJgJYSzT1p6zZpIBa65eBly/T1oVKqVxgupeX1/BWtEmIK3amxMSbm06yeucpasxWpg/tzk8n96ZvRMCVVVRVDPtXGdtcFKaBXzhc+1sY/qCsXxIdVksnSWwCJtRNUvgO2A3MAe5p5XWzgYbbcUZiDCVeFa11EpCUmJj4yNXWJcSlZBVW8sbGk3ywOxur1sy6pgc/ntSLuCtN4HrmEOxaZGxzUVsJPUfD5D9A/9vAw8s+jReinWhpgFJa60ql1MPAK1rrF5RS+67iuruAPkqpWOA0cBdw91XUZzRSttsQdpZ+toL/rE/lk32ncVOKOxIj+fG1vegZ4tvySiy1cPQLY2PAzGRjm4vBP4KRj0C3ofZrvBDtTIsDVN1svnuAh6/kXKXUamASEKqUygae1lovUUo9DnyNMXNvqdb68BW1vAnSgxL2cjyvjNfWp5KUkoOnuxv3jo7m0Wvj6BZ4BTvXluUZ2cT3LDOyiQdFw41/NrI9yGw8IS7S0gD1M+D3wCda64czk/wAACAASURBVMNKqThgfUtO1FrPbeb4GmBNC68vhEMcOl3Ca+tTWXvoDL5e7jwyIY75E+IIC2jhLDqtIWunMenh+8/AWmtkEZ/+kvFfyfYgRLNaFKC01pswnkOdf58G/I+9GtVaMsQnbGV/VjGvfHeC747mE+DtwRPX9WbeuFhC/Fr4XKi2Cg5+aASmMwfAO9AYwhsxH7r0sm/jhXARl9xRVym1EOOZ08EmPvPDmChRrbVeZb8mXjnZUVe01s70Ql5Zd4LNJ84S5OvJw+NiuX9sTMs3CSxMh911a5eqiiB8gBGYBt8pSVuFaEZrd9T9D/BHpdRgjOwRBRgLdfsAnYGlgNMEJ+lBidbQWpOceo6X151gZ3ohof5ePHlzP+4dHY2/dwsGGaxWSFtnTHo4/rWx99L5tUvRY2XtkhCtdMkeVH0hpfyBRH5IdXREa33Mzm1rNelBiZbQWrP+mLF77b5TxUR09ubRib2YOzKKTl4teDZUVQwpq43AVHgS/MKM9EOJ82TtkhBXoLU9KAC01uXABls3SghHsFo133xv7F57OKeUHkGd+MvMQfwoMbJlu9fmHTaC0oH3jLVLkSNh0u9hwG2SgkgIG2oX2220lAzxiUuxWDVfHszltXWpHMsrI6aLLy/cMYRZLdm91lJr7FS7cxFkbqlbu3QHjHgEuie0zQ0I0cG0aIivvZEhPtFQrcXKZ/tz+M/6VNLOVtAn3J/Hr+vNrYO74XG5wFSWB3vfgt3LoCwHgqKMmXjX3Cdrl4Swkasa4mtQiZ/WusJ2zRLCfqrNFj7ac5rXN6aSVVjFgG6def2eYUwZeJnda7WG7F3GFPHDnxprl3pdD9P+BX1ukrVLQrSRlmaDGAssBvyBKKXUUOBRrfVP7Nk4IVrDVGvh3Z2neHNTGrklJob2DOKZ6QO5rt9ldq+trYJDHxmBKTcFvDsbvaUR8yFUho2FaGst7UH9G5gCfA6gtU5RSk20W6taSZ5BdWwV1WZW7chk4aZ0zpZXMzImhBfuGML43qGXDkxFGXX7Lq001i6F9Ydb/wVD5sjaJSEcqMVDfFrrrEZ/yC22b87VkVx8HVNeqYnlWzN4Z8cpSqpqGd87lMevu4bRl9q91mo1tk3fuQiOf1W3dmmaMekhZrysXRLCCbQ0QGXVDfNppZQXRpqjI/ZrlhCXdzinhCWb00k6kIPFqpkysCuPTIy79O61phLYv9rY4uJcqrF2aeKvjfVLgT3arvFCiMtqaYB6DHgJYyfcbOAb4Kf2apQQzbFaNRuPF7BocxpbT57D18ude0ZF89C4WKK6XGLLi7zvjaCU8h7UVkDkCLh9EQyYIWuXhHBSLV2oe5bWb04oxFUz1Vr4eO9plmxJ42RBBV07+/Dkzf2YOzKq+Tx5FjMcq1u7lLEZ3L3r9l2aD92vadsbEEJcsZbO4osFngBiGp6jtb7NPs1qHZkk4XoKyqpZuT2Tt7dnUlhRw6AenXnprgRuGdyt+cW15fmw5y1j36XS0xAYBTf8n7F2ye8Sz6WEEE6lpbn4UoAlwEHAev641nqj/ZrWerJQt/07nlfGks3pfLL/NDVmKzf0D2f+hDhGxYY0PSOvfu3SIjj8ibF2KW6ykbC17xRZuySEE7vahbomrfXLNm6TEBfQWrMl9SyLN6ez8XgBPp5u3JkYyUPjYokLa2a6d3UZHPwAdi2FvIPgFQAjHq5bu9SnbW9ACGFTLQ1QLymlnsaYHFF9/qDWeq9dWiU6lGqzhc/357BkSzpHz5QR6u/Nr2/qy92jopvfIDDvsLF26cD7UFMGEYNh2r+NZ0zeAW17A0IIu2hpgBoM3Adcxw9DfLruvRCtUlRRw6odmby1LZOCsmr6dQ3gH3cM4baE7k1nFa81Gdum714CWTuMSQ+DbofEhyEyUdYuCeFiWhqgZgFxWusaezZGdAxpBeUsTU7nwz3ZmGqtXNs3jPl3xjaf8aEwzUjWuu9tqCqEkDi46S+QcI8kbBXChbU0QKUAQUC+HdsiXJjWmh3phSzenMZ3R/PxdHNj1jU9eHhCLH0jmhiSs5iNDA+7l8DJdaDcod8tRm8p9lpwu0wWciFEu9fSABUBHFVK7eLCZ1AyzVxcUq3FypqDuSzanMah06WE+HnxxHV9uG90NGEBTSyQLc2BvSuMaeJlORDQHSb9Lwy7T3apFaKDaek082ubOi7TzEVzSqpqWb3zFMuTMzhTaiIuzI/54+O4fVgPfDwbPV+yWiF9A+xeCkfXgLYY21skPgR9p4K7S+2rKYRo5Gq3fHfKQCScT1ZhJUu2pPP+7iwqayyM7dWFv94+iEl9wy/eg6myEPavMp4vFZ6ETiEw5qeQOM94ziSE6NAuGaCUUlu01uOVUmUYs/bqPwK01rqzXVsn2o09mUUs3pzG14fP4KYUtw3tzsMTYhnYPfDCgucX1O5aYiyotVRDz9Ew6Unofxt4+jjmBoQQTudyPSg/AK21LCwRFzFbrHx9OI/FW9LYd6qYzj4ePHptLx4YE0PXwEaBprocDr7fYEGtP1xzrzGM13WQY25ACOHULhegLv+ASnQ45dVm3tuVxbLkdLKLqoju4suzMwYye1gkft6NvlKyoFYI0UqXC1DhSqlfNveh1vpfNm6PcGKni6t4a2sGq3ecoqzazMiYEP44bQA39I/AveHzpfoFtUsha7ssqBVCtMrlApQ74I/xzEl0UAeyi1m8OZ0vD+YCcPOgrsyfEEdCz6ALC8qCWiGEDV0uQOVqrZ9tk5Y0QynVH/gZEAp8p7V+3ZHt6SisVs23R/JYvCWdnemF+Ht78NC4GB4YG0NkcIONAWVBrRDCTi4XoK6q56SUWgpMA/K11oMaHJ+KsUOvO7BYa/235urQWh8BHlNKuQGLrqY94vIqa8x8tCebJVvSyThXSY+gTjx1a3/mjOhJgE+DjQFLc40FtXvfMvZckgW1Qggbu1yAuv4q618OvAqsOH9AKeUOvAbciLF9/C6l1OcYwer5Ruc/pLXOV0rdBjxZV5ewg7xSEyu2ZbBqxymKK2sZ2jOIV6fEM3VgVzzObwxotUL6RqO31HBB7c0vyIJaIYTNXfJvFK114dVUrrXepJSKaXR4JJCqtU4DUEq9C8zQWj+P0dtqqp7Pgc+VUl8C7zRVRim1AFgAEBUVdTXN7lC+zyll8ZY0klJyMFs1UwZ0Zf6EWIZHB/+QuFUW1AohHMAR/+TtAWQ1eJ8NjGqusFJqEnA74A2saa6c1nohsBCMVEe2aKirslo1G08UsHhzGsmp5/D1cueeUdHMGxdDdBc/o5DWkLXTmIl36GNZUCuEaHOOCFBNPddqNqBorTcAG1pUsSSLvaTiyhq+PJjLsuQMUvPLiejsze+m9uPukVEE+tY9X5IFtUIIJ+GIAJUN9GzwPhLIsUXFWuskICkxMfERW9TnCspMtXx7JI+klFw2HS/AbNUM6NaZf88Zyq2Du+PlUfd8SRbUCiGcjCMC1C6gj1IqFjgN3AXcbYuKpQdlqKqx8N3RPL5IyWXdsXxqzFZ6BHXi4fGxTBvSnUE9OhvPl2pNkCILaoUQzsmuAUoptRqYBIQqpbKBp7XWS5RSjwNfY8zcW6q1PmyL63XkHlS12cKm42dJSsnh2yN5VNZYCA/w5u6RUUwf2p1hUUE/THqQBbVCiHbArgFKaz23meNruMSEh9bqaD2oWouV5NSzfHEgl68Pn6HMZCbY15OZ1/Rg+pDujIwN+SEFkdUCx76GXYvh5HeyoFYI4fRatGFhe+PKGxZarJqd6YUkHchh7cFciiprCfDxYMrArkwf2p2xvbrg6d4g2JQXGItp9yyHkixjQe3wB2DY/bKgVgjhFK5qw0LhWFarZl9WEUkpuXx5MJeCsmp8vdy5oX8E04d2Z2LfULw9GuxSqzVk7TB6S4c/BWut0Uua8leIv0UW1Aoh2gWX+pvKlYb4tNYcziklKSWHLw7kcrq4Ci8PN66LD2f60O5c1y+cTl6Ntk6vLoeDHxiz8fIOgndnGPGwMYwX1tcxNyKEEK0kQ3xO5nheGUkpOSSl5JBxrhIPN8XEvmFMH9qNG/pHXJgP77yC40b6of3vQHWpMUV85HxjiriXX9vfhBBCXAEZ4nNi6Wcr+CIlh6QDORzPK8dNwdheofx4Ui+mDOxKkK/XxSdZzHDsS2MYL30TuHvBgJkwYj70HClTxIUQ7Z5LBaj2NMSXXVTJlwdySTqQw6HTpQCMjAnhzzMGMnVQN8ICvJs+sewM7Kmb9FCWA4E94fo/wTX3g39Y292AEELYmQzxtaH8UhNfHswlKSWHvaeKARjaM4jpQ7px65BudAvs1PSJWkNmstFbOpIEVrORRXzkI9DnJnBzb/o8IYRoB2SIz0EKK2pYe8gISjvSC9Ea+nfrzG+nxjNtcHeiuvg2f7KpFA68ZwSmgqPgEwSjHjPy4nXp1XY3IYQQDiAByg5Kqmr55vAZkg7kkpx6FotVExfmx8+u78O0Id3pHe5/6QryvjeC0oH3oKYcuiXAjNdg4O3gdYmAJoQQLsSlApQjn0FVVJsvSMpaY7HSM6QTCybGMX1Id/p3C/gh1VBTzDVwNAl2LoZTW+vy4s02Jj1EDm+7GxFCCCchz6CugqnWwoZj+SSl5PLd0TxMtVa6dvbh1iHdmD60O0MjAy8dlABKThsTHva+BeV5EBxjrFu65l7JiyeE6BDkGZSN1JitbEktICkll/9+n0d5tZlQfy9+NLwn04d2JzE6GDe3ywQlrSFtgzGMd2wtaCv0nWL0lnpdL3nxhBACCVAtYrZY2Z5WSFJKDl8dPkNJVS2BnTy5dbDRUxodF4KHewuCSlUxpKw2Mj2cO2FsnT72CWPr9OAYu9+HEEK0Jy4VoGz5DMpq1ezOLCIpJYe1h3I5W16Dv7cHNw6IYPrQbozvHfbDZn+Xk3vA6C0d/ABqKyFyBMx601hYK1unCyFEk1wqQNliP6hT5ypZsS2DLw7kcqbUhI+nG9f3M4LSpPhwfDxbuObIXA3ffwY7F0H2TvDoBIPvMIbxuie0tnlCCNFhuFSAsoWC8mpWbMtkYt8wfn9LP27oH4Gf9xX8mooyYc8y2LsCKs9BSC+Y8jwkzIVOwfZruBBCuBgJUI0Miwpi11M3ENipiaSszbFa4eQ62LUIjn9t5MGLv8XIJB47SSY9CCFEK0iAakQp1fLgVFlobJu+eykUpYNfGEz4lTHpITDSvg0VQggXJwGqNU7vMWbiHfoIzCaIGgPXPQX9bwOPJjKPCyGEuGISoFqqtgoOfWzMxsvZC55+kHC3sai26yBHt04IIVyOSwUou6Q6KkwzhvD2vQ1VRRAaD7e8CEPmgE9n211HCCHEBVwqQNlimjkAVguc+MboLaV+C8od+k+DEY9AzHjZDFAIIdqASwUom8jcCh8/CiWnIKAbTPo9DHsAOndzdMuEEKJDkQDVWFA0hMTCTX+GfreC+xVMNxdCCGEzEqAaC+wBD3zu6FYIIUSHJytIhRBCOCUJUEIIIZySBCghhBBOSQKUEEIIp9QuApRSyk8ptUcpNc3RbRFCCNE27BqglFJLlVL5SqlDjY5PVUodU0qlKqWebEFVvwPet08rhRBCOCN7TzNfDrwKrDh/QCnlDrwG3AhkA7uUUp8D7sDzjc5/CBgCfA/I1rNCCNGB2DVAaa03KaViGh0eCaRqrdMAlFLvAjO01s8DFw3hKaUmA37AAKBKKbVGa21totwCYAFAVFSULW9DCCGEAzhioW4PIKvB+2xgVHOFtdZ/AFBKPQicbSo41ZVbCCwESExM1LZqrBBCCMdwRIBqKtPqZQOK1nr5ZSu2RzZzIYQQDuGIWXzZQM8G7yOBHFtUrLVO0lovCAwMtEV1QgghHMgRAWoX0EcpFauU8gLuAmyS/E4pNV0ptbCkpMQW1QkhhHAge08zXw1sA+KVUtlKqYe11mbgceBr4Ajwvtb6sC2uJz0oIYRwHfaexTe3meNrgDW2vp48gxJCCNfRLjJJtJT0oIQQwnW4VICSZ1BCCOE6XCpASQ9KCCFch0sFKCGEEK5DApQQQgin5FIBSp5BCSGE63CpACXPoIQQom1UW6o5VniMgwUH7XYNR+TiE0II0U5UmavIKMngZMlJThYbP2klaWSVZWHVVgaHDuadW9+xy7VdKkDJQl0hhGiditoK0kvSjSBUcpK04jROFp/kdPlpdF0+bw/lQVTnKPoG92VqzFR6BfWiT1Afu7VJae16O1MkJibq3bt3O7oZQgjhdEprSkkrTiOtJO2CYJRbkVtfxtPNk5jAGHoF9iIuKI5egb3oFdSLqIAoPN09bd4mpdQerXVi4+Mu1YMSQghhKDYV1w/LnQ9GacVp5Ffl15fxdvcmLjCOYRHDLghGkQGReLg5Pjw4vgVCCCFaRWvNOdM5YziuUTAqNBXWl+vk0Ym4wDhGdx9NXGAcvYKMHlF3v+64u7k78A4uTQKUEEI4Oa01BVUFFwSg86+Lq4vry/l7+hMXFMe1kdfSK6hXfTDq6tcVN9X+Jm27VICSSRJCiPZMa82ZijNNDs2V1ZbVl+vs1ZneQb25IfqGC4bmwn3DUaqpTcvbJ5kkIYQQbcyqrZwuP33h0FzdxIVKc2V9uRCfkAt6QueDURefLq4ViGSShBBCtC2rtpJVlkVqcWp9MEorTiO9JB2TxVRfLqxTGHFBcczqM6s+GMUFxhHsE+zA1jueBCghhLCBitoKThSd4FjhMY4VHeNY4TFOFJ+gylxVX6arX1d6BfYisWti/dTt2MBYAr0l+01TJEAJIcQV0FqTU5FTH4iOFx7nWNExssqy6ssEeAUQHxzPrN6ziA+Jp09QH2IDY/H38ndgy9sfCVBCCNEMk9lEanHqhb2iohP1ExYUiqjOUfQL6ceMXjOID4knPjiern5dXeoZkaNIgBJCdHhaa/Ir840eUdHx+oCUWZqJVVsB8PXwpW9wX26Ju4W+wX3re0a+nr4Obr3rcqkAJdPMhRCXU2up5WTJyYuG6BquJ+rh34O+wX2ZEjOF+GCjV9QjoEe7XEvUnsk0cyGEyzpXde6CIHSs6BjpxemYtRkwUv30CepDfEh8fa+ob3BfArwCHNzyjkWmmQshXJbZaiajJKM+CJ0PSGerztaXCfcNJz44nmsjryU+OJ6+IX2JDoh26lQ/HZ0EKCFEu1JSXXLBc6Jjhcc4WXySGmsNAB5uHvQO6s3Y7mON4bm6XlFHX1PUHkmAEkI4Jau2cqr0VH0QOl5k9IrOVJypLxPiE0J8cDx397+7foguNjAWTzfbbwkh2p4EKCGEw1XUVlzQKzpeePyCRa7uyp3YwFiGhQ+rn8odHxJPaKdQB7dc2JMEKCGE3Vi1lUJTIfmV+fU/eZV5F70vq7kwEWp8SDyz+8yu7xX1CuqFt7u3A+9EOIIEKCFEq1SZqyioLGgy4Jx/XVBVgNlqvuA8N+VGqE8o4b7hRAVEkRiRSIRfRP1sugjfCFnkKoB2EKCUUpOAPwOHgXe11hsc2iAhXFxrej3n+Xr4Eu4bToRvBIkRiYT7hte/P/+6S6cuTrFbq3B+dv2WKKWWAtOAfK31oAbHpwIvAe7AYq313y5RjQbKAR8g247NFcLlmcymZgPO5Xo9XXy6EO4bTs+AngyPGH5B0Dn/WnLNCVuy9z9jlgOvAivOH1BKuQOvATdiBJxdSqnPMYLV843OfwjYrLXeqJSKAP4F3GPnNgvR7jTX6ymoLLggAJXWlF50rvR6hLOy6zdOa71JKRXT6PBIIFVrnQaglHoXmKG1fh6jt9WcIkCekgqXp7WmxlqDyWyiylyFyWzCZDFRVlN2Rb0ehSK0U6j0ekS75Yh/EvUAshq8zwZGNVdYKXU7MAUIwuiNNVduAbAAICoqyiYNFaKhhoHjfNCoDyJ1rxseN5lNVFmqmjx+uXM1l05B1smjU32AGR4xvMleT2inUOn1iHbNEd/epqbnNPunUWv9MfDx5SrVWi9USuUC0728vIZfRftcjtYak8VERW0FlbWVWLSl/jPV4H9Hw5lTqtH/pgveq6aPN3d+c2Waq7/Zelpwrtlq/uEv/GYCRbWl+oKeSeOeSlPB5PznlwscTfF088THw4dO7p3w9vCuf+3j4UOYbxg+7j74ePj88F8PHzp5dLrgvY+7D36efvXBx9/TX2a6CZfniACVDfRs8D4SyLFFxVrrJCApMTHxEVvU50hWbaXKXEV5TTkVZiOwlNeW1weZitoKymvLm3xd/2OuoKLG+O/5LQPED84HjgsCQ13gCPUKvWzg8Hb3Nt7Xlal/3eAzb3dv6cUI0UqO+JOzC+ijlIoFTgN3AXfbomJHb7dhsVp+CAqNAsQFgaPRT+PgU15bTqW5skXXdFfu+Hr64ufph7+nf/3rCL8IfD2M141/3JWRHLNhb+CC15fIcN9cuZac36JzW1JPM+eCMdusk0eniwKHt7v3BUFGAocQzs/e08xXA5OAUKVUNvC01nqJUupx4GuMmXtLtdaHbXE9W/Sgik3FHDx7sL7X0qKeSt2PyWJq0TU83Dzw9/S/IGgE+QQRGRCJn6cfvh6++Hv54+fhh6+nb33Zxq/9PP3wcfeRoR4hhEuy9yy+uc0cXwOssee1W+tI4RF+8t1PLjru7e59UU8kzDeMGI8Y/Lz88PO4uKfS3I+Xu5cD7kwIIdoXlxrjsMUQ38DQgay8eeUFAcXX01eyIwshRBuTHXWFEEI4VHM76ro5ojH2opSarpRaWFJS4uimCCGEuEouFaC01kla6wWBgYGObooQQoir5FIBSgghhOtwqQAlQ3xCCOE6XCpAyRCfEEK4DpcKUEIIIVyHBCghhBBOyaUClDyDEkII1+GSC3WVUgVAZguLBwJXGtFaeo6tyoUCZ1tQT3vUmt9/e2qDLeq+mjqu9NwrKd+Ssi0pI9/v9nl9W9UdCARprcMu+kRr3aF/gIX2OsdW5YDdjv49OdPvvz21wRZ1X00dV3rulZRvSdkWlpHvdzu8vq3qvlQ9LjXE10pJdjzH1uVckTPcuz3bYIu6r6aOKz33Ssq3pKwz/P91JEffv7N/ty9Zj0sO8bkapdRu3USeKiFcgXy/RXOkB9U+LHR0A4SwI/l+iyZJD0oIIYRTkh6UEEIIpyQBSgghhFOSACWEEMIpSYASQgjhlCRAtUNKqZlKqUVKqc+UUjc5uj1C2IpSqr9S6g2l1IdKqR87uj3CsSRAOQml1FKlVL5S6lCj41OVUseUUqlKqScBtNafaq0fAR4E5jiguUK02BV+t49orR8D7gRkbVQHJwHKeSwHpjY8oJRyB14DbgYGAHOVUgMaFHmq7nMhnNlyruC7rZS6DdgCfNe2zRTORgKUk9BabwIKGx0eCaRqrdO01jXAu8AMZfg7sFZrvbet2yrElbiS73Zd+c+11mOBe9q2pcLZeDi6AeKSegBZDd5nA6OAJ4AbgEClVG+t9RuOaJwQV6HJ77ZSahJwO+ANrHFAu4QTkQDl3FQTx7TW+mXg5bZujBA21Nx3ewOwoW2bIpyVDPE5t2ygZ4P3kUCOg9oihC3Jd1tclgQo57YL6KOUilVKeQF3AZ87uE1C2IJ8t8VlSYByEkqp1cA2IF4pla2UelhrbQYeB74GjgDva60PO7KdQlwp+W6L1pJs5kIIIZyS9KCEEEI4JQlQQgghnJIEKCGEEE5JApQQQginJAFKCCGEU5IAJYQQwilJgBKijlLq30qpnzd4/7VSanGD9/9USv3SjtfPUEqF2rjOGKXU3Q3eP6iUerWF536olIpr0LaPGnx2h1Jqed3raUqp/7Nlu4UACVBCNLQVGAuglHIDQoGBDT4fCyQ7oF1XIwa4+3KFGlNKDQTctdZpDQ4n1h1v7EvgNqWUb+uaKETTJEAJ8YNk6gIURmA6BJQppYKVUt5Af2CfUspfKfWdUmqvUuqgUmoGgFLq70qpn5yvTCn1jFLqV3Wvf6OU2qWUOtCS3oZS6l6l1E6l1H6l1Jt1+yehlCpXSj2nlEpRSm1XSkXUHe9V936XUupZpVR5XVV/AybU1fOLumPdlVJfKaVOKKVeaKYJ9wCfNTr2IvC/jQtqY7X/BmDa5e5LiCshAUqIOlrrHMCslIrCCFTbgB3AGIzdXQ/U7V1kAmZprYcBk4F/KqUUxp5GDXc4vhP4QCl1E9AHYw+kBGC4Umpic+1QSvWvq2ec1joBsPDD3kh+wHat9VBgE/BI3fGXgJe01iO4MOnqk8BmrXWC1vrfdccS6uofDMxRSjVM2nreOGBPo2PvA8OUUr2bKL8bmNDcPQnRGhKghLjQ+V7U+QC1rcH7rXVlFPBXpdQB4FuMvY0itNb7gHClVHel1FCgSGt9Crip7mcfsBfohxGwmnM9MBzYpZTaX/c+ru6zGuCLutd7MIbwwAiiH9S9fucy9/id1rpEa20CvgeimyjTDShodMwC/AP4fRPl84Hul7muEFdE9oMS4kLnn0MNxhjiywJ+BZQCS+vK3AOEAcO11rVKqQzAp+6zD4E7gK4YPSowAtrzWus3W9gGBbyltW4qENTqHxJoWmjdn+HqBq+bq6OKH+6poZUYAapxYlefunOEsBnpQQlxoWSMZymFWmuL1roQCMLooWyrKxMI5NcFp8lc2AN5F2PriDswghUYGbsfUkr5Ayileiilwi/Rhu+AO86XUUqFKKWa6uU0tB2YXff6rgbHy4CAy5zblCPARUN5Wuta4N/Azxt91BcjoAthMxKghLjQQYzZe9sbHSvRWp+te78KY0bbboze1NHzBeu2jAgATmutc+uOfYMx7LZNKXUQI3A1GzS01t8DTwHf1A0j/hdjyO1Sfg78Uim1s65sSd3xAxjP1VIaTJJoiS+BSc18toSLe12T684RwmZkuw0hXEDdFO8qrbVWSt0FzNVaz7iK+joB6zEmalguUzYCeEdrfX1rrydEUyRACeEClFIT/H3cGQAAAFFJREFUgFcxnl8VAw9prVOvss4pwJG6iR6XKjcC49nY/qu5nhCNSYASQgjhlOQZlBBCCKckAUoIIYRTkgAlhBDCKUmAEkII4ZQkQAkhhHBK/x/XI8KOosEFhQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(ns, bests, label='analyze1')\n", "plt.plot(ns, bests2, label='analyze2')\n", "plt.plot(ns, bests3, label='fftpack.dct')\n", "decorate(xlabel='Wave length (N)', ylabel='Time (s)', **loglog)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2\n", "\n", "One of the major applications of the DCT is compression for both sound and images. In its simplest form, DCT-based compression works like this:\n", "\n", "1. Break a long signal into segments.\n", "2. Compute the DCT of each segment.\n", "3. Identify frequency components with amplitudes so low they are inaudible, and remove them. Store only the frequencies and amplitudes that remain.\n", "4. To play back the signal, load the frequencies and amplitudes for each segment and apply the inverse DCT.\n", "\n", "Implement a version of this algorithm and apply it to a recording of music or speech. How many components can you eliminate before the difference is perceptible?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`thinkdsp` provides a class, `Dct` that is similar to a `Spectrum`, but which uses DCT instead of FFT." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example, I'll use a recording of a saxophone:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "if not os.path.exists('100475__iluppai__saxophone-weep.wav'):\n", " !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/100475__iluppai__saxophone-weep.wav" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from thinkdsp import read_wave\n", "\n", "wave = read_wave('100475__iluppai__saxophone-weep.wav')\n", "wave.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's a short segment:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "segment = wave.segment(start=1.2, duration=0.5)\n", "segment.normalize()\n", "segment.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's the DCT of that segment:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZgc1X3u8e87M9pYtCFhQBKIRYABs3nYwcYGgwAH4efiGDsxso0vCcELdmKzJQY7kHjJjR0STC4GAuRiAyb2RYmxQcYQY2xAYkcIISEWCQk0QgtCQst0//JHnRm1htYMmu6eqtG8n+eZZ7pPV1f/qmZ5u06dPqWIwMzMrGia8i7AzMysGgeUmZkVkgPKzMwKyQFlZmaF5IAyM7NCasm7gLyMGTMmJk6cmHcZZmYD3qOPPro0IsZ2bR+wATVx4kRmzpyZdxlmZgOepJertbuLz8zMCskBZWZmheSAMjOzQnJAmZlZIeUaUJJGSrpD0nOSZks6StJoSdMlzU3fR6VlJekqSfMkPSXp0Ir1TE3Lz5U0Nb8tMjOzesn7COqfgF9FxL7AQcBs4CLg3oiYBNyb7gOcAkxKX+cC1wBIGg1cBhwBHA5c1hFqZmbWf+UWUJKGAx8ArgeIiPURsQKYAtyUFrsJOCPdngLcHJmHgJGSdgZOBqZHxLKIWA5MByb34aaYmVkD5HkEtQfQBvybpMclXSdpW+A9EbEYIH3fMS0/DlhQ8fyFqW1z7e8g6VxJMyXNbGtrq+/WmJlZXeUZUC3AocA1EXEIsJqN3XnVqEpbdNP+zsaIayOiNSJax459x4eWB5znX1/FjQ++mHcZZmZV5RlQC4GFEfFwun8HWWC9nrruSN+XVCw/oeL544FF3bRbD079pwe4/D+fzbsMM7OqcguoiHgNWCBpn9R0AvAsMA3oGIk3Fbgz3Z4GnJ1G8x0JrExdgHcDJ0kalQZHnJTarAftZV9N2cyKK++5+L4I3CJpMDAf+CxZaN4u6RzgFeDjadm7gFOBecCatCwRsUzS3wIz0nLfiohlfbcJZmbWCLkGVEQ8AbRWeeiEKssGcP5m1nMDcEN9qzMzszzl/TkoMzOzqhxQZmZWSA4oMzMrJAeUmZkVkgPKzMwKyQFlZmaF5IAyM7NCckCZmVkhOaDMzKyQHFBmZlZIDigzMyskB5SZmRWSA8rMzArJAWVmZoXkgMpBqRx89bYnmPPaqrxLASC7komZWbE4oHIwd8kqfvb4q3zpJ4/nXYqZWWE5oMzMrJAcUIZ7+MysiBxQZmZWSA4oMzMrJAeUmZkVkgPK8CkoMysiB1QOPCjBzKxnDqgcSXlXYGZWXA4oMzMrJAdUjorS1eepjsysiBxQOXDXnplZz3IPKEnNkh6X9F/p/u6SHpY0V9Jtkgan9iHp/rz0+MSKdVyc2udIOjmfLTEzs3rKPaCALwOzK+5/B/h+REwClgPnpPZzgOURsRfw/bQckvYDzgL2ByYDP5TU3Ee1bxXcwWdmRZRrQEkaD5wGXJfuC/gwcEda5CbgjHR7SrpPevyEtPwU4NaIWBcRLwLzgMP7Zgt6x6d8zMx6lvcR1A+ArwPldH8HYEVEtKf7C4Fx6fY4YAFAenxlWr6zvcpzNiHpXEkzJc1sa2ur53b0is9FmZltXm4BJemjwJKIeLSyucqi0cNj3T1n08aIayOiNSJax44du0X1mplZ32rJ8bWPAU6XdCowFBhOdkQ1UlJLOkoaDyxKyy8EJgALJbUAI4BlFe0dKp9j74K7HM2siHI7goqIiyNifERMJBvk8JuI+BPgPuDMtNhU4M50e1q6T3r8N5F9gGcacFYa5bc7MAl4pI82w8zMGiTPI6jNuRC4VdIVwOPA9an9euDfJc0jO3I6CyAiZkm6HXgWaAfOj4hS35dtZmb1VIiAioj7gfvT7flUGYUXEWuBj2/m+VcCVzauQjMz62t5j+KzAgh/EsrMCsgBZWZmheSAMjOzQnJAmZlZITmgzMyskBxQZmZWSA4oMzMrJAdUDjy1kJlZzxxQZmZWSA4oMzMrJAeUmZkVkgMqB75QoZlZzxxQZmZWSA4oMzMrJAdUDjzM3MysZw4oMzMrJAeUmZkVkgPK3OVoZoXkgDIzs0JyQJmZWSE5oMzMrJAcUGZmVkgOKDMzKyQHlJmZFZIDyszMCskBZWZmheSAMjOzQsotoCRNkHSfpNmSZkn6cmofLWm6pLnp+6jULklXSZon6SlJh1asa2pafq6kqXltk5mZ1U+eR1DtwF9GxHuBI4HzJe0HXATcGxGTgHvTfYBTgEnp61zgGsgCDbgMOAI4HLisI9SKKvDcQmZmPcktoCJicUQ8lm6vAmYD44ApwE1psZuAM9LtKcDNkXkIGClpZ+BkYHpELIuI5cB0YHIfbkqvyZfWNTPbrEKcg5I0ETgEeBh4T0QshizEgB3TYuOABRVPW5jaNtduZmb9WO4BJWk74D+ACyLize4WrdIW3bRXe61zJc2UNLOtrW3Li62z8DTiZmablWtASRpEFk63RMTPUvPrqeuO9H1Jal8ITKh4+nhgUTft7xAR10ZEa0S0jh07tn4bsoVUNVPNzKxSnqP4BFwPzI6If6x4aBrQMRJvKnBnRfvZaTTfkcDK1AV4N3CSpFFpcMRJqc3MzPqxlhxf+xjg08DTkp5IbZcA3wZul3QO8Arw8fTYXcCpwDxgDfBZgIhYJulvgRlpuW9FxLK+2QQzM2uU3AIqIn5H9fNHACdUWT6A8zezrhuAG+pXXWN5mLmZWc9yHyQxkHmYuZnZ5jmgzMyskBxQZmZWSA4oMzMrJAeUmZkVkgPKzMwKyQFlZmaF5ICy3H36+oe57oH5eZdhZgXjgLLcPTB3KVf8YnbeZZhZwTigzMyskBxQZmZWSA6ofqxU9px+HX7z3Ou8uXZD3mWYWR05oPqppxauYM9L7uK3z+d74cXL7nyG0//ld7nWsGjF23zuxpl8+SeP51qHmdWXAyoH9biQ7oyXlgNw35wlPSzZWDf94WWeWrgy1xre3lAC4KU31uRah5nVlwMqR7XMZd6UnuyrxpvZ1soB1U91hFvZCdXJFy8x27o4oHJUS7h0XEuqHvnkjDOzInJA5aDjnM1zr63q9To6u/h8dV4z20o5oHKwYHkdTuanIyiPNDezrZUDqp/yIAkz29o5oHIwYtigmtchOs5BOaE6eE+YbV0cUDnYb+fhABw2cVSv1+EjKDPb2jmgclBKqTJscEuv15FOQXmYuZlttboNKEl/11eFDCTlNLKhuYYP7nR28fXy+e4aNLOi6+kIanKfVDHAdEzy2txUU0IBve/ie+bVN3v/2nVUz6DM+4O6P374FQ678tc5V2G29eipj6lZ0ig287cfEcvqX9LW74kFKwBYsab3s2/f9PuXAJi9uHdB89D8N3r92vW0Nc3IfsnPn867BLOtSk8BtS/wKNUDKoA96l7RAPDD+18AYObLy3u9jlmLsmB6+Y3VvXr+lXcV4wq2G0pbT0CZWX31FFDPRsQhfVLJAFRTF1+yen1pi5+zrn3Ln9MoG8rluq2rKFEXEZ1TUeVhyaq1bDu4hW2H9H4Qjm19Fixbw1MLV3LagTvnXcq7ttWM4pM0WdIcSfMkXZR3PZtT2aV10PgRffa65XKwdkOJf/3vF9jnr39V9/VvKPUuaGYvqt+5sKWr1tVtXbWo5bTa4pVvs2BZbTONHH7lvblfo2tDqcwdjy7sHBCUl/ZSueZu5N+/sJT5bW/VtI5X3ljD+vbevxlb317mnlmv1VTDaVc9wPk/fqymdSxYtoZnXu27y+v09BbrR5LGRsQmV8WTtCPwZkSsbVxp756kZuBq4CPAQmCGpGkR8WwjXu+cG2dw73O1X4fp2L3GbHK/vVRmyap1LF+zntmLV/H6m2v53dyl7Dh8CPPbVvP0Zn4xJl70i5rqKEcQEZQD1qxvRxKr1m5gUHMTg1uauPuZ1zjz/eM7u+MGNesdRwiTLv0lAEftsQMjtxlE68TRDG4Wj7+ygueXrGLSjtvTOnEU989pY/HKt7sdpNFRS5Oyf/YSLFz+NiO2GcTwoYMolYNyBC1NIgKaUi2r1rUzb8lbjNluMBtKwfBhLaxdX+aN1euY37aaBcvXcNCEkey4/RAk0V4qM3LYYF5ftZaH579B21vr+fMP7sGg5ibaS0EQREB7KVjXXuKZRStZumo9H9xnLMMGN9NeCt7eUKKlSTyZzit23Y4IaE//INvLZd5a285b69rZacRQhrQ0s6FUZn2pjMi6O4/6+98AMP/vTqU9befg5ibWtpdoaWpCgpeWrmb7oYPYbmgLQ1qaeGttO4NamhjS0sSaddnR8Qttq2kvlWluEqVyUIqgpamJcmTrbC9l37cZ3EKpHKxe187wYYPYUCqzbkOZz974CJ8+ajcm778zEgxpaWLl2xuQxDaDmylHsGptO9sPbaFZIsjOAzQ3iQ2l4OY/vMQVv5jNvCVvcd4H96S5WbQ0ieam7Lsk1rWXWLu+zJBBTUTQ+Vgp/Wx/9tir/OVPn+Tpy09i2KBmmpvUuU8GpX2xLv3Tb06/C+UIBjU3saFUphzBR/7xt4wYNohffOnYzp9LqRy0NDdRKgdNyqYLW9deYkMpKJWD9nKZ9nR7Q6nMp370cOfPBLLfR0mdP9+u92Hj0byANRtKfOB793HMXjtw42cP76x/QzkLz46vbQa3pG0s09LU1Pk7s6EUHHDZ3QB87eR9OP9De2V/t+VsJs7Kz0R29Eg0STSnv4sN5TLNEm+ubc/ul8qdj60vlRnc3NS5Xe3l2OTn2bFtHY777n0APPmNkxg2uJkmZa/Vdbl6UXejqCRdC/wqIn7Wpf1PgGMj4ry6V9QLko4CLo+Ik9P9iwEi4u8395zW1taYOXNmr16v1kAwG+g6gmFzJH8IvTs97b96kTaG3ZCWJlata6+63K8uOI59dxpew+vo0Yho7dre0xHUsRFxbtfGiLhF0iW9rqb+xgELKu4vBI7oupCkc4FzAXbddddev9ghu47k8Vfe+Y753frm6ftz2bRZnfeHtDR1vhO02n1g77H89vm2nheswaBmdR5RjtluCMOHtTC/beOAlaGDmli7ofuf6Q7bDuaN1esbWudxk8bwwNylAOy143Z8aJ+x/OiBF7doHYdNHNV5BedqjtpjBw4cP4Lla9bzwNylLF65sWPl4Akj+aODdqFcDtrLQalcTt+Df/7NPAAOHD+Cpxau5PPH7k4pguFDB9FeLnP1fdlgoo+/fzwTRm9DKXVTDx3UTAAr1qzn5j+8zOkH7cKkHbfrfBffXgpamsXLb6zm9pkLO2v5ROsE3jNiKESwvhS80PYW40YOY85rq/jgPmNpSUdwLc1Nm3y/4LYnAPjSCZM6j1bKAUTQ1CTKsXEUmVT5GcWsJ+D+OUs6r2Dw1Y/sTZNgfSkY0tJEc1P2z/+XzyzmsVdW8JmjJzJ628GdRzItzWLaE4uY8/rGKx/8r0PHs9OIIQixbM36zrDabkgL2w1pob1UZu6St9hpxFBWr2tnx+2HsnzNem55+JXOdXzxw3t1HnXOX7qaPcZsSzmCJxeu5ODxIzq3KztKy35269vL/NuDLwEwftQw/rh1Qradkf0NNEJPR1CzI+K9W/pYX5P0ceDkiPh8uv9p4PCI+OLmnlPLEVStVq3dwPsuvweAfz/ncI6bNHaLnv/m2g0cmJ4P8NK3T+tVHd0dCV75sQPYZeQwFi5/m+P3HsuE0du8Y5lyOdjjkrtqqqGyjhmXnsjY7bf8F/2JBSs44+oHufpTh/b6BPA37nyGm//wMt88fX+mHj2xV+vo2I67vnQc++3Su3eTHeuoZX8uX72ewS1N/X6QRMe+mHflKbQ09+50eT325/1zljBmuyEcMK5354wXLFvDcd+9jz/7wB5cfGrv/mV+6z+f5YYHX+SvT3svnz+ud4On7571Gpf+/Gl+d+GHGTqouVfr+Mkjr1AqB3965G69ev7m9PYIaomkwyPikS4rOwxo7FvULbMQmFBxfzywKKdatsiWhhPA8KEbJ5v9o4N2qVstW/qPoKkOoxAr9SacIHuX/tjffIRR2/R+Et6O92m1dKPvNHwor725lpE11HHd2a28uLR3Hx3oMGrbwTU9vyjOO35Prrn/hbqMdq3F8fvsWNPzJ4zehl9/9QPstsO2vV5HPebePHn/nTh5/516vwLgk4f3vuepN3oKqK8Bt0u6kezzUACtwNnAWQ2sa0vNACZJ2h14lay2T+VbUt8YPrQ+75Kvn9ra63epRTC6AP+UO+ZFrOUf6on7vade5fR7F07elwsn71vTOv7uY+/j5WW1BX497LXj9jU9f+KYLNx2Hjm0HuX0G93+d4uIRyQdAfwF8JnUPAs4IiJqH8ZWJxHRLukLwN1AM3BDRMzq4WlbhXqdKP1Qje8SbePPIsePQFkXnzqib9/xN8qfHLEru4/ZlqP33CHvUvpUj2+/I+J14DJJY9P9InXtdYqIu4C78q6jr5Xq9EHXenfX9Tf1CJWOI6gmJ5TVmSSO6fKxlIGgp9nMJelySW3Ac8AcSW2SvtE35VlP2j1VUF3UY0hzZxefA8qsLno66XABcAzZiLgdImI02fDtYyR9peHVmfUjZx2WdScNG9y7EVJmtqmeAups4JMR0fnBiYiYD/xpesx6oR6fuD71fdlonAP7cLok696Fk/fh+StO6fUQXjPbVE8BNSgilnZtTOehej+W1mrW8cG4eoTd4Jb+O3qvXk4/OBuuf/Seve/nl+R9aVZHPQ2S6O5j7o39CLy9K/W44N8AHx8BwGETR9f0YU4zq7+eAuogSdVm9RQwsAbkN0DeH0DsoBquRTvtC8d0TobaW/9x3lFsN8QH5Ga2qZ4+B+XO9AaqRz7Vo4tvUHPv13Hg+JE1v/77dxtd8zrMbOvjDvMc5XlRu0pDfFLfzArIAZWjehxB1XIOakg6oX/UHgPr0+lm1j84oHKU94wDB0/IuufOOmxCD0uamfU9B1QOoo5T4tSlm7AYPY1mZptwQOWgo1PugHG9vwJlPeuoZRSfmVmj9O8rmvVTw4cO4sefP4ID8p4FwrNvm1mBOaBycnQBZiaOlFDOJzMrInfxWWGGu5uZVXJADWD1uMSEmVmjOKD6qXqES+cgCR9AmVkBOaD6uVrCpWO4u/PJzIrIAdVPBbUfQvkIysyKzAHVz9WSLRu7CZ1QZlY8Dqh+qp4DHHwEZWZF5IDq72pIFw/iM7Mic0D1U3UJFw+SMLMCc0D1czWdg+pYh/v4zKyAHFD9VF0+B9UxF1/tqzIzqzsHVD9X0+egOubic0KZWQE5oPqtOnwOqvMIygllZsWTS0BJ+p6k5yQ9JennkkZWPHaxpHmS5kg6uaJ9cmqbJ+miivbdJT0saa6k2yQN7uvtyVMt4RK+3IaZFVheR1DTgQMi4kDgeeBiAEn7AWcB+wOTgR9KapbUDFwNnALsB3wyLQvwHeD7ETEJWA6c06dbkhNP9GpmW7tcAioi7omI9nT3IWB8uj0FuDUi1kXEi8A84PD0NS8i5kfEeuBWYIqy4WcfBu5Iz78JOKOvtiNPpx+8CwBH7DE650rMzBqjCBcs/BxwW7o9jiywOixMbQALurQfAewArKgIu8rl30HSucC5ALvuumvNhefp6D3H8NK3T8u7DDOzhmlYQEn6NbBTlYcujYg70zKXAu3ALR1Pq7J8UP1IL7pZvqqIuBa4FqC1tdWdZGZmBdawgIqIE7t7XNJU4KPACRGdZ1QWAhMqFhsPLEq3q7UvBUZKaklHUZXLm5lZP5bXKL7JwIXA6RGxpuKhacBZkoZI2h2YBDwCzAAmpRF7g8kGUkxLwXYfcGZ6/lTgzr7aDjMza5y8zkH9CzAEmJ6m2XkoIv48ImZJuh14lqzr7/yIKAFI+gJwN9AM3BARs9K6LgRulXQF8Dhwfd9uSv/lPk4zK7JcAioi9urmsSuBK6u03wXcVaV9PtkoP+slfw7KzIrIM0mYmVkhOaDMzKyQHFADWHg6CjMrMAeUebJYMyskB5SZmRWSA8rMzArJAWUeZm5mheSAMjOzQnJAmZlZITmgBjCPMjezInNAmc9BmVkhOaDMzKyQHFBmZlZIDigzMyskB5SZmRWSA8rMzArJATWAha+pa2YF5oAyz2ZuZoXkgDIzs0JyQA1gnknCzIrMAWWeScLMCskBZWZmheSAMjOzQnJADWA+BWVmReaAMg8yN7NCckCZmVkhOaAGsPA4czMrsFwDStJfSQpJY9J9SbpK0jxJT0k6tGLZqZLmpq+pFe3vl/R0es5VkgdNbynvMTMrotwCStIE4CPAKxXNpwCT0te5wDVp2dHAZcARwOHAZZJGpedck5bteN7kvqjfzMwaK88jqO8DX2fTwWRTgJsj8xAwUtLOwMnA9IhYFhHLgenA5PTY8Ij4Q2T9VTcDZ/TtZpiZWSPkElCSTgdejYgnuzw0DlhQcX9hauuufWGV9s297rmSZkqa2dbWVsMWmJlZo7U0asWSfg3sVOWhS4FLgJOqPa1KW/SivaqIuBa4FqC1tdUjBMzMCqxhARURJ1Zrl/Q+YHfgyTSeYTzwmKTDyY6AJlQsPh5YlNqP79J+f2ofX2V5MzPr5/q8iy8ino6IHSNiYkRMJAuZQyPiNWAacHYazXcksDIiFgN3AydJGpUGR5wE3J0eWyXpyDR672zgzr7epv7Kh5BmVmQNO4LqpbuAU4F5wBrgswARsUzS3wIz0nLfiohl6fZ5wI3AMOCX6cu2iMeZm1nx5B5Q6Siq43YA529muRuAG6q0zwQOaFR9ZmaWD88kYWZmheSAGsh8EsrMCswBZZ7qyMwKyQFlZmaF5IAawNzDZ2ZF5oAyDzI3s0JyQJmZWSE5oMzMrJAcUAOYr6hrZkXmgDJ8EWIzKyIHlJmZFZIDyszMCskBZWZmheSAMjOzQnJAmZlZITmgBjAPMjezInNAmac6MrNCckCZmVkhOaAGME8kYWZF5oAyX7DQzArJAWVmZoXkgDIzs0JyQA1g4YHmZlZgDihDHmhuZgXkgDIzs0JyQA1gHmZuZkXmgDIPMzezQsotoCR9UdIcSbMkfbei/WJJ89JjJ1e0T05t8yRdVNG+u6SHJc2VdJukwX29LWZmVn+5BJSkDwFTgAMjYn/gH1L7fsBZwP7AZOCHkpolNQNXA6cA+wGfTMsCfAf4fkRMApYD5/TpxpiZWUPkdQR1HvDtiFgHEBFLUvsU4NaIWBcRLwLzgMPT17yImB8R64FbgSmSBHwYuCM9/ybgjD7cjn7te2ceROtuo9hpxNC8SzEze4e8Ampv4LjUNfffkg5L7eOABRXLLUxtm2vfAVgREe1d2quSdK6kmZJmtrW11WlT+q+j9tyBO847mkHNPhVpZsXT0qgVS/o1sFOVhy5NrzsKOBI4DLhd0h5Uv/JDUD1Io5vlq4qIa4FrAVpbWz2GzcyswBoWUBFx4uYek3Qe8LOICOARSWVgDNkR0ISKRccDi9Ltau1LgZGSWtJRVOXyZmbWj+XVt/P/yc4dIWlvYDBZ2EwDzpI0RNLuwCTgEWAGMCmN2BtMNpBiWgq4+4Az03qnAnf26ZaYmVlDNOwIqgc3ADdIegZYD0xNYTNL0u3As0A7cH5ElAAkfQG4G2gGboiIWWldFwK3SroCeBy4vm83xczMGkExQKcTaG1tjZkzZ+ZdhpnZgCfp0Yho7dru4VtmZlZIDigzMyskB5SZmRWSA8rMzAppwA6SkNQGvFzDKsaQDY0vqqLXB8Wvsej1QfFrLHp9UPwaB0J9u0XE2K6NAzagaiVpZrVRJ0VR9Pqg+DUWvT4ofo1Frw+KX+NArs9dfGZmVkgOKDMzKyQHVO9dm3cBPSh6fVD8GoteHxS/xqLXB8WvccDW53NQZmZWSD6CMjOzQnJAmZlZITmgtpCkyZLmSJon6aKca3lJ0tOSnpA0M7WNljRd0tz0fVRql6SrUt1PSTq0AfXcIGlJmqW+o22L65E0NS0/V9LUPqjxckmvpv34hKRTKx67ONU4R9LJFe0N+T2QNEHSfZJmS5ol6cupvRD7sZv6irQPh0p6RNKTqcZvpvbdlV3Fe66k29Kle0iX97kt1fGwpIk91d6g+m6U9GLFPjw4tef1t9Is6XFJ/5Xu9/3+iwh/vcsvskt9vADsQXYNqyeB/XKs5yVgTJe27wIXpdsXAd9Jt08Ffkl2FeIjgYcbUM8HgEOBZ3pbDzAamJ++j0q3RzW4xsuBv6qy7H7pZzwE2D397Jsb+XsA7Awcmm5vDzyf6ijEfuymviLtQwHbpduDgIfTvrkdOCu1/ytwXrr9F8C/pttnAbd1V3sD67sROLPK8nn9rXwV+DHwX+l+n+8/H0FtmcOBeRExPyLWA7cCU3KuqaspwE3p9k3AGRXtN0fmIbIrEe9czxeOiN8Cy2qs52RgekQsi4jlwHRgcoNr3JwpwK0RsS4iXgTmkf0ONOz3ICIWR8Rj6fYqYDYwjoLsx27q25w89mFExFvp7qD0FWQXSb0jtXfdhx379g7gBEnqpvZG1bc5ff63Imk8cBpwXbovcth/DqgtMw5YUHF/Id3/cTZaAPdIelTSuantPRGxGLJ/JsCOqT2v2re0nrzq/ELqPrmho/ss7xpTV8khZO+wC7cfu9QHBdqHqXvqCWAJ2T/uF4AVEdFe5fU6a0mPrwR2aGSNXeuLiI59eGXah9+XNKRrfV3qaOQ+/AHwdaCc7u9ADvvPAbVlVKUtz3H6x0TEocApwPmSPtDNskWrfXP15FHnNcCewMHAYuD/pPbcapS0HfAfwAUR8WZ3i26mlobWWKW+Qu3DiChFxMHAeLJ37e/t5vX6vMau9Uk6ALgY2Bc4jKzb7sI86pP0UWBJRDxa2dzNazWsPgfUllkITKi4Px5YlFMtRMSi9H0J8HOyP8TXO7ru0vclafG8at/Sevq8zoh4Pf3DKAM/YmM3RC41ShpE9s//loj4WWouzH6sVl/R9mGHiFgB3E927makpJYqr9dZS3p8BFk3cMNrrKhvcuo+jYhYB/wb+e3DY4DTJb1E1vX6YbIjqr7ff/U4mTZQvoAWshORu7PxxO7+OVqk+voAAAS1SURBVNWyLbB9xe3fk/U/f49NT6Z/N90+jU1PtD7SoLomsukAhC2qh+yd44tkJ31HpdujG1zjzhW3v0LWbw6wP5ue5J1PdnK/Yb8HaX/cDPygS3sh9mM39RVpH44FRqbbw4AHgI8CP2XTk/x/kW6fz6Yn+W/vrvYG1rdzxT7+AfDtAvytHM/GQRJ9vv/qtiED5YtsRM3zZH3al+ZYxx7ph/8kMKujFrK+33uBuen76NQu4OpU99NAawNq+glZ984GsndP5/SmHuBzZCdU5wGf7YMa/z3V8BQwjU3/2V6aapwDnNLo3wPgWLJukKeAJ9LXqUXZj93UV6R9eCDweKrlGeAbFX8zj6T98VNgSGofmu7PS4/v0VPtDarvN2kfPgP8PzaO9MvlbyWt/3g2BlSf7z9PdWRmZoXkc1BmZlZIDigzMyskB5SZmRWSA8rMzArJAWVmZoXkgDJ7lySVKmaafqJy1uatgaRDJHXMvfYZSf/S5fH7JbV28/xbJU1qdJ02cLT0vIiZJW9HNj1NVZJaYuNcZf3RJcAVNTz/GrL52/53fcqxgc5HUGY1SEcaP5X0n8A9qe1rkmakST+/WbHspem6OL+W9BNJf5XaO49MJI1JU8x0TCj6vYp1/VlqPz495w5Jz0m6Jc0ejaTDJP0+XWvoEUnbS3qg49pCaZkHJR3YZTu2Bw6MiCffxTafXnEUOUfSi+mhB4ATK6bDMauJf5HM3r1haQZqgBcj4mPp9lFk/9yXSToJmEQ2j5qAaWkS39Vk08AcQvZ39xjwKN07B1gZEYelma0flHRPeuwQsqlkFgEPAsdIegS4DfhERMyQNBx4m+ySCZ8BLpC0N9kMAE91ea1WshkMKn1C0rEV9/cCiIhpZLNFIOl24L9Te1nSPOCgd7FtZj1yQJm9e5vr4pseER3XmDopfT2e7m9HFljbAz+PiDUAkqa9i9c7CThQ0pnp/oi0rvVk87EtTOt6gmx+wZXA4oiYARBpFnRJPwX+RtLXyKbGubHKa+0MtHVpuy0ivtBxR9L9lQ9K+jrZPrm6onkJsAsOKKsDB5RZ7VZX3Bbw9xHxfysXkHQBm7/UQDsbu9uHdlnXFyPi7i7rOh5YV9FUIvtbVrXXiIg1kqaTXUDuj8mOlrp6u8trd0vSCcDHya5QXGloWpdZzXwOyqy+7gY+l66XhKRxknYEfgt8TNKwdL7njyqe8xLw/nT7zC7rOi9d3gJJe0vatpvXfg7YRdJhafntK84HXQdcBcyoONqrNJvUhdcTSbsBPwT+OCK6htHeZJMXm9XMR1BmdRQR90h6L/CHNG7hLeBPI+IxSbeRzf79MtmAgg7/ANwu6dNkM1p3uI6s6+6xNAiijY2X2a722uslfQL4Z0nDyI5kTgTeiohHJb1Jdp2has99TtIISdtHdin37nyGbHb1n6dtXBQRp0p6D1mX3+Ienm/2rng2c7McSLqcLDj+oY9ebxeyC+PtG9lFBast8xVgVURc18vX+ArwZkRc3+tCzSq4i89sKyfpbOBhsmsuVQ2n5Bo2Pbe1pVYAN9XwfLNN+AjKzMwKyUdQZmZWSA4oMzMrJAeUmZkVkgPKzMwKyQFlZmaF9D/doqr+A6PWngAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "seg_dct = segment.make_dct()\n", "seg_dct.plot(high=4000)\n", "decorate(xlabel='Frequency (Hz)', ylabel='DCT')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are only a few harmonics with substantial amplitude, and many entries near zero.\n", "\n", "The following function takes a DCT and sets elements below `thresh` to 0." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "def compress(dct, thresh=1):\n", " count = 0\n", " for i, amp in enumerate(dct.amps):\n", " if np.abs(amp) < thresh:\n", " dct.hs[i] = 0\n", " count += 1\n", " \n", " n = len(dct.amps)\n", " print(count, n, 100 * count / n, sep='\\t')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we apply it to the segment, we can eliminate more than 90% of the elements:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20457\t22050\t92.77551020408163\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdFklEQVR4nO3deZhU9Z3v8ffXZhFxYVfCYiOSKKhR7KB5TLxuA4hGfGZ0HpLckcd4L3ciOsksNwGdJ2QSnZjEe10mLsMoiTreqGMWuYoh4DJj7ijQuLBqaEGlZdVmUxRo+nv/qF8diqa6m6pzqk+f7s/refrpql+dU+fbp7vrU7/f75w65u6IiIgAHJF2ASIi0nEoFEREJKJQEBGRiEJBREQiCgUREYl0S7uAuAYMGODV1dVplyEikilLly79wN0HNm/PfChUV1dTW1ubdhkiIpliZu8Wa9fwkYiIRBQKIiISUSiIiEgkkVAwsz5m9qSZvWlmq83si2bWz8wWmNma8L1vWNbM7G4zqzOzZWY2tuB5pobl15jZ1CRqExGRw5dUT+Eu4HfufgrweWA1MAN4zt1HAc+F+wCXAqPC1zTgPgAz6wfMAs4BxgGz8kEiIiLtI3YomNmxwPnAgwDuvtfdtwOTgYfCYg8BV4bbk4GHPecVoI+ZDQYmAAvcvcHdtwELgIlx6xMRkcOXRE/hJGAr8HMze83MHjCz3sDx7r4RIHwfFJYfAqwvWL8+tLXUfggzm2ZmtWZWu3Xr1gR+BBERgWRCoRswFrjP3c8CPubAUFExVqTNW2k/tNF9trvXuHvNwIGHnHvR6X3w0R5+t2Jj2mWISCeURCjUA/Xuvijcf5JcSGwOw0KE71sKlh9WsP5QYEMr7dLMtT9fwl/+66vs+GRf2qWISCcTOxTcfROw3sw+F5ouBlYBc4H8EURTgafC7bnANeEopHOBHWF4aT4w3sz6hgnm8aFNmlm/bTcATU26QJKIJCupj7m4EXjUzHoAa4FryQXOE2Z2HfAecHVYdh4wCagDdodlcfcGM/shsCQs9wN3b0ioPhEROQyJhIK7vw7UFHno4iLLOjC9heeZA8xJoiYRESmdzmgWEZGIQkFERCIKBRERiSgUREQkolAQEZGIQkFERCIKhQxynbMmIhWiUBARkYhCIYOs2EcHiogkQKGQQRo+EpFKUSiIiEhEoZBBGj4SkUpRKIiISEShUIY9jfupnvEM97xQl8r2NacgIpWiUCjDR582AvDgH9alXImISLIUChmkOQURqRSFQgZp+EhEKkWhkGHqMYhI0hQKGaYeg4gkTaGQQeohiEilKBTKkPYbdPUQRKRSFAoxpP2GXT0GEUmaQiHD1GMQkaQpFEREJKJQiCHtN+oaPhKRpCUWCmZWZWavmdnT4f4IM1tkZmvM7HEz6xHae4b7deHx6oLnmBna3zKzCUnVlrSO8lqs4SMRSVqSPYVvAasL7v8YuMPdRwHbgOtC+3XANnc/GbgjLIeZjQamAGOAicC9ZlaVYH0iItKGRELBzIYClwEPhPsGXAQ8GRZ5CLgy3J4c7hMevzgsPxl4zN33uPs6oA4Yl0R9nZWGj0QkaUn1FO4EvgM0hfv9ge3u3hju1wNDwu0hwHqA8PiOsHzUXmSdg5jZNDOrNbParVu3JvQjHD6N2ohIZxU7FMzscmCLuy8tbC6yqLfxWGvrHNzoPtvda9y9ZuDAgSXVm6S036hrTkFEktYtgec4D7jCzCYBRwLHkus59DGzbqE3MBTYEJavB4YB9WbWDTgOaChozytcR0RE2kHsnoK7z3T3oe5eTW6i+Hl3/zrwAnBVWGwq8FS4PTfcJzz+vLt7aJ8Sjk4aAYwCFsetT0REDl8SPYWWfBd4zMxuAV4DHgztDwKPmFkduR7CFAB3X2lmTwCrgEZgurvvr2B9IiLSTKKh4O4vAi+G22spcvSQu38KXN3C+rcCtyZZU2emKQURSZrOaBYRkYhCQUREIgoFERGJKBQyyHWCgohUiEJBREQiCgUREYkoFDJIg0ciUikKBRERiSgUyqB5XhHprBQKMaR2PQOFkohUiEJBREQiCgUREYkoFGLQ3IKIdDYKhTKkfW1kZZGIVIpCQUREIgqFMnSUYSN9BpKIJE2hEEPaw0giIklTKIiISEShICIiEYWCiIhEFAoZpAlmEakUhYKIiEQUCmVwnT4mIp2UQiGDFEkiUimxQ8HMhpnZC2a22sxWmtm3Qns/M1tgZmvC976h3czsbjOrM7NlZja24LmmhuXXmNnUuLVViqETFESkc0qip9AI/K27nwqcC0w3s9HADOA5dx8FPBfuA1wKjApf04D7IBciwCzgHGAcMCsfJCIi0j5ih4K7b3T3V8PtXcBqYAgwGXgoLPYQcGW4PRl42HNeAfqY2WBgArDA3RvcfRuwAJgYt77OSAcfiUilJDqnYGbVwFnAIuB4d98IueAABoXFhgDrC1arD20ttRfbzjQzqzWz2q1btyb5I4iIdGmJhYKZHQ38Cvi2u+9sbdEibd5K+6GN7rPdvcbdawYOHFh6sSIiUlQioWBm3ckFwqPu/uvQvDkMCxG+bwnt9cCwgtWHAhtaae9wdEiqiHRWSRx9ZMCDwGp3/98FD80F8kcQTQWeKmi/JhyFdC6wIwwvzQfGm1nfMME8PrR1YOkchaRQEpFK6ZbAc5wH/AWw3MxeD203AbcBT5jZdcB7wNXhsXnAJKAO2A1cC+DuDWb2Q2BJWO4H7t6QQH2dlqJBRJIWOxTc/Q+0/Jb54iLLOzC9heeaA8yJW1Nnp/MkRKRSdEZzLOm8V9fwkYhUikKhDHqnLiKdlUJBREQiCgUREYkoFMqQ9pi+PuZCRCpFoRBL15xb2NvYROP+prTLEJEKUChIyT77989yxc/+X9pliEgFKBQyqCOMHq3a2NrHW4lIVikUMmhvo4ZuRKQyFAoiIhJRKKRgb2MTN/9mOR98tCeV7b+5aSfzlm9MZdtNTU71jGf4p+fWpLJ9EWmdQiEFz67YyKOL3uOWp1elsv2Jd77E9Y++msq29zXlhr7+6fm6VLYvIq1TKKRof0eYMRYRKaBQSMERlju/oakLn4WW9gmAIlKcQiGGcl/U86HgMUOhC2eKiFSIQqEM7324G4CGj/eWtX7IBJp0ZKmIdDAKhTJs3RXvqKEjQihoCEVEOhqFQhmOOCLuZx7l5xTi1yIikiSFQhkGHN0j1vpRT6ELhkJX/JlFskShUIY+R+VCobr/UWWtbwlNNIuIJE2hUIamMO7TWOb4z4E5BRGRjkWhUIZnV2wCoH7bJ2Wtnz8kdX8XnlRIq5O0Y/c+fr9yUzobF8kAhUIZyj0UNe/97bkw2ba79OdZs3lXrG2nLe0Rsxt++SrTHlnKph2fpluISAelUCjDL/7znVjrz/6PtQAsq99R8rpT5yyOte2ubt0HHwOwT1eOEylKoZCC9xp2l73uhpTf4cY+C7uDXN86/mHFIodvT+N+lpfxJjANHS4UzGyimb1lZnVmNiPtejqKHbv3UT3jmUSfs5w5jd179yey7XIn6ePK9xCqrLxQ+HhPY6zhw9NnzefP73+57PXjemXth9GBElnzzgcfU7+t/DdUH+9pjHVt8RXv7yj7TdH3fruSr/zsD2zYXt485Ed7GqNebqV1a5etHCYzqwLuAf4EqAeWmNlcd0/8M6YfeGkttzyzOtZzfHnUAODQd8+NTc7CVZvZtnsfz7+5mdUbd7Hzk33s2tN4yHPEeaFv7V33+obdHHdUd449svvB6xTUOvKmeZzY/yjGVfdj/JgTGDmwN8f26s7jS9azauNOvj5uOI1NzrsffszLaz/kdys2HXLCXbF/Evdcb+jEIofsFi6+69N9HN3z4D/BLbv2sOL9HRzbqzs1J/Y96LGtH+3BHV57bzsTxhxf9Oducli4ejMjB/Zm5MCjD3nuLeFs9O5V1uI/+Pbd++jdsxvdqw4NjjGz5gOw7keTiq67e+9+HOjdo6ro47v2NLL4nYayX1xue/ZNzhrehwljTih53Zff/pCvPbCI//7lEdw06dSS13+jfgfTH32VZ/7qSxzXq3vbKzTz979dwfL3d/DU9PNKXhfggttfBFre920ZM2s+Xzp5AI9cN67kde9cuIa7nlvDn40dyu1Xn1Hy+ovfaQBg997Gsn73F/z0BT74aC9r/3ESzd/PWJlvcFpiHelYeTP7IvB9d58Q7s8EcPcftbROTU2N19bWlrytpN91i4i0tzd/OJEjuxd/A9IWM1vq7jXN2ztUTwEYAqwvuF8PnNN8ITObBkwDGD58eFkbuursoTy5tL6sdS/43EBefGvrIe3/47+cRK/uVdRv+6Ts5+5KTh9yHBefOuigtjsXHrgi27cvGRXddocnl9ZHR24B/OlZQxjerDfS0vqQ60GseH9ndL9H1RH07H4EZw7rw9kn9uXTfU3c/+9vA/CVz3+GkQN7H1Jza88P8MjL7/Lhx3v5q4tOLjpvkV+/2LoAi9Y20KtHFWcMPe6Qx9zhrudaXr+pyZm/cjOnDD6GEQMOrT3/d9mj6giuv3DkIY9v372PJe808MWT+nP0kYe+NORrb2nftCW//rXnVZfV02hr37Vmb2MT9774dtnrP7Z4PZt2flr2+m393Rzu+s3/rrpVYG6so/UUrgYmuPt/C/f/Ahjn7je2tE65PYU49jc5I2+aB8CCvz6fUccfc9jrujsjZs6L7r9z22Ulb//C218sOr541dlD+cZ5I2hsamLTjk8Z38IQQ76XVM62C9dffPPFDDrmyJLW3bprD1+4dSH3fG0sl50xuORt54f9rj2vmllfGVPy+pPueolVG3dy79fHMun00rcfd999EuZkerUwvNSRXXrXS6zeuJOnb/wSpw05NLTaEnffrdywA3fK2nZTk3PSTfP46rjh/OhPTy95/f/7xgZu/OVrXHb6YO75+tiS13/3w4+Z8avl/PM1Zx8ypHs4Fq7azFubdzH9wpNLXrclWekp1APDCu4PBTakVMthKSUQ4ODxv369432GUl6xf7Izhiby1K0qNRAABh7Tk3U/mhR7HNQob/1jwjvgPkeV/o8JuTcBqzbubHvBFmQxDPK+d/loZv562SFzNe1lzGdKD4O8I44wVvzDBHqVOdQSXQOlzKPnTuzfm19OO7esdQEuGX08l4wuPo+WtI4WCkuAUWY2AngfmAJ8Ld2SKqf5JGs5Xp55UQKVtK+kJ8ZKke8Yl3v00ajjjyn5jUBn8cWR/Xnxf15Y9voL/+Z83tr0UYIVlSbO/9tpQ44F4LLTP5NUOR1WhwoFd280sxuA+UAVMMfdV6ZcVsWU+9pYeMW3wcf1SqiariG/73SeQvs7edAxnDwom4F6Yv/eZQ97ZU2HCgUAd58HzGtzwS6sA00DZU4UCsoEkaI63MlrWZDU60m5L+5pnxWcZf169wSgZ7fsju2LVFKH6yl0BeNG9GPxugbOGt6nrPV1befy3X71GTy7YlNZR7CIdAXqKaTg/HAm9JA+mg8oVf7Il1NOKG9sus9RPfjquPLObRHpCtRTyKCOdG5Je7vwlEE8+60vlx0KItI6hUIK4h6SGTcSHv7GOLaGzwAqx/3/9eyKnEl5uE4dfGxq2xbp7BQKGbRvf7xYOP+zA2OtP/G00j+MTUSyQXMKKYg7/LNnXzIfXy0i0pxCIYNOCh9GFvcdv4hIcwqFFMSdUxjWL/fJoH82dkgS5YiIRBQKGZQffErzM4REpHNSKJQh/1p8/QWHfib94Yh9SGlYXZEgIknT0UdlMLNUPxwr/zEX6iiISNLUU0hBUsM+5V5TQESkJQqFDOrCJzSLSIUpFFIQd04hv7qGj0QkaQqFFJX7oh7NKSRYi4gIKBQyST0FEakUhUIKEjoiFfUVRCRpCoUURCeflfmirp6CiFSKQiFF5b+oa05BRCpDoZCC2MNHUU9BsSAiyVIopCDu0UMHhp9ERJKlUMig/HkO6iiISNIUChmmUBCRpMUKBTP7qZm9aWbLzOw3Ztan4LGZZlZnZm+Z2YSC9omhrc7MZhS0jzCzRWa2xsweN7MecWrryPzAZ1+Xt35ypYiIHCRuT2EBcJq7nwH8EZgJYGajgSnAGGAicK+ZVZlZFXAPcCkwGvhqWBbgx8Ad7j4K2AZcF7O2DivunEA00axZBRFJWKxQcPffu3tjuPsKMDTcngw85u573H0dUAeMC1917r7W3fcCjwGTLXcYzUXAk2H9h4Ar49TWkY0Ml9PMX1azbMoEEUlYktdT+AbweLg9hFxI5NWHNoD1zdrPAfoD2wsCpnD5Q5jZNGAawPDhw2MX3t6u+PxnOLF/bz4/9Li0SxEROUiboWBmC4ETijx0s7s/FZa5GWgEHs2vVmR5p3jPxFtZvih3nw3MBqipqcncELuZceawPm0vKCLSztoMBXe/pLXHzWwqcDlwsR/4TOh6YFjBYkOBDeF2sfYPgD5m1i30FgqXFxGRdhL36KOJwHeBK9x9d8FDc4EpZtbTzEYAo4DFwBJgVDjSqAe5yei5IUxeAK4K608FnopTm4iIlC7unMLPgJ7AgvCRC6+4+1+6+0ozewJYRW5Yabq77wcwsxuA+UAVMMfdV4bn+i7wmJndArwGPBizNhERKVGsUHD3k1t57Fbg1iLt84B5RdrXkjs6SdqQuUkUEckMndGcYToiVUSSplAQEZGIQiGDPO5nb4uItEChkGG6noKIJE2hICIiEYWCiIhEFAoiIhJRKGSYZhREJGkKBRERiSgUMkhHpIpIpSgUMkxHpIpI0hQKGaYeg4gkTaGQQeohiEilKBQySD0EEakUhUKGqccgIklTKIiISEShICIiEYVCBrmuvSYiFaJQyDDTB12ISMIUCiIiElEoZJAOSRWRSlEoZJgOSRWRpCkUREQkolAQEZFIIqFgZn9nZm5mA8J9M7O7zazOzJaZ2diCZaea2ZrwNbWg/WwzWx7Wudt0VfoWaU5BRColdiiY2TDgT4D3CpovBUaFr2nAfWHZfsAs4BxgHDDLzPqGde4Ly+bXmxi3ts5OqSkiSUuip3AH8B046IyqycDDnvMK0MfMBgMTgAXu3uDu24AFwMTw2LHu/rK7O/AwcGUCtYmISAlihYKZXQG87+5vNHtoCLC+4H59aGutvb5Iu4iItKNubS1gZguBE4o8dDNwEzC+2GpF2ryM9pZqmkZuqInhw4e3tJiIiJSozVBw90uKtZvZ6cAI4I0wJzwUeNXMxpF7pz+sYPGhwIbQfkGz9hdD+9Aiy7dU02xgNkBNTY2mXUVEElL28JG7L3f3Qe5e7e7V5F7Yx7r7JmAucE04CulcYIe7bwTmA+PNrG+YYB4PzA+P7TKzc8NRR9cAT8X82UREpERt9hTKNA+YBNQBu4FrAdy9wcx+CCwJy/3A3RvC7W8CvwB6Ac+GLylCn5IqIpWSWCiE3kL+tgPTW1huDjCnSHstcFpS9XQJOiZVRBKmM5pFRCSiUMggndEsIpWiUMgwXWRHRJKmUBARkYhCQUREIgqFDNKUgohUikIhw/Th4iKSNIWCiIhEFAoiIhJRKGSRJhVEpEIUCiIiElEoZJEmmEWkQhQKWaThIxGpEIVChqnDICJJUyiIiEhEoSAiIhGFQgbpymsiUikKhQwzfc6FiCRMoSAiIhGFQgbpymsiUikKhQzT6JGIJE2hICIiEYWCiIhEFAoZpCkFEakUhUKGaUpBRJIWOxTM7EYze8vMVprZTwraZ5pZXXhsQkH7xNBWZ2YzCtpHmNkiM1tjZo+bWY+4tYmISGlihYKZXQhMBs5w9zHA7aF9NDAFGANMBO41syozqwLuAS4FRgNfDcsC/Bi4w91HAduA6+LU1pkd2T33a9PJayKStLg9hW8Ct7n7HgB33xLaJwOPufsed18H1AHjwledu691973AY8Bky726XQQ8GdZ/CLgyZm2d1h1/fiY3XnQyY4f3SbsUEelk4obCZ4Evh2GffzezL4T2IcD6guXqQ1tL7f2B7e7e2Ky9KDObZma1Zla7devWmD9C9gw69kj+dvzn1FMQkcR1a2sBM1sInFDkoZvD+n2Bc4EvAE+Y2UkUnwN1ioeQt7J8Ue4+G5gNUFNTo4NxREQS0mYouPslLT1mZt8Efu3uDiw2syZgALl3+sMKFh0KbAi3i7V/APQxs26ht1C4vIiItJO4w0e/JTcXgJl9FuhB7gV+LjDFzHqa2QhgFLAYWAKMCkca9SA3GT03hMoLwFXheacCT8WsTUREStRmT6ENc4A5ZrYC2AtMDS/wK83sCWAV0AhMd/f9AGZ2AzAfqALmuPvK8FzfBR4zs1uA14AHY9YmIiIlMs/4R27W1NR4bW1t2mWIiGSKmS1195rm7TqjWUREIgoFERGJKBRERCSS+TkFM9sKvFvm6gPIHS3V0aiu0qiu0qiu0nTWuk5094HNGzMfCnGYWW2xiZa0qa7SqK7SqK7SdLW6NHwkIiIRhYKIiES6eijMTruAFqiu0qiu0qiu0nSpurr0nIKIiBysq/cURESkgEJBREQiXTIUWrpOdDtu/x0zW25mr5tZbWjrZ2YLwjWqF5hZ39BuZnZ3qHWZmY1NuJY5ZrYlfKhhvq3kWsxsalh+jZlNrVBd3zez98N+e93MJhU8VtI1wcusaZiZvWBmq8M1yb8V2lPdX63Ulfb+OtLMFpvZG6GufwjtI6zI9djDpyo/Hra9yMyq26o34bp+YWbrCvbXmaG93f7uw3NWmdlrZvZ0uN+++8vdu9QXuU9nfRs4idxHfb8BjG7nGt4BBjRr+wkwI9yeAfw43J4EPEvuQkTnAosSruV8YCywotxagH7A2vC9b7jdtwJ1fR/4uyLLjg6/x57AiPD7rUr6dw0MBsaG28cAfwzbTnV/tVJX2vvLgKPD7e7AorAfngCmhPb7gW+G29cD94fbU4DHW6u3AnX9AriqyPLt9ncfnvdvgP8DPB3ut+v+6oo9haLXiU65JsjV8FC4XXiN6snAw57zCrmLEQ1OaqPu/h9AQ8xaJgAL3L3B3bcBC4CJFairJSVdEzxGTRvd/dVwexewmtxlY1PdX63U1ZL22l/u7h+Fu93Dl9Py9dgL9+OTwMVmZq3Um3RdLWm3v3szGwpcBjwQ7rd2/fqK7K+uGAotXSe6PTnwezNbambTQtvx7r4Rcv/kwKDQnka9pdbSnjXeELrwc/LDNGnUFbrqZ5F7l9lh9lezuiDl/RWGQl4HtpB70Xyblq/HHm0/PL6D3PXbK16Xu+f3161hf91hZj2b19Vs+5X4Pd4JfAdoCvdbu359RfZXVwyFkq4HXSHnuftY4FJgupmd38qyHaHevJZqaa8a7wNGAmcCG4H/lUZdZnY08Cvg2+6+s7VFU64r9f3l7vvd/Uxyl9gdB5zayjZSq8vMTgNmAqeQu958P3IX/mq3uszscmCLuy8tbG5lGxWpqyuGQmvXj24X7r4hfN8C/IbcP8vm/LBQ+L4lLJ5GvaXW0i41uvvm8M/cBPwLB7rE7VaXmXUn98L7qLv/OjSnvr+K1dUR9leeu28HXiQ3Jt/HzPJXfSzcRrT98Phx5IYQ26OuiWEYzt19D/Bz2n9/nQdcYWbvkBu6u4hcz6F991fcSZGsfZG7BOlachMw+cm0Me24/d7AMQW3/5PcOORPOXiy8ifh9mUcPMm1uAI1VXPwhG5JtZB7V7WO3GRb33C7XwXqGlxw+6/JjZsCjOHgibW15CZNE/1dh5/7YeDOZu2p7q9W6kp7fw0E+oTbvYCXgMuBf+PgidPrw+3pHDxx+kRr9VagrsEF+/NO4LY0/u7Dc1/AgYnmdt1fib64ZOWL3NEEfyQ3vnlzO2/7pPALewNYmd8+ubHA54A14Xu/gj/Qe0Kty4GahOv5JbmhhX3k3mFcV04twDfITWjVAddWqK5HwnaXAXM5+EXv5lDXW8CllfhdA18i1w1fBrwevialvb9aqSvt/XUGueutLwNWAN8r+B9YHH72fwN6hvYjw/268PhJbdWbcF3Ph/21AvhXDhyh1G5/9wXPewEHQqFd95c+5kJERCJdcU5BRERaoFAQEZGIQkFERCIKBRERiSgUREQkolAQEZGIQkFERCL/HyCGxkibr55mAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "seg_dct = segment.make_dct()\n", "compress(seg_dct, thresh=10)\n", "seg_dct.plot(high=4000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the result sounds the same (at least to me):" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "seg2 = seg_dct.make_wave()\n", "seg2.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compress a longer segment, we can make a DCT spectrogram. The following function is similar to `wave.make_spectrogram` except that it uses the DCT." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from thinkdsp import Spectrogram\n", "\n", "def make_dct_spectrogram(wave, seg_length):\n", " \"\"\"Computes the DCT spectrogram of the wave.\n", "\n", " seg_length: number of samples in each segment\n", "\n", " returns: Spectrogram\n", " \"\"\"\n", " window = np.hamming(seg_length)\n", " i, j = 0, seg_length\n", " step = seg_length // 2\n", "\n", " # map from time to Spectrum\n", " spec_map = {}\n", "\n", " while j < len(wave.ys):\n", " segment = wave.slice(i, j)\n", " segment.window(window)\n", "\n", " # the nominal time for this segment is the midpoint\n", " t = (segment.start + segment.end) / 2\n", " spec_map[t] = segment.make_dct()\n", "\n", " i += step\n", " j += step\n", "\n", " return Spectrogram(spec_map, seg_length)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can make a DCT spectrogram and apply `compress` to each segment:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1018\t1024\t99.4140625\n", "1016\t1024\t99.21875\n", "1014\t1024\t99.0234375\n", "1017\t1024\t99.31640625\n", "1016\t1024\t99.21875\n", "1017\t1024\t99.31640625\n", "1016\t1024\t99.21875\n", "1020\t1024\t99.609375\n", "1014\t1024\t99.0234375\n", "1005\t1024\t98.14453125\n", "1009\t1024\t98.53515625\n", "1015\t1024\t99.12109375\n", "1015\t1024\t99.12109375\n", "1016\t1024\t99.21875\n", "1016\t1024\t99.21875\n", "1015\t1024\t99.12109375\n", "1017\t1024\t99.31640625\n", "1020\t1024\t99.609375\n", "1013\t1024\t98.92578125\n", "1017\t1024\t99.31640625\n", "1013\t1024\t98.92578125\n", "1017\t1024\t99.31640625\n", "1018\t1024\t99.4140625\n", "1015\t1024\t99.12109375\n", "1013\t1024\t98.92578125\n", "794\t1024\t77.5390625\n", "785\t1024\t76.66015625\n", "955\t1024\t93.26171875\n", "995\t1024\t97.16796875\n", "992\t1024\t96.875\n", "976\t1024\t95.3125\n", "925\t1024\t90.33203125\n", "802\t1024\t78.3203125\n", "836\t1024\t81.640625\n", "850\t1024\t83.0078125\n", "882\t1024\t86.1328125\n", "883\t1024\t86.23046875\n", "891\t1024\t87.01171875\n", "901\t1024\t87.98828125\n", "902\t1024\t88.0859375\n", "900\t1024\t87.890625\n", "900\t1024\t87.890625\n", "894\t1024\t87.3046875\n", "904\t1024\t88.28125\n", "901\t1024\t87.98828125\n", "915\t1024\t89.35546875\n", "913\t1024\t89.16015625\n", "899\t1024\t87.79296875\n", "905\t1024\t88.37890625\n", "905\t1024\t88.37890625\n", "888\t1024\t86.71875\n", "898\t1024\t87.6953125\n", "879\t1024\t85.83984375\n", "893\t1024\t87.20703125\n", "893\t1024\t87.20703125\n", "882\t1024\t86.1328125\n", "874\t1024\t85.3515625\n", "876\t1024\t85.546875\n", "864\t1024\t84.375\n", "879\t1024\t85.83984375\n", "869\t1024\t84.86328125\n", "872\t1024\t85.15625\n", "871\t1024\t85.05859375\n", "878\t1024\t85.7421875\n", "872\t1024\t85.15625\n", "859\t1024\t83.88671875\n", "879\t1024\t85.83984375\n", "889\t1024\t86.81640625\n", "872\t1024\t85.15625\n", "837\t1024\t81.73828125\n", "842\t1024\t82.2265625\n", "825\t1024\t80.56640625\n", "839\t1024\t81.93359375\n", "796\t1024\t77.734375\n", "792\t1024\t77.34375\n", "769\t1024\t75.09765625\n", "836\t1024\t81.640625\n", "919\t1024\t89.74609375\n", "913\t1024\t89.16015625\n", "942\t1024\t91.9921875\n", "837\t1024\t81.73828125\n", "739\t1024\t72.16796875\n", "737\t1024\t71.97265625\n", "726\t1024\t70.8984375\n", "728\t1024\t71.09375\n", "733\t1024\t71.58203125\n", "717\t1024\t70.01953125\n", "716\t1024\t69.921875\n", "676\t1024\t66.015625\n", "712\t1024\t69.53125\n", "697\t1024\t68.06640625\n", "718\t1024\t70.1171875\n", "717\t1024\t70.01953125\n", "718\t1024\t70.1171875\n", "681\t1024\t66.50390625\n", "707\t1024\t69.04296875\n", "691\t1024\t67.48046875\n", "681\t1024\t66.50390625\n", "709\t1024\t69.23828125\n", "684\t1024\t66.796875\n", "743\t1024\t72.55859375\n", "710\t1024\t69.3359375\n", "712\t1024\t69.53125\n", "714\t1024\t69.7265625\n", "719\t1024\t70.21484375\n", "708\t1024\t69.140625\n", "725\t1024\t70.80078125\n", "700\t1024\t68.359375\n", "726\t1024\t70.8984375\n", "716\t1024\t69.921875\n", "725\t1024\t70.80078125\n", "692\t1024\t67.578125\n", "675\t1024\t65.91796875\n", "747\t1024\t72.94921875\n", "741\t1024\t72.36328125\n", "730\t1024\t71.2890625\n", "701\t1024\t68.45703125\n", "721\t1024\t70.41015625\n", "747\t1024\t72.94921875\n", "725\t1024\t70.80078125\n", "744\t1024\t72.65625\n", "720\t1024\t70.3125\n", "716\t1024\t69.921875\n", "723\t1024\t70.60546875\n", "721\t1024\t70.41015625\n", "734\t1024\t71.6796875\n", "730\t1024\t71.2890625\n", "718\t1024\t70.1171875\n", "730\t1024\t71.2890625\n", "723\t1024\t70.60546875\n", "749\t1024\t73.14453125\n", "727\t1024\t70.99609375\n", "728\t1024\t71.09375\n", "746\t1024\t72.8515625\n", "739\t1024\t72.16796875\n", "757\t1024\t73.92578125\n", "741\t1024\t72.36328125\n", "751\t1024\t73.33984375\n", "775\t1024\t75.68359375\n", "749\t1024\t73.14453125\n", "768\t1024\t75.0\n", "763\t1024\t74.51171875\n", "771\t1024\t75.29296875\n", "758\t1024\t74.0234375\n", "745\t1024\t72.75390625\n", "756\t1024\t73.828125\n", "744\t1024\t72.65625\n", "743\t1024\t72.55859375\n", "757\t1024\t73.92578125\n", "779\t1024\t76.07421875\n", "760\t1024\t74.21875\n", "770\t1024\t75.1953125\n", "759\t1024\t74.12109375\n", "737\t1024\t71.97265625\n", "739\t1024\t72.16796875\n", "751\t1024\t73.33984375\n", "762\t1024\t74.4140625\n", "754\t1024\t73.6328125\n", "811\t1024\t79.19921875\n", "899\t1024\t87.79296875\n", "832\t1024\t81.25\n", "800\t1024\t78.125\n", "756\t1024\t73.828125\n", "748\t1024\t73.046875\n", "727\t1024\t70.99609375\n", "744\t1024\t72.65625\n", "725\t1024\t70.80078125\n", "720\t1024\t70.3125\n", "755\t1024\t73.73046875\n", "737\t1024\t71.97265625\n", "766\t1024\t74.8046875\n", "747\t1024\t72.94921875\n", "743\t1024\t72.55859375\n", "727\t1024\t70.99609375\n", "726\t1024\t70.8984375\n", "746\t1024\t72.8515625\n", "764\t1024\t74.609375\n", "751\t1024\t73.33984375\n", "734\t1024\t71.6796875\n", "741\t1024\t72.36328125\n", "760\t1024\t74.21875\n", "750\t1024\t73.2421875\n", "784\t1024\t76.5625\n", "730\t1024\t71.2890625\n", "757\t1024\t73.92578125\n", "761\t1024\t74.31640625\n", "734\t1024\t71.6796875\n", "744\t1024\t72.65625\n", "757\t1024\t73.92578125\n", "714\t1024\t69.7265625\n", "740\t1024\t72.265625\n", "738\t1024\t72.0703125\n", "763\t1024\t74.51171875\n", "766\t1024\t74.8046875\n", "745\t1024\t72.75390625\n", "751\t1024\t73.33984375\n", "759\t1024\t74.12109375\n", "756\t1024\t73.828125\n", "756\t1024\t73.828125\n", "756\t1024\t73.828125\n", "755\t1024\t73.73046875\n", "746\t1024\t72.8515625\n", "756\t1024\t73.828125\n", "738\t1024\t72.0703125\n", "757\t1024\t73.92578125\n", "764\t1024\t74.609375\n", "765\t1024\t74.70703125\n", "762\t1024\t74.4140625\n", "768\t1024\t75.0\n", "773\t1024\t75.48828125\n", "782\t1024\t76.3671875\n", "773\t1024\t75.48828125\n", "766\t1024\t74.8046875\n", "755\t1024\t73.73046875\n", "766\t1024\t74.8046875\n", "772\t1024\t75.390625\n", "810\t1024\t79.1015625\n", "739\t1024\t72.16796875\n", "717\t1024\t70.01953125\n", "722\t1024\t70.5078125\n", "739\t1024\t72.16796875\n", "725\t1024\t70.80078125\n", "736\t1024\t71.875\n", "759\t1024\t74.12109375\n", "769\t1024\t75.09765625\n", "749\t1024\t73.14453125\n", "710\t1024\t69.3359375\n", "748\t1024\t73.046875\n", "720\t1024\t70.3125\n", "732\t1024\t71.484375\n", "721\t1024\t70.41015625\n", "734\t1024\t71.6796875\n", "763\t1024\t74.51171875\n", "747\t1024\t72.94921875\n", "754\t1024\t73.6328125\n", "755\t1024\t73.73046875\n", "764\t1024\t74.609375\n", "801\t1024\t78.22265625\n", "768\t1024\t75.0\n", "780\t1024\t76.171875\n", "773\t1024\t75.48828125\n", "764\t1024\t74.609375\n", "775\t1024\t75.68359375\n", "740\t1024\t72.265625\n", "794\t1024\t77.5390625\n", "796\t1024\t77.734375\n", "769\t1024\t75.09765625\n", "751\t1024\t73.33984375\n", "782\t1024\t76.3671875\n", "758\t1024\t74.0234375\n", "777\t1024\t75.87890625\n", "794\t1024\t77.5390625\n", "784\t1024\t76.5625\n", "788\t1024\t76.953125\n", "773\t1024\t75.48828125\n", "783\t1024\t76.46484375\n", "784\t1024\t76.5625\n", "785\t1024\t76.66015625\n", "806\t1024\t78.7109375\n", "807\t1024\t78.80859375\n", "797\t1024\t77.83203125\n", "785\t1024\t76.66015625\n", "794\t1024\t77.5390625\n", "766\t1024\t74.8046875\n", "790\t1024\t77.1484375\n", "746\t1024\t72.8515625\n", "762\t1024\t74.4140625\n", "813\t1024\t79.39453125\n", "801\t1024\t78.22265625\n", "782\t1024\t76.3671875\n", "776\t1024\t75.78125\n", "755\t1024\t73.73046875\n", "780\t1024\t76.171875\n", "784\t1024\t76.5625\n", "805\t1024\t78.61328125\n", "791\t1024\t77.24609375\n", "803\t1024\t78.41796875\n", "799\t1024\t78.02734375\n", "795\t1024\t77.63671875\n", "797\t1024\t77.83203125\n", "806\t1024\t78.7109375\n", "781\t1024\t76.26953125\n", "795\t1024\t77.63671875\n", "797\t1024\t77.83203125\n", "893\t1024\t87.20703125\n", "775\t1024\t75.68359375\n", "787\t1024\t76.85546875\n", "746\t1024\t72.8515625\n", "767\t1024\t74.90234375\n", "749\t1024\t73.14453125\n", "749\t1024\t73.14453125\n", "738\t1024\t72.0703125\n", "736\t1024\t71.875\n", "747\t1024\t72.94921875\n", "760\t1024\t74.21875\n", "737\t1024\t71.97265625\n", "752\t1024\t73.4375\n", "756\t1024\t73.828125\n", "772\t1024\t75.390625\n", "740\t1024\t72.265625\n", "737\t1024\t71.97265625\n", "766\t1024\t74.8046875\n", "791\t1024\t77.24609375\n", "765\t1024\t74.70703125\n", "771\t1024\t75.29296875\n", "786\t1024\t76.7578125\n", "770\t1024\t75.1953125\n", "761\t1024\t74.31640625\n", "765\t1024\t74.70703125\n", "756\t1024\t73.828125\n", "758\t1024\t74.0234375\n", "765\t1024\t74.70703125\n", "785\t1024\t76.66015625\n", "769\t1024\t75.09765625\n", "781\t1024\t76.26953125\n", "792\t1024\t77.34375\n", "798\t1024\t77.9296875\n", "809\t1024\t79.00390625\n", "778\t1024\t75.9765625\n", "782\t1024\t76.3671875\n", "776\t1024\t75.78125\n", "791\t1024\t77.24609375\n", "794\t1024\t77.5390625\n", "783\t1024\t76.46484375\n", "771\t1024\t75.29296875\n", "792\t1024\t77.34375\n", "785\t1024\t76.66015625\n", "812\t1024\t79.296875\n", "809\t1024\t79.00390625\n", "799\t1024\t78.02734375\n", "798\t1024\t77.9296875\n", "803\t1024\t78.41796875\n", "800\t1024\t78.125\n", "805\t1024\t78.61328125\n", "803\t1024\t78.41796875\n", "799\t1024\t78.02734375\n", "802\t1024\t78.3203125\n", "804\t1024\t78.515625\n", "809\t1024\t79.00390625\n", "784\t1024\t76.5625\n", "791\t1024\t77.24609375\n", "814\t1024\t79.4921875\n", "788\t1024\t76.953125\n", "816\t1024\t79.6875\n", "810\t1024\t79.1015625\n", "820\t1024\t80.078125\n", "823\t1024\t80.37109375\n", "813\t1024\t79.39453125\n", "799\t1024\t78.02734375\n", "807\t1024\t78.80859375\n", "799\t1024\t78.02734375\n", "789\t1024\t77.05078125\n", "813\t1024\t79.39453125\n", "819\t1024\t79.98046875\n", "809\t1024\t79.00390625\n", "784\t1024\t76.5625\n", "809\t1024\t79.00390625\n", "810\t1024\t79.1015625\n", "785\t1024\t76.66015625\n", "838\t1024\t81.8359375\n", "821\t1024\t80.17578125\n", "822\t1024\t80.2734375\n", "800\t1024\t78.125\n", "815\t1024\t79.58984375\n", "827\t1024\t80.76171875\n", "820\t1024\t80.078125\n", "792\t1024\t77.34375\n", "818\t1024\t79.8828125\n", "813\t1024\t79.39453125\n", "824\t1024\t80.46875\n", "795\t1024\t77.63671875\n", "788\t1024\t76.953125\n", "796\t1024\t77.734375\n", "802\t1024\t78.3203125\n", "800\t1024\t78.125\n", "796\t1024\t77.734375\n", "823\t1024\t80.37109375\n", "804\t1024\t78.515625\n", "811\t1024\t79.19921875\n", "808\t1024\t78.90625\n", "815\t1024\t79.58984375\n", "812\t1024\t79.296875\n", "822\t1024\t80.2734375\n", "793\t1024\t77.44140625\n", "803\t1024\t78.41796875\n", "806\t1024\t78.7109375\n", "812\t1024\t79.296875\n", "796\t1024\t77.734375\n", "804\t1024\t78.515625\n", "807\t1024\t78.80859375\n", "821\t1024\t80.17578125\n", "793\t1024\t77.44140625\n", "799\t1024\t78.02734375\n", "810\t1024\t79.1015625\n", "818\t1024\t79.8828125\n", "813\t1024\t79.39453125\n", "825\t1024\t80.56640625\n", "804\t1024\t78.515625\n", "821\t1024\t80.17578125\n", "809\t1024\t79.00390625\n", "828\t1024\t80.859375\n", "813\t1024\t79.39453125\n", "838\t1024\t81.8359375\n", "836\t1024\t81.640625\n", "818\t1024\t79.8828125\n", "808\t1024\t78.90625\n", "819\t1024\t79.98046875\n", "820\t1024\t80.078125\n", "814\t1024\t79.4921875\n", "901\t1024\t87.98828125\n", "894\t1024\t87.3046875\n", "888\t1024\t86.71875\n", "780\t1024\t76.171875\n", "773\t1024\t75.48828125\n", "750\t1024\t73.2421875\n", "750\t1024\t73.2421875\n", "730\t1024\t71.2890625\n", "761\t1024\t74.31640625\n", "775\t1024\t75.68359375\n", "782\t1024\t76.3671875\n", "788\t1024\t76.953125\n", "748\t1024\t73.046875\n", "752\t1024\t73.4375\n", "771\t1024\t75.29296875\n", "746\t1024\t72.8515625\n", "778\t1024\t75.9765625\n", "777\t1024\t75.87890625\n", "760\t1024\t74.21875\n", "734\t1024\t71.6796875\n", "711\t1024\t69.43359375\n", "754\t1024\t73.6328125\n", "745\t1024\t72.75390625\n", "758\t1024\t74.0234375\n", "744\t1024\t72.65625\n", "755\t1024\t73.73046875\n", "749\t1024\t73.14453125\n", "723\t1024\t70.60546875\n", "784\t1024\t76.5625\n", "761\t1024\t74.31640625\n", "758\t1024\t74.0234375\n", "709\t1024\t69.23828125\n", "769\t1024\t75.09765625\n", "773\t1024\t75.48828125\n", "769\t1024\t75.09765625\n", "756\t1024\t73.828125\n", "747\t1024\t72.94921875\n", "787\t1024\t76.85546875\n", "770\t1024\t75.1953125\n", "749\t1024\t73.14453125\n", "769\t1024\t75.09765625\n", "748\t1024\t73.046875\n", "761\t1024\t74.31640625\n", "759\t1024\t74.12109375\n", "775\t1024\t75.68359375\n", "756\t1024\t73.828125\n", "774\t1024\t75.5859375\n", "776\t1024\t75.78125\n", "760\t1024\t74.21875\n", "783\t1024\t76.46484375\n", "744\t1024\t72.65625\n", "766\t1024\t74.8046875\n", "761\t1024\t74.31640625\n", "788\t1024\t76.953125\n", "774\t1024\t75.5859375\n", "753\t1024\t73.53515625\n", "754\t1024\t73.6328125\n", "765\t1024\t74.70703125\n", "736\t1024\t71.875\n", "782\t1024\t76.3671875\n", "768\t1024\t75.0\n", "778\t1024\t75.9765625\n", "767\t1024\t74.90234375\n", "774\t1024\t75.5859375\n", "772\t1024\t75.390625\n", "769\t1024\t75.09765625\n", "774\t1024\t75.5859375\n", "776\t1024\t75.78125\n", "796\t1024\t77.734375\n", "762\t1024\t74.4140625\n", "766\t1024\t74.8046875\n", "765\t1024\t74.70703125\n", "783\t1024\t76.46484375\n", "770\t1024\t75.1953125\n", "799\t1024\t78.02734375\n", "779\t1024\t76.07421875\n", "774\t1024\t75.5859375\n", "791\t1024\t77.24609375\n", "797\t1024\t77.83203125\n", "781\t1024\t76.26953125\n", "754\t1024\t73.6328125\n", "790\t1024\t77.1484375\n", "790\t1024\t77.1484375\n", "801\t1024\t78.22265625\n", "783\t1024\t76.46484375\n", "787\t1024\t76.85546875\n", "805\t1024\t78.61328125\n", "758\t1024\t74.0234375\n", "785\t1024\t76.66015625\n", "788\t1024\t76.953125\n", "806\t1024\t78.7109375\n", "818\t1024\t79.8828125\n", "776\t1024\t75.78125\n", "807\t1024\t78.80859375\n", "802\t1024\t78.3203125\n", "782\t1024\t76.3671875\n", "812\t1024\t79.296875\n", "803\t1024\t78.41796875\n", "803\t1024\t78.41796875\n", "787\t1024\t76.85546875\n", "799\t1024\t78.02734375\n", "786\t1024\t76.7578125\n", "813\t1024\t79.39453125\n", "813\t1024\t79.39453125\n", "813\t1024\t79.39453125\n", "803\t1024\t78.41796875\n", "815\t1024\t79.58984375\n", "792\t1024\t77.34375\n", "807\t1024\t78.80859375\n", "829\t1024\t80.95703125\n", "797\t1024\t77.83203125\n", "814\t1024\t79.4921875\n", "793\t1024\t77.44140625\n", "802\t1024\t78.3203125\n", "775\t1024\t75.68359375\n", "816\t1024\t79.6875\n", "804\t1024\t78.515625\n", "808\t1024\t78.90625\n", "809\t1024\t79.00390625\n", "814\t1024\t79.4921875\n", "808\t1024\t78.90625\n", "823\t1024\t80.37109375\n", "811\t1024\t79.19921875\n", "806\t1024\t78.7109375\n", "819\t1024\t79.98046875\n", "805\t1024\t78.61328125\n", "826\t1024\t80.6640625\n", "826\t1024\t80.6640625\n", "807\t1024\t78.80859375\n", "818\t1024\t79.8828125\n", "818\t1024\t79.8828125\n", "812\t1024\t79.296875\n", "816\t1024\t79.6875\n", "815\t1024\t79.58984375\n", "827\t1024\t80.76171875\n", "830\t1024\t81.0546875\n", "852\t1024\t83.203125\n", "827\t1024\t80.76171875\n", "834\t1024\t81.4453125\n", "835\t1024\t81.54296875\n", "835\t1024\t81.54296875\n", "829\t1024\t80.95703125\n", "822\t1024\t80.2734375\n", "818\t1024\t79.8828125\n", "827\t1024\t80.76171875\n", "834\t1024\t81.4453125\n", "829\t1024\t80.95703125\n", "846\t1024\t82.6171875\n", "829\t1024\t80.95703125\n", "829\t1024\t80.95703125\n", "833\t1024\t81.34765625\n", "837\t1024\t81.73828125\n", "837\t1024\t81.73828125\n", "815\t1024\t79.58984375\n", "834\t1024\t81.4453125\n", "833\t1024\t81.34765625\n", "840\t1024\t82.03125\n", "855\t1024\t83.49609375\n", "853\t1024\t83.30078125\n", "853\t1024\t83.30078125\n", "846\t1024\t82.6171875\n", "852\t1024\t83.203125\n", "856\t1024\t83.59375\n", "859\t1024\t83.88671875\n", "851\t1024\t83.10546875\n", "845\t1024\t82.51953125\n", "874\t1024\t85.3515625\n", "861\t1024\t84.08203125\n", "877\t1024\t85.64453125\n", "853\t1024\t83.30078125\n", "861\t1024\t84.08203125\n", "859\t1024\t83.88671875\n", "866\t1024\t84.5703125\n", "868\t1024\t84.765625\n", "870\t1024\t84.9609375\n", "856\t1024\t83.59375\n", "859\t1024\t83.88671875\n", "864\t1024\t84.375\n", "864\t1024\t84.375\n", "876\t1024\t85.546875\n", "872\t1024\t85.15625\n", "872\t1024\t85.15625\n", "863\t1024\t84.27734375\n", "859\t1024\t83.88671875\n", "878\t1024\t85.7421875\n", "860\t1024\t83.984375\n", "864\t1024\t84.375\n", "875\t1024\t85.44921875\n", "862\t1024\t84.1796875\n", "867\t1024\t84.66796875\n", "867\t1024\t84.66796875\n", "864\t1024\t84.375\n", "864\t1024\t84.375\n", "876\t1024\t85.546875\n", "875\t1024\t85.44921875\n", "860\t1024\t83.984375\n", "865\t1024\t84.47265625\n", "881\t1024\t86.03515625\n", "867\t1024\t84.66796875\n", "869\t1024\t84.86328125\n", "873\t1024\t85.25390625\n", "869\t1024\t84.86328125\n", "873\t1024\t85.25390625\n", "873\t1024\t85.25390625\n", "862\t1024\t84.1796875\n", "865\t1024\t84.47265625\n", "871\t1024\t85.05859375\n", "869\t1024\t84.86328125\n", "871\t1024\t85.05859375\n", "866\t1024\t84.5703125\n", "877\t1024\t85.64453125\n", "861\t1024\t84.08203125\n", "881\t1024\t86.03515625\n", "882\t1024\t86.1328125\n", "874\t1024\t85.3515625\n", "875\t1024\t85.44921875\n", "866\t1024\t84.5703125\n", "870\t1024\t84.9609375\n", "883\t1024\t86.23046875\n", "870\t1024\t84.9609375\n", "871\t1024\t85.05859375\n", "877\t1024\t85.64453125\n", "866\t1024\t84.5703125\n", "877\t1024\t85.64453125\n", "863\t1024\t84.27734375\n", "873\t1024\t85.25390625\n", "871\t1024\t85.05859375\n", "883\t1024\t86.23046875\n", "862\t1024\t84.1796875\n", "853\t1024\t83.30078125\n", "858\t1024\t83.7890625\n", "857\t1024\t83.69140625\n", "855\t1024\t83.49609375\n", "847\t1024\t82.71484375\n", "837\t1024\t81.73828125\n", "850\t1024\t83.0078125\n", "864\t1024\t84.375\n", "879\t1024\t85.83984375\n", "883\t1024\t86.23046875\n", "871\t1024\t85.05859375\n", "888\t1024\t86.71875\n", "881\t1024\t86.03515625\n", "830\t1024\t81.0546875\n", "870\t1024\t84.9609375\n", "877\t1024\t85.64453125\n", "886\t1024\t86.5234375\n", "863\t1024\t84.27734375\n", "871\t1024\t85.05859375\n", "886\t1024\t86.5234375\n", "871\t1024\t85.05859375\n", "896\t1024\t87.5\n", "872\t1024\t85.15625\n", "870\t1024\t84.9609375\n", "877\t1024\t85.64453125\n", "863\t1024\t84.27734375\n", "886\t1024\t86.5234375\n", "898\t1024\t87.6953125\n", "884\t1024\t86.328125\n", "908\t1024\t88.671875\n", "878\t1024\t85.7421875\n", "865\t1024\t84.47265625\n", "864\t1024\t84.375\n", "888\t1024\t86.71875\n", "870\t1024\t84.9609375\n", "862\t1024\t84.1796875\n", "866\t1024\t84.5703125\n", "889\t1024\t86.81640625\n", "879\t1024\t85.83984375\n", "884\t1024\t86.328125\n", "880\t1024\t85.9375\n", "876\t1024\t85.546875\n", "864\t1024\t84.375\n", "877\t1024\t85.64453125\n", "858\t1024\t83.7890625\n", "894\t1024\t87.3046875\n", "890\t1024\t86.9140625\n", "893\t1024\t87.20703125\n", "891\t1024\t87.01171875\n", "896\t1024\t87.5\n", "892\t1024\t87.109375\n", "906\t1024\t88.4765625\n", "878\t1024\t85.7421875\n", "893\t1024\t87.20703125\n", "898\t1024\t87.6953125\n", "888\t1024\t86.71875\n", "903\t1024\t88.18359375\n", "911\t1024\t88.96484375\n", "911\t1024\t88.96484375\n", "901\t1024\t87.98828125\n", "909\t1024\t88.76953125\n", "911\t1024\t88.96484375\n", "921\t1024\t89.94140625\n", "922\t1024\t90.0390625\n", "916\t1024\t89.453125\n", "923\t1024\t90.13671875\n", "928\t1024\t90.625\n", "920\t1024\t89.84375\n", "922\t1024\t90.0390625\n", "915\t1024\t89.35546875\n", "930\t1024\t90.8203125\n", "914\t1024\t89.2578125\n", "917\t1024\t89.55078125\n", "918\t1024\t89.6484375\n", "921\t1024\t89.94140625\n", "921\t1024\t89.94140625\n", "937\t1024\t91.50390625\n", "931\t1024\t90.91796875\n", "923\t1024\t90.13671875\n", "921\t1024\t89.94140625\n", "934\t1024\t91.2109375\n", "930\t1024\t90.8203125\n", "933\t1024\t91.11328125\n", "932\t1024\t91.015625\n", "930\t1024\t90.8203125\n", "930\t1024\t90.8203125\n", "933\t1024\t91.11328125\n", "933\t1024\t91.11328125\n", "949\t1024\t92.67578125\n", "941\t1024\t91.89453125\n", "945\t1024\t92.28515625\n", "936\t1024\t91.40625\n", "956\t1024\t93.359375\n", "948\t1024\t92.578125\n", "936\t1024\t91.40625\n", "941\t1024\t91.89453125\n", "949\t1024\t92.67578125\n", "941\t1024\t91.89453125\n", "940\t1024\t91.796875\n", "951\t1024\t92.87109375\n", "941\t1024\t91.89453125\n", "941\t1024\t91.89453125\n", "930\t1024\t90.8203125\n", "930\t1024\t90.8203125\n", "924\t1024\t90.234375\n", "919\t1024\t89.74609375\n", "911\t1024\t88.96484375\n", "934\t1024\t91.2109375\n", "892\t1024\t87.109375\n", "929\t1024\t90.72265625\n", "922\t1024\t90.0390625\n", "927\t1024\t90.52734375\n", "917\t1024\t89.55078125\n", "856\t1024\t83.59375\n", "835\t1024\t81.54296875\n", "852\t1024\t83.203125\n", "870\t1024\t84.9609375\n", "878\t1024\t85.7421875\n", "872\t1024\t85.15625\n", "894\t1024\t87.3046875\n", "865\t1024\t84.47265625\n", "889\t1024\t86.81640625\n", "871\t1024\t85.05859375\n", "873\t1024\t85.25390625\n", "864\t1024\t84.375\n", "859\t1024\t83.88671875\n", "867\t1024\t84.66796875\n", "833\t1024\t81.34765625\n", "853\t1024\t83.30078125\n", "874\t1024\t85.3515625\n", "843\t1024\t82.32421875\n", "848\t1024\t82.8125\n", "844\t1024\t82.421875\n", "817\t1024\t79.78515625\n", "865\t1024\t84.47265625\n", "807\t1024\t78.80859375\n", "752\t1024\t73.4375\n", "775\t1024\t75.68359375\n", "772\t1024\t75.390625\n", "778\t1024\t75.9765625\n", "767\t1024\t74.90234375\n", "784\t1024\t76.5625\n", "800\t1024\t78.125\n", "807\t1024\t78.80859375\n", "826\t1024\t80.6640625\n", "805\t1024\t78.61328125\n", "788\t1024\t76.953125\n", "820\t1024\t80.078125\n", "809\t1024\t79.00390625\n", "803\t1024\t78.41796875\n", "799\t1024\t78.02734375\n", "806\t1024\t78.7109375\n", "839\t1024\t81.93359375\n", "846\t1024\t82.6171875\n", "914\t1024\t89.2578125\n", "888\t1024\t86.71875\n", "839\t1024\t81.93359375\n", "836\t1024\t81.640625\n", "830\t1024\t81.0546875\n", "845\t1024\t82.51953125\n", "828\t1024\t80.859375\n", "834\t1024\t81.4453125\n", "854\t1024\t83.3984375\n", "847\t1024\t82.71484375\n", "846\t1024\t82.6171875\n", "845\t1024\t82.51953125\n", "863\t1024\t84.27734375\n", "867\t1024\t84.66796875\n", "855\t1024\t83.49609375\n", "844\t1024\t82.421875\n", "864\t1024\t84.375\n", "865\t1024\t84.47265625\n", "860\t1024\t83.984375\n", "868\t1024\t84.765625\n", "871\t1024\t85.05859375\n", "868\t1024\t84.765625\n", "857\t1024\t83.69140625\n", "885\t1024\t86.42578125\n", "908\t1024\t88.671875\n", "872\t1024\t85.15625\n", "848\t1024\t82.8125\n", "813\t1024\t79.39453125\n", "763\t1024\t74.51171875\n", "761\t1024\t74.31640625\n", "779\t1024\t76.07421875\n", "783\t1024\t76.46484375\n", "776\t1024\t75.78125\n", "784\t1024\t76.5625\n", "811\t1024\t79.19921875\n", "812\t1024\t79.296875\n", "777\t1024\t75.87890625\n", "794\t1024\t77.5390625\n", "794\t1024\t77.5390625\n", "813\t1024\t79.39453125\n", "805\t1024\t78.61328125\n", "828\t1024\t80.859375\n", "796\t1024\t77.734375\n", "806\t1024\t78.7109375\n", "817\t1024\t79.78515625\n", "840\t1024\t82.03125\n", "786\t1024\t76.7578125\n", "818\t1024\t79.8828125\n", "832\t1024\t81.25\n", "835\t1024\t81.54296875\n", "768\t1024\t75.0\n", "841\t1024\t82.12890625\n", "833\t1024\t81.34765625\n", "836\t1024\t81.640625\n", "824\t1024\t80.46875\n", "830\t1024\t81.0546875\n", "837\t1024\t81.73828125\n", "837\t1024\t81.73828125\n", "858\t1024\t83.7890625\n", "847\t1024\t82.71484375\n", "870\t1024\t84.9609375\n", "866\t1024\t84.5703125\n", "843\t1024\t82.32421875\n", "867\t1024\t84.66796875\n", "849\t1024\t82.91015625\n", "869\t1024\t84.86328125\n", "860\t1024\t83.984375\n", "862\t1024\t84.1796875\n", "846\t1024\t82.6171875\n", "854\t1024\t83.3984375\n", "871\t1024\t85.05859375\n", "861\t1024\t84.08203125\n", "882\t1024\t86.1328125\n", "892\t1024\t87.109375\n", "877\t1024\t85.64453125\n", "895\t1024\t87.40234375\n", "887\t1024\t86.62109375\n", "873\t1024\t85.25390625\n", "894\t1024\t87.3046875\n", "890\t1024\t86.9140625\n", "878\t1024\t85.7421875\n", "894\t1024\t87.3046875\n", "879\t1024\t85.83984375\n", "890\t1024\t86.9140625\n", "895\t1024\t87.40234375\n", "889\t1024\t86.81640625\n", "896\t1024\t87.5\n", "898\t1024\t87.6953125\n", "901\t1024\t87.98828125\n", "879\t1024\t85.83984375\n", "890\t1024\t86.9140625\n", "888\t1024\t86.71875\n", "917\t1024\t89.55078125\n", "902\t1024\t88.0859375\n", "921\t1024\t89.94140625\n", "915\t1024\t89.35546875\n", "927\t1024\t90.52734375\n", "923\t1024\t90.13671875\n", "928\t1024\t90.625\n", "923\t1024\t90.13671875\n", "914\t1024\t89.2578125\n", "918\t1024\t89.6484375\n", "927\t1024\t90.52734375\n", "926\t1024\t90.4296875\n", "919\t1024\t89.74609375\n", "916\t1024\t89.453125\n", "928\t1024\t90.625\n", "916\t1024\t89.453125\n", "933\t1024\t91.11328125\n", "925\t1024\t90.33203125\n", "930\t1024\t90.8203125\n", "930\t1024\t90.8203125\n", "934\t1024\t91.2109375\n", "933\t1024\t91.11328125\n", "935\t1024\t91.30859375\n", "939\t1024\t91.69921875\n", "934\t1024\t91.2109375\n", "938\t1024\t91.6015625\n", "944\t1024\t92.1875\n", "937\t1024\t91.50390625\n", "937\t1024\t91.50390625\n", "935\t1024\t91.30859375\n", "937\t1024\t91.50390625\n", "937\t1024\t91.50390625\n", "954\t1024\t93.1640625\n", "940\t1024\t91.796875\n", "942\t1024\t91.9921875\n", "955\t1024\t93.26171875\n", "949\t1024\t92.67578125\n", "941\t1024\t91.89453125\n", "947\t1024\t92.48046875\n", "940\t1024\t91.796875\n", "943\t1024\t92.08984375\n", "946\t1024\t92.3828125\n", "962\t1024\t93.9453125\n", "954\t1024\t93.1640625\n", "956\t1024\t93.359375\n", "957\t1024\t93.45703125\n", "962\t1024\t93.9453125\n", "960\t1024\t93.75\n", "944\t1024\t92.1875\n", "969\t1024\t94.62890625\n", "969\t1024\t94.62890625\n", "969\t1024\t94.62890625\n", "968\t1024\t94.53125\n", "970\t1024\t94.7265625\n", "969\t1024\t94.62890625\n", "975\t1024\t95.21484375\n", "963\t1024\t94.04296875\n", "965\t1024\t94.23828125\n", "975\t1024\t95.21484375\n", "969\t1024\t94.62890625\n", "966\t1024\t94.3359375\n", "969\t1024\t94.62890625\n", "984\t1024\t96.09375\n", "981\t1024\t95.80078125\n", "977\t1024\t95.41015625\n", "982\t1024\t95.8984375\n", "980\t1024\t95.703125\n", "980\t1024\t95.703125\n", "981\t1024\t95.80078125\n", "982\t1024\t95.8984375\n", "981\t1024\t95.80078125\n", "987\t1024\t96.38671875\n", "982\t1024\t95.8984375\n", "982\t1024\t95.8984375\n", "977\t1024\t95.41015625\n", "982\t1024\t95.8984375\n", "986\t1024\t96.2890625\n", "984\t1024\t96.09375\n", "986\t1024\t96.2890625\n", "987\t1024\t96.38671875\n", "996\t1024\t97.265625\n", "995\t1024\t97.16796875\n", "997\t1024\t97.36328125\n", "999\t1024\t97.55859375\n", "1002\t1024\t97.8515625\n", "999\t1024\t97.55859375\n", "1004\t1024\t98.046875\n", "1002\t1024\t97.8515625\n", "999\t1024\t97.55859375\n", "1000\t1024\t97.65625\n", "1007\t1024\t98.33984375\n", "1010\t1024\t98.6328125\n", "1012\t1024\t98.828125\n", "1004\t1024\t98.046875\n", "1000\t1024\t97.65625\n", "1008\t1024\t98.4375\n", "1005\t1024\t98.14453125\n", "1011\t1024\t98.73046875\n", "1011\t1024\t98.73046875\n", "1009\t1024\t98.53515625\n", "1005\t1024\t98.14453125\n", "1007\t1024\t98.33984375\n", "1010\t1024\t98.6328125\n", "1010\t1024\t98.6328125\n", "1005\t1024\t98.14453125\n", "1008\t1024\t98.4375\n", "1012\t1024\t98.828125\n", "1009\t1024\t98.53515625\n", "1012\t1024\t98.828125\n", "1013\t1024\t98.92578125\n", "1010\t1024\t98.6328125\n", "1012\t1024\t98.828125\n", "1014\t1024\t99.0234375\n", "1016\t1024\t99.21875\n", "1010\t1024\t98.6328125\n", "1014\t1024\t99.0234375\n", "1015\t1024\t99.12109375\n", "1012\t1024\t98.828125\n", "1019\t1024\t99.51171875\n", "1015\t1024\t99.12109375\n", "1016\t1024\t99.21875\n", "1019\t1024\t99.51171875\n", "1019\t1024\t99.51171875\n", "1016\t1024\t99.21875\n", "1018\t1024\t99.4140625\n", "1018\t1024\t99.4140625\n" ] } ], "source": [ "spectro = make_dct_spectrogram(wave, seg_length=1024)\n", "for t, dct in sorted(spectro.spec_map.items()):\n", " compress(dct, thresh=0.2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In most segments, the compression is 75-85%.\n", "\n", "To hear what it sounds like, we can convert the spectrogram back to a wave and play it." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wave2 = spectro.make_wave()\n", "wave2.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's the original again for comparison." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wave.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an experiment, you might try increasing `thresh` to see when the effect of compression becomes audible (to you).\n", "\n", "Also, you might try compressing a signal with some noisy elements, like cymbals." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "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 }