music-synthesizer-for-android/lab/Second order sections in matrix form.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from scipy import signal"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook contains a detailed description of implementing classic second order filter sections using matrix arithmetic.\n",
      "The matrix formulation lends itself to robust, fast evaluation on SIMD architectures. Second order filter sections are\n",
      "extremely versatile, as a single section can implement many different types of filter, including lowpass, highpass,\n",
      "bandpass, notch, bell, low shelf, and high shelf. Being second order, rolloff is 12dB/octave. For sharper frequency\n",
      "response, higher order IIR filters can be decomposed into multiple second order sections.\n",
      "\n",
      "All second order filters can be described using two state variables. The evolution of the state is a 2x2 matrix. In addition,\n",
      "there is a 2-vector describing the contribution of the input to the state, a 2-vector describing \"reading out\" the state to\n",
      "the output, and a scalar for the zero-delay contribution of the input to the output. For mathematical convenience, we group\n",
      "the last two into a single 3-vector. Thus, a filter is described as the tuple $(\\mathbf{A}, B, C)$, and each step is as follows:\n",
      "\n",
      "\\begin{align}\n",
      "    out_n & = C \\cdot concat([[x_n], y_n]) \\\\\n",
      "    y_{n+1} & = B \\cdot x_n + \\mathbf{A} \\cdot y_n\n",
      "\\end{align}\n",
      "\n",
      "Where the $\\cdot$ notation represents dot product, vector times scalar, or matrix times vector multiplication, respectively.\n",
      "\n",
      "Question: would it be more convenient to have a single 3x3 matrix?\n",
      "\n",
      "There is some redundancy in this representation. Obviously, multiplying A and B times some scalar, and dividing C by the same,\n",
      "is a filter with the same results. More subtly, the y vector can be transformed by any non-degenerate linear transformation.\n",
      "However, these transforms do matter when the filter parameters are modulated, which is especially common in synthesizers, and\n",
      "also matter for roundoff behavior (particularly important when evaluating the filters using 32 bit floating point arithmetic).\n",
      "\n",
      "The reference code for evaluating a filter in this form is straightforward:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def eval(filt, inp):\n",
      "    A, B, C = filt\n",
      "    result = []\n",
      "    y = zeros(2)\n",
      "    for x in inp:\n",
      "        result.append(dot(C, concatenate([[x], y])))\n",
      "        y = B * x + dot(A, y)\n",
      "    return array(result)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now how do we get a filter in this form? Starting from standard biquad coefficients is, again, pretty easy. Here, the\n",
      "y vector represents the two state variables in transposed direct form II. For review, the transfer function of this\n",
      "filter is:\n",
      "\n",
      "\\begin{equation} H(z) = \\frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1  + a_1 z^{-1} + a_2 z^{-2}} \\end{equation}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def from_biquad(coefs):\n",
      "    b0, b1, b2, a1, a2 = coefs\n",
      "    A = array([[-a1, 1], [-a2, 0]])\n",
      "    B = array([b1 - a1 * b0, b2 - a2 * b0])\n",
      "    C = array([b0, 1, 0])\n",
      "    return A, B, C"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can test this against a straightforward implementation (in this case, direct form I)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def biquad(coefs, inp):\n",
      "    b0, b1, b2, a1, a2 = coefs\n",
      "    out = zeros(len(inp))\n",
      "    x1 = 0\n",
      "    x2 = 0\n",
      "    y1 = 0\n",
      "    y2 = 0\n",
      "    for i in range(len(inp)):\n",
      "        x = inp[i]\n",
      "        y = b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2\n",
      "        out[i] = y\n",
      "        x2, x1 = x1, x\n",
      "        y2, y1 = y1, y\n",
      "    return out\n",
      "\n",
      "impulse = zeros(100)\n",
      "impulse[10] = 1\n",
      "coefs = 1.0207, -1.7719, .9376, -1.7719, 0.9583  # peaking\n",
      "plot(biquad(coefs, impulse))\n",
      "plot(eval(from_biquad(coefs), impulse))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHXCP/DPmZlzuGmIN1RAUUBAUdRFrVxd0szLo1Tm\n02qvp1w1lywr23az9tJK+5Raz65r2fN6bE2tLNNfN9vSycjQTPGSF0rM0CS5qIWIijBzzpzz/f2B\nogh4m4Gx+X7erxevFzPz5ZwvZ2M+8/l+Z1ZFCCFARETSsfl7AkRE5B8MACIiSTEAiIgkxQAgIpIU\nA4CISFIMACIiSXkdAFOmTEFkZCR69erV4ONvvvkmUlNT0bt3bwwaNAh5eXnenpKIiHzA6wCYPHky\nnE5no49369YNGzduRF5eHv7yl7/gt7/9rbenJCIiH/A6AAYPHoyIiIhGH7/pppsQHh4OABg4cCCK\ni4u9PSUREflAs+4BvPrqqxg9enRznpKIiBrhaK4Tff7551iyZAm+/PLL5jolERFdQrMEQF5eHqZN\nmwan09ngclF8fDwOHjzYHFMhIgoYcXFxOHDgwDX/fJMvAR0+fBjjxo3D8uXLER8f3+CYgwcPQgjB\nLyHw17/+1e9zuF6+eC14LXgtLv3l7QtnrxvAxIkTsWHDBpSVlSEmJgZZWVkwDAMAkJmZiWeeeQYn\nTpzA9OnTAQCqqmLbtm3enpaIiLzkdQCsWLHiko8vXrwYixcv9vY0RETkY/wk8HUmPT3d31O4bvBa\nnMdrcR6vhe8oQgi//4MwiqLgOpgGEdHPirfPnWwARESSYgAQEUmKAUBEJCkGABGRpBgARESSYgAQ\nEUmKAUBEJCkGABGRpBgARESSYgAQEUmKAUBEJCkGABGRpBgARESSYgAQEUmKAUBEJCkGABGRpBgA\nRESSYgAQEUmKAUBEJCkGABGRpBgARESS8ioApkyZgsjISPTq1avRMY888ggSEhKQmpqKXbt2eXO6\nJrWzoBQ/HKvw9zSIiJqNVwEwefJkOJ3ORh9fs2YNDhw4gIKCArzyyiuYPn26N6drUpNfnYvfvfa6\nv6dBRNRsvAqAwYMHIyIiotHHP/zwQ0yaNAkAMHDgQFRUVODYsWPenLLJ6JYLLo/b39MgImo2TboH\nUFJSgpiYmNrb0dHRKC4ubspTXjOPZcAwDX9Pg4io2Tia+gRCiDq3FUVpcNzs2bNrv09PT0d6enoT\nzqo+jzCgm3qznpOI6Grk5OQgJyfHZ8dr0gCIiopCUVFR7e3i4mJERUU1OPbCAPAHj9BhWGwARHT9\nuvjFcVZWllfHa9IloIyMDLz+es3Gam5uLlq1aoXIyMimPOU1M9kAiEgyXjWAiRMnYsOGDSgrK0NM\nTAyysrJgGDWvojMzMzF69GisWbMG8fHxCAsLw9KlS30y6aZgwoCHDYCIJKKIixfp/TEJRam3V9Dc\n2s4cgQ5Bcfhm3v/6dR5ERFfK2+dOfhL4LBM6GwARSYUBcJYJA4bgHgARyYMBcJalGDAFGwARyYMB\ncJYFAx42ACKSCAPgLEvR2QCISCoMgLMsxYAHbABEJA8GwFmCewBEJBkGwFmWosNkAyAiiTAAzhI2\nA5bCBkBE8mAAnGMzYIEBQETyYACcVdMAuARERPJgAJxj09kAiEgqDIBz7AYsGxsAEcmDAQBAN0xA\nERDcBCYiiTAAAFS5a574BRsAEUmEAQCgsrrmiZ8NgIhkwgAAUO02ANMBYWcDICJ5MAAAVOsGYIQB\nNjYAIpIHAwA1DcBmhgJsAEQkEQYAgDNuHTYrGMDZdwQREUmAAQDApRuwCQ0wtdoNYSKiQMcAQM0e\ngCJUwFJr3xJKRBToGAA41wBUKJaKMy42ACKSg9cB4HQ6kZSUhISEBMybN6/e42VlZRg5ciT69OmD\nlJQULFu2zNtT+ly1rsMGFYql1bwjiIhIAl4FgGmamDFjBpxOJ/Lz87FixQrs27evzpiFCxeib9++\n2L17N3JycvD444/D4/F4NWlfcxkG7EKDYqk1nwkgIpKAVwGwbds2xMfHIzY2FqqqYsKECVi9enWd\nMR07dsSpU6cAAKdOnUKbNm3gcDi8Oa3PuQ2jpgEIDVVuLgERkRy8eiYuKSlBTExM7e3o6Ghs3bq1\nzphp06Zh6NCh6NSpE06fPo1Vq1Z5c8omUa3rsEOFTahcAiIiaXgVAIqiXHbMc889hz59+iAnJwcH\nDx7E8OHDsWfPHrRs2bLOuNmzZ9d+n56ejvT0dG+mdlXchgG7osLGBkBE17GcnBzk5OT47HheBUBU\nVBSKiopqbxcVFSE6OrrOmM2bN+NPf/oTACAuLg5du3bF/v37kZaWVmfchQHQ3NweA3ZoUIQKFxsA\nEV2nLn5xnJWV5dXxvNoDSEtLQ0FBAQoLC6HrOlauXImMjIw6Y5KSkpCdnQ0AOHbsGPbv349u3bp5\nc1qfc3sMOBQVdmio1tkAiEgOXjUAh8OBhQsXYsSIETBNE1OnTkVycjIWLVoEAMjMzMQf//hHTJ48\nGampqbAsC88//zxat27tk8n7itvQa5aALBUugw2AiOSgCCGE3yehKPDnNCa/uASbDn+B454iPDFo\nFp78z+F+mwsR0ZXy9rmTnwQGoJsGHDYNdkWFmw2AiCTBAACgewyoNhUOaHAZ3AMgIjkwAADopl6z\nCayocHvYAIhIDgwA1CwBqXYNDoUNgIjkwQAAYJg1S0B2RYXOBkBEkmAAADAsA6pdhUNR4fKwARCR\nHBgAqNkDUG0qVJvGBkBE0mAAoKYBaA4NDpsKw2QAEJEcGAAAPJYBzV7TANxcAiIiSTAAABiWfjYA\nVBgWGwARyYEBgPMNQLNr0E02ACKSAwMAgEcYCFI1qHbuARCRPBgAAExxQQOw2ACISA4MAAAeoSNI\nVaHZVXi4B0BEkmAAoKYBBDlUBDk0GGwARCQJBgAAEwaCVY0NgIikwgDA2QagqghSNXgEGwARyYEB\nAMCEjuBzewCCDYCI5MAAQM0SUJCqIpgNgIgkwgAAYCkGQjUNQaoKkw2AiCTBAABgwUCwpiLIocID\nNgAikgMDAICl1OwBhGgaLDYAIpIEAwCAUAyEBp1dAgIDgIjk4HUAOJ1OJCUlISEhAfPmzWtwTE5O\nDvr27YuUlBSkp6d7e0qfsxQDIUE1DcBUuARERHJwePPDpmlixowZyM7ORlRUFPr374+MjAwkJyfX\njqmoqMBDDz2ETz75BNHR0SgrK/N60r4mbDV7AMGqCosNgIgk4VUD2LZtG+Lj4xEbGwtVVTFhwgSs\nXr26zpi33noLd911F6KjowEAbdu29eaUTULYdIRoKkKDNFhsAEQkCa8CoKSkBDExMbW3o6OjUVJS\nUmdMQUEBysvLccsttyAtLQ1vvPGGN6dsEkIx0CJEQ7CmQihsAEQkB6+WgBRFuewYwzCwc+dOfPbZ\nZ6iqqsJNN92EG2+8EQkJCXXGzZ49u/b79PT05t0rsBlsAER03cvJyUFOTo7PjudVAERFRaGoqKj2\ndlFRUe1SzzkxMTFo27YtQkJCEBISgiFDhmDPnj2XDIBmZ9cRGqwihA2AiK5jF784zsrK8up4Xi0B\npaWloaCgAIWFhdB1HStXrkRGRkadMbfffjs2bdoE0zRRVVWFrVu3okePHl5N2pcsSwB2D0KDVIQF\na7BsbABEJAevGoDD4cDChQsxYsQImKaJqVOnIjk5GYsWLQIAZGZmIikpCSNHjkTv3r1hs9kwbdq0\n6yoAXLoHMB2w2RSEaCpgYwMgIjkoQgjh90koCvw1jbKTVWj3QluI/65C6fHTiPpHR4hnK/0yFyKi\nq+Htc6f0nwQ+49IBSwUAhAaxARCRPKQPgCq3AeXCAHDoNfsCREQBTvoAqHYbUCwNAKCpdsCyQfeY\nfp4VEVHTYwDoBhShnr/DUlHl4jIQEQU+6QOgyqXXLgEBAEwNlS6+FZSIAp/0AVBtGLBd0AAUS0W1\nmw2AiAKf9AHg0g3YhFZ7W7G0
       "text": [
        "<matplotlib.figure.Figure at 0x109465510>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "State Variable Filter design"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "However, biquads are not necessarily your best bet for designing second order filters, for two reasons. First, they don't\n",
      "modulate well, because the state space for one set of filter parameters probably doesn't make sense with a different filter.\n",
      "Second, roundoff behavior can be bad. We'll analyze that in more detail later.\n",
      "\n",
      "A much better design paradigm is to discretize the State Variable Filter. Historically, many people used the naive\n",
      "discretization, originally popularized by Hal Chamberlin's book, _Musical Applications of Microprocessors_ (see\n",
      "<https://ccrma.stanford.edu/~jos/svf/svf.pdf>), but this only works well at low frequencies and becomes\n",
      "unstable for frequencies around one sixth the sampling rate, or even worse for higher resonance values. Fortunately,\n",
      "Andrew Simper of Cytomic has done a thorough development of discretization using trapezoidal integration. This is equivalent\n",
      "to the bilinear transform, which is not exactly trivial because the trapezoidal integration (unlike the naive discretization)\n",
      "contains a delay-free loop. But fortunately Simper has done all the work for us:\n",
      "\n",
      "[Solving the continuous SVF equations using trapezoidal integration and equivalent currents](http://www.cytomic.com/files/dsp/SvfLinearTrapOptimised2.pdf)\n",
      "\n",
      "For simplicity, we'll start with just a lowpass filter. As we'll see later, the full derivation is very general and can\n",
      "do the entire RBJ suite of filter designs."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def svf_lp(f, res):\n",
      "    g = tan(pi * f)\n",
      "    k = 2 - 2 * res # 1/Q\n",
      "    a1 = 1/(1 + g * (g + k))\n",
      "    a2 = g * a1\n",
      "    a3 = g * a2\n",
      "    A = array([[2*a1 - 1, -2*a2],\n",
      "               [2*a2, 1 - 2*a3]])\n",
      "    B = array([2 * a2, 2 * a3])\n",
      "    C = array([a3, a2, 1 - a3])\n",
      "    return A, B, C"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Simper prefers a resonance control that ranges from 0 to 1. If you prefer the classic Q, just use this conversion:\n",
      "\n",
      "\\begin{equation} res = 1 - .5/Q \\end{equation}\n",
      "\n",
      "Let's test the filter out:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def db(x):\n",
      "    return 20 * log(x)/log(10)\n",
      "def plot_fr(filter):\n",
      "    impulse = concatenate([[1], zeros(5000)])\n",
      "    w, h = signal.freqz(eval(filter, impulse))\n",
      "    semilogx(w, db(abs(h)))\n",
      "axis([3e-3, pi, -60, 10])\n",
      "for f in [0.00316, 0.01, 0.0316, 0.1, 0.316, 0.49]:\n",
      "    plot_fr(svf_lp(f, 0.293))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdYFFcXB+DfbKN3BBVUEFBE0BhrEgsWVFDRoKJorLFE\nY4waEzWJJVGDJeYzsSRGY4sFYiFYEmILphrsDeyoKIIKio26O98fV7EAC7s7y+yO532efYDdKQdG\nz969c++5HM/zPAghhJg9mdgBEEIIEQYldEIIkQhK6IQQIhGU0AkhRCIooRNCiERQQieEEIkwOKEP\nHToU7u7uCAoKKn4uOzsbISEhqFOnDjp27Ii7d+8aehpCCCHlMDihDxkyBAkJCc89N2fOHISEhODc\nuXNo37495syZY+hpCCGElIMTYmLR5cuX0a1bN5w8eRIA4O/vj/3798Pd3R0ZGRkIDg7GmTNnDA6W\nEEJI2YzSh56ZmQl3d3cAgLu7OzIzM41xGkIIIc8w+k1RjuPAcZyxT0MIIS89hTEO+qSrpWrVqrhx\n4wbc3NxKbOPr64uLFy8a4/SEECJZDRs2xLFjx0p9zSgt9PDwcKxZswYAsGbNGvTo0aPENhcvXgTP\n8zo/pk+fLui22rbR9bWKPifGQ8g49D0WXTu6doZcO12el/K1O378eJm51+CEHhUVhddffx1nz55F\njRo1sGrVKkyePBm7d+9GnTp1sG/fPkyePNnQ0xQLDg4WdFtt2+j6mi6xVTYhY9P3WHTt9EPXTr/n\nTUFlXztBRrnog+M4iHTqSjVjxgzMmDFD7DCIHujamS8pXzttuZNmihqZKbceiHZ07czXy3rtqIVO\nCCFmhFrohBDyEqCETgghEkEJnRBCJIISOiGESAQldEIIkQhK6IQQIhGU0AkhRCIooRNCiERQQieE\nEImghE4IIRJBCZ0QQiSCEjohhEgEJXRCCJEISuiEECIRlNAJIUQiKKETQohEUEInhBCJoIROCCES\nQQmdEEIkghI6IYRIBCV0QgiRCKMl9ISEBPj7+8PPzw9z58411mkIIYQ8xvE8zwt9ULVajbp162LP\nnj3w8PBA06ZNsXHjRtSrV+/piTkORjg1IYRImrbcaZQWelJSEnx9feHl5QWlUom+ffsiPj7eGKci\nhBDymFES+vXr11GjRo3inz09PXH9+nVjnIoQQshjCmMclOM4YxzWKNbN6I5q/ge1bsO/+OuUeKIc\nOvcs6Xp83bZ/MZwSv96zz5dxbL6U/djz+sRe2h+olOc5Toe/ZXlxlPI6h1J/X+1/L67UbfSP69kT\nc89v/cKuJY7EcwD39PmSX5++yD35WRf8Mzvrsn0pSn9F17zBPflFwMnA/n08+fmZ79nrHDgZB8g5\ncHL2PScH+1nBQaZkrz13dJ3zmC7b63ZsKytf1Kr1Cayt62jdzigJ3cPDA2lpacU/p6WlwdPTs8R2\nM2bMKP4+ODgYwcHBxghHK7dXJiD1YhoKC4C8fKCgACjI51FQAOQVAIV5QH4BexTm8ygs4JFfwF7n\nAKhUPFQWgEr1+GEBWCgBlSWgUgIWT15T8lCqAAvV022Vysf7qwCZnH/aL8bz4MF+5nn+cbLgAfDg\n+cfJg3/y8+NtAXA8oAHA8xqA56Hh1VDzPHhooNFooNEUQQNAw6uh0fDQ8JrHDzV4AEWaIqjVRSjU\nFKJQXYRCdQEK+UIUqQtRqClEkboIBY+/5qvzkVeUB0u5BaxVNrBWWsFGaQNrpQ1sFNawUlnDXmUH\nJytn2KlsIXvxHzCPx1GXngpLdhHygIYHz6vZLho1oNEAGp595TXgNTygVhdvC42a/axWA4VFQGEh\nUFTIvi8qAl/4+Hv149cK8oG8PCAvD3xuHtvf0gqwtgJs7QC7xw97e/bV2RFwdgGsrIqvUOm/C//0\nd3jhVyr1ewAa8NDwgJrnoQEPNf/kwa5xEa+Bhn+8zePXNY+3KeR55Gs0KHj8tfh7tQaFPI9cjRoP\n1exhKZPBVi6HjVwOO7kcTnIlnJUKOCsUcFYq4aJUoIpK9fSjfAXeubRd17KOwWv725R6DJ5d7yIe\nmiIe/OMHigBNoebx98+8VqiBJlcD9SM11I800DxSQ/OQ/Vx0T43CrELwGh5KFyUULkpY1FDBspYF\nrLysYOljBVU1VTkJXpdWm24tvD//PIw9e+Lw8OHrqF79Xa3bGuWmaFFREerWrYu9e/eievXqaNas\nmeRuivI8kJ8P3L8P3LvHvj77vS5f791jOaJzZ6BLF/bV2Vns37B8ao0aWblZyHyQiZsPbyLzYeZz\n31/NuYqU2ynIyctBXde6qOdajz2q1ENAlQDUcakDGWfCI2fz84G7d4GsLODGDfZIT2dfr18HLl4E\nzp9n78516rDHq68CjRsDr7wC2NiI/RtopeF55BQVIauwEFlFRbhVUIDrBQVIy8vD1fx8pOXn43Je\nHjIKCuBnZYUAa2vUs7FBY1tbtLC3h6tKJfavIKiiB0UoyChAQXoBcs/n4mHKQzw68wgPjj4AeMCh\npQOcOjjB9U1XqKpU7u/O8zz++88HQUHbYWsbWGbuNEpCB4Bff/0V48aNg1qtxttvv40pU6Y8f2Iz\nT+hC4nmWH375BdixA0hMZPmgSxega1cgIIB9mjRXOXk5OHP7DFJupyDlVgpSbqfg1M1TuF9wH229\n2qKtV1u0826HOi51zKq7DgC7eJmZLLGfOQMcOQIcOgScPg3Urg20bAmEhABt25rHu3QpHqrVOPvo\nEVIePcKphw9x6P59JN27BzeVCi3s7dHe0RGdnZ1R1cJC7FCNgud55F3OQ86fOcj6JQvZv2bDvoU9\nPN7zgEuYC+vKqQSnTkXAzS0K7u6RlZ/Qy0MJvWy5uSyp79jBHjLZ0+QeHAxYWoodoTDSctLw++Xf\nsS91H/al7oOaV6OtV1t08euC7v7dYa20FjtE/RUUACdPAvv3A3v2AH/9Bfj7s4vYpw9Qt67YERpE\nzfNIefgQ/9y7h9137mB3djZ8razQxcUFUW5u8DfxTyeGUD9S43bcbaR9lQbNIw1qz6sNl64uRm+M\nXLz4IZTKKqhVaxIldHPF86yx9yS5nzzJGntdurBH9epiRygMnudx6c4l7E3di60pW5F0PQkR9SIw\nqOEgtKzZ0vxa7i/Kzwf+/ReIiwM2bQLc3IDISGDwYElcxEKNBv/cu4f427cRc/MmPCwsMMDdHf3d\n3eGiVIodnlHwPI/s37JxYdwFWPlYoe4PdWFR1XifUq5enYfCwtvw9Z1PCV0qsrKAhASW3H/7DfD2\nftp6b9KEteal4Pq961h/cj3WHF+D3MJcDGw4ECMbj0Q1u2pih2Y4tRr4+29gwwYgNhbo0AEYPZp9\n/DL3Ny6w1vveO3ewNiMDO7Oz0adKFYzz9JRsq11ToMGVWVdwY8UN1NtQD07BTkY5z40bPyAn52/U\nq7eKEroUFRUB//zDkvvOnSzZh4ay5N61KxthY+54nseRG0ew8uhKbDy1Ef2D+mNSy0nwtC85asos\n3bsHrFsHLFkCKBTA1KlARIRk3pkz8vPxbXo6vktPx2v29pjp7Y0gW1uxwzKKO3vvIDkqGX5L/eDW\ny03w42dmxuD27Z8RGBhLCf1lcOkSS+xbtgC3bgErVwLNm4sdlXAyHmRgwT8L8MPRHxBZPxKTW06G\nl6OX2GEJg+fZxZs5E3jwAJg2jXXJSKDFDgC5ajWWpadjztWr6ODkhM+9vVH78XBPKbl/7D5Ohp2E\n32I/VImoIuixb9/ehhs3VqBBg+1l505eJCKeWvI0Gp6PieH5qlV5fvx4nn/4UOyIhHXzwU1+yp4p\nvPNcZ/79X9/n7+beFTsk4Wg0PP/bbzzfuDHPt2jB8//8I3ZEgrpXWMh/lprKu/z5Jz8jNZXPLSoS\nOyTB3Tt8j/+ryl/8nT/vCHrcrKxd/LFjHbTmTml8riPP4Tg2kOLkSeDmTSAoCPj9d7GjEk4Vmyr4\nov0XODvmLB4UPEC9JfWw/sR6
       "text": [
        "<matplotlib.figure.Figure at 0x1094af490>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the classic BLT lowpass filter response, with complete attenuation at Nyquist and the characteristic frequency warping."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "axis([3e-3, pi, -60, 10])\n",
      "for f in [0.00316, 0.01, 0.0316, 0.1, 0.316, 0.49]:\n",
      "    plot_fr(svf_lp(f, 0.8))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdYVMfXx793G0uvUmSVIihiwRKxJCpqiL0raOwtUaOx\nxRKNkdiwRBMTNeZNbMnPRCxR1ERjJbbYu1gRBUQQpUmH3Xn/GCUWypa7e3fX+TzPPiy3zHzJmrPn\nnjlzDkcIIWAwGAyGySMSWgCDwWAw+IEZdAaDwTATmEFnMBgMM4EZdAaDwTATmEFnMBgMM4EZdAaD\nwTATdDbow4cPh5ubG+rVq1d6LD09HaGhoahZsyY++OADZGZm6joNg8FgMCpBZ4M+bNgw7Nu375Vj\nixYtQmhoKG7fvo127dph0aJFuk7DYDAYjErg+NhYdP/+fXTt2hVXr14FAAQEBOCff/6Bm5sbUlJS\nEBISgps3b+oslsFgMBjlo5cYempqKtzc3AAAbm5uSE1N1cc0DAaDwXgJvS+KchwHjuP0PQ2DwWC8\n9Uj0MeiLUIu7uzsePXoEV1fXN67x8/NDXFycPqZnMBgMsyUoKAiXLl0q85xePPRu3bph48aNAICN\nGzeiR48eb1wTFxcHQojGrzlz5vB6bUXXaHpO3WNCvPjUoe1Y7LNjn50un50mx835s7t8+XK5tldn\ng96/f3+0aNECt27dQrVq1bB+/XrMmDEDBw4cQM2aNXH48GHMmDFD12lKCQkJ4fXaiq7R9Jwm2gwN\nn9q0HYt9dtrBPjvtjhsDhv7seMly0QaO4yDQ1AYlIiICERERQstgaAH77EwXc/7sKrKdbKeonjFm\n74FRMeyzM13e1s+OeehmBiFAdDTwzjuAQiG0Gv1QWFKI/XH7sfv2bnxY70OEeIcILYkf8vOBDRuA\nI0eA1FSgfn1g9GigTh2hlfFCXH4+1iQn43R2NuQiEbq7uGCUhwdkIvPxK5UFStybcQ95sXmosbQG\nbIJseBu7qCgVWVnH4erah3nobwPPngGDBgFjxgCdOwM5OUIr4o8SVQn2x+3H8OjhqLq8KpaeXAqF\nnQJhW8Ow5foWoeXpzrZtQM2awN9/A127ArNnA66uQEgIsHix0Op0olilwvz799HswgWIAUR4e2Oc\npyd2P3mCdy9exOOiIqEl8gIhBLH9YlGYWAjnbs642v0qVIUq3sZPTFyG69f7VHiNXtIWGYbn/Hmg\nXz+gbVsgLg749FNq3LdvB0zdAbqRdgNdf+8KZytn9KvTD3PbzIXCjj5+9AjogU6bOiElJwWfNv1U\nYKVaUFgITJgAHDwIREUBLVr8d+7994Fhw4D27emjF4/JBYYis7gYfa5fh4jjcKFxY1STy0vPdXV2\nxpf37yP08mWcbNQI1mKxgEp159HPj1CYVIhG/zaCSCpC+r50pGxMQdWPqvIyPsfJKr+ICISAU5sV\nKhUh33xDSJUqhERF/Xe8sJCQli0J+fxz4bTxwaF7h4jrUley8dLGcq+5n3Gf1Pq+Fpl+YDpRqVQG\nVKcj6emEtGpFSM+ehGRllX/dw4eE+PgQsm2b4bTxwP38fFLn9Gny6e3bpKScz0WlUpGBsbFk1M2b\nBlbHL/kJ+eS4y3Hy7Oqz0mNpu9PIhfcu8DbH3bvTyJEjqNB2MoNuwqSlEdKlCyHBwYTcu/fm+ceP\nqR349VfDa+OD9RfXE9elruRI/JFKr32S+4Q0+7kZGbVrlP6F8UF8PCG1axMyaRIhSmXl1586RYir\nKzXuJsDF7GzieeIE+TYxsdJrs4uLid+pUyQ6Lc0AyvhHpVKRy50vk/i58a8cVxYqyTGnYyQ/MZ+X\neW7fnsAMurkSE0OIQkHI1KnUGy+Pq1ep9/7vv4bTpisqlYp8cegL4rvCl9xIu6H2fblFucR3hS85\nGHdQj+p44Nw5QqpWJeTbbzW7LyKCkG7d9KOJR2IyMkiV48fJ1tRUte/Z//Qp8f33X5JfUqJHZfrh\n8bbH5HTt00RZ+OYX87XwayR5XTIv89y6NYYZdHOjpISQOXMIcXcnZO9e9e7ZvZvajwcP9CqNF/KL\n80n/bf1Js5+bkdQc9Q3CC7Zc20IarGlAlCo1vF4h+PNPQlxcCNm+XfN7CwoI8fcn5K+/+NfFEzvT\n0kiV48fJofR0je/tcfUqWXj/vh5U6Y/irGJywvMEyTiaUeb5h2sekthBsbzMdePGiEoNuokvl71d\nJCXRRc9jx4ALF4AOHdS7r0sXYNIkoHt3IDdXvxp1ITErEa03tEaJqgSHBx+Gq/WbNYAqo09gH1hK\nLPG/K//Tg0IdWbMGGD4c2LUL6NVL8/stLIAVK+gianEx//p0ZN2jRxh9+zb+qlcPbR0dNb5/sa8v\nliclIaukRA/q9EP8rHg4dXCCQ0uHMs87tHVA5pFMXlK0Can8M2cG3UTYs4fmlrdvD+zfD3h4aHb/\nlClAw4Y080XFXyYVbxy6dwjBPwejd+3eiOoTBUuppVbjcByHZR8sw6zDs5BXnMezSi0pKgI+/hj4\n7jv6bdy8ufZjdewIeHnRfHUjQUUI5sTHY+79+4hp0ADv2NlpNU5NKyt0cnLCt0lJPCvUD+kH0pG2\nIw01ltQo9xpLP0sQJUHBgwKd5yNEjfROXp4FtEDAqU2KggJCJk4kpHp1Qo4f132s994jZNYsfrTx\ngUqlIouOLSLuX7uTQ/cO8TZu3y19yYKjC3gbT2uSkwlp0YKQHj0qzmTRhFOnCKlWjZB8fhbbdOFZ\ncTHpdfUqaXH+PEmpaDFHTe7k5hLnY8dIelERD+r0R2FaITnheYKkH6w8tHSl2xWSukXz8OHrXL3a\ni4VcTJk7d2ha8v37wMWLwLvv6jaehQXNS9+0CfjtN14k6kR2YTZ6b+mNP27+gTMjz6CtT1vexo5s\nF4nl/y5Hao6AzVVOnQKCg+lj1fbtgJae6xs0bQo0aAD8+CM/42lJfH4+3r14EfYSCQ43aAA3mRp5\n0pXgZ2WFzs7O+CE5mQeF+kFVosKNATfg1t8Nju0qDy3ZBdvh2Zlnus+rKqz0GmbQjZRNm6gxHz4c\n+OMPwMmJn3FdXWkId+JE4PRpfsbUhuuPr6PJT03gbuOOo0OPopp9NV7Hr+FUA4ODBiMiJoLXcdWi\npASYP58uWqxcCXz5Jf+7u+bNAxYtEmxR5LfUVDS9cAEjPDywtlYtWPD4902tVg3fP3yIAqWStzH5\nJG4K7ePgE+mj1vW2TWyRfTZb53lVKjU+a52fA7REwKmNmmfPCBk6lJBatQi5eFF/8+zaRTNfEhL0\nN0dZPCt8RmYdmkWcFzuTDRc36HWup3lPicsSF3L7yW29zvMK164R0qwZIe3aEaJGDrZOhIcTEhmp\n3zle42lRERkUG0tqnTpFLmRn622ejpcvk5+MMOf+/sL75HTt06QoQ/2QUNHTInLU9ihRlei26e3c\nuWAWcjElLl+mC58AcO4cfarWF1270mQJQ2W+qIgK6y+uR62VtfAg6wEujb6EIQ2G6HVOJ0snjA8e\njwXHFuh1HgBAXh6tvxISAgwZQleu9V0dbc4c4JtvDFK0hxCCTampqHP2LOzEYpx/5x00tLXV23xT\nq1XD14mJUBlRAb+EJQlIWZ+CoENBkDpI1b5P6iSF1FWKvFu6LdKrVGrcr9NXhg4IOLXRUVxMyMqV\nND3ZkLs6VSpCBg8mpHdv9TYraktMfAxpuKYhaf5zc3I66bT+JiqDjPwM4rLEhdx5ekc/EyiVhPzy\nC93lFR5OSFKSfuYpj/BwQhYv1usUF7KzSbuLF0n9M2fIKb4WditBpVKRRmfPkl1GsHtUWawkt8fd\nJqcDT2u96/N6v+vk0cZHOun4919ftrHIGFGpCLl+nZAVK+jGPwcHQt59l5DbBowMvKCggJDmzQmZ\nPZv/se8+vUt6RfUi1b+pTn6/+rtgdVYijkSQoTuH8juoUklrq9SrR2svnDjB7/jqcu0aIW5uhOTk\n8D703bw80u/6deJ+4gRZmZRE
       "text": [
        "<matplotlib.figure.Figure at 0x109738650>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As promised, the rest of the RBJ designs. Note that in this version, for all types other then the shelf filters, the calculation\n",
      "of a1, a2, a3 has been reworked a bit to be done directly in terms of sine and cosine of the frequency, rather than tangent in\n",
      "Simper's derivation (but should be mathematically equivalent). This version should be slightly faster to compute and also more\n",
      "numerically robust."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def svf_core(f, res, m, shelfA = 1):\n",
      "    k = 2 - 2 * res\n",
      "    if shelfA == 1:\n",
      "        s = sin(2 * pi * f)\n",
      "        rs = 1 / (2 + k * s)\n",
      "        rc = rs * cos(2 * pi * f)\n",
      "        a1 = rs + rc\n",
      "        a2 = s * rs\n",
      "        a3 = rs - rc\n",
      "    else:\n",
      "        g = tan(pi * f) * shelfA\n",
      "        a1 = 1/(1 + g * (g + k))\n",
      "        a2 = g * a1\n",
      "        a3 = g * a2\n",
      "    A = array([[2*a1 - 1, -2*a2],\n",
      "               [2*a2, 1 - 2*a3]])\n",
      "    B = array([2 * a2, 2 * a3])\n",
      "    C_v0 = array([1, 0, 0])\n",
      "    C_v1 = array([a2, a1, -a2])\n",
      "    C_v2 = array([a3, a2, 1 - a3])\n",
      "    C = m[0] * C_v0 + m[1] * C_v1 + m[2] * C_v2\n",
      "    return A, B, C\n",
      "\n",
      "def svf_lp(f, res):\n",
      "    return svf_core(f, res, [0, 0, 1])\n",
      "\n",
      "def svf_bp(f, res):\n",
      "    return svf_core(f, res, [0, 1, 0])\n",
      "\n",
      "def svf_hp(f, res):\n",
      "    k = 2 - 2 * res\n",
      "    return svf_core(f, res, [1, -k, -1])\n",
      "\n",
      "def svf_notch(f, res):\n",
      "    k = 2 - 2 * res\n",
      "    return svf_core(f, res, [1, -k, 0])\n",
      "\n",
      "def svf_peak(f, res):\n",
      "    k = 2 - 2 * res\n",
      "    return svf_core(f, res, [1, -k, -2])\n",
      "\n",
      "def svf_bell(f, Q, gaindB):\n",
      "    A = 10 ** (gaindB / 40.)\n",
      "    k = 1/(Q * A)\n",
      "    return svf_core(f, 1 - 0.5 * k, [1, k * (A * A - 1), 0])\n",
      "\n",
      "def svf_lowshelf(f, Q, gaindB):\n",
      "    A = 10 ** (gaindB / 40.)\n",
      "    k = 1./Q\n",
      "    return svf_core(f, 1 - 0.5 * k, [1, k * (A - 1), A * A - 1], 1/sqrt(A))\n",
      "\n",
      "def svf_highshelf(f, Q, gaindB):\n",
      "    A = 10 ** (gaindB / 40.)\n",
      "    A2 = A * A\n",
      "    k = 1./Q\n",
      "    return svf_core(f, 1 - 0.5 * k, [A2, k * (A - A2), 1 - A2], sqrt(A))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "axis([3e-3, pi, -60, 10])\n",
      "f0 = 0.013\n",
      "res = 0.293\n",
      "plot_fr(svf_lp(f0, res))\n",
      "plot_fr(svf_bp(f0, res))\n",
      "plot_fr(svf_hp(f0, res))\n",
      "plot_fr(svf_notch(f0, res))\n",
      "plot_fr(svf_peak(f0, res))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FFXbwOHfZrOp1FACJEKAhBRBQBQkikZ6EUQUBCwI\n+IIUURCFFz4VVAyiiAovAlKkiKCixgIRlCJKB6mhNyGhphBI393z/TEQAoSQsruTbJ77uuaa3dnZ\nmSeZ5NmzZ04xKKUUQgghSjwXvQMQQghhG5LQhRDCSUhCF0IIJyEJXQghnIQkdCGEcBKS0IUQwkkU\nOaH369cPX19fGjRokL0tISGBNm3aUK9ePdq2bUtSUlJRTyOEEOIOipzQ+/btS3R09A3bJk6cSJs2\nbTh06BCtWrVi4sSJRT2NEEKIOzDYomPRiRMn6Ny5M3v27AEgJCSEdevW4evry9mzZ4mIiODAgQNF\nDlYIIcTt2aUO/dy5c/j6+gLg6+vLuXPn7HEaIYQQOdj9pqjBYMBgMNj7NEIIUeq52uOg16paqlWr\nxpkzZ6hateot+wQGBnL06FF7nF4IIZxWw4YN2blzZ66v2aWE3qVLF+bPnw/A/Pnz6dq16y37HD16\nFKVUgZe3337bpvvmtU9BX8vvNj0WW8ZR2GPJtZNrV5RrV5Dtznztdu3addvcW+SE3qtXL8LDwzl4\n8CB33XUX8+bNY/To0axatYp69eqxevVqRo8eXdTTZIuIiLDpvnntU9DXChKbo9kytsIeS65d4ci1\nK9z24sDR184mrVwKw2AwoNOpHWrcuHGMGzdO7zBEIci1K7mc+drllTulp6idFefSg8ibXLuSq7Re\nOymhCyFECSIldCGEKAUkoQshhJOQhC6EEE5CEroQQjgJSehCCOEkJKELIYSTkIQuhBBOQhK6EEI4\nCUnoQgjhJCShCyGEk5CELoQQTkISuhBCOAlJ6EII4SQkoQshhJOQhC6EEE7CLpNEC6EXZVFkJWaR\ndeHqEp+FJcWCNc2qLelWcAGDiwFcwOhpxLWCK64Vry4Vrq9dXKW8I0oWSeiixMmKzyJlfwqp+1NJ\nP55O+kltyTiZQebZTIzljJgqmzBVMWGqZMJYxoiLpwtGLyMu7i7a5ABWLflb06yYE82Yk8xkJWZl\nP7ZctuDm64Z7LXc8anlcX+p64BXihbufOwaDQe9fhRA3kBmLRLGklCLjdAap+1NJ3Z+ancBT96di\nzbDiFeqFd6g3HnWuJ1v3Wu64+7njYip6ydqaZSUjNoOMkxnZHxjpJ9NJO5KmxZBmxSvEC69Qr+y1\nd6g3HnU9pGQv7Cqv3CkJXehKWRXpJ9JJiUkhNSY1e526PxUXbxe8Q721pBnmlf3Yrbqb7qXjrIQs\nUg/k+LA5kEpqTCqZZzLxDPLU4r3bG+8wb7zCvPAM9LTJB40QktCF7pRFkXYs7YaknbIvhdSDqZgq\nmbQEGOaN991XE3ioF6aKJr3DLjBLqoXUA1d/xn3Xf9aM0xl41PW4/jNe/Xk9gzxxcZNEL/JPErpw\nCKUU5kQzacfSSD+aTurh1OwEnnYoDTdft+yS67WE5hXihWs557+VY0mzkHow9YYkn7IvhfR/0/Gs\nc2OJ3jPQE486HiXyA03Yny4JPTo6mldffRWLxcKLL77IqFGj8h2UKJ6UUlguW8iIyyDjdAbpx9JJ\nO5qmrY+lkXY0DRR41tUSkmeg5/VqhxAvjN5GvX+EYseSbiHtUJr2beXah99R7QPR4GrQfo91tN+n\nR+2rj2t74O7njtFLfp+lkcMTusViITg4mN9//x0/Pz/uv/9+vv76a0JDQ/MVlHAsZVGYL5nJis8i\n82wmmXGZZMRlaOvYjOuP4zIAcPfTbj5eSzbZCbyOJ64+rrrXbzsDpRRZ8VnZH5bpx9JJP379cUZc\nBkYvI25+brjX0K5HzscmX5PW0qeyCdfyck2cSV650y7fdbds2UJgYCABAQEA9OzZk6ioqBsSurAN\nZdWa3llSLTeszZfMmBPMmBPNZCVk3bA2J1xtonf1dXOyGWNZI6aKJtyqu2nJoYaWHLwbeuNe4/rz\n0lA9UhwYDAbcKrvhVtmNck3L3fL6tYR/7UP32vrK7iskrEgg83wmWRezyLqYhTXVimsl1+wEn734\nmDCWM+JazvW2a9fyrri4Sx1/SWGX/87Y2Fjuuuuu7Of+/v5s3rzZHqcqMmVRWLOsYAFlViiLynVt\n79dV1o2J2ZpqxZJmwZpqzTVhX3tdZSpcPFxw8XLB6GnExcsFF08XXMtrHWRMPqbstbu/O64+rpgq\natuyH1dwxWCUElxJkjPhl7mnTJ77WjOtZMVnZSf47CVe+4DPOJmBOdmMJdlyfX35+nMUWht+T5fs\nxeiZ93MXDxdc3FwwmAwYTAZcTC4Y3G7z2GTA4Jbj8bXFxaD9Xbpww/pap7B8bTNe7URmoFR8S7FL\nQi9Jv7iPW+/knvWXsBjBajBgcQGrC1gMBm3tYkAZcjw3GLAawOpydW0waO/L8diK9poyGLAYDKjs\nbTeus7dfXbIMrpiVEbNyx6yMWJSJLIORLBeXq4sRs6cLWd4uZLoYMRtdMBtctD9YIOev3ZAJnAPD\n+RzbDLmvcz52dQWjUVvnXHLb5uEBnp65rz08wNsbKlSA8uW1dYUK4OV143lLLaUgNRUSE7UlKQku\nX4a0tBuX9PQbn2dmgtmsLRbLjevbPHaxWHBXCneltPNe+7qe8/nNizfgpcBXYbUYsVhNWC2u2jrF\nhPWyCatVe56FO2kuHmQaPMg0uJOJB2bcMbuYMOOKxeCKBVcsmLBixIorVlxRyqj9sykjSrlgsLqA\n1aitlQsoAyi0xxgwZD+/9lhbG65uc7FeXed4fu11gzJgNSiUATCAMmi/A4X2POdjddNzuBpKjr/b\n7ONc3fPa82vvy97XkMtxyPHaTc+zj5+9Td34wh3YJaH7+flx6tSp7OenTp3C39//lv3GjRuX/Tgi\nIoKIiAh7hJOn+xcEceBSOhalMCuFBYVFgTnnc64+zn791teskL2Pmev7aq9d33b9sRULZB8vCytp\nykKaspKGhTRlIQMrbrjggQueBiOeGPEwGPHEBQ+MeBmMeHD9ucfVfTwxUg5XKuKGj8ENH9zwUsar\ndW/az52zCi7nttxyRM7l2vbMTMjIuJ5vLl68/jg9XVuuXIFLl7RcdW3JytISfMWKUK0a+PndutSu\nDTVqlNDErxScPQsnT8Lp07cuFy5cT+CurtqnXMWK2lK2rPaJmPPT8dpSvrz2C3Nzu/0n7e0eu7ho\nv8ybFqtSJBsMJMD1Ranra6W4AlzJuc6xpFxdWwAvgwH3awvgcdNz95zPDQY8AJPBgCvgenVtvPY8\n52ODAePVbdmPb9rmgjYolYvBwLW86pLb2srVJH/1uQID196jfRAYlLr+HqUVTrXtwA3vu3a9tfde\n+ycy5Pi8NKjr+2Q/z/FY209px7q2Pcef0rX3b9m1lS27t14/1qLb//nZ5aao2WwmODiYP/74gxo1\natC0aVO5KVoIVqVIs1pJsViuL/l4fsViIT4ri3OZmdqSlUWW1Yqvm9v1xWTC182Najm2VXNzI8DD\nA3cX+9WZZmRoST4xEc6cgdjYW5djx7QPhLCwG5emTaFSJbuFVnCXL8OWLbBjB+zfDzExcOCAlnRr\n1wZ//+vLtU8rX9/rSdzd3W6hpVosHEtL42RGBnEZGcRlZt6yvpCZiZfRiI+rK5VMJnxMJnxcXW9Y\nlzUaKXPT4n3TczeDoUR9Ky/pdGm2uGLFiuxmi/379+e///1vvoMStpdqsXAuM5OzOZJ8dsK/uv1M\nZiaxGRkEeHgQ6u1NmJcXYd7ehHp5EeLlhZfRcc3kLl68niP374e9e2HrVi0ntmihLQ89BFfvuzvG\n6dOwdi1s2KAtR45A48Zw333aJ05oqLY46FPHqhQn09PZnZLC3pQUjqSlcfTqEp+VRYCHBwEeHvi5\nu1PD3Z0abm7Z6+publR1c8PN
       "text": [
        "<matplotlib.figure.Figure at 0x1097dcbd0>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "axis([3e-3, pi, -40, 40])\n",
      "for gain in range(-30, 31, 10):\n",
      "    plot_fr(svf_bell(f0, .5, gain))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYVMfeB/AvvfeyVMVCEY2AKNYoiDUGxJhojD3GxGiK\nvknEe280Jl4VW2wxibHGaKJ4ExHssWCJir2BgihKXeqylF3Y9nv/GMUYI0H3LM35PM95ljozwPLd\nOXNm5ugREYHjOI5r8vQbugEcx3GcMHigcxzHNRM80DmO45oJHugcx3HNBA90juO4ZoIHOsdxXDMh\nSKCr1WoEBQUhIiICAFBSUoL+/fvDx8cHAwYMQGlpqRDVcBzHcbUQJNBXrlwJf39/6OnpAQBiYmLQ\nv39/pKWlITw8HDExMUJUw3Ecx9VC60DPzs7Gvn378M477+DhGqX4+HiMHz8eADB+/HjExcVpWw3H\ncRz3D7QO9BkzZmDJkiXQ139UVH5+PkQiEQBAJBIhPz9f22o4juO4f6BVoO/ZswfOzs4ICgrC03YQ\n0NPTqxmK4TiO43THUJtvPn36NOLj47Fv3z5UVVWhrKwMY8eOhUgkglgshouLC/Ly8uDs7PzE97Zt\n2xZ37tzRpnqO47gXTkBAAK5cufL3nySBJCYm0quvvkpERJ999hnFxMQQEdHChQspOjr6ia9/3qq/\n+OILQb+2tq951s/V9WMNQch2PG9Z/G/3fPjf7tk/3pz/drVlp6Dz0B8OrcyaNQu///47fHx8cPTo\nUcyaNUuwOkJDQwX92tq+5lk/9yxtq29Ctu15y+J/u+fD/3bP9/HGoL7/dnoPEr/e6enpPXXcvTmZ\nO3cu5s6d29DN4J4D/9s1Xc35b1dbdvKVojrWmHsPXO34367pelH/dryHznEc14TwHjrHcdwLgAc6\nx3FcM8EDneM4rpnggc5xHNdM8EDnOI5rJnigcxzHNRM80DmO45oJHugcx3HNBA90juO4ZoIHOsdx\nXDPBA53jOK6Z4IHOcRzXTPBA5ziOayZ4oHMcxzUTPNA5juOaCR7oHMdxzQQPdI7juGZCq0CvqqpC\n165dERgYCH9/f/zrX/8CAJSUlKB///7w8fHBgAEDUFpaKkhjOY7juKfT+hZ0MpkM5ubmUKlU6NWr\nF5YuXYr4+Hg4Ojpi5syZWLRoESQSCWJiYh6vmN+CjuM47pnp9BZ05ubmAACFQgG1Wg07OzvEx8dj\n/PjxAIDx48cjLi5O22o4juO4f6B1oGs0GgQGBkIkEiEsLAzt27dHfn4+RCIRAEAkEiE/P1/rhnIc\nx3G1M9S2AH19fVy5cgVSqRQDBw7EsWPHHvu8np4e9PT0tK2G47SWmZmJZcuWoW3btujcuTMCAgJq\nzjA5rjnQOtAfsrGxwZAhQ3Dx4kWIRCKIxWK4uLggLy8Pzs7Of/s9c+fOrXk7NDQUoaGhQjWH4x5z\n6dIlREZGYsSIEUhOTsaPP/6IlJQUvP/++1i6dCnvdHCNVmJiIhITE+v0tVpdFC0qKoKhoSFsbW0h\nl8sxcOBAfPHFFzh48CAcHBwQHR2NmJgYlJaW8ouiXIPZt28fxo8fj++//x7Dhw+v+XhpaSnCwsIw\ndOjQxzoXHNeY1ZadWvXQ8/LyMH78eGg0Gmg0GowdOxbh4eEICgrCiBEjsGHDBnh5eSE2Nlabajju\nua1duxZz585FfHw8unfv/tjnbG1tceDAAbz88suwt7fHRx991ECt5DhhaD1t8bkr5j10TscOHDiA\nd955B4mJiWjbtu1Tv+7evXt4+eWXERMTg9GjR9djCznu2dWWnTzQuWYpPz8fQUFB2LZtG8LCwv7x\n65OTkxEaGoojR46gY8eO9dBCjns+Op2HznGNjUajwYQJEzBx4sQ6hTkAtG/fHsuWLcOoUaMgk8l0\n3EKO0w3eQ+eaneXLlyM2NhYnTpyAkZFRnb+PiDB69GjY2Njgu+++02ELOe758SEX7oVx6dIlDBw4\nEOfOnUOrVq2e+fulUimCgoKwbNkyDBs2TAct5Djt8CEX7oVQUVGBUaNGYdWqVc8V5gBbT/Hzzz9j\nypQpyM7OFriFHKdbvIfONRuTJk2CWq3G5s2btS5r/vz5OHz4MA4fPgwDAwPtG8dxAuE9dK7Zezhm\nvnr1akHKmzVrFogIixYtEqQ8jqsPvIfONXn37t1DSEgI9u3bh86dOwtWblZWFjp37ozdu3ejW7du\ngpXLcdrgPXSu2VKpVBg9ejRmzpwpaJgDgKenJ77//nu89dZbkEqlgpbNcbrAe+hckzZnzhwkJSVh\n//790NfXTf9kypQpqKiowNatW3VSPsc9Cz5tkWuWTpw4gZEjR+Ly5ctwcXHRWT0ymQxdunTBrFmz\nMHbsWJ3Vw3F1wQOda3Zyc3MREhKCH374Aa+88orO67t27RrCw8Nx5syZWveF4Thd42PoXLMil8sR\nFRWFqVOn1kuYA0DHjh3xxRdf4PXXX0dlZWW91Mlxz4r30LkmhYgwZswYaDQa/Pzzz/V6YwoiwoQJ\nE1BVVYXt27fzm2JwDYL30LlmY/HixUhNTcWGDRvqPVD19PSwdu1aZGRkPHHDFo5rDAS7BR3H6VpC\nQgJWr16NpKSkBrsXqKmpKXbt2oWQkBC0a9cOUVFRDdIOjvs7fMiFaxKSk5MRFhaGhIQEdO3ataGb\ngwsXLmDw4MFISEjgi464esWHXLgmLSsrC6+++iqWLVvWKMIcADp37ozNmzdj2LBhuHPnTkM3h+MA\n8EDnGrm8vDz07dsXH330UaObAz5kyBDMnTsXAwYM4Dszco2CVoGelZWFsLAwtG/fHh06dMCqVasA\nACUlJejfvz98fHwwYMAAlJaWCtJY7sVSUFCA8PBwTJw4ETNmzGjo5vyt9957D++//z769u2LvLy8\nhm4O94LTagxdLBZDLBYjMDAQFRUVCA4ORlxcHDZt2gRHR0fMnDkTixYtgkQieWJWAB9D52pTUlKC\nsLAwDB06FF999VVDN+cfzZ8/H9u2bUNiYiKcnZ0bujlcM1ZvK0WjoqLwwQcf4IMPPsDx48chEokg\nFosRGhqKW7du1blR3ItNKpWiX79+CA0NxeLFi5vMfO85c+YgLi4Ox44dg4ODQ0M3h2um6iXQ7927\nhz59+uDGjRto0aIFJBIJALYYw97evub9ujSKe3GVl5dj4MCBCA4OxqpVq5pMmAPsuT5r1iwcPnwY\nv//+O+zt7Ru6SVwzVFt2CjIPvaKiAsOHD8fKlSthZWX1ROVN6Z+SaziFhYWIiIhAQEAAVq5cKdjz\nRqUqQ0nJQRQXJ6Cs7Az09c1haGgHIyM7WFl1hovLBJiYuGtdj56eHmJiYjBz5kz06dMHBw4cgLu7\n9uVyXF1pHehKpRLDhw/H2LFjaxZZPBxqcXFxQV5e3lPHFOfOnVvzdmhoKEJDQ7VtDtdEpaenY9Cg\nQRg1ahS++uorQcK8ouIq7t79N6TSk7Cx6QkHh0i0aDELGo0CKpUEKlUJSkp+x/nzL8HGphfc3afB\n3n6gVnXq6elh8eLFcHR0RK9evXDw4EH4+Pho/bNwL67ExEQkJibW7YtJCxqNhsaOHUvTp09/7OOf\nffYZxcTEEBHRwoULKTo6+onv1bJqrhk5e/Ysubi40Nq1awUpT6kspbS0j+jUKSfKzv6WlMqyWr9e\npaqg3NyNdPasD127NpTk8vuCtGP9+vXk6upKFy5cEKQ8jiOqPTu1StWTJ0+Snp4eBQQEUGBgIAUG\nBtL+/fupuLiYwsPDydvbm/r3708SieSZGsW9OHbv3k2Ojo6UkJCgdVkajYby8n6iP/5wpVu33qHq\n6sJn+n61uooyMr6ikycdKDNzGWk0aq3btGvXLnJycqJDhw5pXRbHEdWenXzpP9dgvvvuO8ybNw+7\nd+9Gly5dtCpLpSpHauo7kMluwsfnB9jYPP9yfJnsNm7dmggDA0u0a7cFxsbaTUM8ceIERowYgdmz\nZ2PatGlalcVx/AYXXKOiVqvx
       "text": [
        "<matplotlib.figure.Figure at 0x109450510>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "axis([3e-3, pi, -40, 40])\n",
      "for gain in range(-30, 31, 10):\n",
      "    plot_fr(svf_lowshelf(f0, .707, gain))\n",
      "    plot_fr(svf_highshelf(f0, .707, gain))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8FGX+wPHPtuxuet80IAmQ0LtBLAjnAVJEbNwpKiqc\neOpZzoJdPE/BfvaC/hTLKWBBBBRQTgXpHUJJIISENNI3yZbs7jy/PwZCIknYNBLged9rbmZn5pnn\n2ez65dlnnnkejRBCIEmSJJ3xtO1dAEmSJKl1yIAuSZJ0lpABXZIk6SwhA7okSdJZQgZ0SZKks4QM\n6JIkSWeJVgnoHo+HgQMHcvnllwNQUlLCqFGjSEpKYvTo0ZSVlbVGNpIkSVIjWiWgv/baa/Tq1QuN\nRgPAnDlzGDVqFGlpaVx66aXMmTOnNbKRJEmSGtHigH7kyBGWLVvG9OnTOf6M0uLFi5k6dSoAU6dO\nZdGiRS3NRpIkSTqFFgf0++67jxdffBGt9sSlCgoKsFgsAFgsFgoKClqajSRJknQKLQroS5YsITIy\nkoEDB9LQCAIajaamKUaSJElqO/qWJF67di2LFy9m2bJlOBwOrFYrN954IxaLhfz8fKKiosjLyyMy\nMvKktN26dePgwYMtyV6SJOmc079/f7Zv317/QdFKfvnlFzFhwgQhhBAPPvigmDNnjhBCiNmzZ4uZ\nM2eedH5zs37qqada9dzGzmnqMW/3tYfWLEdzryU/u+aRn13T95/Nn11jsbNV+6Efb1p5+OGHWbly\nJUlJSaxatYqHH3641fIYMWJEq57b2DlNPdaUsp1urVm25l5LfnbNIz+75u3vCE73Z6c5FvFPO41G\n02C7+9lk1qxZzJo1q72LITWD/OzOXGfzZ9dY7JRPiraxjlx7kBonP7sz17n62ckauiRJ0hlE1tAl\nSZLOATKgS5IknSVkQJckSTpLyIAuSZJ0lpABXZIk6SwhA7okSdJZQgZ0SZKks4QM6JIkSWcJGdAl\nSZLOEjKgS5IknSVkQJckSTpLyIAuSZJ0lpABXZIk6SwhA7okSdJZQgZ0SZKks4QM6JIkSWcJGdAl\nSZLOEi0K6A6Hg6FDhzJgwAB69erFI488AkBJSQmjRo0iKSmJ0aNHU1ZW1iqFlSRJkhrW4inobDYb\nvr6+uN1uLrroIl566SUWL15MeHg4Dz30EM8//zylpaXMmTOnbsZyCjpJkqQma9Mp6Hx9fQGorq7G\n4/EQEhLC4sWLmTp1KgBTp05l0aJFLc1GkiRJOoUWB3RFURgwYAAWi4WRI0fSu3dvCgoKsFgsAFgs\nFgoKClpcUEmSJKlx+pZeQKvVsn37dsrLyxkzZgz/+9//6hzXaDRoNJqWZtNmMjP/xdGjXwIaQING\no62zrrutPfZeam8fP7+h7ZOvU/8+HVqtLzqdb71rrdbc4LET5xg79N9aOnMoikJ5eTmlpaWUlJRQ\nWVmJzWbDbrdjs9nqLHa7HYfDgdvtxu1243K56qzr26coCkKImuWPrxva19i5TdHU5t6mnN+W1z6V\nFgf044KCghg/fjxbtmzBYrGQn59PVFQUeXl5REZG1ptm1qxZNdsjRoxgxIgRrVUcr0VH/42IiGsQ\nQgEEIGptK8f+2LW3lWPn1N6uL+2Jbe+u7cbjsaMoNhTFjsdjw+U6isNhQ1FseDwNre01r4Worgn8\nPj4xmEwJmM2JmEyJNWuTKR6dznRa/8ZSx+HxeMjMzOTgwYPk5OTULLm5ueTm5lJcXExpaSkVFRX4\n+/sTGhpKSEgIAQEB+Pr6Yjab8fX1PWk7KCgIg8GAXq9Hr9fXbNe3T6/Xo9Vqayp7tbcb2+fNuU3R\nlue35rU3bNjAxo0ba16np6c3fJ2W3BQtKipCr9cTHByM3W5nzJgxPPXUUyxfvpywsDBmzpzJnDlz\nKCsrkzdFTwMhPMcCfBVOZy4ORwZ2e8ax9SEcjgwcjiwMhrBaAb520O+K0Rjd3m9DagVCCLKysti4\ncSPbtm1j37597N+/n4yMDCIjI+nWrRtxcXHExsYSExNTsw4LCyM0NJSgoCB0Ol17vw2pHo3FzhYF\n9F27djF16lQURUFRFG688UYefPBBSkpKmDx5MllZWcTHx7NgwQKCg4O9LpTUdoTwNBjsbbY0DIYQ\nQkPHExY2juDgS9Bqje1dZMkLiqKwdetWVq5cydq1a9m4cSMajYahQ4cyaNAgevbsSXJyMt27d6/p\nyCCdmdosoLeEDOgdjxCCysodlJQspbh4GVVVuwkOHkFY2HhCQ8diMnVq7yJKtVRWVrJkyRIWLVrE\nTz/9RGRkJKNGjWL48OGkpKQQFxcn76mchWRAl5rF5SqmpGQ5xcVLKSlZjtEYS1jYOEJDxxMYeD5a\nbavdgpG85HQ6Wbx4MfPnz2flypVceOGFXH311YwZM4a4uLj2Lp50GsiALrWYEB6s1g0UFy+jpGQp\nTmcucXF3ExNzJwZD8KkvILXI4cOHee+99/jwww/p06cPU6ZMYdKkSYSGhrZ30aTTTAZ0qdVVVe0l\nK2sOxcVLiImZQVzcvfj41N+bSWq+1atX8+KLL/L7779z0003cfvtt5OcnNzexZLakQzoUpux2zPJ\nzn6Bo0e/xGK5gU6dHsBk6tzexTqjCSFYtWoVzzzzDNnZ2Tz00EPccMMN+Pn5tXfRpA5ABnSpzTmd\neRw58ip5eR8SHj6Jzp1n4uub1N7FOuNs27aNBx54gOzsbB5//HGuv/569Hp5r0I6QQZ06bRxuUrI\nyXmDnJw3CQ0dS9euL+PjE9HexerwcnJyeOyxx1i+fDlPPfUU06dPl4FcqlebDs4lSbUZDKHExz/F\n0KEZ+PhY2LSpD/n5n8p/vBvg8Xh4/fXX6d+/PzExMezfv5/bb79dBnOpWWQNXWpTFRVb2LdvGj4+\nkSQlvYfZnNDeReowtmzZwowZMwgICOCdd96hR48e7V0k6Qwga+hSuwkIGMzgwZsICbmULVvOIzv7\nFYTwtHex2lV1dTVPPPEE48aN4+6772bVqlUymEutQtbQpdPGZksnLe02PJ4qkpPn4u/fv72LdNql\npqZy0003ERUVxQcffEB0tBw7R2oaWUOXOgRf3+7077+KmJgZ7Ngxiuzsl8+Zf9SFELz22mtccskl\n3H777SxZskQGc6nVyRq61C7s9kz27JmMj08MPXp8fFY/bVpaWsqtt97KkSNH+PLLL+natWt7F0k6\ng8kautThmM3xDBy4GpOpM1u2DKaiYkt7F6lNbN68mcGDB9OpUyfWrFkjg7nUpmQNXWp3R48uID39\nTuLjnyEmZsZZMUKgEIJ33nmHWbNm8fbbb3PNNde0d5Gks4R8sEjq8Gy2NFJTr8XPrw9JSe+h1/u3\nd5Gazel0cuedd7Jhwwa+/fZbunXr1t5Fks4isslF6vB8fZMYNGg9Wq2RbduGYbdntneRmiU/P5+R\nI0dSUlLCunXrZDCXTisZ0KUOQ6czk5z8IdHR09m27QLKy9e2d5GaZNOmTZx33nlcdtllfPXVV/j7\nn7m/MqQzk2xykTqk4uIf2LdvKl27vkJU1A3tXZxTWrhwIXfeeSfvv/8+kyZNau/iSGcx2YYunZGq\nqlLZtetyIiOvIyHhGTSajveDUgjBCy+8wJtvvsn333/PgAED2rtI0lmuzdrQs7OzGTlyJL1796ZP\nnz68/vrrAJSUlDBq1CiSkpIYPXo0ZWVlLclGOkf5+fVm0KANlJX9SmrqZDweW3sXqQ6Xy8WMGTP4\n4osvWLdunQzmUrtrUQ09Pz+f/Px8BgwYQGVlJYMHD2bRokV89NFHhIeH89BDD/H8889TWlrKnDlz\n6mYsa+iSlxTFyf7907HbD9K37/cYDGHtXSTKy8uZPHkyOp2O+fPnExAQ0N5Fks4RbVZDj4qKqqmV\n+Pv707NnT3Jycli8eDFTp04FYOrUqSxatKgl2UjnOK3WSI8e8wgKuoht2y7G4chu1/Lk5eUxfPhw\nunbtyuLFi2UwlzqMVmuUzMzM
       "text": [
        "<matplotlib.figure.Figure at 0x1097aa810>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Fast matrix evaluation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The matrix form above is reasonably efficient to implement. However, to maximally take advantage of the power of modern\n",
      "SIMD, we can construct a single 4x4 matrix that computes _two_ samples in a single iteration. Thus, the iteration of the\n",
      "state is represented by the square of the $\\mathbf{A}$ matrix above."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mk4mat(filter):\n",
      "    A, B, C = filter\n",
      "    AA = dot(A, A)\n",
      "    AB = dot(A, B)\n",
      "    CA = dot(C[1:], A)\n",
      "    cb = dot(C[1:], B)\n",
      "    return array([[C[0], 0, C[1], C[2]],\n",
      "                 [cb, C[0], CA[0], CA[1]],\n",
      "                 [AB[0], B[0], AA[0][0], AA[0][1]],\n",
      "                 [AB[1], B[1], AA[1][0], AA[1][1]]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then, evaluation is again quite simple, each iteration processing two samples. The central 4x4 matrix-vector multiply is especially\n",
      "well suited for implementation in NEON, see cpp/src/neon_iir.s for NEON assembly language."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def eval4(A, inp):\n",
      "    out = zeros(len(inp))\n",
      "    y = zeros(4)\n",
      "    for i in range(0, len(inp), 2):\n",
      "        y[0:2] = inp[i:i+2]\n",
      "        y = dot(A, y)\n",
      "        out[i:i+2] = y[0:2]\n",
      "    return out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's test it works the same:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filter = svf_bp(f0, res)\n",
      "plot(eval(filter, impulse))\n",
      "plot(eval4(mk4mat(filter), impulse))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVGX+B/DPmQt3ECFFZTAURrkZqChltWKGJimWuhta\nSWWuP11XbdutvbSbtZuXbXe7WZu25dpl0cptJUVKM9JSJBdSE1JgobgoXnCEmWGYmTPn94dJTXIZ\nmBnOAJ/36zUv5hyec853HnU+PueZc0aQJEkCERHRtxRyF0BERJ6FwUBERHYYDEREZIfBQEREdhgM\nRERkh8FARER2nA6GvLw8xMTEQKvVYv369W22WbFiBbRaLRITE1FcXNy6fu3atYiPj8eYMWOwYMEC\ntLS0OFsOERE5yalgEEURy5cvR15eHkpKSpCdnY3S0lK7Nrm5uSgvL0dZWRk2bdqEpUuXAgCqqqrw\nyiuvoKioCMePH4coiti6dasz5RARkQs4FQyFhYWIjo5GZGQk1Go1MjMzsWPHDrs2OTk5yMrKAgCk\npKRAp9Ohvr4eQUFBUKvVMBqNsFqtMBqNCA8Pd6YcIiJyAaeCoba2FhEREa3LGo0GtbW1DrUJCQnB\nww8/jOHDh2PYsGEIDg7Grbfe6kw5RETkAk4FgyAIDrVr664bFRUVePbZZ1FVVYW6ujro9Xq89dZb\nzpRDREQuoHJm4/DwcFRXV7cuV1dXQ6PRdNimpqYG4eHhyM/Px6RJkxAaGgoAmDNnDg4ePIi7777b\nbvvo6GhUVFQ4UyYRUb8TFRWF8vLybm3r1IghOTkZZWVlqKqqgtlsxrZt25CRkWHXJiMjA6+//joA\noKCgAMHBwQgLC8Po0aNRUFCA5uZmSJKEvXv3Ii4u7qpjVFRUQJIkPiQJjz/+uOw1eMqDfcG+YF90\n/HDmP9ROjRhUKhU2bNiA6dOnQxRFLFq0CLGxsdi4cSMAYMmSJUhPT0dubi6io6Ph7++PzZs3AwCS\nkpKwcOFCJCcnQ6FQYNy4cfjpT3/qTDlEROQCTgUDAMyYMQMzZsywW7dkyRK75Q0bNrS57SOPPIJH\nHnnE2RKIiMiFeOVzL5Kamip3CR6DffEd9sV32BeuIUiS5NFf1CMIAjy8RCIij+PMeydHDEREZIfB\nQEREdhgMRERkh8FARER2GAxERGSHwUBERHYYDEREZIfBQEREdhgMRERkh8HQRfc9/yrOXjTIXQYR\nkdswGLrozeonsWVfgdxlEBG5DYOhi2xKAw5VfCl3GUREbsNg6CJJZcCJc8flLoOIyG0YDF1gtoiA\n2oRaK0cMRNR3MRi64KzOAIgqGPxOwCra5C6HiMgtGAxdcL7RAEVLKJTmYBw88bXc5RARuQWDoQsu\nNBqgEP0x0JqAD45ynoGI+iYGQxdcaDJAZfPHCL8x+LyK8wxE1DcxGLrgQpMeKlsAkoYl4KSOIwYi\n6psYDF2gMxjgBX9MiR+DeokjBiLqmxgMXaAzGuAl+GP6uBi0+JVD32yWuyQiIpdjMHRBY7MB3gp/\nhAT5wqv5WuwpOiV3SURELsdg6IJLzXr4KgMAAIOkBOz7kvMMRNT3MBi6oKnFAD+VPwBAO2AMimo5\nz0BEfQ+DoQsMZgP81JeDITkiARVNHDEQUd/DYOgCg8UA/2+D4dbrxuCCkiMGIup7GAxdYLToEeh9\neY5hSmIUrD5nUHehSeaqiIhci8HQBc2iAYE+l0cMXmolfA2x2H2kROaqiIhci8HQBSabAQN8/VuX\nhyoTsP8rzjMQUd/CYOgCs2QfDLGhY/DF6WMyVkRE5HoMhi5okfQYGBDQuvwjbRK+Nh2VsSIiItdj\nMHSBRTBgoP93I4ZZExNxyfcobDZJxqqIiFyLwdAFosKAkMDvnUoaPggK0R8HS/ilPUTUdzAYusCq\n0CP0e8EAAKGWROQW8XQSEfUdTgdDXl4eYmJioNVqsX79+jbbrFixAlqtFomJiSguLm5dr9PpMG/e\nPMTGxiIuLg4FBQXOluNWksqAwcEBduuiA5JwqPILmSoiInI9p4JBFEUsX74ceXl5KCkpQXZ2NkpL\nS+3a5Obmory8HGVlZdi0aROWLl3a+ruVK1ciPT0dpaWlOHbsGGJjY50px+0ktQGDBtiPGCYMT8RJ\nHUcMRNR3OBUMhYWFiI6ORmRkJNRqNTIzM7Fjxw67Njk5OcjKygIApKSkQKfTob6+HpcuXcKBAwfw\nwAMPAABUKhUGDBjgTDlupW82A4INAb5edutvS0rCOSWDgYj6DqeCoba2FhEREa3LGo0GtbW1nbap\nqalBZWUlBg0ahPvvvx/jxo3D4sWLYTQanSnHrc7pDIDFHwqFYLd+alI0rN71qDnXKFNlRESupXJm\nY0EQOm8EQJLsP84pCAKsViuKioqwYcMGTJgwAatWrcK6devw5JNPXrX96tWrW5+npqYiNTXVmbK7\n5XyjAQprwFXrvdRK+BsTkHP4GJbNvKnH6yIiAoD8/Hzk5+e7ZF9OBUN4eDiqq6tbl6urq6HRaDps\nU1NTg/DwcEiSBI1GgwkTJgAA5s2bh3Xr1rV5nO8Hg1zONxqgFP3b/F2EOhEfl37BYCAi2fzwP81P\nPPFEt/fl1Kmk5ORklJWVoaqqCmazGdu2bUNGRoZdm4yMDLz++usAgIKCAgQHByMsLAxDhgxBREQE\nTp26/PWYe/fuRXx8vDPluFVDkwEqW9vBkBSWhONnOc9ARH2DUyMGlUqFDRs2YPr06RBFEYsWLUJs\nbCw2btwIAFiyZAnS09ORm5uL6Oho+Pv7Y/Pmza3bv/DCC7j77rthNpsRFRVl9ztP06DXQyW1HQyp\nsYnIqXmthysiInIPQfrhBICHEQThqjkKOTyZvRvPFz6H88/kXfW7ugtNCP/bEDQ/fgk+Xk5lLRGR\nSzjz3skrnx10yWiAt9D2iGFYaCDUpmH4sOhUD1dFROR6DAYHNZoM8FG0HQwAMFhKwp6jnGcgot6P\nweCgS816+CjbD4aY4EQUVvPWGETU+zEYHKRvMcBPdfV1DFdMGpmECgNHDETU+zEYHGSwGOCvbn/E\nMGtCEhq8ivndDETU6zEYHGSwGODv1X4wjNeGA5Dw37LadtsQEfUGDAYHGS16BHq3HwwKhYBQ83j8\np7CoB6siInI9BoODTKIBQT7tzzEAwOjA8fi04r89VBERkXswGBxkshkwwLf9EQMATBoxDl9dYjAQ\nUe/GYHCQWTJggF/HwTB74nicVzMYiKh3YzA4yAw9gv07DoYbYodDUlhQVFbXQ1UREbkeg8FBFsGA\n0MCO5xgUCgEhpvH4z2FOQBNR78VgcJCoMCAkoOMRAwBoA8fhACegiagXYzA4SFQaEBrUeTBMihyP\nUh2DgYh6LwaDg2wqPQYN6DwYZk8Yj3OcgCaiXozB4ACbTQLUBoeC4aaESEjKZhz735keqIyIyPUY\nDA5oNLYANhX8fNSdtlUoBAQ3j8N7BZyAJqLeicHggHM6AwRr56OFK7QB47G/nKeTiKh3YjA44Owl\nPRRdCIZJkeNRcpHBQES9E4PBAReaDFCKHV/D8H2zksfjnIrBQES9E4PBARebDFBJjo8YUq8bCVHV\nhBNVZ91YFRGRezAYHHDRYIC6C8GgUAgY2Dwe7x484saqiIjcg8HggAa9Hl6C48EAADGBKdh36rCb\nKiIich8GgwMuNRvgLTg+xwAAk6MnouRSoZsqIiJyHwaDAxqbDfBRdG3EcNdNKbjgXcjvgCaiXofB\n4ICmFgN8VV0LhqSooVCI/vj4aIWbqiIicg8GgwOaWvRdDgYAGCJOxPbDnGcgot6FweAAg9mAAK+u\nzTEAQOI1KTj4NecZiKh3YTA4wGgxIMCr6yOGafETUW7iiIGIehcGgwOaRQMCvLseDHfdPB4Gv+PQ\nN5vdUBURkXswGBzQLOoR5NONOYaQAPg0R2H7p0fdUBURkXswGBzQYjMg2K/rcwwAcK0yBbuOcp6B\niHoPBoMDWiQDBvh1fcQAABPD
       "text": [
        "<matplotlib.figure.Figure at 0x109489d90>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With the 4x4 matrix formulation, each iteration performs 16 multiplies and 12 adds to compute two samples, or 8 multiplies and\n",
      "6 adds per sample. This is slightly better than the 9 multiplies and 8 adds for the most general form of Simper's code,\n",
      "and roughly the same as most of the specialized forms. However, raw number of multiply/add operations tells only a small part\n",
      "of the speed story. All these operations are grouped for easy computation in 4x wide SIMD instructions. On NEON, this yields\n",
      "1.47ns per IIR on a Nexus 5, which is 6x a straightforward scalar implementation of direct form I biquad evaluation."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Error measurement"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to evaluate these filters in floating point arithmetic for speed. IIR filters are notoriously finicky about\n",
      "roundoff error.\n",
      "\n",
      "For these tests, we'll run a filter in double math, then the same filter in float, and compare."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def eval4_f32(A, inp):\n",
      "    A = array(A, dtype=float32)\n",
      "    inp = array(inp, dtype=float32)\n",
      "    out = zeros(len(inp), dtype=float32)\n",
      "    y = zeros(4, dtype=float32)\n",
      "    for i in range(0, len(inp), 2):\n",
      "        y[0:2] = inp[i:i+2]\n",
      "        y = dot(A, y)\n",
      "        out[i:i+2] = y[0:2]\n",
      "    return out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "impulse = concatenate([[1], zeros(99)])\n",
      "filter = mk4mat(svf_lp(.1, .75))\n",
      "plot(eval4(filter, impulse) - eval4_f32(filter, impulse))\n",
      "grid(True)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEGCAYAAABxfL6kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VFWeB/BvhQQQWcIiYUkkyJZVKjTLNC1SdLrAAInN\n4ii2OAHHnh6OqHSP3W0zDj3OQKAdFFCPbbsAag8wOtI42KRFTIHKJoQCWQSMRLIQQCDR7NubP66V\nW9krVfe9eqn6fs7JgVf13ru3foFfbn7vvvssmqZpICKigBHi7w4QEZFaTOxERAGGiZ2IKMAwsRMR\nBRgmdiKiAMPETkQUYHRP7IsXL0ZERAQSExOVnO/RRx9FfHw84uLi8Nhjjyk5JxFRINE9sS9atAiZ\nmZlKzuVwOJCdnY2TJ0/i5MmT+Oyzz7B3714l5yYiChS6J/YpU6agb9++jV7LyclBSkoKxo8fjzvv\nvBNnz5716FwRERGorq5GVVUVKioqUFNTg0GDBunRbSKiTivUH43+/Oc/x8svv4yRI0fi0KFDWLJk\nCfbs2dPucbGxsZg+fToGDx4MTdOwdOlSjBkzxoAeExF1HoYn9tLSUhw4cAD33HNPw2vV1dUAgHff\nfRcrVqxodkxkZCR27dqFffv2ISsrCwUFBdA0DXa7HTNmzMAdd9xhWP+JiMzO8MReX1+P8PBwHDt2\nrNl7c+fOxdy5c1s99uDBg0hJSUGPHj0AACkpKThw4AATOxGRG59q7JWVlZg0aRKsVivi4uLw5JNP\ntntM7969MXz4cLzzzjsAAE3TcOLECY/ai4mJwd69e1FXV4eamhrs3bsXcXFxvnwEIqKA41Ni7969\nO7KysuB0OnHixAlkZWXhk08+abTPggULMHnyZJw9exZRUVHYuHEj/vznP+O1116D1WpFQkIC3nvv\nPY/aS0tLQ0JCAsaOHQur1Qqr1YpZs2b58hGIiAKORdWyveXl5Zg6dSo2b97MUTQRkR/5PN2xvr4e\nVqsVERERmDZtGpM6EZGf+ZzYQ0JC4HQ6kZ+fj3379sHhcCjoFhEReUvZrJg+ffpg1qxZOHLkCGw2\nW8PrQ4cORWFhoapmiIiCwogRI/Dll196daxPI/ZvvvkGxcXFAICKigrs3r0bSUlJjfYpLCyEpmn8\n0jSsWLHC730wyxdjwVgwFm1/5eTkeJ2bfRqxX7p0Cf/wD/+A+vp61NfXY+HChUhOTvbllAEtNzfX\n310wDcZCYiwkxkINnxJ7YmIisrOzVfWFiIgU4HrsBkpPT/d3F0yDsZAYC4mxUEPZPPZWG7BYoHMT\nREQBx5fcyRG7gTgVVGIsJMZCYizUYGInIgowLMUQEZkQSzFERNSAid1ArB9KjIXEWEiMhRpM7ERE\nAYY1diIiE2KNnYiIGjCxG4j1Q4mxkBgLibFQg4mdiCjAsMZORGRCrLETEVEDJnYDsX4oMRYSYyEx\nFmowsRMRBRjW2ImITIg1diIiasDEbiDWDyXGQmIsJMZCDSZ2IqIAwxo7EZEJscZOREQNmNgNxPqh\nxFhIjIXEWKjhc2LPy8vDtGnTEB8fj4SEBGzYsMGj44qKgNJSX1snIqKmfK6xFxUVoaioCFarFaWl\npfjBD36Av/zlL4iNjRUNtFInevhhYNIk4B//0ZfWiYgCk19r7IMGDYLVagUA9OzZE7GxsSgsLGz3\nuLIyoKrK19aJiKgppTX23NxcHDt2DJMmTWp338pKoKZGZevmx/qhxFhIjIXEWKgRqupEpaWlmD9/\nPtavX4+ePXs2ei89PR3R0dEAgPDwcFitVlRV2VBTI7+RNpsNALeDZdvFLP3x57bT6TRVf/y57XQ6\nTdUfI7cdDgc2bdoEAA350ltK5rHX1NRg9uzZSElJweOPP964gVbqRMnJ4ut3v/O1dSKiwOPXGrum\naXjooYcQFxfXLKm3JRhLMURERvA5sX/66ad46623kJWVhaSkJCQlJSEzM7Pd44IxsTctQwQzxkJi\nLCTGQg2fa+x33HEH6uvrO3xcMCZ2IiIj+G2tmJEjgbQ04Nln9WydiKhz6pRrxXDETkSkDyZ2A7F+\nKDEWEmMhMRZqMLETEQUYv9XYQ0OB++8H3nhDz9aJiDqnTldjr60F6uo4Yici0oNfErtr8a9gS+ys\nH0qMhcRYSIyFGn5J7JWV4s/aWn+0TkQU2PxSYy8oACIjgZkzgfff17N1IqLOqdPV2F0j9mArxRAR\nGYGJ3UCsH0qMhcRYSIyFGrx4SkQUYPxSY9+/H5gyBRg/Hjh0SM/WiYg6p05ZY+/ZkyN2IiI9+C2x\n9+oVfImd9UOJsZAYC4mxUIOJnYgowPilxr5li1iH/fp1ICdHz9ab+/RToLwcsNuNbZeIqCM6ZY3d\nHyP2Q4eAOXOABx4Atm0ztm0iIqMETWI/dUo8sWnjRuDDD4Fly4BNm4xrH2D90B1jITEWEmOhhs/P\nPPWG0Yn9wgXgrrtE+WfWLPFaVhbwk5+IOfX/9E/G9IOIyAh+qbFnZABffy1q7SUlerYu/OhHwD33\nAI8/3vj1s2eBSZOAGzcAi0X/fhARecqXGrtfRuxVVcaN2OvrgWPHgL/9rfl7Y8YAXbsCly4BQ4bo\n3xciIiME/A1K+flA376ivZbExQGnT+vfD4D1Q3eMhcRYSIyFGj4n9sWLFyMiIgKJiYkeH+OqsdfW\nAvoWgoBz54DRo1t/Pz5eXFglIgoUPif2RYsWITMzs0PHVFYCN90EdOkiHpGnp/YSu5EjdpvNZkxD\nnQBjITEWEmOhhs+JfcqUKejbt2+HjqmsBLp3B8LC9C/HmCmxExEZwW/L9nbrZp7EfuqU/iUhgPVD\nd4yFxFhIjIUahsyKSU9PR3R0NAAgPDwceXlWdO9uQ1gYkJXlQJ8+8lcw1zdW1fbx4w7cuAEALb9/\n+rQDtbXAlSs2RESob5/bLW+7mKU//tx2Op2m6o8/t51Op6n6Y+S2w+HApu/vmnTlS28pmceem5uL\n1NRUfP75580baGEuZkoKsHQpsHgx4HQCgwb52oOWVVcDvXsD334rpjW2ZsoU4OmngWnT9OkHEVFH\ndcq1YoyosX/1FRAV1XZSBzgzhogCi8+JfcGCBZg8eTLOnTuHqKgobNy4sd1jjErs7dXXXYy6gNq0\nDAGI2v7SpcDevfq3byYtxSJYMRYSY6GGzzX2LVu2dPgYMyb2d9/Vrx9teekl4I9/FAl+6lT/9IGI\nAktAz4ox24jddcHEJTsbWLECeOON4BuxN41FMGMsJMZCjYCusXua2AcPFhdar17Vry9NlZQAf//3\nwIsvigXKLl4EvvnGuPaJKHAxsUOs7BgXB5w5o19fgMb1w4cfBmbMEMk9NBSYPBn4+GN92zcT1lIl\nxkJiLNTwa2IPDdUvsX/3HVBcDAwd6tn+Rs6MycsT68GvXStfmzo1+MoxRKSPgB2xnz8PjBoFhHj4\nCY2os7vqh1lZgM0mYuASbImdtVSJsZAYCzUMT+yaZszFU0/LMC7GLt8rEru78ePFg73FXbJERN4z\nPLHX1oqadmho8CV2V/0wK6v5Xa5hYeJpTp98om8fzIK1VImxkBgLNQxP7K4yDCCSWW2tPu2cPdux\nxB4ZCZSVAdev69Mfl9xcoKICiI1t/l6wlWOISB9+T+xmGbFbLIDVCuzZo09/AFE/dNXXW3rGajAl\ndtZSJcZCYizUCMjErmkdT+wA8NRTwG9+I0bUemmpDOMycSLwxRdi0TIiIm8ZntirqvRP7FeuiBp+\n//4dO85uF6N292mIKmVlORpG7C3p1k1cRP30U33aNxPWUiXGQmIs1PDLiL1bN/F3vRJ7YaHn89eb\nevZZYN06cSeoaoWF4lGAbf0mMXWqmDVDROStgCzFlJUBPXt6d2x0NPDII8C//IvSLgEAKipsmDat\n5fq6y9y5wOuvA//+72KZg0DFWqrEWEiMhRoBmdjLy4Gbb/b++F//Gjh8WNTDVWqrvu5y++3AsWPA\nZ5+JssyRI2r7QESBLyATe1mZ
       "text": [
        "<matplotlib.figure.Figure at 0x109b1fd50>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For comparison, here's a similar test with biquad:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Adapted directly from Audio EQ Cookbook, by Robert Bristow-Johnson\n",
      "def rbj_lpf(f, res):\n",
      "    w0 = 2 * pi * f\n",
      "    Q = 1/(2 - 2 * res)\n",
      "    cw0 = cos(w0)\n",
      "    b1 = 1 - cw0\n",
      "    b0 = .5 * b1\n",
      "    b2 = b0\n",
      "    alpha = sin(w0)/(2*Q)\n",
      "    a0 = 1 + alpha\n",
      "    a1 = -2 * cw0\n",
      "    a2 = 1 - alpha\n",
      "    ra0 = 1/a0\n",
      "    return b0 * ra0, b1 * ra0, b2 * ra0, a1 * ra0, a2 * ra0\n",
      "\n",
      "def biquad_f32(coefs, inp):\n",
      "    b0, b1, b2, a1, a2 = map(float32, coefs)\n",
      "    out = zeros(len(inp), dtype = float32)\n",
      "    x1 = float32(0)\n",
      "    x2 = float32(0)\n",
      "    y1 = float32(0)\n",
      "    y2 = float32(0)\n",
      "    for i in range(len(inp)):\n",
      "        x = float32(inp[i])\n",
      "        y = b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2\n",
      "        out[i] = y\n",
      "        x2, x1 = x1, x\n",
      "        y2, y1 = y1, y\n",
      "    return out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll design the same filter in both biquad and SVF forms, then check that the impulse response is the same."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coefs = rbj_lpf(.1, .75)\n",
      "plot(biquad(coefs, impulse))\n",
      "filter = mk4mat(svf_lp(.1, .75))\n",
      "plot(eval4(filter, impulse))\n",
      "grid(True)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FNXdP/DPbHZzzyabOySBQBJy4RY0SNWKWIoxVCIq\nldjnsaDUUq2tWlvQ2seiVUPaWtRi/eXRSmn1AWprgZfKqhSjVo1oMYJyC5Al90BuJCGX3ezO7w+a\nHEIuJDvD7jD7eb9evF47u2dnTj7AfnPO2ZmRZFmWQUREPsfg7Q4QEZF3sAAQEfkoFgAiIh/FAkBE\n5KNYAIiIfBQLABGRj1JcAKxWKzIyMpCWloaioqJBr5eUlCA8PByzZs3CrFmz8Pjjjys9JBERqcCo\n5M1OpxP33HMPdu7ciYSEBMyePRv5+fnIzMwc0O7qq6/G9u3bFXWUiIjUpWgEsHv3bqSmpiI5ORkm\nkwkFBQXYtm3boHY814yISHsUFYCamhokJSX1bycmJqKmpmZAG0mS8NFHH2HmzJlYuHAh9u/fr+SQ\nRESkEkVTQJIknbfNJZdcgqqqKgQHB2PHjh1YvHgxDh8+rOSwRESkAkUFICEhAVVVVf3bVVVVSExM\nHNAmLCys/3FeXh7uvvtuNDc3IzIyctC+amtrlXSHiMinpKSk4MiRI+7vQFbA4XDIkydPlisqKuSe\nnh555syZ8v79+we0qa+vl10ulyzLsvzJJ5/IEydOHHJfCruiG7/85S+93QXNYBYCsxCYhaD0c1PR\nCMBoNGL9+vXIzc2F0+nEihUrkJmZieLiYgDAypUr8be//Q3PP/88jEYjgoODsXnzZiWH1D2bzebt\nLmgGsxCYhcAs1KOoAABnpnXy8vIGPLdy5cr+xz/84Q/xwx/+UOlhiIhIZTwTWGOWL1/u7S5oBrMQ\nmIXALNQj/WceyeskSeL5AkREY6D0c5MjAI0pKSnxdhc0g1kIzEJgFuphASAi8lGcAiIiukhxCoiI\niNzCAqAxnN8UmIXALARmoR4WACIiH8U1ACKiixTXAIiIyC0sABrD+U2BWQjMQmAW6mEBICLyUVwD\nICK6SHENgIiI3MICoDGc3xSYhcAsBGahHhYAIiIfxTUAIqKLFNcAiIjILSwAGsP5TYFZCMxCYBbq\nYQEgIvJRml4D+N8dH+OXb69F3bptXuoVEZF26XoNoOz4MdSH7UB9c4e3u0JEpDuaLgD1bY2AnwPP\nvfGut7viMZzfFJiFwCwEZqEexQXAarUiIyMDaWlpKCoqGrbdp59+CqPRiNdee23U+2483QSpOwLb\nvnpLaTeJiOgcitYAnE4n0tPTsXPnTiQkJGD27NnYtGkTMjMzB7VbsGABgoODcfvtt+Pmm28e3JEh\n5rKmrb4bTrkXR53vwv5UubvdJCLSJa+uAezevRupqalITk6GyWRCQUEBtm0bvGD7+9//HkuWLEFM\nTMyY9t/maMI3U66B09iOki+OKekqERGdQ1EBqKmpQVJSUv92YmIiampqBrXZtm0b7rrrLgBnKtZo\ndbgakRQVg2TntSje6RvTQJzfFJiFwCwEZqEeo5I3j+bD/L777sPatWv7hyojDVeWL1+O5ORkAEBE\nRAQ6qo8jKToKuam52LLzeZSUZGLevHkAxD8CvW330Up/vLldVlamqf54c7usrExT/eG2d7b7Htts\nNqhB0RpAaWkp1qxZA6vVCgAoLCyEwWDA6tWr+9tMnjy5/0O/sbERwcHBeOGFF5Cfnz+wI0PMZRl/\nloQP7/gIoUEBmFY8Bad/eRLBgSZ3u0tEpCteXQPIyclBeXk5bDYb7HY7tmzZMuiD/dixY6ioqEBF\nRQWWLFmC559/flCb4TgDGpEyPgpTk2MR1J2Cl94pVdJdIiI6i6ICYDQasX79euTm5iIrKwtLly5F\nZmYmiouLUVxcrKhjjac6AUiIDg8GAGSH5mLzp/pfBzh7qOfrmIXALARmoR5FawAAkJeXh7y8vAHP\nrVy5csi2GzZsGPV+y2sa4dcT1b+9NCcXPy/5KYDH3eonERENpNlrAW0q+Rwrtt2BznWfAwA6uuwI\nezwSNT+pw/ioMG91k4hIM3R7LaDjJxsRKIsRQGiQP4w9sThcfdKLvSIi0g/NFoDaliaEGqIHPGdy\nWlB5stlLPfIMzm8KzEJgFgKzUI9mC0DdqUaEm6IGPBcoW1Db0uKlHhER6YtmC0BjZxMigwaOAIIN\nFjSc0ncB6Dvxg5jF2ZiFwCzUo9kC0NTViOjggSOAMKMFJ9r1XQCIiDxFswWgzdGEePPAEUC4fySa\nTuu7AHB+U2AWArMQmIV6NFsA2p2NGG8ZOAKwBFnQ3KXvRWAiIk/RbAHokpowMXrgCCAqxII2u75H\nAJzfFJiFwCwEZqEezRYAu7ERyXEDRwCxYRZ0OPVdAIiIPEWzBcDp34TU8QNHAOMiItHp0ncB4Pym\nwCwEZiEwC/VosgC0dnQDBgdiI0IGPD/eYkG3xDUAIiI1aLIAHKltgqEnCgbDwBvOJEVbYPfT9wiA\n85sCsxCYhcAs1KPJAnCsvhH+vdGDnp8YZ4HTX98FgIjIUzRZAKoamxDoihr0fEK0GTCdht3h9EKv\nPIPzmwKzEJiFwCzUo8kCUN3ciBDD4BGA0c8AqSccxxtavdArIiJ90WQBqD/VBLNx8AgAAPwcFtga\n9LsQzPlNgVkIzEJgFurRZAE4cboRkYGDRwAA4O+yoLqJ6wBEREppsgA0dzUNuhBcH71fEprzmwKz\nEJiFwCzUo8kCcMo++EJwfUJ84JLQRESeoMkC0O5qxLiIoUcAYaZInGjnGoAvYBYCsxCYhXo0WQC6\n0IQJ0UOPAML9LWjq5AiAiEgpTRaAHr9GTIwdegRgCbKgpVu/BYDzmwKzEJiFwCzUo7gAWK1WZGRk\nIC0tDUVFRYNe37ZtG2bOnIlZs2bh0ksvxa5du867z94hLgTXJ9oHLglNROQJkizLsrtvdjqdSE9P\nx86dO5GQkIDZs2dj06ZNyMzM7G9z+vRphIScuajbvn37cOONN+LIkSODOyJJkGUZHV12hBWGwLnG\nPuhaQACwasNreHnvy6hd95q73SYi0oW+z013KRoB7N69G6mpqUhOTobJZEJBQQG2bds2oE3fhz8A\ndHR0IHqYuf0+x+qah7wQXJ/4CAs6Zf0uAhMReYqiAlBTU4OkpKT+7cTERNTU1Axqt3XrVmRmZiIv\nLw/PPvvsiPs8WtcIk2Po+X8ASIi0oEfS7xQQ5zcFZiEwC4FZqMeo5M2SNPRv6edavHgxFi9ejA8+\n+AC33XYbDh06NGS75cuXo74bgK0DTz/9NLKzs/u/8tX3l54YMwkOY0v/9rmvX+zbfbTSH29ul5WV\naao/3twuKyvTVH+47Z3tvsc2mw1qULQGUFpaijVr1sBqtQIACgsLYTAYsHr16mHfk5KSgt27dyMq\nauBv+X1zWT976e94Zd8rw87xV59sQ9LTCZCfaHe320REuuDVNYCcnByUl5fDZrPBbrdjy5YtyM/P\nH9Dm6NGj/R3cs2cPAAz68D9bXVsTwk3DrxOMjwoD/LrR2e1Q0nUiIp+nqAAYjUasX78eubm5yMrK\nwtKlS5GZmYni4mIUFxcDAP7+979j+vTpmDVrFu69915s3rx5xH2eaG9ERMDwBcJgkCDZI1BRr891\ngLOHer6OWQjMQmAW6lG0BgAAeXl5yMvLG/DcypUr+x+vWrUKq1atGvX+mruaEB82bsQ2RocFlSdb\nMDU5dmydJSKifpo7E7jV3oTY0OFHAMCZS0LX6PSS0H2LPsQszsYsBGahHs0VgE5nG6JCzCO2CUKk\nri8JTUTkCZorAD1yO6LCRi4AIQYL6k/p82Qwzm8KzEJgFgKzUI/2CoDUhuiwsBHbhBktONnOEQAR\nkRKaKwAOqR2xESOPAMIDLGju0mcB4PymwCwEZiEwC/VorgA4/doRFz7yCMASZEGrji8JTUTkCdor\nAKY2xEeOXABiQiN1e0lozm8K
       "text": [
        "<matplotlib.figure.Figure at 0x10949fa10>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And we can see that the error is much less with the matrix state variable approach."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(biquad(coefs, impulse) - biquad_f32(coefs, impulse))\n",
      "plot(eval4(filter, impulse) - eval4_f32(filter, impulse))\n",
      "grid(True)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VdW1wPHfzUAIYxIgCSSRBBIyMQRFqSgYxIBMQREF\ntCjKUypFSp+vH2trH1qVobaKilq1Dqh9TCoFESJQCCpTyjwTCAQySwgBwpBxvz+ONychN8kdk3tz\n1/fzyceee849Z2dVz8pe+5y9DUophRBCCLfj0dwNEEII0TwkAQghhJuSBCCEEG5KEoAQQrgpSQBC\nCOGmJAEIIYSbcpoE8MQTTxAUFESfPn1sPtfmzZvp379/9Y+vry+rV6+2QyuFEKLlMDjLewA//PAD\n7dq149FHH+XgwYN2O++FCxeIjIwkJyeH1q1b2+28Qgjh6pymBzB48GD8/f1rfZaRkcHIkSMZMGAA\nQ4YM4fjx4xafd8WKFYwaNUpu/kIIcQOnSQCmPPXUU7z99tvs2rWL1157jRkzZlh8jqVLlzJ58mQH\ntE4IIVybV3M3oD4lJSVs376dBx98sPqzsrIyAL7++mvmzJlT5zuhoaGsW7euejsvL49Dhw4xYsQI\nxzdYCCFcjNMmgKqqKvz8/Ni7d2+dfePHj2f8+PGNnmP58uWMHz8eT09PRzRRCCFcms0lIHOe3pk1\naxZRUVH069fP5A3dlA4dOhAREcGXX34JgFKKAwcOWNS2JUuWSPlHCCHqYXMCePzxx0lJSal3/9q1\nazl58iQnTpzggw8+4OmnnzZ53OTJkxk0aBDHjx8nLCyMTz75hH/+85989NFHJCQk0Lt3b4se5czM\nzCQnJ4e77rrL4t9JCCHcgV0eA83MzGTs2LEmH9/81a9+xdChQ5k4cSIAMTExbNmyhaCgIFsvK4QQ\nwgYOfwooJyeHsLCw6u3Q0FCys7MdfVkhhBCNaJLHQG/sZBgMhqa4rBBCiAY4/CmgkJAQsrKyqrez\ns7MJCQkxeVxubq6jmyOEEC1Kz549OXnypFXfdXgPIDk5mc8++wyAHTt24OfnZ7L+n5ubi1JKfpRi\nzpw5zd4GZ/mRWEgsJBYN/2RkZFh9f7a5BzB58mS2bNlCYWEhYWFhvPTSS5SXlwMwffp0Ro0axdq1\na4mMjKRt27Z88skntl6yxcvMzGzuJjgNiYVOYqGTWNiHzQlgyZIljR6zaNEiWy8jhBDCzpx6LiB3\nNXXq1OZugtOQWOgkFjqJhX04zXTQBoMBJ2mKEEK4DFvundIDcEKpqanN3QSnIbHQSSx0Egv7kAQg\nhBBuSkpAQgjhwqQEJIQQwmKSAJyQ1Dd1EgudxEInsbAPSQA/++orqGemaiGEaJFkDAAoK4OYGKiq\nAnnBUAjhSmQMwEYffghRUVBYCJcuNXdrhBCiabh9AigpgVdegfnzIS4ODh1q7hZJfbMmiYVOYqGT\nWNiH2yeAN9+ExETo3x/69AELlx0WQgiX5dZjAOfPQ3Q07NgBkZGwcCGcOAHvvNOkzRBCCKvJGICV\nXn8dJkzQbv6g9QBMLGsshBAtklsngN27Ydw4fbtvXy0BNHefSOqbOomFTmKhk1jYh1sngPx8CA7W\nt7t0AR8fkDXrhRDuwK3HAIKCYN8+6NpV/ywpCX77Wxg1qkmbIoQQVpExACtUVEBREQQG1v5cxgGE\nEO7CbRPAuXPQqRN4etb+3DgO0JykvqmTWOgkFjqJhX24bQK4sf5vJO8CCCHchduOAaxdC2+/DevW\n1f782jUICNCmhPD2brLmCCGEVWQMwAr19QB8feGmm+D48aZvkxBCNCVJACY090Cw1Dd1EgudxEIn\nsbAPt00AeXn1J4C+fWUcQAjR8rntGMCDD2o/Dz1Ud9/KlfDRR7BmTZM1RwghrCJjAFZw5hKQEEI0\nBbdNAA2VgHr00N4TKClp2jYZSX1TJ7HQSSx0Egv7cNsE0FAPwMNDmx4iN7dp2ySEEE3JLccASkq0\neYBKSsBgMH3MnXfCq6/CXXc1SZOEEMIqMgZgIeNf//Xd/EHrAeTlNV2bhBCiqbllAmio/m/UnAlA\n6ps6iYVOYqGTWNiHWyaA/PzaU0Cb0q2b9ACEEC2b2yYAZ+4BJCYmNs+FnZDEQiex0Eks7MMtE4Cz\nl4CEEKIpuGUCcPYegNQ3dRILncRCJ7GwD7dNAI2NAUgPQAjR0tmcAFJSUoiJiSEqKooFCxbU2Z+a\nmkrHjh3p378//fv355VXXrH1kjYzpwTUqRNcuQLXrzdNm2qS+qZOYqGTWOgkFvbhZcuXKysrmTlz\nJhs3biQkJIRbb72V5ORkYmNjax131113sXr1apsaak/mlIAMBu1lsbw8iIhomnYJIURTsqkHkJaW\nRmRkJOHh4Xh7ezNp0iRWrVpV5zgnedkYgMpKKCysuxi8Kc1VBpL6pk5ioZNY6CQW9mFTAsjJySEs\nLKx6OzQ0lJycnFrHGAwGtm3bRr9+/Rg1ahRHjhyx5ZI2KywEf3/zlnuUcQAhREtmUwnI0NBcCj+7\n+eabycrKok2bNqxbt4777ruP9PR0k8dOnTqV8PBwAPz8/EhISKiu9Rkzvq3bfn6JBAebd3xVFeTl\n2ff65mwnJiY26fVk23W2jZylPc21bfzMWdrTlNupqal8+umnANX3S2vZNBncjh07ePHFF0lJSQFg\n3rx5eHh48Nxzz9X7nYiICHbv3k1AQEDthjTRZHApKfD667B+fePHvvyyNgj86qsOb5YQQlil2SaD\nGzBgACdOnCAzM5OysjKWLVtGcnJyrWMKCgqqG5eWloZSqs7NvymZ8wiokYwBND+JhU5ioZNY2IdN\nJSAvLy8WLVrEiBEjqKysZNq0acTGxvL+++8DMH36dL788kvee+89vLy8aNOmDUuXLrVLw61lziOg\nRjIGIIRoydxuPYDf/AbCw+G3v2382D174IknYN8+hzdLCCGsIusBWMCcdwCMpAcghGjJ3C4B5OWZ\nPwYQGAhFRVBe7tg23UjqmzqJhU5ioZNY2IfbJQBLegCentC5MxQUOLZNQgjRHNxuDKBDBzh7Fvz8\nzDv+5pvh/ffh1lsd2y4hhLCGjAGYqbQUrl2Djh3N/46MAwghWiq3SgAXL2o3fzNeYK7WHAlA6ps6\niYVOYqGTWNiHWyYAS8jawEKIlkoSQCOaowdQc74Tdyex0EksdBIL+3CrBFBcbP7gr5GMAQghWiq3\nSgCu0gOQ+qZOYqGTWOgkFvYhCaAR0gMQQrRUbvUewBtvQGYmvPmm+d8pK4N27bRpoT3cKl0KIVyB\nvAdgJmt6AK1aaS+PnTvnmDbd6Fe/gs8/B+dIy0KIlkwSgBmaqgx06BD861/w5z+nMnw4ZGQ4/prO\nTmq9OomFTmJhH5IAzNBUCeCTT7Tpp//+d7j3Xhg4EFaudPx1hRDuyaYFYVzNxYuWPwYKTfMyWFkZ\nfPEF/PgjREUlMmwY9OoFCxbA/fc79trOTJ731kksdBIL+3CrHkBxsXU9gG7dIDfX/u2pac0aiImB\nqCj9s+HD4fBheQpJCOEYbpUArC0BdesGOTn2b09NH3+slX9Ar2/6+MCoUdq4gLuSWq9OYqGTWNiH\nJAAzhIQ4tgeQmwtbt8KECXX3PfAAfP21464thHBfbvUeQGAgHDhg/oIwRjt3wq9/Dbt2OaZd8+dr\nT/x8+GHdfVeuaD2QU6egUyfHXF8I4bpa7HsASim2nt1qt/M5Yw9AKa38M22a6f1t28I998A33zjm\n+kII9+XUCeDUhVOM/r/RdjnX9evazbZ1a8u/GxysvQhWUWGXptSyb5/WroED9c9urG+OH+++ZSCp\n9eokFjqJhX04dQLIvpTNxdKLVFZV2nwuaxaDMfLyctzawDt3wuDBDbdrzBhITYXLl+1/fSGE+3L6\nBABQfL3Y5nNZMxV0TSEhjnkSaNcuGDCg9mc3PuPcsSPceSesXWv/6zs7ed5bJ7HQSSzswyUSQNG1\nIpvPZW3938hR7wKYSgCmPPAA
       "text": [
        "<matplotlib.figure.Figure at 0x1094a7450>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The difference is even more dramatic at lower frequencies, where the error of direct form biquad evaluation\n",
      "is large."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "impulse = concatenate([[1], zeros(499)])\n",
      "coefs = rbj_lpf(0.01, .75)\n",
      "filter = mk4mat(svf_lp(0.01, .75))\n",
      "plot(biquad(coefs, impulse) - biquad_f32(coefs, impulse))\n",
      "plot((eval4(filter, impulse) - eval4_f32(filter, impulse)))\n",
      "grid(True)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXZ+PHvTBa2QBKWBLJAICSEQEgCSIQKBhEQkMhi\nX0BbpWJFcSlqFe2vFXyrIK3Lq0X7IlZQa8GqRaJCXkAcRCBGkSCyBkkkCSFsIRCWrOf3x2MyQhaS\nmXNmvT/Xlas9M5NzHm7h3HPuZzNpmqYhhBDC65id3QAhhBDOIQlACCG8lCQAIYTwUpIAhBDCS0kC\nEEIILyUJQAghvJTLJIC77rqL0NBQEhIS7D7X559/TnJyct1PmzZtSE9P16GVQgjhOUyuMg9gy5Yt\nBAQEcMcdd7B7927dzltSUkLv3r0pLCykdevWup1XCCHcncs8AQwfPpzg4ODLXvvhhx8YN24cgwcP\nZsSIERw4cKDF533//fcZP3683PyFEOIKLpMAGnLPPffwt7/9jW+++Ya//vWvzJkzp8XnWLVqFTNm\nzDCgdUII4d58nd2AxpSVlbF9+3Z++ctf1r1WUVEBwH/+8x/mz59f73ciIiJYt25d3XFRURHff/89\nY8eONb7BQgjhZlw2AdTU1BAUFMTOnTvrvTdlyhSmTJly1XP8+9//ZsqUKfj4+BjRRCGEcGt2lYDy\n8/MZOXIk/fr1o3///rzyyisNfu6hhx4iJiaGxMTEBm/oDenQoQM9e/bkgw8+AEDTNL777rsWtW/l\nypVS/hFCiEbYlQD8/Px46aWX2LNnD5mZmbz66qvs27fvss+sXbuWQ4cOkZOTw+uvv859993X4Llm\nzJjBsGHDOHDgAJGRkSxfvpx3332Xf/zjHyQlJdG/f/8WDeXMy8ujsLCQ66+/3p4/ohBCeCxdh4FO\nmjSJBx98kFGjRtW9du+99zJy5EimTZsGQFxcHJs3byY0NFSvywohhLCBbqOA8vLy2LlzJykpKZe9\nXlhYSGRkZN1xREQEBQUFel1WCCGEjXRJAGVlZdx66628/PLLBAQE1Hv/yocMk8mkx2WFEELYwe5R\nQJWVlUydOpVf/epXTJo0qd774eHh5Ofn1x0XFBQQHh7e4OeOHj1qb3OEEMKrREdHc+jQIZt+164n\nAE3TmDVrFvHx8cydO7fBz6SlpfH2228DkJmZSVBQUIP1/6NHj6JpmvxoGvPnz3d6G1zlR2IhsZBY\nNP3zww8/2HwPt+sJYOvWrfzzn/9kwIABJCcnA7Bw4UKOHDkCwOzZsxk/fjxr166ld+/etGvXjuXL\nl9tzSa+Ql5fn7Ca4DImFlcTCSmKhD7sSwHXXXUdNTc1VP7dkyRJ7LiOEEMIALr0WkLeaOXOms5vg\nMiQWVhILK4mFPlxmOWiTyYSLNEUIIdyGPfdOeQJwQRaLxdlNcBkSCyuJhZXEQh+SAIQQwktJCUgI\nIdyYlICEEEK0mCQAFyT1TSuJhZXEwkpioQ9JAEII4aWkD0AIIdyY9AEIIYRoMUkALkjqm1YSCyuJ\nhZXEQh+SAIQQwktJH4AQQrgx6QMQQgjRYpIAXJDUN60kFlYSCyuJhT4kAQghhJeSPgAhhHBj0gcg\nhBCixSQBuCCpb1pJLKwkFlYSC33YnQDuuusuQkNDSUhIaPB9i8VCYGAgycnJJCcn88wzz9h7SSGE\nEDqwuw9gy5YtBAQEcMcdd7B79+5671ssFl588UXS09Obboj0AQghRIs5tQ9g+PDhBAcHN/kZubEL\nIYTrMbwPwGQysW3bNhITExk/fjx79+41+pJ2OXIELl1ybhukvmklsbCSWFhJLPTha/QFBg4cSH5+\nPm3btmXdunVMmjSJgwcPNvjZmTNnEhUVBUBQUBBJSUmkpqYC1v/gRh2vX29hwQI4eDCVVq3gySct\n9O9v3PXkuHnHtVylPc48zs7Odqn2OPM4OzvbpdrjyGOLxcKKFSsA6u6XttJlHkBeXh4TJ05ssA/g\nSj179mTHjh107Njx8oY4uQ/g97+HQ4fggw8gIwN++1vYuRO6dnVak4QQ4qpceh5AcXFxXeOysrLQ\nNK3ezd/ZduyAd96BZcvA1xduvhnuvBPmzXN2y4QQwjh2J4AZM2YwbNgwDhw4QGRkJG+++SZLly5l\n6dKlAHzwwQckJCSQlJTE3LlzWbVqld2N1pOmwe9+B4sWQZcu1teffBI++UT1CTjaleUPbyaxsJJY\nWEks9GF3H8DKlSubfP/+++/n/vvvt/cyhtmyBU6cUN/4fy4wEG6/HZYvh/nzndM2IYQwktevBTRn\nDvTo0XC5Z/t2uPtu2LPH4c0SQohmsefe6dUJoKYGwsNh82aIjW38/S+/hOhohzZNCCGaxaU7gV1Z\nZiZ07tzwzR/AbIaRI1WCcCSpb1pJLKwkFlYSC314dQL4z39gypSmP5OaCvJ3TQjhiby2BKRpqqyz\nejUkJjb+uYMH4cYb4ccfwWRyWPOEEKJZpARkg1271A19wICmPxcTA1VVkJvrmHYJIYSjeG0CqC3/\nXO1bvcnk+DKQ1DetJBZWEgsriYU+vDIBaBp8+OHV6/+1pB9ACOGJvLIP4PXX4dVX1Vo/5makQOkH\nEEK4KnvunYavBupqnn5aze795JPm3fzh8n6AXr2MbZ8QQjiKR5eAiorUGP+0NDWpa+dOteBbVhb0\n79/88zi6H0Dqm1YSCyuJhZXEQh8enQD+3/9TN/8TJ1TJ57//Gx56CEJCWn4u6QcQQngaj+0DOHsW\nundXa/yXlEBKiirfbNkCbdq0/HzSDyCEcEXSB9CAtWvhuuvUUg+dO8Px4+rG7eNj2/mkH0AI4Wk8\ntgS0ezdcc4312NfX9ps/OLYfQOqbVhILK4mFlcRCHx6bAPbvh7599T2n9AMIITyJx/YBxMfDqlVX\nX+qhJaQfQAjhamQtoCtUVcHhw6pur6eYGLh0CY4e1fe8QgjhDB6ZAI4eVR2/toz2aYrJBP36wd69\n+p73SlLftJJYWEksrCQW+rA7Adx1112EhoaSkJDQ6GceeughYmJiSExMZOfOnfZe8qqOHFFDQI0Q\nH298AhBCCEewOwH85je/ISMjo9H3165dy6FDh8jJyeH111/nvvvus/eSV2VkAnDEE0BqaqqxF2gG\niwUeewxeeknNqXAWV4iFq5BYWEks9GF3Ahg+fDjBwcGNvp+ens6dd94JQEpKCmfOnKG4uNjeyzZJ\nngDss2ABzJwJQUFq2YzkZNWnIoTwLIb3ARQWFhIZGVl3HBERQUFBgaHXNDoB7NmjlpQ2ijPrm8uX\nq9FTWVlqKY2VK+GRR2D8eDh3zvHtkVqvlcTCSmKhD4fMBL5yiJKpkTGUM2fOJCoqCoCgoCCSkpLq\nHvVq/4M35/jIEQgLs2CxNO/zLTm+/vpUzGZYvdpCx476n9+Zx2fOwOOPp2KxwN69FvbuVe/ffz+k\np1uYNQv+/W/Htq+WK8TH2cfZ2dku1R5nHmdnZ7tUexx5bLFYWLFiBUDd/dJWuswDyMvLY+LEieze\nvbvee/feey+pqalMnz4dgLi4ODZv3kxoaOjlDdFxHsCgQfC//3v5TGA9jRihyiQ33GDM+Z3l0Ueh\nvByWLKn/XnGx6v/Yvl3/4bVCCNu59DyAtLQ03n77bQAyMzMJCgqqd/PX27Fj0LWrcefv108tNeFJ\niotV+efJJxt+PzQUfv97+OMfHdsuIYRx7E4AM2bMYNiwYRw4cIDIyEjefPNNli5dytKlSwEYP348\nvXr1onfv3syePZvXXnvN7kY3pbpaLfxmZI4ZMgS++sq4819Z/nCEF16A22+H8PDGP3P//bBxo+pj\ncRRnxMJVSSysJBb6sLsPYOXKlVf9zJKGagoGOXUKAgPB39+4a1x7Lfz5z8ad39EqK2HFClXeaUr7\n9nDHHfDaa/Dccw5pmhDCQB63FtCuXfCrXxlboqmpgY4dIScHunQx7jqOsm6d2iznagkA1P4KQ4dC\nfj60bm1824QQTXPpPgBHM7r+
       "text": [
        "<matplotlib.figure.Figure at 0x10949fc50>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By visualizing the state of the filter, it's possible to see intuitively where the difference in numerical robustness\n",
      "comes from. In the SVF, where the state corresponds to physical quantities, the norm of the state vector is close to\n",
      "the system's energy."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def eval_showstate(filt, inp):\n",
      "    A, B, C = filt\n",
      "    result = []\n",
      "    y = zeros(2)\n",
      "    for x in inp:\n",
      "        result.append(y)\n",
      "        y = B * x + dot(A, y)\n",
      "    return result\n",
      "\n",
      "impulse = concatenate([[1], zeros(49)])\n",
      "filt = svf_lp(0.02, 0.5)\n",
      "xs, ys = zip(*eval_showstate(filt, impulse))\n",
      "axes().set_aspect(1)\n",
      "grid(True)\n",
      "plot(xs, ys, '-o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "[<matplotlib.lines.Line2D at 0x1094aa250>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAADdCAYAAAChd+XWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXtcVVX6/98oJokZZYoplNMBBW9g4ZdsCimDo1JM9e1X\ndFEsbegmWd/pJjpSmuI09U2jGZ35zkjOhZxb6QAiOXmgJpFKqWlsEg1KjkJl4TVEDuv3x5YTBw5y\nOftw9obn/Xrtl6191t77wybWc9bzrGc9fkophSAIgiC0op+vBQiCIAjGRAyEIAiC4BYxEIIgCIJb\nxEAIgiAIbhEDIQiCILhFDIQgCILgFo8NRGFhIREREYSHh7Nq1Sq3fdLT0wkPDycqKordu3cD8Omn\nnzJ58mTncf7557NmzRpP5QiCIAg64edJHoTD4WDs2LFs27aNUaNGMWXKFHJzc4mMjHT2KSgoIDs7\nm4KCAnbu3MkjjzxCaWmpy32ampoYNWoUZWVlhIaGdv+nEQRBEHTDoxlEWVkZYWFhjB49mgEDBpCS\nksKmTZtc+mzevJnU1FQAYmNjqauro7a21qXPtm3bsFgsYhwEQRAMhEcGwm63uwzqISEh2O32DvtU\nV1e79Hnttde48847PZEiCIIg6Iy/Jxf7+fl1ql9rL1bL6xoaGvj73//ebvzioosu4vDhw90XKQiC\n0AexWCzs27fPo3t4NIMYNWoUBw4ccLYPHDhASEjIWftUV1czatQoZ3vLli1cccUVDBs2zO0zDh8+\njFLK8EdqaqrXn5GXV0xiYgbTpi0lMTGDvLziLl8/ePBEQDkPi2VRl+/TG96l6DTeITr1Pfbv3+/J\n8A54OIOIiYmhoqKCqqoqRo4cycaNG8nNzXXpk5ycTHZ2NikpKZSWlhIUFERwcLDz89zcXO644w5P\nZBiC0aNHe/X++fklPPLIVvbvf855bv/+DACSkuLcXnP6NBw//v2xbFkRx4/f4tJn//7neOmlJe3e\noyNNa9YUceqUPwMHNpKentit+7TG2+9SL0SnvohO4+GRgfD39yc7Oxur1YrD4WDevHlERkaybt06\nANLS0pg1axYFBQWEhYURGBjI+vXrndefOHGCbdu28etf/9qzn6KXoxSsWlXkYhxAG9zvuWcJEybE\nuRiC48fhxAlwOGDwYO0IDAS73f2ve9u2/gQGwkUXtX8MHera3rmzhMcf75rBEgTBXHhkIABmzpzJ\nzJkzXc6lpaW5tLOzs91eGxgYyNdff+2pBEMQFBTUqX4dfes+fRo++QTKy7Vj927t3xMn3P+qhg/v\nT0bG94ag2RgMHgwDB0LLMJHV2khR0UVt7mG1OvjLX+Drr12Pw4e1fz/6qO35Q4eKgLYGa9myJUyb\nFsfgwZ16HW7p7Lv0NaJTX0Sn8fDYQAga0dHRHfZx5yb6+OMMbrgBGhvjKC/XjMOll0J0tHY8/TRE\nRcGcOY0UFbW9Z0iIg+nTO6cxPT2Rjz/+FQcPfn/OYlnEggUznMals7PnadP8KSlpe/7f/+7P8OFw\nySXf/wyTJ2v/tvAsnpXOvEsjIDr1RXQaD48S5XoCPz8/DC6x0yQkLGbbtuVtzoeELCEjYxmTJ8OE\nCdoMoDXujIvFsojVq2d0yaWTn1/Cyy+/SX19fwICHCxYkNAtl5DVupiiorY/i9W6hL//fRn/+Y/r\nDKi8HAICXA1GdDRYLNCv1VIJb8U2BKEvocfYKTMIL9PUBO++C7//Pdhs7l+3xdKf++8/+32aB8iX\nX17SYnDvmnFovo8eg216eiL792e0MVgLFsxgwACYOFE7Zs/WPlMKvvjie2Pxhz/A44/DN9/ApEnf\nG41jx0p4+eWtfPaZxDYEwecog2NkiXl5xSoxMUNNm7ZUxcTcrfLyip2f7dmj1KJFSo0erdT48Uqt\nXKnUNddkKG2odD2s1sU9pnn79u263Ssvr1hZrYvVtGlLldW62OXn7yyHDyv11ltKvfiiUrNnKzV4\ncIaC7T59R51Fz3fpTUSnvphFpx5jp8wguklbl4+Nhx7ayuuvw65dcdTWwp13whtvaN+Q/fxg4sRE\nDh50/63bjOgxG7nwQrj2Wu0AiI/3p7i4bb+dO/vz0kswaxaMGePRIwVB6CQSg+gm7fngR45cwu9+\nt4xp06B//7bX6RUD6K20914nT17CFVcso6AABg2CpCTtiIvTVmsJguCKxCB8SH29+1cXHt6f665r\n/zq9YgC9lfZiG8uWzSApSXM4ffgh5OdDZiZ8/LE2+0hK0mYXLZL0BUHwECkY1A2KiuCDDxpbnbUB\nEBDg6HE9XcFms/lawllJSopj3rzhWK1LmDYtE6t1ictKLT8/LZidkQH//Cfs3w//7//B9u2aK6/5\ns3ff1RIFvYnR32UzolNfzKJTD2QG0QX+8x/4yU/g00/h0UcTyc3tPfEEIzF1ahRPPx3fqb4XXQR3\n3aUdjY2wc6c2u3jgAbDbwWrVZhczZmjxDpBltILQWSQG0Qp3g8eVV8bxzDOQmwtPPQUPP6z5vSWe\nYGyqq6GgQDMYNpu27Payy0rYvn0r1dUtDXsGq1db5Xcn9Cr0GDvFQLTAXTLaRRdl0NBg5e6748jM\nhHY2nRUMTn09lJTAffct5osv3Cf4FRYu84EyQfAOeoydEoNowZo1bTfE+/rr55g48U1eeeXsxsEs\nfkkz6PSGxoAASEyEH/zAvVf1s8/6c+RI1+5phncJolNvzKJTDzw2EIWFhURERBAeHt5u0Z/09HTC\nw8OJiopi9+7dzvN1dXXceuutREZGMm7cuDa1qnuaU6fcDx7+/m7WqwqmZODA1osLNI4edTB6NMyd\nC++8o62WEoQ+jydZdo2NjcpisajKykrV0NCgoqKi1J49e1z65Ofnq5kzZyqllCotLVWxsbHOz+bM\nmaN+85vfKKWUOn36tKqrq2vzDA8ldolrr/V9prPgXfLyipXFssjl92uxPK3y8orVl18q9cILSkVG\nKjV2rFI/+5lSNTW+ViwI3UOPsdOjVUxlZWWEhYU5C2ikpKSwadMmIiMjnX02b95MamoqALGxsdTV\n1VFbW0tAQABvv/02r776KqDVljj//PM9keMR+/bB3r2JnH9+BkeOyMqk3kpHe1o99hg8+ijs2AH/\n938QEaHlWcyfr62Icpf8KAi9FY9cTHa7ndDQUGc7JCQEu93eYZ/q6moqKysZNmwY99xzD5dffjn3\n3XcfJ0+e9EROtykuhquvhsWL4/jDH6ztrsE/G2bxS5pBp7c1JiXFUVi4DJstk8LCZW1+v35+cNVV\n8Nvfwuefw8yZ8Mwz2lboP/0pVFb2jE69EJ36YhadeuDRDMKvZTWas6BaOXT9/PxobGxk165dZGdn\nM2XKFBYuXEhWVhbPPvtsm+vnzp3rnKUEBQURHR1NfHw88P0vq7vtJ5+0sW4d/OUv8Vx/PdhsTTz1\n1HSX/jabrcP7NeOpHm+3y8vLDaXHXbu8vNwwenbtshEeDjt3xvPRR/DMMzaio+G//iueqVOhocHG\nOecY6/21bhvpffaGtlHfp81mIycnB9CvLKpHy1xLS0vJzMyksLAQgJUrV9KvXz+efPJJZ5/777+f\n+Ph4UlJSAIiIiKC4uBilFFOnTqXyzNexd955h6ysLPLy8lwF6rzMtTnPob7eny++aKS+PpHt2+OI\niNDtEUIvp75e24Tx//5P2/bjrrs0F9SECb5WJgjf4/NlrjExMVRUVFBVVUVDQwMbN24kOTnZpU9y\ncjIbNmwANIMSFBREcHAwI0aMIDQ0lL179wKwbds2xo8f74mcDmnOcygqWk5JSSZVVcsJCNjK/v1u\nSqMJQjsEBEBKCmzbpmVun3eelql95ZWa0Th2TOuXn1+C1bqY+PhMrNbF5OfL/2eCyfA0yl1QUKDG\njBmjLBaLWrFihVJKqbVr16q1a9c6+zz00EPKYrGoSZMmqQ8++MB5vry8XMXExKhJkyapm2++2eur\nmBITvbdKySx7xJtBpxk0KuWq
       "text": [
        "<matplotlib.figure.Figure at 0x1094aa090>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By contrast, with biquad, the state vector is a fairly extreme transform. Thus, even a relatively small error in the\n",
      "state variables can lead to a significant delta in energy, resulting in much larger output errors. The transform\n",
      "becomes more extreme as the frequency decreases. Note that this is the biquad in matrix form, which is most similar\n",
      "to transposed direct form II. However, error behavior is broadly similar, at least in floating point, for all the\n",
      "biquad evaluation schemes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filt = from_biquad(rbj_lpf(0.02, 0.5))\n",
      "xs, ys = zip(*eval_showstate(filt, impulse))\n",
      "axes().set_aspect(1)\n",
      "grid(True)\n",
      "plot(xs, ys, '-o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10974ccd0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAATQAAAEACAYAAAA9aookAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXtcVHX+/18DeMtKJXVUBkVHbhoOtBZtu8KUDQNS5K00\ntxJt2zLLbusNJWlLhNzaVay1Wi/0rRCrNclBQMph3J+4bqm1LWZEUDBcthA0Tbn5/v1xZGSYAWHO\nmSvv5+MxDz2Hzzmv9wHOm8/n/fm8P28ZEREYhmE8AC9nG8AwDCMV7NAYhvEY2KExDOMxsENjGMZj\nYIfGMIzHwA6NYRiPQbRDy8vLQ0hICAIDA5Genm61zbJlyxAYGAiVSoXjx4+bzi9evBhyuRxhYWFi\nzWAYhhHn0Nra2vDEE08gLy8PJSUlyMrKwsmTJ83a5Obm4ttvv0VpaSnefPNNLFmyxPS1RYsWIS8v\nT4wJDMMwJkQ5tKNHj2LixIkICAhAv379MH/+fOzdu9esTU5ODhYuXAgAiIyMRGNjI2prawEA06ZN\nw7Bhw8SYwDAMY0KUQzMajfD39zcdKxQKGI3GXrdhGIaRAlEOTSaT9ahd5+yqnl7HMAzTG3zEXOzn\n54fKykrTcWVlJRQKRbdtqqqq4Ofn1yuN6upqMWYyDOMGKJVKfPvtt6LuIaqHNnXqVJSWlqKiogLN\nzc3Izs5GQkKCWZuEhAS8/fbbAIAjR45g6NChkMvlPdaorq4GEUn+AeIArAGwDgBd/j8BSASQdPn/\nwkepTMK+fUV2sWPdunV2uS/rsZ47aRERysrKxLgjACIdmo+PD7Zs2QKtVotJkyZh3rx5CA0NxRtv\nvIE33ngDADBjxgxMmDABEydOxKOPPorXX3/ddP3999+P2267Dd988w38/f2xY8cOcU/TC+TyNgBf\nA2iflW3vrH4OYD0AA4C1AFJQViZDcvLbdrGjoqLCLvdlPdZzJy2pEDXkBIC4uDjExcWZnXv00UfN\njrds2WL12qysLLHyNrNt2xrcffc6EDUDWASgfRjsBcGZ5UNwbAInTy6BTmdAfHyUw21lGKZn9NlM\ngfj4KISHKwHcAaA/gOMA/gAgDEABrjgzoad28aIcCxe+Bp3OIKkdiYmJkt6P9VjPHbWkQkZELr3B\no0wmg71M1OkMeOqpfJSVdXReaQCGAMiCZU/NgEGDXoNSORpjxlyLZctiuMfGMBIhxbveZ3togNBL\n27RJixtumAcgBcABAFoA7Yt9O/fU8nHhQja++uqvKCh4CU89lS+6x6bX60Vdz3qs5wlaUtGnHRog\nOLXMzKVQKlsAvAhABWA+gMdgHmK0HIaWlfWzyzCUYRjb6NNDzo7odAZkZBzAkSOVOHPGH8BoAEUA\nsi+3SLn8sZwwUCrXYNMmLQ8/GUYEUrzr7NA6odWuRUHBS5ePOjqvtQBe6vBv+9cLAPjghhtOIjNz\nKTs1hrERjqHZAbV6BJTKNZePoiDE1OYBKIf5MLTd2b0EIAX19dk2xdQ8OQbDeu6txzE0D+DXv1Zh\n0yYttNpkDBmSCGGiYCmAdwEswJWFuB1jagJlZeuRkXHAgdYyDNMRHnJ2g/nwsx0DBg16DxcujIIQ\nUzMnOjoFer3leYZhuoeHnHZm2bKYDsNPgeuv34PRo3+Ej8/nVq8ZOLDNEaYxDGMFdmid6Bg3aF+n\nptUmIzo6BZMnL8XPP5/Fd999iNbW5RAS2q+gVCbhySc1Nus5AtZjPVfUkgrRuZyeTnx8lGnmUqtd\ni//+t30I2j6bmYxhw37ALbeMxZNPxvIsJ8M4EY6h9QK1OgVFRSkW5zluxjDi4RiagxkwoNXqeY6b\nMYxrwA6tE93FDaxNEnh5JcHPbzS02rVQq1Og1a7t1Vo0T47BsJ5763EMzcNpj49lZCTj4kVvDBzY\nBl9fBXbuNOLSpStr0srK1pi1ZxjGMXAMTSTW16oBWm0y8vJedIJFDOOecAzNBWhqst7JvXjR28GW\nMAzDDq0TvY0biJ0o8OQYDOu5t547xtBEO7S8vDyEhIQgMDAQ6enpVtssW7YMgYGBUKlUOH78eK+u\ndXWsTRR4eydh7tzeLbB1NXQ6g80THQzjNEgEra2tpFQqqby8nJqbm0mlUlFJSYlZG51OR3FxcURE\ndOTIEYqMjOzxtZfje2JMdAj79hWRVruWoqPXkVa7lh5/vIhGjCiiadPWUHT0OoqJWUP79hU528we\ns29fESmVSQSQ6aNUJrnVMzDuhxTvuqhZzqNHj2LixIkICAgAAMyfPx979+5FaGioqU1OTg4WLlwI\nAIiMjERjYyNqa2tRXl5+1WvdhY7ZBIDQu9m9Ox+HDnU986nTGbB5cwGamnwwYECr0+oTEAEtLUBT\nE3DxovDZsKGgQ52FdvvXIyMjmWduGZdGlEMzGo3w9/c3HSsUCvzrX/+6ahuj0Yjq6uqrXusM9Ho9\n1Gq1qHts3lyAn37q2iGYF2fRA1Djm2/W4Mcfgd/+NsrMubR/pDrX0KCHTKY2O+fjAwwcCAwYIPz7\n00/STXRI8f1kPefoOfrZpECUQ5PJZD1qRy687MIeXG3mc/Nmyx5QRcV6PPZYMvz8ojBwIEyfdifT\n+dPxvK9vz9ueOAGo1VfODRgAeHfyU1ptKwoKLO3njAjG1RHl0Pz8/FBZWWk6rqyshEKh6LZNVVUV\nFAoFWlparnptO4mJiaah6dChQxEeHm76y9E+EyPVcfs5Mfc7f75jSXv95X/VANqg1+tRV1fV6Qn1\nANS49VZvpKRI+zx6vR5NTcCvfy0c19UBp051/3xq9QiUla3p0IME+vUrwO9+F+uU7yfrOUdPrVbb\n9f56vR47d+4EANP7LRoxAbiWlhaaMGEClZeXU1NT01UnBYqLi02TAj25lsg9JgU6Yy2ofsMNq2nU\nqCKqrCSKiVlj9rX2j1a71tmmm+g80fG73xXRpElEdXXOtozxVKR410XfITc3l4KCgkipVFJqaioR\nEW3dupW2bt1qarN06VJSKpU0ZcoU+vzzz7u91sJABzu0gwcPSnKfzg5h374i2riRSKkkevrp12jQ\noPsIWEfAAwQUkVK52iGziGKe7/nnicLCiH780TF6tsB67qlF5AKznAAQFxeHuLg4s3OPPvqo2fGW\nLVt6fK2n0HnmUzgH/Pe/BmzZYkRra3t5PD0GDXoPDzwwxeVnEFNSgNZW4M47gU8/FWJ3DONKcC6n\ng3H33E8iYNUqoLBQ+AwbdvVrGKYncC6nG+LuuZ8yGZCWBkRHAzExQGOjsy1imCuwQ+uEvfPXLHM/\nBT1HLYmQ4vlkMuCVV4Bf/xqIjQXOnrWvXm9gPffUkgp2aA7GWu6nj08S7rvPvXI/ZTJg0ybgppuA\nuDjg55+dbRHDcAzNKeh0BmRkHDBtEjlunAb5+VH49FNgwgRnW9c7Ll0CliwBTp4E9u8HBg92tkWM\nuyLFu84OzUV4/XUgPR345BNg4kRnW9M7Ll0CHnkE+O47QKcDrrnG2RYx7ghPCtgBZ8VEHn8cWLMG\nuP124M037bd1jz2ez8sLeOstYOxYICEBuHDBvnrdwXruqSUVXFPAhfjDH4CvvjLg8cfz0dbmXjUK\nvLyA7duBhx4CZs4E9u4VckUZxpHwkNPFcPd1aq2twAMPCDOfe/YIye8M0xN4yOmBuPs6NR8f4J13\nhMmBuXOB5mZnW8T0JdihdcLZMRF7FzN2xPP5+ADvvQf06wfccYceLS12lzTh7J+fJ+m5YwyNHZqL\n0VUx45gY91qn1q8fsGuXMAN6//1wqFNj+i4cQ3NBOq9Ti4zU4G9/i8KePcBvfuNs63pHUxMwezZw\n3XXCUNSHp6GYLuB1aH2IggIh2L57t7DjrDtx8SJwzz3A8OHA229b7pDLMABPCtgFV42JxMQIQ7h7\n7wVeesn2dWrOeL6BA4GPPhJ2y120CGizY9qqq/783FHPHWNoPABwI+64A1i+3IDVq/Nx6ZJ7rVMb\nNAjIyRH2hHvkEeDvfxfWrjGM
       "text": [
        "<matplotlib.figure.Figure at 0x10974cd10>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Modulation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The main reason to prefer one state space over another for the same transfer function is performance under modulation,\n",
      "that is when the filter parameters don't stay constant. Biquad filters are notoriously bad at this. State variable filters\n",
      "perform well under modulation and are highly regarded for these use cases.\n",
      "\n",
      "Let's test it and see how it does."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def saw(freq, n): return 1-2 * (linspace(0, freq * n, num=n, endpoint=False) % 1)\n",
      "def sweep(n, maxf = 0.5):\n",
      "    out = zeros(n)\n",
      "    ph = 0\n",
      "    for i in range(n):\n",
      "        out[i] = sin(ph)\n",
      "        f = exp(5 * (float(i) / n - 1)) * pi * 2 * maxf\n",
      "        ph += f\n",
      "    return out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def eval_mod(freq, res, inp, filtf):\n",
      "    A, B, C = filt\n",
      "    result = []\n",
      "    y = zeros(2)\n",
      "    for i in range(len(inp)):\n",
      "        x = inp[i]\n",
      "        A, B, C = filtf(freq[i], res)\n",
      "        result.append(dot(C, concatenate([[x], y])))\n",
      "        y = B * x + dot(A, y)\n",
      "    return result\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The spectrogram looks good."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 10000\n",
      "res = 0.9\n",
      "specgram(eval_mod(0.25 + 0.2 * sweep(n, 0.1), res, saw(0.05, n), svf_lp))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvcmrdUt65veLbnW7Oe33fbfJm5k4M0WpsJDBYNkFAmHP\nBBYeSnhQSBOhicBjTyT8D3iQI00EnmiugRB4YAQ1kMGgabmUslN58zZfd7rdrCa6GkTE2uucmyKN\n8+YtkE5AEKfZezWxIp73eZ/3jVgixhh5Ls/luTyX5/LPusj/1BfwXJ7Lc3kuz+WXX57B/rk8l+fy\nXP4FlGewfy7P5bk8l38B5Rnsn8tzeS7P5V9AeQb75/Jcnstz+RdQnsH+uTyX5/Jc/gWUnwv2f/AH\nf8CrV6/4tV/7tX/yM3/8x3/MD37wA37913+dv/u7v/taL/C5PJfn8lyeyy9efi7Y//7v/z5//dd/\n/U/+/6/+6q/40Y9+xN///d/zZ3/2Z/zRH/3R13qBz+W5PJfn8lx+8fJzwf43f/M3ubi4+Cf//5d/\n+Zf823/7bwH4jd/4De7u7nj9+vXXd4XP5bk8l+fyXH7h8gtr9p999hmffPLJ/Pu3vvUtfvrTn/6i\nh30uz+W5PJfn8jWWryVA+3THBSHE13HY5/JcnstzeS5fU9G/6AE+/vhjPv300/n3n/70p3z88cdf\n+Zw4/z7c/8Mverrn8lyey3P5F1W+973v8aMf/egXPs4vDPa/8zu/ww9/+EN+93d/l7/927/l/Pyc\nV69effWD9/8A/y4idUBvRsx6Su1qQoiAJCJIrRsrxr5hOjZMfYO/1/CGVN8CryN81sNnw6m6A7A/\nVV1DdQXVdarnZ/CrwL9mbrv1nm69Z7Xa0a32OBTHqeM4rThOHeNDA59J+EzA56UdTvWzAV5p+EEL\nv9LADxqqb3vWLx9Yv3hg8+Ke9vLASD3XgZr7/+mHmP/+f8b+Q830o4bgFbwCXgKvQF87tptbzja3\nbLepPcaOQ1hziCv2Yc10V+NfV/jXBvfaoKbA9sM7Nh/cp/bVPUNsGGgYaRioCVFCZK4aRyMGWtnT\niB4TLfe7Cx72qd7vzlm3O843N5yvbzjf3FCZES8UnlQDEkFAxJieYQxMsWKKNWNuBRGDTTVatPAI\nERAi8On/8r/x3T/5H4kI4mIUKDwal2q0uKjZ+S0PYcPOb+lDy7m850zccSbvORd3GGnRwqGFRwmH\nFAEAIZLnKYjEICAIYoAYBFFKghAEKQlSgACZx6EkEAEbK6ZQMUWDi4YzcZfPnc6rooeYjxkFloo9\nKw6sOIjU9qKlJ1fRsmXHpbjhkvdcihtaev6PP/l3/Ld/8m8QgMRTM1HNI2eip+VIx4EVRzoGGiaq\nuQL5Ez0dRzqOQCT52amHlyUiUPj5CDUjGsdEhV0c16GxmPI08Kh8tPSsAAx2fl4Gy4oDG3as2bNm\nj0exY5N/W3Okm49V2txjrDjwv//J/8n/8Cf/+XxOFzUFJSIQhUAQ2fKQ644tD9SMVEyYfAeC+Ogc\n5Xh2vlLNLRe854r3XHHDFRrHJTdz3bCbe5E8fab5TOkJlf+LRX+XUvp9OcbTp372s/kd/lfO+Nbp\n+1+TUvJzwf73fu/3+Ju/+RvevXvHJ598wp/+6Z9irQXgD//wD/nt3/5t/uqv/orvf//7rFYr/vzP\n//yfPlgFUQm80ODBjwIrTL710w3HmCYhbUDXA1JJ4qCIe0kwiigECAMIwABN/tkBPRAhOLABQgSb\nP2ZjEq6aCNuI6DyycUjj0NIihSNUEhTIylOFCV8ZHBV+NPidhmOEMYDzEF06rfKgA1QRqogyHqMn\najXSyh6ZobEMsyOWyo3EUWH7CJ50jQGQIHTEaEujBtbiwBkP4AWTbYlWYW3LdN8QbiTxnSK+lsjJ\nUzUTq82ec3vHhXjPXqSJBWDRc/+Wvq6Y6DjMU7COA6qNCA2hE9hzyVrfc17dcF295Uq9RQnPgY4D\n63mwt/TJaNDTih47Q1Rq07ky2DNhcLlPAvfigU/ET78yIQURhU+fEx6HRimPiZaGkZGac+5SFXdc\ncIvBzv2s8ARkBi7DRIVHUckp1Txd3fzp9HQCIo+lNAUlAS0sRhVosKj5eaZWioAQESnD3K+aiY4j\nW+6xVE9gVlAz0JL6rKFHEGkY6OiZqAhINI6OfgbMt7zgnjNuuOQ1r4gI1uxnUO04LiA61XIfqeUR\nuIhs0PQJUud+O1XFgW4+y47NDPDlGSk8Wx7YZMDdsHt0lp52NuAbdqw44NAn40fDQMuG3fxMtzzw\nMZ8xkAjLIJo0BvI1lifW0tMtjlQzYuLJVMkYIIjZGHs093KLFVsOcsWD3OKRrLN5/oRP812n5ysI\nDDSc6E0ak8kgh9lkCOKjfpTERR+K+eeYf4/IhblMbZmfCv/zYPn/V/m5YP8Xf/EXP/cgP/zhD/+/\nnS1CDBCsIgiBiApcSA8i5gcSQWqPNg5lHMZMEME3Bi8NBIGfNFgDwZwMqAhAD1EAFkKAYMF6IIAC\nhggu01oTQQWijERBmuQCtHbUekTGgLGOyTRMESYr8UcNYzr8P/U8ElSUaZIe/AlA0vBX0WOcw04B\n0ccE8pbTbFQRJT1GWBox0HFkCC3KBRgVbqxxu5p4J+CdIL4GMQXMdqK76tnYB865BSIupkkVY7oq\nER8zjkpYOnq24oGVOBAria8ktlOMUbMRdwlexDteijfzgxxJk2+gyX7DxDZPVovJHkXipJIwM63C\n8Muk+Ak7PuLzJxCVWOpyqgQkFZZWDKzZYzGcc8cZ95xxzzl3M9iX6tAzVyyTrWZ8xCBPJigZhaes\nVePYsJvrmj0D9ewtFVZ3gsY0MComIocZaJsMaaWK/LBjHhUjNQ0jLf18B6XfkhE4IjPw3HHO53xE\nzZgNygOf8CnXvJuvobSPgVvOAH+CnDCPhXI18sm37jlD4/BIBpp8bSegqpi44JaXvOElb3jBG46s\nFr22yc9vmscAwJ41O9Zo1ggia3accc8lN3QcueSGAys0DkHEoWc+vWTwy6pjAt8qTlRxQkeH9BHp\nI8qn+TCqlnst6em4lRfZyNxyEe8455aATNcmNuzZcKSb6YteeAwnNxkM7tHz1bgnREKlno2F3cu5\nr4VIfVmK/E8F9l/32aSO6HZCtxOmHVGNxY0Vbqqwo8HZChkDWvvEwMwAFUyqYZKCGBU+aLJpTYy9\nJrH2GCD4BPIxkhA5MD8UEREyJBZuAk5KelqcV/RTCwJ8VPigCFHie413ycj4SsFawF5D1YCSICqS\nFTFAlS8ozNPlxOWrGRhGGuR//ZuMosJJRZTidB8VyetoIVaCaAReKrxQie3Lkbbas/YGM7a4twYn\nDG6sED2oIaCto/ITDQPSBZxXDK7hwW+xrgIvwEvwgpU4UNees3qPqRyr6kAvWo6xoxHDzNbX7NmK\nB865I+QJv2OTZBEkijCD6Bn3jNQz4EACkiU8Lifpr/7WCy64neWJwsYLa2uzeCEJjAsDYjGPZII1\n+9mIFB8hIPJZTkB+ct5T+9jI1LBgaMvPlZ8lgYYRg6XjQEAhCKgnIPv05yUrLIao8E+LweD4L3+r\n42M+mw3E8jgOzYoD3+KnbNjxbX6CxmWDl9jwmj16hpbUByfpIoHOiY27bJjELNPY/O2lh2SwnHNH\nw8gVN/Q0OMwj1qoIMyM/444tOzp6zrhfCConQWqknu9nzX4+0mmGDPw3v6XTGMbTMLDhgYiYe6yM\noWXvDjSzQVP4ZFTStEfmakJgLfeEmLyrc25pQ88qHujCkVU4MIkaJwy97IhSEITE4FhxYMuODQ+z\nEXMYRhpC7kFByNRhwjxi8gLlPToEdPBo73FSYaXGKoOVGi+SxySf+IFfI/x+g6WOCO0xzUTdHmm6\nI3UzMIiWgS5JO1bPrLaWA43qESoiJESh8MKkY4l89ZoE9oITa48RCOkpiyzdVCCqiKwConKIyhMl\n2GBwk2GYuoXIkTm4FYSoiEoT
       "text": [
        "<matplotlib.figure.Figure at 0x1094d7410>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By contrast, the spectrum of the modulated biquad is a lot dirtier and has high energy artifacts at high modulation\n",
      "rates. The problem gets worse at higher resonance settings."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rbj_lpf_mat(f, res):\n",
      "    return from_biquad(rbj_lpf(f, res))\n",
      "specgram(eval_mod(0.25 + 0.2 * sweep(n, 0.1), res, saw(0.05, n), rbj_lpf_mat))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvUmvZUt23/eLZnenu21mvi4fy2RRlgQIAgwYBAwQIDyw\nAQ00FkeENBE04VdQEf4GRo000ZBDAwJMcMoZR+bEECyyxKZenzdvc/rdRYQHsdfZcXZmsWTUe2mg\neAMI7Jt5dhM7dsR/rfVfK1aoEELguTyX5/JcnsuvddH/fzfguTyX5/JcnssPX57B/rk8l+fyXP4B\nlGewfy7P5bk8l38A5Rnsn8tzeS7P5R9AeQb75/Jcnstz+QdQnsH+uTyX5/Jc/gGUXwr2/+bf/Bte\nvXrFP/tn/+wXnvOHf/iH/PZv/zb//J//c/7iL/7ie23gc3kuz+W5PJdfvfxSsP/X//pf86d/+qe/\n8Pc/+ZM/4Wc/+xl/9Vd/xX/4D/+Bf/fv/t332sDn8lyey3N5Lr96+aVg/7u/+7tcXV39wt//03/6\nT/zBH/wBAL/zO7/D09MT33333ffXwufyXJ7Lc3kuv3L5lTn7r776itevX5/+/dlnn/Hll1/+qrd9\nLs/luTyX5/I9lu/FQTvNuKCU+j5u+1yey3N5Ls/leyr2V73Bp59+yhdffHH695dffsmnn376znnq\n5Y/h7r/+qo97Ls/luTyXf1Dlt37rt/jZz372K9/nVwb7f/kv/yU//elP+Vf/6l/x53/+51xeXvLq\n1at3T7z7r/DHIdoSFsiG+r4WhKH64bgHvkvqGwdv13CfVNcBHdAPtQRuxjpfwT8G/omCfzIcfXK6\nI7YtS+oe+H+G+l+Avwzg34B7E4/+O/hkAT9+OdaXi/joYqg58UVUADW80P/+v8F//+/hLzX8pYKg\n4HPgM+LxRfLuciyBaqizob01cIzHzLRcfPJ4qsuPnti0Kzbdim0Xj/1dTvjWEL42hG8Mlp7Z612s\nn++oPj4QlCIohqPilre84js+4ls+4ls6Mr7ks1PdsOIzvuQ1P+c1X/KaL/j5+kf89dOP+eunH/M3\n69/CW8Xl5cOpXsyeWLFhyZb/+yf/B//LT34HhzmrBQ0zDszZM+PAgi0XbFix5oINc3a05DQUtOS0\n5MzdkVW/O1VLhzMaZxTOKrxS2DZgWh+PnYcWVAsMtc0t9aLguMw5Lgp8pin3LdWhoTo0lPsWf9RD\nVfijRnUB7QK692jnUWEY51LNMM7T2gzfT75hDz/5P+En/+vwzRmuk2vN8M1dMl5DUmXe+EntzmsI\nCq80nlgBNB4VPDoE9OnhgJrcV47p71Kn6RTl3dVw7Id3ltq+p63JPX/yl/CT35q0P3COHdkwF9K6\nAe6Bt8OxB+aTc+QectwB2+HaLYQDuBZ8A76FYEF/Bvr1UD+L56nhfLbDO7mxdjVs17Bbw3YDxx0s\n5lIVizmYEigCqiDObzP0wU//M7z+J2M3f09MyS8F+9///d/nz/7sz3j79i2vX7/mj/7oj+i6DoB/\n+2//Lf/iX/wL/uRP/oQf//jHzOdz/uN//I+/+GYZ8YXKEF+uDFAE6BT0avio6nxyGCJoHokfIwOU\nBlWBUsQbXQw/boD1cHJDHDky4g+gczAFmByyIn5k+WAbIjBfAJfDsQRWwO3wew0cNRwMHDM4lBAK\n8Bk4E9+hJw5kiIPTSXsBNXxNb8DYCNwXxMlwzUkuqZeBrGixRYstOrKipTeWTmf0OqMzllCb4ZUV\nKIXGU9iGhd5xpR645i2VPTLTB+b2wLzYU4cZjaloZhXNTYXGUd0cubhZc7m8Z2k37wBoTsuMAys2\nXPFIR8aaCyqOZHQYHAU1c/ZcsOaGe7bliourB5azR+Y3a7xWzIttrNmOOXsqjlQcsTgKGgJqwAtF\nQL0D9nP2LNgx40BJTU6LJmBw5HR0NFS+oexrsrZDtxG8sB5lQbsoxHQHug1oAfmeM3DVJqCNw6qe\nHIX3mqzrMLVD7QJsiPfufBxWAoK/KHeszFMBZzeMDwEGmRfZMN5eyRgnKhuH5CiKkswPuZcAaAAW\nwHI4zoGnpO5ANQHjHKZ37z5fqh/uPQVhOcrvoiiFyfU554JB7mOGd8wm18vfIhjlvMXw3n1SSc4L\nnKY2PVBD2IB/BP8EYR1/V0fQZayqZFTG5O8DETJkPNjhXAPaQqhA34C6ATUbYMcO1/qhT6S/hrGg\nayhtPLdw0HaD/tdA7gL6CCof+iof2iNgL+/5PZdfCvZ//Md//Etv8tOf/vS/7WkZYAPMAmruYRFQ\nM0/YazgogtPgFOiAyoBZgCrEjtgqqBTBQuxdUXNllr0h9vaG8evLqGxB1aAWYBaQqQj2jgj4b4Bv\nGSebyJBseMQV8PHwqEcLjxWwhNpBKMHPwOXQ6zj4NGACysUjXhFOGowaP6bIqZwzI0S98OSXNbOL\nHdXljtnljrqvOLax0s5wW0vQhhAModNo5ymyhoXZca0eecUb5vrAQu9Z2h07FuzMim21YnezItQK\nExyzcs9l8cSL8o5LHtmxYM+cHQt6LBkdFUcWbLnkiY6MJVsqjuS0WDrKAexXrLnigXWx4rK45oJH\nVjzgMCxYs2B7Au+KAwU1hp6MdtAzA2o4ltRUHJlxmNQjBQ1FaLHKkWFwtNEa8B1l35J1PboJKA86\ni99Bu6htq3astIxArQADygSMcmQKIOC9Iut67NGhdz7qEqlGbTnXcBn+PwUuzQh+AgpSZLILeNwS\nx7smaqcCQKKMCEAJ2PVEoNoN95sNY/YF8BL4hnGct8PU6Ia/u8lUEmGSgrkIBNHSBewTawg/PFfe\nNU/6Qar0h5Sp1TG1qkui0BILVqxwyzjlVdLGQbcLW/BrcDtww7uaNt5D2QFg50MVQdUwkgJ+OG9B\nFDYiPOWaWdLWIukbAf2hmBoqH4E+yL0ZDXzVQBisHSUCXDyojh+k/Mo0zv+nUoDKAmbWYxYdetVj\nFj3OWpzJcFh8yFC5R8979MJh5j1k4OcGV1p8ZvHaJINPelhGQ0Mc0cPsZU6cQddgc8hzqPL4ERtG\nI6AlDrBLoja/HC6XySPaQJaDXYDW8T9VBnoOZgbWosqAunSoCxePC09wmuA0fjjyP/1uHAnrYfaI\n1lPFZ+srT7U4clGtucruueKeg5mxL2bsszn7ak6jKtqmpD2UdFmJJlCohqXacqUeecmbAVSjBl1S\nU9qaSh2YZXvmsy0ax62+49bccWPecsETZhhpHZaagpyGigNLdiewX7FmyYbFQKUs2LFky4otF6y5\nZM01j9xwzxOX9FiWbFkOYB/bVFNR849+7yNK6hPYC8FQUp+1v+JI5Y+UoaYMNQVNpCOUxqtITdjg\nyX2P7RyqDSfWTA3gEDSoAVyUaMQpADtQRcDiUNpjTE/QCtN6TBfQXRgBJ7U8U4oiFQLpOQKw/fC3\nZaT5ivj83/vHRJpS6B0REDnw0eSelhGwpwzmljg/GiKV8TSM81TISBHQhdFSEAG0GNrWJe9XT54n\nwCRzRyxUOXegGU9zLBU0cg85JlTl733EKCSEblGM7IAAY0qVGVA9mA60B6vic89kcRjOzYdnLYb7\nDDRMOMZ/qm5o+xZ4ZKRlpaaWjvRjaqkM76QL4CWEq3i/cIjP8McB9AdF4EOEtHxYsK8GsK96sllL\nPquxs5bOFXTO0zqF9xade+ysI5u12FkDAfoyp8sLeqvwxpxzohoIgZP6HBwnW1edEDR+/UJBpWHh\nI+3iGOSCigPgijhoF5xrPRVQKigs2BKMjmqCMpEe0hkYjSo82aIlu67JXjTYqw7nDN4ZvI9H9z//\nj7j/XOMLiwuW4BXKBCgDahmwlx2zYs9l+cgr+x0v1bcc1IytFlhdsG+X7PMlexvwOkMrT6kaFmrH\nNQ+84rsR5IngWOqaSh+psgMzdhgcN9xzw1tuuWfJBlD0GBpyjlSUA52yHIC8x3LFE1c8cckTHs3l\nAPAXA59+Ofx2xSNXPNCRMWd/
       "text": [
        "<matplotlib.figure.Figure at 0x1094cb390>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It's visually even easier to see when the modulation has sharp steps in it."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "res = 0.1\n",
      "specgram(eval_mod(0.15 + 0.1 * sign(sweep(n, 0.1)), res, saw(0.05, n), svf_lp))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvUmvLct15/eLJrvdnXNu81pSrBKpIgRDluAGgg0IpfLM\nNRA8MSCNBMkDQQPrGxigvoLNkQBDgCeaGhrQgmeaCa6yNXKVGqpUksjHx3e70+wmu4jwYMXKzL3v\nJWiIj08AeQKIm+fuc3Y2kRH/tdZ/NWFSSonH9tge22N7bD/Rzf5T38Bje2yP7bE9th9/ewT7x/bY\nHttj+yloj2D/2B7bY3tsPwXtEewf22N7bI/tp6A9gv1je2yP7bH9FLRHsH9sj+2xPbafgvZDwf63\nf/u3ef/99/mFX/iFH/g3v/d7v8fP/dzP8Yu/+Iv8+Z//+ed6g4/tsT22x/bYfvT2Q8H+t37rt/iT\nP/mTH/j7b33rW3z729/mr//6r/mDP/gDfvd3f/dzvcHH9tge22N7bD96+6Fg/yu/8ivc3Nz8wN//\n8R//Mb/5m78JwC//8i9ze3vL97///c/vDh/bY3tsj+2x/cjtR+bsv/vd7/LlL395+v+XvvQlvvOd\n7/yop31sj+2xPbbH9jm2z8VBe1lxwRjzeZz2sT22x/bYHtvn1PyPeoKPP/6Yf/iHf5j+/53vfIeP\nP/74rb8z5ddg+Jsf9XKP7bE9tsf2U9W++tWv8u1vf/tHPs+PDPa/9mu/xje/+U1+/dd/nT/7sz/j\n+vqa999//+0/HP4G/mWCATgs+vEdJ22AFbAB1kAJmNwtEEe4fQm3L+DuhfwcTX6cQo5fWsMvXcN/\ndgO/dI35F2tsNeCqEVcN2Gog/t8F4d944r8pCP+XJ9UWvo70fwF8CagX3Uf480/h//ke/HnuH2/g\nv3gf/vP35fhsC7cW3uR+MHBF7kmO/8s34H/4BnQGOqjKls31HdurezbXd9SbE22q6Jj7dbzjJr7J\n/ZYulbwyT3mF9IPZYG3EmoC1EZMiw1/UDH9ZMfxFTf+XFcWzgerrJ6qvH6m/fsI8TbRtQ9c2tF3D\n0Ff4uqeoO3zTU9Q9xiQgW24GUjTEaInBEoMjJYOxCWMixiasTVSmpaalNtKTMfSUDBT0lITkIIFJ\nifvf/5+5+p/+RyrTU5qOynaUpscSMQkgYRIEHB0lnZHxGPFcmTuuuGPHHVfmjoCfxqulYkyF3GeU\nPo4lD6cd+3bHw+mKfbtj1ezZrO7ZrB5Yr+4pXI81CUPEkAid5+HFNQ8vd3J8tWX37I7d01t2z27Z\nPrtj8J593LCPa/ZxQ8DRuBO1b2lcS2l7TmPDaVhxHBtOY8PKndgWD2z9Pdvigdq2/P03/je+/o3/\nHs+IIbFnc9abdGLFkVU+pmTO5kjAUTDgGShyl4GWdyeHiJ16YsRzYM2eDQez5kTDmsPc00GuaY5y\nzAv2mFacUsOJFW2qCcYSjSXgiMaSMKRkSNh8BIwewTOy454tD9NxOUf+7Tf+D37xG7+GZ8SnEc+I\nJYLORmMYk+P7x4/4/ulDPj19xPePH7Irb3m//pT3m094v/mUtd/P2LG4tiPIuRl5YMs9u+k4UGDk\nCTAkPCMrDqw5suIwjcFlS/lCKf/rCPk6crR5Tun4H1nlJ5feUQHwv/KrfIXtdN7Piyn5oWD/G7/x\nG/zpn/4pL1++5Mtf/jK///u/zzAMAPzO7/wO//pf/2u+9a1v8bWvfY31es0f/uEf/uCTVYADAhBz\nVwZoeawQcK1yXyECoM5HZ+B7BryFwcGdR96mzxfwYBx4A2WCOmBWA0XVUlUtVdVRli19U9PVDX2V\niJUj6fnXwBbMVcKuArYZ5ViOxJsDcdMSq5FoEs4PFPURv72jeOIxz3oGUzFQMsSK0ZXYq4C9Cdjr\ngLkODM2Ive5JJ0dsHdZHqrpnXR64cneszIGTaTjS4AkYEmu7Z8cdT3jNc17QppqEISRHn0oArIky\nqUzApki7jbTvQ+w8xpW4q5HqSy2bp3s2zT3GRfZuh3GG0ZYMFrwdqVxHbY/U9pjBfm4BTzfWdH3F\n2BcMY5mnd57qKeHKBGWPLwN12RKxhOjoYkUfC0Y8te2obUtnW94rP6Okp2CYjiOe1tQZuGvGPFUN\niYYTjsCOe6645SofBwpONBkYAqPxMhXyvY+hgGAIfUGbVjDKXDRJhJV1EedH/AIIQukIN46xtAxX\njuEjw2Z9x9X6DTfrV9w0r4nWcoqN9NQwJo+1AWcj1gaMTk2TwEdMjKztga2948qKoGo48YIjN7zG\nEybAiFiGLMRqc2LLwyTkRjz37EjsRLjh8Aw4AiU9NSciTsYfS8RlfSllAJJnjHmEHIGK7lzEmP0E\n8tojlsIMWBMzbBk8YxY0ckwYkhG4xJDhzk1XjVgMkZ6SVzzlNU8WQihO9+MZKYwILkskTeLK4IGr\n8g2FGbgu3vDl5u/YuAduitfcFK+5Ll7j7YioHg0nGgYKWSsZeB1hmnclPRXd9O7lTmUslwJwxeEt\naDMLEPtB0KxzW693ZMVmEnd7+gz2nvEHnOFHaz8U7P/oj/7oh57km9/85v+/q1X5ikuwB0gXvbro\nK0Qjvs7HKn+vM6JFG4sMscvdgjHgwBQR6hHb9JRly6o85H7kWG5wZSJVjqGqiWpRKNjvIm41UKw7\n/LrDlx3jds/QnBiLjmhGvO1pyiNN7VivE2w7jsOG47ghDY5gC+xuxO/6qaeqx21agitJzmBtoKx6\nGn9iZx/YcH+myQQsK05s2LPjnmtuZyA0NQfWBOykQagWYdeW+NwzmAqzAbcK1E9b1td7rspbrIlg\nDMEUtLbGuorC9tSmzZN6P91D1snoY4UdIfWesQU6R8pCO+k7bRx2DYUZqYuWEUcXSwiJMDpCcjgf\naNyJNnU858VbYHFgTU9JR8UdV/SppIpd7j116mjsibU9sbEHtnbPQDEBtWMkZHBTDW3EM6aSPtUc\n45oi9Pg44NKITQGTx67IOmZJT3KGceMZVp7+mWeIlpW9Z2PuuLZveGJfY4kMrmBw8s1xggt9Ik+y\nZvKQGdLZu7zhDQ0nalp23GMREB3xWU+Xc23Ys+WeK+64zsIt4Gipz96TAlhFP92LgCTMUKSaa6Bk\nIHHCEqnoJl1TrvfwFtjr+cYshEb8hQ2yz2cXSBVh4yZ7Q7+nerLqzE2G5ZoTNj8HkAVeMYG9NmMS\n18UtT4tX+Qpy/3Uey5qWEc8dV1nArOiosMQJUPVnEY4tAUfCTKBcZAGw5sBmAfic3cn8XrWnaTa5\nabzq6fmkr7OVsOHAjnt6RGkr/qnA/nNtK0SbErVlpmUSs5YfOadOVMtXWmeX/78xUFsorGjx0wnl\nFRibsD5gywFb9/gmUbsTjT+xskfWHEjWMfqSrhyhSdJXwDrBOmFWAd/0lFVL5VtKe6KzLcb0JAIj\nMlkKBhpa1jiMMaTCMdYl3TqAAbcOlM0g5yla+Jf/JZQ9XbCEWGBNwnvRqBtzmrQnXSB9nnANLSvk\n3j0jezZnC1A1EceIIxKKgn5V0V4FMAlbRfx2pK46Vu6EMyNDUTLYkr4siAFqf6D2R2ovC8Zly0I1\nrugObOsHhqJgXBWE6IRuSXnyJ1j5I407sPJHVvbAiGflBJR3/oExeVZWqIHmX31tep6ZghhJmEmr\nL+kZh4LuruF0tyXeO+wh8cHVJ5idpb5q4SpR2Q7PwIoDcaFJau9tSVw5XDlSbw/shjdURUvjjzTF\nicYeKemne/BZYCg4HVhDfuflAgQcYQHwfqIiLimWWZcO+btHrrjjKa+p6Pjnv/ozWJJYhRQ4Ijvu\nWXPkPT6jWIiPjpKEYc2Bio73+IyAO9OOFXSWfUniBCxDvk+F8p4SS6DhhCVQ5zm34pR141MmgMxk\nCXRUrDgipI6c6e3rmresDJ2pSnaoXpww2F99zoojD2x4ww0PbEgYrrjjint23LHjfqEpy6i76exy\npZGChJ3eWcvD4k2FSakq6dkg
       "text": [
        "<matplotlib.figure.Figure at 0x109489510>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "specgram(eval_mod(0.15 + 0.1 * sign(sweep(n, 0.1)), res, saw(0.05, n), rbj_lpf_mat))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvVuobklW7/mLiDnnd1tr7bUvmVmZWVlVWOXx2HYhco54\n5xRC03Ae6rmEAi8oIogtjSDNobv1yfemoMF+8FGwn/ShqAZfBNuj9mnso+1ptcoLZmVV7sx9WZdv\nfd83LxGjHyLGN8c319qZZWbuLVgrIJjrMi8xY0b8xxj/MWKEExHhttyW23Jbbsu/6OL/uRtwW27L\nbbktt+X5l1uwvy235bbclm+Dcgv2t+W23Jbb8m1QbsH+ttyW23Jbvg3KLdjflttyW27Lt0G5Bfvb\ncltuy235NijvC/Y/8zM/wyuvvMJnP/vZZ57zS7/0S3znd34n3/u938uf/dmffaQNvC235bbcltvy\n4cv7gv1P//RP85WvfOWZ///yl7/M1772Nb761a/ym7/5m/zCL/zCR9rA23JbbsttuS0fvrwv2P/Y\nj/0Yd+/efeb/f+/3fo+f/MmfBOAHfuAHODs74+HDhx9dC2/Lbbktt+W2fOjyoTn7t956izfeeGP/\n+8c//nG+/vWvf9jb3pbbcltuy235CMtH4qCdZlxwzn0Ut70tt+W23Jbb8hGV6sPe4PXXX+fNN9/c\n//71r3+d119//dp57ugzcPW3H/Zxt+W23Jbb8m1VPv3pT/O1r33tQ9/nQ4P95z//eb70pS/xhS98\ngT/+4z/m9PSUV1555fqJV38LPy+wA94B3i310Q03fQC8NtZ7H3/Md37yr/lXn/wrvvMTf8Vrp//I\n3/zhPf7mj+7xN394j6/+n/fodzOgBhqg5hPf8Zh/8yNf5d/+yF/zb3/kq/zr736T2WXPfN2VY89/\n/Msf4g//4kf5o7/8Yf7wL36U5StrPvtD/5nP/tD/w2d/8D/zHd/xd5x061z7NYvtlj/6vz7NH/7p\nd+yPn3zwhB/9nr/jh7/nb/nR//rveP0T51zdXXB1Z8H67oL2qGG127JqS91t+Z//F/gf/kNNPwv0\n84BzQjVE6iFRDxFJjovqiIvqmPPqmIvqiIQn4REcCU9Dx5LNvs6kpZKBKg1UEglpIPSC7yH0QhgE\ntuC2wCbX3ldsTudcnS64ujunPW6Yxx3ztGNRjr4T3M7h2lx7X7FZzNku5mzmc7pZTSASiPhybIaB\nZhiYDT3NMBC6iG8Fv0v4NuGi5E9Vwa/9r/Br/z3QT2oAaRzSAI3LNmgvuB7cIPmcWOpgfrbV5ftQ\nlap2rJijnpvK0ZdhVIPU+TS3BbTfdtC/7Ble8vQvBYaXPH4rVI8T9aNI/SjhN5Kf6xmPYiqwOZlz\ndW/J+t6C9b0lLgj/2//0kP/w39WstlsW2x204FpgJ/lYHVaZgSwhLRyyBAngN7m9fiO53W1uMy3I\nrrxPKO3SWpt7OpCnnvTUIU8c8tTjEFwQfJVwQXAzgQWwLLUZv5v0QFeeo/NZyOffAU7KcVH+nszx\nnbH+2v8Ov/bfAHeBe/mYTh39cUV3HOhPAsM8UL850PxjpH4zV7eQg2sAOCv1HLi4PtZkDixzH7Jw\nuf3ngjsHdyb5u/vyQq78PPkW+/fU96nLe+o7H5k2nJefk7mvHgH+x/8CH/tu7b2PjCl5X7D/iZ/4\nCf7gD/6AR48e8cYbb/Drv/7r9H0PwM///M/z7//9v+fLX/4yn/nMZ1itVvzWb/3Ws2+mE68BZsCc\n/NHhcAJmvB4nSulkh5QfBXA4PO5gNptRWy52uP01DgEBpx/FgkFTagXOg3Oyv8YhOAGXAJEyaWXS\naPu//Dcn9vpyD2H8fXp/yX/XNl47pzxn+rdnV/b3Jpn3NsCjQiRKxUBFdIHkQxYsztE2NX2Y0c0a\nujhDoiOkhN9ETtZX+JSK+JH9Ufsz1o5tUxN8oEqRuh+oEEKSEVgVZJNpW+lSJ5L7PMp43nvVgTyB\ny1FqSHOHHDnSUZ7Ufgt+l4HQ7+RwDNgJ7cDFSfty5xM2gns3Ea6E9K7DIfhB8CK4O8Ax14FMwaXL\nx9oPLMOWyg0sZQcVzHctq4uBpu3xreRzu/Gafd+Udrge2GRhLtvS/ghuKM+084fy/WtGkF6W/6sw\nW2dB4a4Sfk0WEpE8T49KPS7zBNNXA3uh4nblOvt/xzivZ+V9gukT/W6X7AUFAeQY0j1HetkhL0O6\n45BaCHXEk5i1A36XCDshbBNsZOynDrgq37Ur7zEjA+/GnNOVed0LbgPUks87Iiuds3LuI0YF9V3z\nLnqcgn9RFNiWZ1+W9mjdmXfXqmOs57mU9wX73/7t337fm3zpS1/61p6mmoSC65w84GRSC+gSzLU3\nJmJWkag3nv6sI+2wiL1cNRv9cPrc9xSmOoMtQr24TNEW+D9MUXgeqOip6Wio3EDFQO0D0XnaMGPN\niiuOWHNEaBMn6zUnV2uO11csN9u9YKEIst1xk+tJQzur8ZUw7zskODxCkHQdDPV1jPYkDnCuHPWT\nFOElk2vtJyngL7Uj1Y545Ij3suYWLoTgMkDuAXE6Qe19tF2qzSWKlSLwuJy3IIPgMVmbrLhuqThG\nBUOgbgfq9cBKthkoAyw2sDpjtFRuuocOcfI5LpItNv2oU419YJwOarUsQE7IwFfe163JgPSUEXyc\nebe7ptbkNmtVoaCgrYCt31P7b5fbugfQHQeWB2qJuNyHsnKkE0c8dcT7DjmGMAjVUAC+laypr9lb\nXbSlDVdk7dka/KpgutLeXfkmKkwVMu4C9xnZhQb4m3Lu49JHi3Kvefn7zFyvOmdkb0WTyvOmdWuO\nQ+mrfy6w/8ifph2jZmDHNROXBbmD1ey+hqd21FuN/n1R+rDoBG4YP54C/nP0Mf+7H/pgN7cAP9X2\nP2jRe3jV8SURUtzX5bBh1vec9mtiXyHiwQscCZvjhp2rqONAnWI+xoGqHljUidp3LLuA64XQpzxJ\nJdqH87kfIH82OTym2jEsPLFUcVBtE2GbCAghFoGh30prIH/DCM4LYQB/LlQ7oHK4XnAduF7G4ZIY\nAVbLVJhg7q9amFoUDftJLnNgVrRuW5cc6gcKck+Bb+a/fe414C1GcI6T56hQWpj7bU2N5e+rclww\nAsfAAbA4BULHCNoquKZCWBWhZO7RlWus1aFz52jS9sRoDQxk8OvNPfQ+eq8In/s0uI0Q3gZ/IVRv\nATOH87K3vIGscT9h1Na1b4fSzqnRb4VOZBxzMGrWl+U7XABfL330Lpli2jIKDKUSdxwqsvq8Mhb2\n/acCfMeo3VuhqePvOemNLxbsa8bBs2A0rybUwp4HtFzYtQ6w5FkhWW/U6q8D6/4vquloe1TI6C1u\nLFMz5J+u1f+7H3bfsvB25f5TUP/ogT7z7UEy4HtJ+CS4GAkRZIjI0OcznWdwnlh7ovf45AkxIQWU\nfEg4nzJVnbKZ7NIhLamf5nM/eEPDBMQ58A4JjlRlzj71Dh9cFjbKg9ui37P8Xamw/eRGDodFzaEW\nn8z/VBPWBiuQqPZlNWXVHqdDT8+xOogeE1kjXbPXpj/3Khlo5oxaudZo2hvKM612qoBp+V/HIaB0\n5lzrT5hOpenQVhCTcq0wgrXOYdsPNYe0VWeeo0LHArxW866f+458juskc9wOqORQo24YQVK/p1pj\namWo8NexEjikDRXgbVFKSYWxMFoOFXBqvolSWFYhoJynVJm2WYWyCku9h/AsqPpIy4sHe8gvrnzV\n1AQXcucoD/ZM0H2WZu85nKnv0Ytq6s7J2pA+t36v5/KMRn904tj6GFzRuPOY9OXpH92o8CQCAzUD\nEU9wA94lcJCCo/c1u9mcljk75vgO5puWxdWO1aZl1nZ49W84wbms3aalIy0g1eAq8OVb+qKZ7Xlc\nO+kMuHmXsuNaIhILh9MVrTwqmpuqgGrBS++pwKPasR06+ncFJgXwxtxTz1GqYTGpxhp0G0at0dY0\nOW7IQK8a9ZRSrM25+lwFedWk
       "text": [
        "<matplotlib.figure.Figure at 0x10974e750>"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With very high modulation, the biquad can actually become unstable, while the state variable formulation performs fine."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "r = 0.9\n",
      "specgram(eval_mod(0.25 + 0.185 * sign(sweep(n, 0.1)), res, saw(0.05, n), svf_lp))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvbmvPUl25/c5Ebnde9/6W6qqu7rZ3VyEkTAUZQgiBYkC\nPUEEZhw5JDAAQToEHf4Nzf+hLTq0BDqSQ4MgJEcEKA4hh5YcsYdbV3ctv+0td8vMiDhjRERm3Hyv\n2AN2VQ3QfAEE8r17b0ZGxvI953zPiQhRVeUpPaWn9JSe0k91Mv+5K/CUntJTekpP6ctPT2D/lJ7S\nU3pK/wLSE9g/paf0lJ7Sv4D0BPZP6Sk9paf0LyA9gf1TekpP6Sn9C0hPYP+UntJTekr/AtKPBfvf\n+Z3f4f333+cXf/EXP/c3v//7v88v/MIv8Eu/9Ev89V//9Rdawaf0lJ7SU3pKP3n6sWD/27/92/zZ\nn/3Z537/p3/6p3z/+9/nb/7mb/jDP/xDfu/3fu8LreBTekpP6Sk9pZ88/Viw/9Vf/VWur68/9/s/\n+ZM/4bd+67cA+OVf/mVubm749NNPv7gaPqWn9JSe0lP6idNPzNn/8Ic/5Jvf/Ob0/ze+8Q0++uij\nn7TYp/SUntJTekpfYPpCHLTLHRdE5Iso9ik9paf0lJ7SF5Sqn7SADz/8kB/84AfT/x999BEffvjh\ng9/Jt34W/vHvftLHPaWn9JSe0r+o9HM/93N8//vf/4nL+YnB/t/+23/L9773PX7jN36Dv/qrv+Lq\n6or333//4Q//8e/Y3L1CnSHsa/y+wu9rwuGRKphUswqwsGm2vNd9wgfdx7zXfsyVf8Wnf2n59N8b\nPv1Lw6f/3uCPDdAALdBw8bOeD/6HfcoHnv/rAy09DQMNAy09P3j9Lf7h1Xf4x9ff5h9efQfbep69\nfMXzF694/vIVl+c3tPRTrtzAx/+P5Ud/YfnRX1b86C8sZ9+ED35Fef9XlA/+e6X7GcuRFUfpOLJi\nkJq12bO2c/6r7/5f/Dff/TeM1IzUGAIdR1qOdByxePasT/JaD6x0z0YPrHXPnT/nY/d1PnZf40fu\na9zpJVXtqKqRqnbU1cgz3nLNO57xlmfyloOueMf1lD2WC+64kDsu5JZz7mkZaOhpUxtZPBUOi8fi\nARi1xmmF04qRmoGakYaRhkFqWno6jnQS36eWcbrf4vFYbrnkhiv+4rv/N9/+7r/jPT47yRUjAUvA\nLLIQMDhqbrjiLdfccM07rlhx4JK7+E7cseIwPbPCoQj3nJ/kDTvOdMsZW865xxBwVHgsTiocFQNN\nes94bdIb50/ir+asCH3qzTx6yl+MVFMbrTiwIo7N//O7/y//63f/q6l85dRC9ietOLeNIgQ1CDrd\n20i83nLJbWqVWy5p6dMnN1xyi8fysX6NH/IhH+vXecMzviN/z3f4u5jl76gZpzoIypazk9664epB\nGwg6/T5fY+8pglLh4hgp8tymA//bd/+Of/fdbxezr8OrxaqnCh6rHhs81jgq46jEYc383PzkCj/N\n+dw2eSRZPIbAp7zPx3xtygNNcU9PRz/105o9Kw7TGMnZpufEmTBMcyUnTeM2YKa+M4RpbpVt9kv8\nH6z52bnNvyCm5MeC/W/+5m/y53/+57x+/ZpvfvOb/MEf/AHjGDv/d3/3d/n1X/91/vRP/5Sf//mf\nZ7PZ8Ed/9EefW1bbHVFnGYNCENRbQkhfanG1QKVgFSowzciq2XHR3PCifcUL/wljvWZXrXhn1kgC\neKjJUkLIMiNQMxYTa886Xe/sFZtmR9MNmI0iTaBuHV3dszZ7zthOA3HFgZqeLS0rWmpahJbGjJw3\nB152R765OXB2LuxYs2fDjjU9XSrjkMqJ5V1wOw0ni2dTQHvNyJYNuwhFtPRsZMeZbDknAtPGXuBr\ny6gVfWjiJBKNLy0gEif+GVue8Zb3+ZRBmqlm59zjqIqaxme3HNPUSgIOR10AmiI4qRglfuqLIZSn\ndSlQGwYCJomPJn1S47Dp94rFUzPSpkm1YYfHcsd6AmVPRcchtV+cdAGDEGL7pL6aYXw7Cas8sUFY\ns+ec+2kUtBxZyWGazI6aHRu2qUdGajoObNhOffcYsOkEY4LHFvASBYRfCK6G4UTAN4ycseUFr6hS\nncMEj+YU2B/5PEj8TtI3AWEgCt7nvOaadwRMautZfA00nMs9L3gNAhu2fJ0f8ZLPuOIdZ2yxuAd1\n37DjmncArNkTxb6brstUgn3ujxnGo2icFYsoCK55l0qLYiCIUEkC+EJxiAAaS210oA09XehpQx+F\nt6lwxuJMxUG6qZ5CwKKccw/p+nV+xJGueJM4vuNMnGflUvBm4ZXfyeJPxkQowD6Lmip46jDSBEcd\nHFajgKhb96WsgPqxYP/Hf/zHP7aQ733ve/9JD2vrniAWGkGdxbsKAhElymwVqgCVIlax9UhrD2zM\nPVdywzPecotnjaFmxQzy+WoBjwFsAvuGgRUHzhJYnrHltbljVR+ouwHZBEwdqFpHU/Ws5MiG3QQC\nKw40HNmwpiVQIQgVNX0C1Fs+4I5LwvSEe84SoPQT+OV8zv30X4VLmmW8c/5tnpQjm/R9BrNaRg6S\n4Md2wKnmFzB0HDljyyW3POMtHktHfK8L7nBUD4B5lYRSVwBbqd0HzGSRjNQEzCQQylwV0/RIxx3n\n3HFBT4OjQpFp2LcF6GXBumNDT8c7nvEp7zNSTzZJxZi0rB11AskIyj3rSdTuqRagowhVEvwbdgw0\nVMkmyRplVgmigIp1veCOa254zmue8XaClwzfHnvy9iVQ5LxMeUzmXOFYceCKm6ThydTyHju1mab3\nKAVLqTHm3+a/s3DLwhw4qVfAsGGP5waL55I7PuATXvCaK27YsEXQwiqpqRjZsINkSVxxU/T5WFgl\ns16/BPsKT6NpXmhsgxIaO+250jt8fi+JyoEhYHUWd32GWGkZqGjDwMr1XLo7Lv1dtETshm21Yawa\njrYlMBCS1SQoK/ZsdIfwMQblwGqyiG7lkp6Wc+5n7NDtA4uzwsW+lDiW8lzxxW9E0yhUEFVq72nc\nSOscjRuxWfOt/X8esP8iU8cRLxVqDKGyuKaC4BPIS8ogNiC1RyqPqTxNfaStD6yrPRsTQW9FS8MZ\nFSAPgN4iCBalxtEUmv2aAxv2nLFjZQ809UDVeWRUpFZs66lrR2Oy1nGczLeWAx1CQ5U0RU0aypEz\ndlxzw1UCuKxRVrgE2rNZ+F/+2kvWHJKeRpo80ZK44I6WfhomWVPLAy2DPcAV7yZd1xbDShIB8pJX\nvOA1L3g9mexZg16zx1FN0z5r76VAaumnvlMER01IAKNpAmet1qXPXBILDbaYVMqKAxVuEnIZVPtf\n2/A1vs81b3mWwDz/Nk/E57yhTxZQ1JgMO87SBIsAX6V+mAXGEUEfiJ8jbYK9qLlnIdqme9ZJIXiP\nzxho0MlSmdsjt29uu6wx5wnvsAVJEK+n0KdTP2RBKyj/469F4M2gGpKFkHtjBvlYiqTSM4Uw0yS5\nx/RBf0Zh3UztkgVL1rRHal7wmmu94YJ7NrqPapMInoogBoMm22eH8grgxOYwqddn0TQLqPwLo0o3\n9nSuZ5Wu0SpVxMD//Cue68MdB9txtC1HY1ERmuBofdTc6+DY2g3WnBGMobctEkA8mBHsAFY8q7oH\ngcqOrGgmkeOp2FPThIE2DLR+jFdx1MbTmZ4zs8WZmpXuWYUDK40ZBUl4RRDUAAbUgFoIMrcCCKJQ\njY569FM2PmBDwISADZqEAXD5hUMv8BWDfUuPF0ewhlBbXLA4sagKBAE1qArGekztsLXDVI6mOtKZ\nAyu7Z212rNnRcU5DwGLJPP1jYB9BYDjRGrOWs7JHmqanWo2Ieoz12C7y3o0ZHoB9x56OipqGig5B\nsQkgN2y54JbrNOmACRBLDa6l55d+7YoDh+k39aSpztQRUDC2cmKRnLEF4JJNAvs2McGuAG+XmOwb\nrrjhnHsCZtIoO464JLRm3SOb
       "text": [
        "<matplotlib.figure.Figure at 0x1094a4550>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "specgram(eval_mod(0.25 + 0.185 * sign(sweep(n, 0.1)), res, saw(0.05, n), rbj_lpf_mat))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvU2sLUtW5/dbEZEfe+9z7sd7r+oVPAqsgrKFWohBY9Nq\nQGLiHjBgDCNETxAThj20Go/akoc1YmCGjDxhUGKKDG08QrZsI5mCAuqB+1Hv4957ztk7PyJiebAi\nMnPvc2691/U+Sq46IcXNfe7emRkZmflfa/3XR4iqKo/tsT22x/bYfqib+0EP4LE9tsf22B7b598e\nwf6xPbbH9th+BNoj2D+2x/bYHtuPQHsE+8f22B7bY/sRaI9g/9ge22N7bD8C7RHsH9tje2yP7Ueg\nfSzY/9t/+295++23+bmf+7nX/ub3fu/3+PrXv87P//zP85d/+Zef6QAf22N7bI/tsX369rFg/9u/\n/dv8yZ/8yWu//+Y3v8m3vvUt/vqv/5o/+IM/4Hd/93c/0wE+tsf22B7bY/v07WPB/ld+5Vd4/vz5\na7//4z/+Y37rt34LgF/8xV/kxYsXvPfee5/dCB/bY3tsj+2xfer2qTn7f/zHf+SrX/3q8vdP/MRP\n8O67737awz62x/bYHttj+wzbZ+Kgvay4ICKfxWEf22N7bI/tsX1GLXzaA7zzzjt85zvfWf5+9913\neeedd+79zn3tv0C//fef9nSP7bE9tsf2I9V++qd/mm9961uf+jifGux//dd/nW984xv8xm/8Bn/x\nF3/Bs2fPePvtt+/9Tr/997yT/5p85xn/fsf0Dz3T3/dM7/b3D5qAuPbdkyNvfe2feeunrV+/+QHv\n/1ni/f+Yef/PEh/8eSINAWhKb/nxr73i537p3aV//V+8x1W65Trdch3vuE63/G9/96/482//Mv/r\nt/81f/7tXyZdO77yM//Ej339n/jKz/wTb779Xa645cAdB+7YxSPf+vM9f/NnV3zrP+751p9d8fUf\ne49f/oW/41//wrf55V/4O378p2647ffcdAfbtgeO7Dmy5w77/D//+7/i3/z7/4aJlomWhpknvOIJ\nr3jKS3oGTuyWPtCz57iM48AdHSOeRCBu+kxTtp7Eh7zJB7zBB7zJh7zJ9XTHl8fv8vbwXd4e36fT\n8Wysd+3u3q3wpKU7Mg1z6ZNtdcbljM8ZlzIuK++6d/gH95N8x/8k/+C+ShK/GfstT/SGt+IHfCl+\nyP/037/L//DvBrITkheyg+yF0bWcZMfJ7RikRxF6PbHLAzs90acBP4EbFT/aVhyoAzxmswqgIGpb\nFcihdiEHcAlcUttGSN4xtYG58UxtIIvQHDPNKRGOieaU0U6ttwqdIoBkcKVLtvMtPYNmQZOgyT4T\ngEahUbRRcnD8/n+Af/ffdUTvSc7RpEhIyXqOiIodD0FUyALZCVmE5Ow7HzM+KSFmfFR0EuuzoBOM\nXcvpquN41XE69Ixdu3mCGnJ2PHl1w9NXtzy5ueXpqxucZpvPMqex9wxXDeNVw3DdEntPF+fSJ7o4\nQxSIoFHQKCCgjUJDuV5h9g1TCMzeencb6W5nutuZ//A/Tvz+7wLtpvuCJWUOKNN72WTZqmHIsOlj\nwZe82XagPdAL9OX/jwpHkCMwA3tgt+kRdC7fzetApP6TC37Nm+3xvM//CNM/rD29tGN86f/8v2l+\n9mfX6/mMmJKPBfvf/M3f5E//9E95//33+epXv8rv//7vM88zAL/zO7/Dr/3ar/HNb36Tn/mZn+Fw\nOPCHf/iHrz2WiFp3ingtD/xrRtVgE5aBa2yyOwVfb6azlwzBngJfdvRQvqnfWlMERZSyVRAQb+eS\nDmhBGsAbQggP9XItm092eDl7waV8qPvZL/V8LJ+w130/7v/k4tiU/xO49zeqD471oaYIueyZcZSp\nw5HxmhHV0kFU8WpCoWVkx4mIp2OgY6QtwiJowmk56nJPWP/W8qxovRfgVHFqAsVlRRJIUkgrwMo6\nxetbT71O+95lICuS67nKfgIqQsYzS8NES0IQN+O94ppISBF1oNleds22n5SJzfWZ0jrPNjaZFcm6\ngENsPZNvmJ1nbgOpcUz+xBBacOWeqJgQieBjudaMHUdBnSe2DWPTMIaG2Qc0OFQFk5pCdxrp3Uiv\nE/080uUJN2WaU2SvI3Hw5OzI2ZOTQ5PQjRP9MNLFCXGKBsi9kDtBOxPGIpk2zoQXCUWKUMr4lHCp\nXGtiAdTUOOadZw6B2XtS65ZnrkkzbZoIU6aZE37OkCGrkL2dN+8FPLhJcbMu2wdhcCtoK9BW0K2Y\nEjkbn8zAyTApBs/UNIzPWqYvNWTnaNNMlybaNNPeGf6JsGBV9kJuSg+GBe5O157Kg1nPPUGeIUaY\nEgwJYl6H/3m0jwX7P/qjP/rYg3zjG9/45Gc0lFix+SGwv/z9FSZNGzZehgrytV2qdJ+g1V0aoCs9\nXB7iP1Oq/v/JXfEJxlpFRy4Tn6tYWcBzBfr6mDoygUjLRMdAwBegN2sgEIu1kF8rYGx4GVesChBc\nERAuZ+QCrNmC/OUh9eJz3afuX+ZCHKg4ogQmaTlJTxKHeIdvlKzRflf2YwZGsya0AQ2lu82Y6rnL\nSeo0zXhGaRl8yxBaUvBMLnJ0e0SyadPMuAxNUtwMUsGpWL658UzSchf23MmOU+iIEszeE5vlp+6G\nZ9zwNN3QzjNBzFrYjSNlWs+1z3gxZy2knZAOawfwd4nmGPFHxZ9WxYFLq6b01DvGtuUkHaemI7WO\nLplA6dNIl0ZkBDdh15pNkCaE5B2xEdQLIWZCysiopqkv2kvZpoueL8a1fUYq8CdgWr9LB8/wRsft\nkz13b+yIvefw8sTVyxPySmlfzau1UZRE7co89Y7YC6pCeJkJIZtQiGpWRcWvBJogZwP5Ko8uH9fP\nsn1qGuf7aq6cuU7W5U3YmIw4TKvvMVBe8F02P5LNTp77wP8aVHNlDH05x76M50yofPbtv/zVr3x+\nB/+Mm6A4smnxBaBdUZ+rEHCSUdFi7SgZRyQw0TDSkfE0RBITl7rYL/1ywN6C8+bItMyY8pQAirUw\n48rfZ7c8AFXT8qDlOZHtS20XZODsTPl1FfyFYrEonkRThFVGaGTCS0KcWZZnQBFNSCDrU6ZyAczZ\n6KG0c6QrRwpmHfkcOcyJqw9PqBf+zb9MPDu+QlsbfyARXIKQSeWYSw8QXOTAiW6eeKo3THPL0fec\nfM8x2LZhLmPPdrGXhnChus4A0V30DntHtr93m+M03AfdfH5cbR2xC0yh4+R2xByQUWgGRU4zYQCZ\nrDPBr/48uJPSxIy/VZqCVDKDTCb8iFTWdmVws+3PWLb5geupxr/f3MeNUAgaOYwnuhubV20c4S4S\n7hLhLtqxZ0zYlGO6RmnajG+VprG5EjU2wz0FDmWcWsYm61DacvoKbZ+XvvjFg32d5Ar2Pfc1gu2D\n5FkoHGmwF45Ko1wC/KWUcOgZ4F9MYx1HxypQqna//PQhOfuA6lL7JxDL/9WvfmVRJH6g7XuMdQvK\n259Vyqe+9CpnXy6/1414OO8mOhSHIvzSrzSrNvnAGM4IKimasci5jK8vbfms9b5i2lOlU7YXohX0\ncwF6VrrLkfAay0sohBzxOeFSNq5Wxfh/b9qmiBp141aqsD7TUrXLABoKJdEJbsqEIRGmjB8zKvDf\n/kvPPCmzd2R/oW0ozK5o7ME8MyJqbLtEekb6NNLoRJsn2jjRu4n9eGI3DTQxGrWyBfct1bWdy8tt\ntXYvhiSVsgplLoLNRw6CLFRboXTKe6b1WMUykiO4W/C3m/uj8Kv/dZm/SXGDPqyl1/FXnIAzmoRT\n+U0VCFXJ3B7nAUvENxkvE13l++ux6uft81rmTxyIV1zQBcHzXtC9kPYm7CQqclLcyyIEBJoyJhcg\nF5bj89Izf3BgX6mTHfcn/4KCl16RFiQYj+vICA5Z7t6We1mfWi09L7qpOwMP3GYcfenVEfTaGd+q\nQdsn8PMyvr7YtmX+AdJmIuzd
       "text": [
        "<matplotlib.figure.Figure at 0x109484bd0>"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Conclusion"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Evaluating IIR filters in matrix form yields an especially fast implementation, especially on modern SIMD processors.\n",
      "The technique is general to many different filter design strategies.\n",
      "\n",
      "Given freedom of filter design strategy, a trapezoidal integrated State Variable Filter has many advantages over\n",
      "the traditional biquad formulas. It is both more numerically robust and far more conducive to modulation, even up to\n",
      "audio rate. It is as versatile as the traditional biquad element, and yields the same menu of \"EQ cookbook\" responses\n",
      "without especially difficult math."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    }
   ],
   "metadata": {}
  }
 ]
}