You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
music-synthesizer-for-android/lab/biquad in two.ipynb

305 lines
79 KiB

{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fast matrix biquad\n",
"==================\n",
"\n",
"This notebook shows the math behind a very fast NEON implementation of biquad IIR filters.\n",
"One of the main primitives NEON is very highly optimized for is multiplying a matrix by\n",
"a vector, especially of size 4.\n",
"\n",
"A biquad filter has two state variables (this is most apparent in transposed direct form II).\n",
"For this implementation, the input vector consists of two input values and the two state\n",
"variables. The output vector is two output values and the new state vector. The relationship\n",
"is linear, ie the transformation from input to output vector is simply a matrix-vector\n",
"multiplication. Each iteration computes two samples, using a single 4x4 matrix-vector multiply."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%pylab inline"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"prompt_number": 117
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start with a simple implementation of biquad for reference; we'll compare the error of the matrix approach against this scalar implementation. This is in 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"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 118
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a starting step, here's the same biquad in a 2x2 matrix formulation. This version is derived from\n",
"transposed direct form II, with a little expansion of the effect _this_ input sample has on the _next_\n",
"state vector (this is the B array below). The A matrix describes the evolution of the IIR state."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def bqmat(coefs, inp):\n",
" b0, b1, b2, a1, a2 = coefs\n",
" A = array([[-a1, 1], [-a2, 0]])\n",
" B = array([b1 - a1 * b0, b2 - a2 * b0])\n",
" y = array([0, 0])\n",
" out = np.zeros(len(inp))\n",
" for i in range(len(inp)):\n",
" x = inp[i]\n",
" out[i] = y[0] + x * b0\n",
" y = dot(A, y) + B * x\n",
" return out"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 119
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's test it out."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"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(bqmat(coefs, impulse))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 120,
"text": [
"[<matplotlib.lines.Line2D at 0x10cecce10>]"
]
},
{
"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 0x10ced4110>"
]
}
],
"prompt_number": 120
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(biquad(coefs, impulse) - bqmat(coefs, impulse))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 121,
"text": [
"[<matplotlib.lines.Line2D at 0x10d3a9e50>]"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtclGXeP/DPJJjHRFNBGQzjIEcBUylLwxQRFDO1x0Nb\nrvoqNjOzbbdfu+3zpJvHdVtTKTvqY7WhZZpmSpo62qqEBh5BQZIcQFAxUjwL1++P7zPIaYaZuY/M\nfN+vFy9l5p7rvriV+3tf3+tkEEIIMMYYczt3aV0Bxhhj2uAAwBhjbooDAGOMuSkOAIwx5qY4ADDG\nmJviAMAYY25KNwFg6tSp8Pb2RmRkpCzlDR8+HB07dkRycnKD915//XX06tULYWFhWL58uSznY4yx\n5kY3AWDKlClIT0+XrbxXX30Vn376aYPXV61aheLiYpw8eRI5OTmYMGGCbOdkjLHmRDcBYODAgejY\nsWOd1woKCpCYmIi+ffti0KBBOHnypN3lPfbYY2jXrl2D19977z38z//8T833Xbp0cb7SjDHWjOkm\nADTmueeew/Lly3Hw4EEsXrwY06dPl1xmQUEB1qxZg379+iEpKQmnTp2SoaaMMdb8eGhdAWsqKyux\nf/9+PPnkkzWv3bx5EwCwfv16vPHGGw0+YzQasXXrVpvl3rhxA61bt8aBAwewYcMGTJ06FXv27JG3\n8owx1gzoNgBUV1fDy8sL2dnZDd4bM2YMxowZ02QZBoOhwWtGo7Hms6NHj8aUKVOkV5YxxpohSSkg\ns9mMwYMHIzw8HBEREVi2bFmDY0wmEzp06ICYmBjExMRg7ty5dpV9zz33oGfPnli3bh0AQAiBI0eO\nOFS/xta5Gz16NHbu3AkA2L17N3r16uVQmYwx5jKEBGfPnhXZ2dlCCCEuX74sgoODRU5OTp1jdu3a\nJZKTk5ssa8KECaJbt27C09NTGI1GsXLlSnH69GkxfPhwERUVJcLCwsSbb75pd90eeeQR0aVLF9G6\ndWthNBrFtm3bhBBCVFRUiBEjRojIyEgxYMAAceTIEQd+YsYYcx2SUkA+Pj7w8fEBALRr1w6hoaEo\nKSlBaGho/SDTZFlpaWmNvt5UTt+aH374odHXO3TogM2bNztVJmOMuRLZRgEVFhYiOzsbsbGxdV43\nGAzYt28foqKikJSUhJycHLlOyRhjTAJZOoErKysxbtw4LF26tMHY+z59+sBsNqNNmzbYunUrRo8e\njby8PDlOyxhjTAqpOaSbN2+KYcOGiSVLlth1vL+/vygvL2/wekBAgADAX/zFX/zFXw58BQQEOH3/\nlpQCEkJg2rRpCAsLw6xZsxo9pqysrKYPIDMzE0IIdOrUqcFxBQUFEELwlxB44403NK+DXr74WvC1\n4Gth+6ugoMDpe7ikFNDevXvx2WefoXfv3oiJiQEAzJ8/H2fOnAEApKSkYN26dVixYgU8PDzQpk0b\nrFmzRsopGWOMyURSAHjkkUdQXV1t85gXXngBL7zwgpTTMMYYU4Cu1wJyV3FxcVpXQTf4WtzB1+IO\nvhbyMAghhNaVAGi4qE6qwhhjzYaUeye3ABhjzE1xAGCMMTfFAYAxxtwUBwDGGHNTHAAYY8xNcQBg\njDE3xQGAMcbcFAcAxhhzUxwAGGPMTXEAYIwxN8UBQKIVK4Dycq1rwRhjjuO1gCS4fh1o2xaIjgZ2\n7AC8vLSuEWPM3fBaQBrJywN69QIGDgSGDwcuXdK6RowxZj8OABLk5ABhYcCSJUCfPsCIEUBlpda1\nYowx+0gKAGazGYMHD0Z4eDgiIiKwbNmyRo+bOXMmgoKCEBUVhezsbCmn1JXcXCA0FDAYgNRUoEsX\n6hNgjLHmQFIA8PT0xJIlS3D8+HFkZGTgnXfeQW5ubp1jtmzZglOnTiE/Px8ffPABnn/+eUkV1hNL\nCwAA7rqLWgDHjmlbJ8YYs5ekAODj44Po6GgAQLt27RAaGoqSkpI6x2zatAmTJ08GAMTGxqKiogJl\nZWVSTqsblhaARVgYvcYYY82BbH0AhYWFyM7ORmxsbJ3Xi4uL4efnV/O90WhEUVGRXKfVzO3bQEEB\ndQJbhIZSAGhmg5kYY25K0qbwFpWVlRg3bhyWLl2Kdu3aNXi//hAlg8HQaDmzZ8+u+XtcXJyu9/0s\nKAC6dwdat77zmpcX0L49UFQE1Ip5jDEmG5PJBJPJJEtZkgPArVu3MHbsWPzud7/D6NGjG7zv6+sL\ns9lc831RURF8fX0bLat2ANC72vn/2kJD6T0OAIwxJdR/OJ4zZ47TZUlKAQkhMG3aNISFhWHWrFmN\nHjNq1Ch88sknAICMjAx4eXnB29tbyml1oX7+34L7ARhjzYWkFsDevXvx2WefoXfv3oiJiQEAzJ8/\nH2fOnAEApKSkICkpCVu2bEFgYCDatm2LVatWSa+1DuTkAEOHNnw9NBQ4dEj9+jDGmKN4KQgnPfAA\njfnv37/u6yYT8N//DfzwgybVYoy5GSn3Tg4ATqiups7es2eBe+6p+15ZGbUCystpghhjjCmJ1wJS\n2ZkzQKdODW/+ANC1K934z51Tv16MMeYIDgBOyMlpvAMYoJu/ZT4AY4zpGQcAJ+TmNj4E1CIsjIIE\nY4zpGQcAJ9hqAQDcAmCMNQ8cAJzALQDGmCvgAOAgIbgFwBhzDRwAHFRaCrRsCXTubP0YPz/g8mWg\nokK9ejHGmKM4ADjoxIm6K4A2xmAAQkK4FcAY0zcOAA46exawspZdHdwPwBjTOw4ADjp3DrBnLTvu\nB2CM6R0HAAeVldFs36ZwC4AxpnccABxkbwvg/vuBwkLFq8MYY07jAOCgsjL7AkC3btRfwBhjesUB\nwEH2poA6dQKuXgWuXVO+Towx5gwOAA6yNwVkMAA+PjRvgDHG9EhyAJg6dSq8vb0RGRnZ6Psmkwkd\nOnRATEwMYmJiMHfuXKmn1IwQ9rcAAAoAnAZijOmV5E3hp0yZghdffBHPPPOM1WMeffRRbNq0Seqp\nNHfpEuDpCbRpY9/x3bpxC4Axpl+SWwADBw5Ex44dbR7TXHb6aoq96R8L7ghmjOmZ4n0ABoMB+/bt\nQ1RUFJKSkpDTjAfHO5L+ATgAMMb0TXIKqCl9+vSB2WxGmzZtsHXrVowePRp5eXmNHjt79uyav8fF\nxSEuLk7p6jnEmRZARoZy9WGMuR+TyQSTySRLWbJsCl9YWIjk5GQcPXq0yWN79uyJn376CZ06dapb\nkWawKfyKFcDhw8B779l3/ObNwLvvAlu2KFsvxpj70vWm8GVlZTWVy8zMhBCiwc2/ueAUEGPMlUhO\nAU2cOBG7d+/GhQsX4Ofnhzlz5uDWrVsAgJSUFKxbtw4rVqyAh4cH2rRpgzVr1kiutFbOnQPCw+0/\nngMAY0zPZEkByaE5pIDGjgUmTACefNK+42/fBlq3ptnAHor3tjDG3JGuU0CuxNFOYA8PWhLi3Dnl\n6sQYY87iAOAAexeCq40ngzHG9IoDgAMc7QQGuB+AMaZfnJm20/Xr9OXl5djnmksAyM8Hak/PiI62\nb+tLxljzxQHATufO0dO/weDY55pLAJg4EWjXDmjbFrh8Gbh1C9i/X+taMcaUxAHATs6kfwAKAHpf\n/eL8eWoBnD8PtGwJ3LxJ9S4qAoxGrWvHGFMK9wHYydERQBbNoQWwfTsweDDd/AH6c+RIYMMGbevF\nGFMWBwA7SWkB6D0ApKcDCQl1Xxs7FvjqK23qwxhTBwcAO7lqC6C6Gti2rWEAGDYMOHSI5zAw5so4\nANjJmTkAwJ1tIfU6yfnIEaB9e+D+++u+3qoVMHw48PXX2tSLMaY8DgB2cjYF1Lo1ff36q/x1ksN3\n3zV8+rfgNBBjro0DgJ2cTQEB+k4D2QoAiYm0n4FegxdjTBoOAHZyNgUE6DcAVFYCBw7QCKDGtGsH\nPPYY4ALbOTPGGsEBwE7OpoAA/QYAkwno149u9NZwGogx18UBwA5VVZQG6dzZuc/rNQA0NvyzvuRk\nYOdOWgaDMeZaOADY4cIFWgPI
"text": [
"<matplotlib.figure.Figure at 0x10d390750>"
]
}
],
"prompt_number": 121
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the scale of the error -- it's floating point rounding.\n",
"\n",
"Now for the 4x4 matrix approach. It's really four 2x2 matrices tiled. The top\n",
"left corner is the influence of the two input samples on the two output samples;\n",
"it's the first two values of the biquad's impulse response, banded. The top\n",
"right corner is the contribution of the state vector to the output. The bottom\n",
"left is the effect of the input on the state, and the bottom right is the\n",
"evolution of the IIR. Since each iteration is two samples, this matrix is\n",
"simply the square of the A matrix above.\n",
"\n",
"Note that each iteration processes _two_ samples. The matrix-vector product\n",
"is 16 multiply-accumulates, so it amortizes to 8 per sample. This compares\n",
"well to the 5 multiply-accumulates of the direct form, but organized for extremely\n",
"fast evaluation in SIMD, and NEON in particular."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def biquadm(coefs, inp):\n",
" b0, b1, b2, a1, a2 = coefs\n",
" c1 = b1 - a1 * b0\n",
" c2 = b2 - a2 * b0\n",
" A = array([[b0, 0, 1, 0],\n",
" [c1, b0, -a1, 1],\n",
" [-a1 * c1 + c2, c1, -a2 + a1*a1, -a1],\n",
" [-a2 * c1, c2, a1*a2, -a2]])\n",
" out = np.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": 122
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, let's test it out."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(biquad(coefs, impulse))\n",
"plot(biquadm(coefs, impulse))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 123,
"text": [
"[<matplotlib.lines.Line2D at 0x10d3b9510>]"
]
},
{
"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 0x10d3b94d0>"
]
}
],
"prompt_number": 123
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The error is similar to that of the 2x2 matrix version, but not quite the same."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(biquadm(coefs, delta) - biquad(coefs, delta))\n",
"plot(biquadm(coefs, delta) - bqmat(coefs, delta))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 124,
"text": [
"[<matplotlib.lines.Line2D at 0x10d70b550>]"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVNX7B/DPCLibgguyubEIKCCKkjsuaKKQW4lWllpZ\nZtm3srIstRTlZ1kupWVibqlpuaXydR0zN1RwhdxRFsUFUXFjO78/zncQZLa7zNwL87xfL17KzJ1z\nD1e5zz3bczSMMQZCCCE2p5LSFSCEEKIMCgCEEGKjKAAQQoiNogBACCE2igIAIYTYKAoAhBBio1QT\nAEaOHAlnZ2cEBATIUt5zzz0HR0dHREZGlnnv888/R/PmzeHv74+5c+fKcj5CCClvVBMARowYgfj4\neNnK+/jjj7Fs2bIyry9evBgZGRk4c+YMkpOTER0dLds5CSGkPFFNAOjcuTMcHR1LvXbhwgX06dMH\nISEh6NKlC86cOWN2ed27d0fNmjXLvL5gwQJ8+eWXxd/Xr19ffKUJIaQcU00A0OfNN9/E3LlzceTI\nEcycORNjxoyRXOaFCxewatUqtG3bFhERETh//rwMNSWEkPLHXukKGJKbm4sDBw7ghRdeKH4tLy8P\nAPDnn39i0qRJZT7j7u6OrVu3Gi338ePHqFatGg4fPox169Zh5MiR+Pvvv+WtPCGElAOqDQBFRUWo\nU6cOkpKSyrw3cOBADBw40GQZGo2mzGvu7u7Fn+3fvz9GjBghvbKEEFIOSeoCSktLQ7du3dCiRQu0\nbNkSc+bM0Xvce++9B29vbwQFBem9oevzzDPPoGnTpli7di0AgDGGEydOCKqfvjx3/fv3x65duwAA\ne/bsQfPmzQWVSQghFQaT4OrVqywpKYkxxti9e/eYj48PS05OLnXM5s2bWZ8+fRhjjB08eJCFhobq\nLSs6Opq5uLgwBwcH5u7uzuLi4tilS5fYc889x4KCgpi/vz/7+uuvza5bp06dWP369Vm1atWYu7s7\n27ZtG2OMsZycHNa3b18WEBDAOnTowE6cOCHmRyeEkHJPw5h86aD79++Pd999Fz169Ch+7a233kK3\nbt0wZMgQAICvry/27NkDZ2dnuU5LCCFEBNlmAaWmpiIpKQmhoaGlXs/IyICHh0fx9+7u7khPT5fr\ntIQQQkSSJQDk5uZi8ODBmD17tt659083MvQNzhJCCLEuybOA8vPzMWjQILz88svo379/mffd3NyQ\nlpZW/H16ejrc3NzKHOfl5YULFy5IrQ4hhNgUT09P0euZJLUAGGMYNWoU/P398f777+s9JioqCkuX\nLgUAHDx4EHXq1NHb/3/hwgUwxuiLMUyaNEnxOqjli64FXQu6Fsa/pDw4S2oB7Nu3D8uXL0dgYCCC\ng4MBADExMbhy5QoAYPTo0YiIiMCWLVvg5eWFGjVqYPHixVJOSQghRCaSAkCnTp1QVFRk8rh58+ZJ\nOQ0hhBALUHUuIFsVFhamdBVUg67FE3QtnqBrIQ9Z1wFIodFooJKqEEJIuSHl3kktAEIIsVEUAAgh\nxEZRACCEEBtFAYAQQmwUBQBCCLFRFAAIIcRGUQAghBAbRQGAEEJsFAUAQgixURQACCHERlEAIIQQ\nG0UBgBBCbBQFAEIIsVEUAMqJvXuBF14AKGEqIUQukgPAyJEj4ezsjICAAL3va7Va1K5dG8HBwQgO\nDsbUqVOlntLmHDgADBwI7N4NnDqldG0IIRWF5E3hR4wYgXfffRfDhw83eEzXrl2xceNGqaeySUeO\nAM8/DyxZAuzYAfzxB2Ag1hJCiCCSWwCdO3eGo6Oj0WNooxdxTp8G+vUDFi4EIiKAQYN4ACCEEDlY\nfAxAo9Fg//79CAoKQkREBJKTky19ygrjp5+AMWN4CwAA2rcHbt0Czp5Vtl6EkIpBcheQKa1bt0Za\nWhqqV6+OrVu3on///jhr4A42efLk4r+HhYXZ/L6fiYnA118/+b5SJWDAAN4KmDBBuXoRQpSj1Wqh\n1WplKUuWPYFTU1MRGRmJkydPmjy2adOmOHr0KJycnEpXhPYELqWwEKhTB0hL43/q7NoFfPwxHxsg\nhBBV7wmclZVVXLmEhAQwxsrc/ElZ584B9euXvvkDQJcuwOXLQGqqItUihFQgkruAhg4dij179uDm\nzZvw8PDAlClTkJ+fDwAYPXo01q5di/nz58Pe3h7Vq1fHqlWrJFfaFiQlAa1bl33d3p6PCfz5J/DB\nB9avFyGk4pClC0gO1AVU2vjxgKMj8NlnZd/buhWYOhXYt8/69SKEqIuqu4CIOImJ+lsAANCjB5CS\nAmRmWrdOhJCKhQKACjHGA0BwsP73K1cGOnQAEhKsWy9CSMVCAUCFUlOB6tUBZ2fDx/j5AbSkghAi\nBQUAFTI0AFySvz/vBiKEELEoAKiQsf5/HWoBEEKkogCgQsb6/3X8/IB//wWKiqxTJ0JIxUMBQGUY\nA44eNd0CqF2bLxK7csU69SKEVDwUAFTm6lWeBsLDw/SxNA5ACJGCAoDK6AaANRrTx9I4ACFECgoA\nKmNO/78OtQAIIVJQAFAZc2YA6fj5lY8AUFREq5YJUSMKACpz+rT5Wz76+/MuILWnUPr9dz6m8dpr\nwMWLSteGEKJDAUBFGAPS080bAAZ4umg7OyAry7L1kmrNGmDWLKBpU6BdO+Dtt4GCAqVrRQihAKAi\nt28DDg5ArVrmf0bXClCr+/f5ZvYvvwxMmgScOQPs3UuZTAlRAwoAKpKeDri7C/uM2scB4uP5U3/d\nuvz7unX5fgb//a+y9SKEUABQlYwM4QFA7S2AP/4ABg0q/dpzz1EAIEQNJAeAkSNHwtnZGQFGRi7f\ne+89eHt7IygoCElJSVJPWWGlpwNubsI+o+YWwOPHfPOa/v1Lv/7ss8CFC+ofuyCkopMcAEaMGIH4\n+HiD72/ZsgXnz5/HuXPn8PPPP+Ptt9+WesoKS0wXkJpbANu38xlNDRuWft3BAejWjb9PCFGO5ADQ\nuXNnODo6Gnx/48aNePXVVwEAoaGhyMnJQRY9+uklpgvIzQ148ADIzrZMnaTQ1/2j07s3dQMRojSL\njwFkZGTAo8S8Rnd3d6Snp1v6tOWSmBaARqPObqD8fGDTJmDgQP3v9+4NbNtG2UwJUZK9NU7y9IbF\nGgOJbiZPnlz897CwMISFhVmwVuojZgwAeBIAOnaUv05iabWAl5fhNQ1Nm/KMpsePm5/6ghACaLVa\naLVaWcqyeABwc3NDWlpa8ffp6elwM3CXKxkAbJGYLiBAnS2AjRuBAQOMH6PrBqIAQIj5nn44njJl\niuiyLN4FFBUVhaVLlwIADh48iDp16sDZ2Ga3Nio3F3j0CHByEv7Zxo2BEjFWFY4fB9q2NX4MTQcl\nRFmSWwBDhw7Fnj17cPPmTXh4eGDKlCnIz88HAIwePRoRERHYsmULvLy8UKNGDSxevFhypSuijAze\n/WNOGuinubvz7iM1SUnhM5SMCQsDoqOBe/eErX4mhMhDcgBYuXKlyWPmzZsn9TQVntjuH4AHjowM\neesjxY0bPNePqYZejRp8lfDu3UBUlHXqRgh5glYCq4SYGUA6rq5PdhJTA93Tvzmtmc6dgYQEy9eJ\nEFIWBQCVEDsDCACqVAEcHYHr1+Wtk1jJyXxg2hxqXshGSEVHAUAlpHQBAerqBjKn/19HjTOYCLEV\nFABUQkoXEKCugWAhLQAfH+DSJSAvz7J1kkturtI1IEQ+FABUQkoXEKC+AGBuC6BKFaBRI+D8ecvW\nSQ5xcTydNe1lQCoKCgAqUVG6gO7c4V/m7moGlI9xgOXLgS++AGbOBIYMoUympGKwSioIYlxeHk/m\nJmV9nLs733lLaSkpgK8vUEnAo4XaxwF+/x0YPx7YuZMHq5s3+fqF7dsBe/oNIuUYtQBUIDOTp0y2\nsxNfhlq6gFJSzO//11FzC2DvXuC99/iKZV231qRJQOXKwGefKVs3QqSiAKACulXAUri5qSMACOn/\n11FzC+Cnn3jXT2Dgk9fs7IDf
"text": [
"<matplotlib.figure.Figure at 0x10d70b5d0>"
]
}
],
"prompt_number": 124
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 124
}
],
"metadata": {}
}
]
}