diff --git a/config b/config index f7f2910..53a80d6 100644 --- a/config +++ b/config @@ -33,4 +33,5 @@ source "scripts/paths" # Project-specific variables below. [ -z "${CC}" ] && CC="gcc" +[ -z "${CXX}" ] && CXX="g++" [ -z "${CFLAGS}" ] && CFLAGS="-g -O2 -W -Wall" diff --git a/python/.gitignore b/python/.gitignore new file mode 100644 index 0000000..7c28a64 --- /dev/null +++ b/python/.gitignore @@ -0,0 +1,5 @@ +.ipynb_checkpoints/ +__pycache__/ +build/ +iir.py +pyiir_wrap.cxx diff --git a/python/README b/python/README new file mode 100644 index 0000000..0b933c6 --- /dev/null +++ b/python/README @@ -0,0 +1,25 @@ +Experimental Python bindings +============================ + +These bindings are currently experimental, as it is not yet clear how best +to create and maintain them. They are built upon the C++ interface as the +object oriented nature of it sits much better with Python. + +At present, it is necessary to compile and install the C++ library before +working with the bindings. + +To generate the bindings, swig 3.0.2 has been used, but earlier versions +may work too: + + swig -c++ -python -py3 pyiir.i + +To build the extension, distutils is used: + + python3 setup.py build_ext + +Note that the -L ${LIBDIR} option may be handy if libiir++ is not installed +into the system library path. + +To install the extension: + + python3 setup.py install diff --git a/python/filter_analysis.ipynb b/python/filter_analysis.ipynb new file mode 100644 index 0000000..bf64d79 --- /dev/null +++ b/python/filter_analysis.ipynb @@ -0,0 +1,325 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:1cb8a128cf0baef93386ff39d25370890fe1a95c2ff36083b1a83c2cad364c1c" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "heading", + "level": 1, + "metadata": {}, + "source": [ + "Filter analysis laboratory" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First configure the filter as desired. The filter is independent of sample rate, but is specified as a ratio to it. We use a specific value here simply to aid comprehension and label all the graph axes as desired. Setting it to 1 would be the same as specifying everything as a ratio." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import iir\n", + "\n", + "# change this to alter the sample rate we are displaying at\n", + "sample_rate = 100 # Hz\n", + "nyquist_freq = sample_rate / 2\n", + "\n", + "# build the filter object, fi, here\n", + "corner_freq = 10 # Hz\n", + "order = 4\n", + "gain = 1.0\n", + "filter_desc = \"butterworth_lowpass(%d, %f, %f)\" % (order, gain, corner_freq / nyquist_freq)\n", + "print(filter_desc)\n", + "fi = iir.Filter(filter_desc)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "butterworth_lowpass(4, 1.000000, 0.200000)\n" + ] + } + ], + "prompt_number": 65 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we set up the plotting library and helper functions." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pylab\n", + "\n", + "# number of points to plot on graphs (arbitrary)\n", + "num_points = 1000\n", + "\n", + "# helper function: calls the Filter.response() method on each member of the array\n", + "respfunc = pylab.vectorize(lambda f: fi.response(f/sample_rate))" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 66 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Filter characteristics" + ] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "-3dB points" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "halfgain = np.array(fi.responseMinus3dB(gain)) * sample_rate\n", + "print(halfgain)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[ 9.99999995]\n" + ] + } + ], + "prompt_number": 67 + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Impulse response" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# start at t=0, impulse at t=tau, end at t=10\u00d7tau\n", + "tau = sample_rate * (1 / corner_freq) # number of samples\n", + "\n", + "x = pylab.hstack([pylab.zeros(tau), pylab.ones(9*tau)])\n", + "t = pylab.linspace(0, tau / 10, x.size)\n", + "fi2 = iir.Filter(fi)\n", + "\n", + "y = np.array(x.size, dtype=np.double)\n", + "y.resize(x.size)\n", + "for i in range(0, y.size):\n", + " y[i] = fi2.f(x[i])\n", + "\n", + "fig, axes = plt.subplots(figsize=(10,6), dpi=100)\n", + "axes.plot(t, y, 'r')\n", + "axes.plot(t, x)\n", + "axes.set_xlabel('Time (s)')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 68, + "text": [ + "" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF/CAYAAACL7od8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9lJREFUeJzt3XuQXOV55/FvayQkEBIIBOImkMGAuaxiMGDADgzB2AIX\nJmW2oGScxJckSlVIqFBrQ+LaeFwpdssp7HgDLpZlweVgA0nFToKNgCSQCWALAwYErBAlCYRGEhch\nZCFssHo0vX+8PahpZqa7z3nPpbu/n6qu6Zk+OufBBzE/P+/lgCRJkiRJkiRJkiRJkiRJkiRJkiRJ\nkiRJkprcArwCPD3J55cBK4GngJ8Ai3OqS5Ikqev8JnASkwerM4B96u+XAA/nUZQkSVK3WsTkwarR\nPGBjtqVIkiSV07TI5/sisDzyOSVJknrKIlp3rM4BVhG6VpIkSX1neqTzLAZuIsyx2jbRAUcddVRt\n3bp1kS4nSZKUqXXA+zv9QzGC1eHAD4HPAmsnO2jdunXUarUIl1MRhoaGGBoaKroMJeC9627ev+7l\nvetulUrlqCR/rp1gdTtwNjAfGAG+Csyof3Yj8JeE4b8b6j+rAqclKUaSJKmbtROslrb4/PfrL0mS\npL4We1WgetTg4GDRJSgh71138/51L+9df6rkeK2ac6wkSVI3qFQqkCAn2bGSJEmKxGAlSZIUicFK\nkiQpEoOVJElSJAYrSZKkSAxWkiRJkRisJEmSIjFYSZIkRWKwkiRJisRgJUmSFInBSpIkKRKDlSRJ\nUiQGK0mSpEgMVpIkSZEYrCRJkiIxWEmSJEVisJIkSYrEYCVJkhSJwUqSJCkSg5UkSVIkBitJkqRI\nDFaSJEmRGKwkSZIiMVhJkiRFYrCSJEmKxGAlSZIUicFKkiQpEoOVJElSJAYrSZKkSKYXXYBK5rnn\n4MEH4YQTwmvu3KIrkiSpa1RyvFatVqvleDl17NVX4ZRT4PTTYf16WLUK9tsPTjwRFi+GK6+EAw8s\nukpJkjJXqVQgQU4yWCmoVuG88+AjH4Frrgk/GxsLAeuZZ+Bf/xWWL4e77oLjjiu0VEmSsmawUjp/\n9mewejX8+McwMDDxMX/3d/ClL8Htt8Nv/Va+9UmSlKOkwcrJ64LbboM77wxfJwtVAL/7u/D3fw9L\nl8Itt+RXnyRJXcKOVb978skwBHjffWEeVTueew4++Um49FL4q7+CaeZzSVJvsWOlzm3dCp/+NFx3\nXfuhCuDYY2HFChgehs9+NszPkiRJdqz61tgYnH9+WPH3jW8kO8fbb8PFF8Ps2fD978OMGXFrlCSp\nIE5eV2dWrIDPfz6s+JueYjuzX/86dL0MV5KkHuJQoDqzfDn89m+nC1UAM2fCD38Iv/wlXHaZw4KS\npL5msOpXy5fDBRfEOZfhSpIkwKHA/vTyy2GTz1dfjTt057CgJKlHOBSo9t1zT9hiIXbwGe9c/epX\ncNFFsH173PNLklRyBqt+tHx5WBGYhZkz4Z/+CRYtgjPOgHXrsrmOJEkl5FBgv6lWw4OUV62Cgw/O\n9lrf/nbYQPSOO2BwMNtrSZIUkUOBas+KFXDkkdmHKoA//mO49Va45BK46absrydJUsHaCVa3AK8A\nT09xzN8Ca4CVwEkR6lJWYq4GbMd558FDD8G118IVV4RNRSVJ6lHtBKvvAEum+PwC4P3A0cAfAjdE\nqEtZyTtYARxzDDz8MGzcGB6H873vhZ3fJUnqMe0EqweBbVN8/ingu/X3PwP2BRakrEtZGBmBzZvh\ntNPyv/a8efCDH4Shweuug1NPhfvvz78OSZIyFGOO1aHASMP3G4HDIpxXsd1zD3ziEzAwUFwNZ50V\nuldXXQV/8Aehe/bYY+DCBklSD0j5PJN3NM+a76vfklu2wLJlXbDh+KMfhEPOhwuLLqQCXAIfuBjW\nvwhnr4dd/w4HHgAHHAgHzIcZe8S/7Ngu2LkTdlZhtArV0Yav9deuXeG4XWO734/VoDYWwt9YLXyt\n1YBaw7/pzf/KVxq+VHZ/bfwZQGWi4yY7doLzt/hRYn31N1iS4okRrDYBCxu+P6z+s/cYGhp65/3g\n4CCDPbIEf/16eOop+Ju/KbqSKVSrcN+18I2bYG7RxYwbAI4Mr82b4fHH4fE74YFn4PAjwl5Yc+fC\nvvvCPvuE15w5YX5Wtfru19tvw44d8MYb8Oabu9/v2LH7fbUa/vycOWF3+L32gr3rX2fPhlmzYOYe\nYS+uPWbCHnuE14wZMG1a6PSNv6ZNawhDld0BaVxz+KrVJn5B+OdpfP+uP9/0PvzgvW8bP5+q+9dc\n51Q6OVaSutzT63/O0y8+vvsHryY7T7v/5VwE/Aj4LxN8dgFwef3r6cC36l+b9ew+Vj/9KVx5ZRjh\nKq377oOvfKXkRda9/XZYSbh2bXjsTuPr9dfDg6Nnznz3a/Zs2H//d7/22w/mzw/v588PgcqwIElq\nQ9J9rNrpWN0OnA3MJ8yl+iow/iyUG4HlhFC1Fvgl8PlOi+h21WoXPBaviNWASc2aBR/7WHhJktRF\n2glWS9s45vK0hXSzajWMGJXa8uVhRZ4kScqMO69HUPqO1fPPhyG0k08uuhJJknqawSqC0geru+8O\nD12e5u2WJClL/qaNYOfOkgere++FJVNtni9JkmIwWEVQ+o7VypVwyilFVyFJUs8zWEVQ6snrb74Z\ndjB93/uKrkSSpJ5nsIqg1B2r1avh6KOLfYyNJEl9wmAVQamD1bPPwvHHF12FJEl9wWAVQemD1XHH\nFV2FJEl9wWAVQalXBa5aZbCSJCknBqsISt+xcihQkqRcGKwiKO2qwF//Gl58MUxelyRJmTNYRVDa\njtWaNbBoUUlTnyRJvcdgFUFpg5XzqyRJypXBKoLSBivnV0mSlCuDVQSlXRVox0qSpFwZrCKwYyVJ\nksBgFUUpVwWOjobJ68ceW3QlkiT1DYNVBKXsWL3wAixYALNnF12JJEl9w2AVQSmDlcOAkiTlzmAV\nQSknrztxXZKk3BmsIrBjJUmSwGAVRSmDlR0rSZJyZ7CKoHSrAms1WL3aYCVJUs4MVhGUrmO1cSPs\nvTfMm1d0JZIk9RWDVQSlC1arVjm/SpKkAhisIijdqsBnn3UYUJKkAhisIrBjJUmSwGAVRekmr9ux\nkiSpEAarCErVsarV3GpBkqSCGKwiKFWw2rIlhKsFC4quRJKkvmOwiqBUwWq8W1WpFF2JJEl9x2AV\nQalWBfooG0mSCmOwiqBUHSsnrkuSVBiDVQSlWhXoVguSJBXGYBWBHStJkgQGqyhKE6y2bw+vhQuL\nrkSSpL5ksEpp167wdWCg2DqA0K36wAdgmrdVkqQi+Bs4pVKtCFy9OgQrSZJUCINVSqUZBgTYsAGO\nOKLoKiRJ6lsGq5RKtSJwZMT5VZIkFchglVKpOlYjI3D44UVXIUlS3zJYpVSqYLVhgx0rSZIKZLBK\nqTTBqlZzKFCSpIIZrFIqzarA7dvDg5f32afoSiRJ6lsGq5RK07Ea71ZVKkVXIklS3zJYpVSaVYHO\nr5IkqXAGq5RK17GSJEmFaSdYLQFWA2uAqyb4fD5wD/Ak8AzwuVjFdYNSBSu3WpAkqVCtgtUAcD0h\nXB0PLAWOazrmcuAJ4IPAIPANYHrUKkusNMHKoUBJkgrXKlidBqwF1gNV4A7goqZjXgLm1t/PBbYC\no/FKLLfSrAp0KFCSpMK16iwdCow0fL8R+HDTMTcB9wObgTnAJdGq6wKlmbzuUKAkSYVr1bGqtXGO\nvyDMrzqEMBz4bULA6gulGAocG4NNm+CwwwouRJKk/taqY7UJaBxfWkjoWjU6E7im/n4d8AJwLPBY\n88mGhobeeT84OMjg4GBHxZZRKYLVq6/CnDmw554FFyJJUncaHh5meHg49Xla7SY5HXgOOJcw1PcI\nYQL7sw3HfBPYDnwNWAD8HFgMvN50rlqt1k4DrLt8//tw111w220FFvHoo7BsGTz+eIFFSJLUOyph\nw+2Od91u1bEaJaz6u5ewQvBmQqhaVv/8RuB/AN8BVhKGFr/Me0NVzypFx8r5VZIklUI72yLcXX81\nurHh/WvAhdEq6jKlWBXoikBJkkrBnddTKsWqQPewkiSpFAxWKTkUKEmSxhmsUipNsLJjJUlS4QxW\nKZUiWDkUKElSKRisUio8WFWrsGULHHJIgUVIkiQwWKVW+KrAzZthwQKY3jfPvZYkqbQMVikVvirQ\n+VWSJJWGwSqlwocCnV8lSVJpGKxSKjxYudWCJEmlYbBKqRTByo6VJEmlYLBKqfDJ6wYrSZJKw2CV\nUuEdqw0bHAqUJKkkDFYpuSpQkiSNM1ilVGjH6le/gjffhAMOKKgASZLUyGCVUqHBauNGOOwwqFQK\nKkCSJDUyWKVUaLByfpUkSaVisEqp0FWBzq+SJKlUDFYpFTp53WAlSVKpGKxScihQkiSNM1ilVGiw\nsmMlSVKpGKxSMlhJkqRxBquUCgtWtZrBSpKkkjFYpVTYqsBf/AKmTYN99ing4pIkaSIGq5QKWxVo\nt0qSpNIxWKVU2FCgwUqSpNIxWKVUaLByqwVJkkrFYJVSYcFqwwY7VpIklYzBKiWHAiVJ0jiDVUqF\nrQo0WEmSVDoGqxRqNdi1q8ChQOdYSZJUKpUcr1Wr1Wo5Xi57O3fC7NlhODBXY2Ow556wfTvMmpXz\nxSVJ6n2VSgUS5CQ7VikUNr9q61bYe29DlSRJJWOwSqGwYPXSS3DwwQVcWJIkTcVglYLBSpIkNTJY\npVDYikCDlSRJpWSwSqGw5wQarCRJKiWDVQoOBUqSpEYGqxQMVpIkqZHBKgWDlSRJamSwSsFgJUmS\nGhmsUti5s4DJ67WawUqSpJIyWKVQSMdqxw4YGAg7r0uSpFIxWKVQSLCyWyVJUmkZrFIwWEmSpEYG\nqxQMVpIkqZHBKgWDlSRJamSwSqGQVYEGK0mSSqudYLUEWA2sAa6a5JhB4AngGWA4RmHdwI6VJElq\nNL3F5wPA9cDHgE3Ao8CdwLMNx+wLfBv4BLARmB+/zHIyWEmSpEatOlanAWuB9UAVuAO4qOmYzwA/\nIIQqgNci1ldqBitJktSoVbA6FBhp+H5j/WeNjgb2A/4DeAz4nWjVlZzBSpIkNWo1FFhr4xwzgJOB\nc4G9gBXAw4Q5We8yNDT0zvvBwUEGBwfbLLOccg9Wb70VXvPm5XhRSZJ63/DwMMPDw6nP0ypYbQIW\nNny/kN1DfuNGCMN/b9VfDwC/QYtg1QtyXxX48stw0EFQqeR4UUmSel9zw+drX/taovO0Ggp8jDDU\ntwjYA7iUMHm90b8AHyVMdN8L+DCwKlE1XSb3jpXDgJIklVqrjtUocDlwLyE43UxYEbis/vmNhK0Y\n7gGeAsaAmzBYZcNgJUlSqbUKVgB311+Nbmz6/tr6q69UqzBrVo4XNFhJklRq7ryeQiEdq4MOyvGC\nkiSpEwarFHbudChQkiTtZrBKoVrNeVWgwUqSpFIzWKXg5HVJktTIYJWCwUqSJDUyWKWQa7AaHYXX\nX4cDD8zpgpIkqVMGqxRyDVavvALz58PAQE4XlCRJnTJYpZDrI20cBpQkqfQMVink2rEyWEmSVHoG\nqxQMVpIkqZHBKgWDlSRJamSwSsFgJUmSGhmsUjBYSZKkRgarFFwVKEmSGhmsUrBjJUmSGhmsUsgt\nWI2NhQ1CDzooh4tJkqSkDFYp5Bastm6FOXNg5swcLiZJkpIyWKWQW7ByGFCSpK5gsErBYCVJkhoZ\nrFLIbVXgyy8brCRJ6gIGqxTsWEmSpEYGqxQMVpIkqZHBKgWDlSRJamSwSmjXrvB1YCCHixmsJEnq\nCgarhNx1XZIkNTNYJZTbisBazWAlSVKXMFgllFvHascOmDYN9t47h4tJkqQ0DFYJOXFdkiQ1M1gl\nZLCSJEnNDFYJGawkSVIzg1VC1WpOk9cNVpIkdQ2DVUI7d9qxkiRJ72awSsihQEmS1MxglVBuwWrz\nZjjkkBwuJEmS0jJYJZRrx8pgJUlSVzBYJWTHSpIkNTNYJZTLqsA334TRUZg7N+MLSZKkGAxWCeWy\nKnB8GLBSyfhCkiQpBoNVQrkMBW7e7IpASZK6iMEqodyClfOrJEnqGgarhAxWkiSpmcEqIYOVJElq\nZrBKaOfOHFYFGqwkSeoqBquE7FhJkqRmBquEDFaSJKmZwSqhzINVreZ2C5IkdZl2gtUSYDWwBrhq\niuNOBUaBT0eoq/QyD1ZvvAEDAzBnToYXkSRJMbUKVgPA9YRwdTywFDhukuO+DtwD9MU24ZkHK4cB\nJUnqOq2C1WnAWmA9UAXuAC6a4Lg/Af4R2BKzuDLLfFXg+ONsJElS12gVrA4FRhq+31j/WfMxFwE3\n1L+vxSmt3OxYSZKkZq2CVTsh6VvA1fVjKzgUGIcT1yVJ6jrTW3y+CVjY8P1CQteq0YcIQ4QA84Hz\nCcOGdzafbGho6J33g4ODDA4OdlRsmeQSrI44IsMLSJKkccPDwwwPD6c+T6vu0nTgOeBcYDPwCGEC\n+7OTHP8d4EfADyf4rFar9c4o4bJlcNJJ8Ed/lNEFLrkELr4YLr00owtIkqTJVCoVSDAK16pjNQpc\nDtxLWPl3MyFULat/fmOnF+wV1WrGk9edYyVJUtdpFawA7q6/Gk0WqD6frpzusXOnk9clSdK7ufN6\nQpnOsXLXdUmSupLBKqFMg9W2bbDnnrDXXhldQJIkZcFglVCmwcphQEmSupLBKiGDlSRJamawSijT\nVYEGK0mSupLBKqFMVwUarCRJ6koGq4QcCpQkSc0MVgllGqxeeslgJUlSFzJYJZR5x8o9rCRJ6joG\nq4QcCpQkSc0MVglltipwbAxeftmOlSRJXchglVBmqwK3boU5c2DmzAxOLkmSsmSwSiizoUCHASVJ\n6loGq4QMVpIkqZnBKiGDlSRJamawSshgJUmSmhmsEspsVaDBSpKkrmWwSqBWg9FRmD49g5MbrCRJ\n6loGqwSq1RCqKpUMTm6wkiSpaxmsEnDXdUmSNBGDVQKZBatdu2DLFliwIIOTS5KkrBmsEshs4vqW\nLTBvXobtMEmSlCWDVQJutSBJkiZisEogs+cEGqwkSepqBqsE7FhJkqSJGKwSMFhJkqSJGKwSMFhJ\nkqSJGKwS8HE2kiRpIgarBOxYSZKkiRisEnBVoCRJmojBKoFMOlbVKrz+Ohx4YOQTS5KkvBisEsgk\nWL3yChxwAAwMRD6xJEnKi8EqgUyClcOAkiR1PYNVApmsChwZgUMPjXxSSZKUJ4NVAplMXn/+eTjy\nyMgnlSRJeTJYJZDJUODzz8NRR0U+qSRJypPBKoHMgpUdK0mSuprBKgGDlSRJmojBKoHowWp0FDZs\ngCOOiHhSSZKUN4NVAtFXBW7cGDYGnTUr4kklSVLeDFYJRF8V6MR1SZJ6gsEqgehDgc6vkiSpJxis\nEjBYSZKkiRisEogerNatM1hJktQDDFYJRJ+8bsdKkqSeYLBKIJOhQCevS5LU9QxWCURdFfiLX4QT\nzp8f6YSSJKko7QarJcBqYA1w1QSfXwasBJ4CfgIsjlJdSUXtWL3wQhgGrFQinVCSJBVlehvHDADX\nAx8DNgGPAncCzzYc8zxwFrCdEML+D3B61EpLJGqwcuK6JEk9o52O1WnAWmA9UAXuAC5qOmYFIVQB\n/Aw4LFJ9pRQ1WDlxXZKkntFOsDoUGGn4fmP9Z5P5IrA8TVFlF3VVoBPXJUnqGe0Eq1oH5zsH+AIT\nz8PqGXasJEnSRNqZY7UJWNjw/UJC16rZYuAmwhyrbROdaGho6J33g4ODDA4OtllmuURdFWiwkiSp\ncMPDwwwPD6c+TztL0aYDzwHnApuBR4ClvHvy+uHA/cBngYcnOU+tVuuk+VVe550HX/oSfPzjKU80\nOgqzZ8Mbb8DMmVFqkyRJ6VXCav2Ol+y307EaBS4H7iWsELyZEKqW1T+/EfhLYB5wQ/1nVcKk954U\nbShwZAQWLDBUSZLUI9oJVgB311+Nbmx4//v1V1+IFqycuC5JUk9x5/UEoq0KdH6VJEk9xWCVQLSO\nlZuDSpLUUwxWCURbFWjHSpKknmKwSiDqHCuDlSRJPcNglYCT1yVJ0kQMVglECVbbtoV9rPbfP0pN\nkiSpeAarBKKsChwfBqx0vPeYJEkqKYNVAlE6Vs6vkiSp5xisEoiyKtBgJUlSzzFYJRCtY+XEdUmS\neorBKgGHAiVJ0kQMVh3atSvMNx8YSHkid12XJKnn5LkkrVar1XK8XDbefhv23Td8Taxahb33hh07\nIj10UJIkxVQJq/Y7zkl2rDoUZRhwZAQOOshQJUlSjzFYdSjaikAnrkuS1HMMVh1y4rokSZqMwapD\nUYKVE9clSepJBqsORX2cjSRJ6ikGqw45FChJkiZjsOpQ6mA1Ogpr1sDRR0erSZIklYPBqkOpVwWu\nXAmHHw7z5kWrSZIklYPBqkOpO1YPPABnnRWtHkmSVB4Gqw6lDlb/+Z8GK0mSepTBqkOpVgWOjcGD\nDxqsJEnqUQarDqXqWK1aBfvtB4ccErUmSZJUDgarDqWavO4woCRJPc1g1aFUHSsnrkuS1NMMVh1K\nHKxqtRCszj47ek2SJKkcDFYdShys1q4Nf/CII6LXJEmSysFg1aHEqwLH51dVKtFrkiRJ5WCw6lDi\njpXzqyRJ6nkGqw4lXhXo/CpJknqewapDiTpWL74Ib70FxxyTSU2SJKkcDFYdShSsxocBnV8lSVJP\nM1h1KNHkdYcBJUnqCwarDiXqWLnjuiRJfcFg1aGOg9VLL8Frr8GJJ2ZWkyRJKgeDVYc6XhX44IPw\n0Y/CNP+nliSp1/nbvkMdd6ycXyVJUt8wWHWo42Dl/CpJkvqGwapDHa0K3Lo17GF10kmZ1iRJksrB\nYNWhjjpWDz0EZ54J06dnWpMkSSoHg1WHOgpWt98O55yTaT2SJKk8bKV0qO1Vgf/wD/D443DzzZnX\nJEmSysFg1aG2OlYjI3D55XDXXTB7di51SZKk4jkU2KGWwWpsDH7v9+CKK+DUU3OrS5IkFc9g1aGW\nqwK/+c1w0NVX51aTJEkqh3aC1RJgNbAGuGqSY/62/vlKoKf3FpiyY/Xkk/DXfw233goDA7nWJUmS\nitcqWA0A1xPC1fHAUuC4pmMuAN4PHA38IXBD5BpLZdJg9dZb8JnPhI7VokV5l5W54eHhoktQQt67\n7ub9617eu/7UKlidBqwF1gNV4A7goqZjPgV8t/7+Z8C+wIJ4JZbLpKsCv/xlWLwYLrss95ry4H8g\nupf3rrt5/7qX964/tVoVeCgw0vD9RuDDbRxzGPBK6upKqFqFGdNrsP5F+OlPYcWK8HXbNvj5z6FS\nKbpESZJUkFbBqtbmeZrTxIR/7sIFj7R5uvJas/UEZl5wIUx/NuyqfuaZsHQpnHwyzJpVdHmSJKlA\nrdorpwNDhDlWAH8OjAFfbzjmfwPDhGFCCBPdz+a9Hau1wFHJS5UkScrNOsIc8qim10+8CNgDeJKJ\nJ68vr78/HXg4dhGSJEm94nzgOULH6c/rP1tWf427vv75SuDkXKuTJEmSJEmSWnFD0e7V6t5dRrhn\nTwE/ARbnV5ra0M7fPYBTgVHg03kUpba0c+8GgSeAZwjzWlUere7ffOAewnSaZ4DP5VaZWrmFMCf8\n6SmOKTSzDBCGBBcBM2g9J+vDOCerLNq5d2cA+9TfL8F7Vybt3L/x4+4HfgxcnFdxmlI7925f4P8R\ntrKB8Ita5dDO/RsC/mf9/XxgK61X5Ssfv0kIS5MFq44zS+xnBbqhaPdq596tALbX3/+M3f+RV/Ha\nuX8AfwL8I7Alt8rUSjv37jPADwj7BAK8lldxaqmd+/cSMLf+fi4hWI3mVJ+m9iCwbYrPO84ssYPV\nRJuFHtrGMf6CLl47967RF9md4lW8dv/uXcTux061u0+dstXOvTsa2A/4D+Ax4HfyKU1taOf+3QSc\nAGwmDCddkU9piqDjzBK7FRl1Q1HlqpN7cA7wBeAjGdWizrVz/74FXF0/tkLrfeyUj3bu3QzCiutz\ngb0I3eOHCfM+VKx27t9fEIYIBwn7Of4b8BvAjuzKUkQdZZbYwWoTsLDh+4Xsbl1Pdsxh9Z+pWO3c\nOwgT1m8izLGaqn2qfLVz/z7E7o185xO2UqkCd2ZenabSzr0bIQz/vVV/PUD4xWywKl479+9M4Jr6\n+3XAC8CxhO6jyq3wzOKGot2rnXt3OGEuwem5VqZ2tHP/Gn0HVwWWRTv37gPAvxMmSu9FmGh7fH4l\nagrt3L9vAl+tv19ACF775VSfWltEe5PXC8ssbijavVrdu/9LmHT5RP3V/Q9/7C3t/N0bZ7Aql3bu\n3X8jrAx8GvjTXKtTK63u33zgR4TfeU8TFiOoHG4nzH3bSegMfwEziyRJkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRJkiSV2f7s3gvtJcJmiU8QHu9xfUbXvBz43BSffwr47xldW5IkKRdfBa7M+BoVQnCb6hFe\nFcKO2TMyrkVSj5tWdAGS+t74A04HCbtTAwwB3yU8E289YZf4a4GngLvZHZI+BAwTnrl2D3DQBOf/\nCLAaGK1//6eEHcxXEnZdhvBQ1RXAx9P+w0jqbwYrSWX1PuAcwjDd94B/IzwE/C3gk4Tu0nXAxcAp\nhMf0XDPBeT7Kux92exXwQcJDjBsfW/EIcFbUfwJJfWeq1rgkFaVG6EztAp4h/J/Ae+ufPU14aOox\nwAmEhxNDeEDx5gnOdTjwUMP3TwG3Af9cf43bDCyJUr2kvmWwklRWO+tfx4Bqw8/HCP/tqhCG9M5s\n41yVhvefJHSmLgS+ApxYP+c0QqCTpMQcCpRURpXWh/AccABwev37GcDxExz3IrvnXlUIHaxh4Gpg\nH2Dv+mcH14+VpMQMVpKKVmv4OtF7eG8nqUboYv1X4OuEFX1PAGdMcP6HCHOwIHS6biUMBz4O/C/g\njfpnpxEmy0uSJGkS49st7DHFMdMI4czpEZJSGSi6AEnKwQBwHCE8TeRCYCthiFCSJEmSJEmSJEmS\nJEmSJEmSJEmSJEmSJEmS1J7/D3fX0KOwLllLAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 68 + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Impulse response" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "x = pylab.hstack([pylab.zeros(tau), 1, pylab.zeros(9*tau - 1)])\n", + "t = pylab.linspace(0, tau / 10, x.size)\n", + "fi2 = iir.Filter(fi)\n", + "\n", + "y = np.array(x.size, dtype=np.double)\n", + "y.resize(x.size)\n", + "for i in range(0, y.size):\n", + " y[i] = fi2.f(x[i])\n", + "\n", + "fig, axes = plt.subplots(figsize=(10,6), dpi=100)\n", + "axes.plot(t, y, 'r')\n", + "axes.plot(t, x)\n", + "axes.set_xlabel('Time (s)')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 69, + "text": [ + "" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAF/CAYAAABzDnyLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH39JREFUeJzt3X2UXXV97/H3yUyChQIBYqFCaJBnuJVgFIIKHQRpQCDW\npjyIIASVrprC6tUlUpcS/8CWLkx5SrkQeYhcAVvtcgGSRLgwSOUhqCSAEiCh6UoIhAchIKCZmZz7\nx++czJ7DmZlzZu/fPufMfr/W2mv2mbPn7F/YJPnk+/t99wZJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ\nkiQVzA3ARuCJEY65EngWWAkcmsegJEmSxpsjCUFquNB1AnBXZf9w4OE8BiVJkjQeTWP40PV/gFMT\nr1cBu8YekCRJUruZEPnzdwfWJV6vB/aIfE5JkqS2Ezt0AZRqXpdzOKckSVJb6Y78+c8DUxOv96h8\nb4i99967vGbNmshDkSRJysQaYJ9mfyh2pet24KzK/kzgdUK34xBr1qyhXC6P2+2448rcemvrxxFr\nu/jii1s+BjevXxE3r11nb16/zt2AvccSitJWum4F/gKYQli7dTEwsfLetYTOxROA1cBbwDkpz9eR\n+vrCJkmSiitt6Dq9gWPmpTxHxzN0SZKkPBbSF954D109PT2tHoJS8Pp1Lq9dZ/P6FU9tZ2GrlCtz\npOPSjBlwzjkwr/A1P0mSOl+pVIIxZCgrXTnYvDlskiSpuAxdORjv04uSJGl0hq4cGLokSZKhKweG\nLkmSZOjKgaFLkiQZunJg6JIkSYauHNi9KEmSDF05sNIlSZIMXTkwdEmSJENXZOUyDAwYuiRJKjpD\nV2TVsGXokiSp2AxdkRm6JEkSGLqiq3Yt2r0oSVKxGbois9IlSZLA0BWdoUuSJIGhKzpDlyRJAkNX\ndIYuSZIEhq7oqmHLhfSSJBWboSuyatiy0iVJUrEZuiLr64NttjF0SZJUdIauyPr6YNttDV2SJBWd\noSsyQ5ckSQJDV3SGLkmSBIau6Pr6YLvt7F6UJKnoDF2Rbd4cQpeVLkmSis3QFZnTi5IkCQxd0Rm6\nJEkSGLqiM3RJkiQwdEXX1wfveQ9s2QIDA60ejSRJahVDV2R9fTBpUtisdkmSVFyGrsg2b4aJE8Nm\n6JIkqbgMXZH19Rm6JEmSoSs6Q5ckSQJDV3SGLkmSBIau6FxIL0mSwNAVXbLS5fMXJUkqLkNXZHYv\nSpIkMHRF55ouSZIEhq7oDF2SJAkMXdEZuiRJEhi6okt2L7qQXpKk4jJ0ReZCekmSBIau6JxelCRJ\nYOiKztAlSZLA0BWdoUuSJEE2oWsWsAp4FriwzvtTgKXACuBJ4OwMztkxDF2SJAnSh64u4GpC8DoI\nOB04sOaYecBjwHSgB/gO0J3yvB3D7kVJkgTpQ9dhwGpgLdAH3AbMrjnmBWCHyv4OwKtAf8rzdgy7\nFyVJEqSvOO0OrEu8Xg8cXnPMIuBeYAOwPXBKynN2FKcXJUkSpK90lRs45h8J67neR5hiXEgIX4Vg\n6JIkSZC+0vU8MDXxeiqh2pX0EeCSyv4a4L+B/YFfJA+aP3/+1v2enh56enpSDq09GLokSepsvb29\n9Pb2pv6cUsqf7waeBo4hTB8uJyymfypxzAJgE/AtYFfgl8AHgN8mjimXy40UzTrPfvvBHXfAjTfC\njjvCRRe1ekSSJCmNUqkEY8hQaStd/YTuxGWETsbrCYHrvMr71wLfBm4EVhKmM7/K0MA1rtm9KEmS\nIJtbNyypbEnXJvZfAU7K4DwdKdm9+Ic/tHo0kiSpVbwjfWSu6ZIkSWDois7QJUmSwNAVnaFLkiSB\noSs6Q5ckSQJDV1Tl8mDosntRkqRiM3RFNDAApRJ0dVnpkiSp6AxdEVWrXGDokiSp6AxdERm6JElS\nlaErIkOXJEmqMnRFVH0EELiQXpKkojN0RWSlS5IkVRm6Iqo+dxEMXZIkFZ2hKyIrXZIkqcrQFZGh\nS5IkVRm6IjJ0SZKkKkNXRHYvSpKkKkNXRFa6JElSlaErIrsXJUlSlaErIitdkiSpytAVkaFLkiRV\nGboiMnRJkqQqQ1dEdi9KkqQqQ1dEVrokSVKVoSui2u7F/n4ol1s7JkmS1BqGroiSla5SCbq6QvCS\nJEnFY+iKKBm6wClGSZKKzNAVkaFLkiRVGboiSnYvgh2MkiQVmaErouRCerDSJUlSkRm6InJ6UZIk\nVRm6IjJ0SZKkKkNXRIYuSZJUZeiKqF7ociG9JEnFZOiKqF73opUuSZKKydAVkd2LkiSpytAVkWu6\nJElSlaErIkOXJEmqMnRFZOiSJElVhq6I7F6UJElVhq6I7F6UJElVhq6I7F6UJElVhq6IXNMlSZKq\nDF0RGbokSVKVoSsiQ5ckSaoydEVUbyG93YuSJBWToSsiK12SJKnK0BWR3YuSJKkqi9A1C1gFPAtc\nOMwxPcBjwJNAbwbn7AhWuiRJUlV3yp/vAq4GjgWeBx4FbgeeShwzGVgI/CWwHpiS8pwdw9AlSZKq\n0la6DgNWA2uBPuA2YHbNMZ8BfkQIXACvpDxnxzB0SZKkqrSha3dgXeL1+sr3kvYFdgbuA34BnJny\nnB3D7kVJklSVdnqx3MAxE4EPAscA2wIPAQ8T1oCNa1a6JElSVdrQ9TwwNfF6KoPTiFXrCFOK71S2\nnwGHUBO65s+fv3W/p6eHnp6elENrPbsXJUnqfL29vfT29qb+nFLKn+8GniZUsTYAy4HTGbqQ/gDC\nYvu/BLYBHgFOBX6TOKZcLjdSNOss3d3wzjuDweu66+DRR2HRotaOS5IkjV2pVIIxZKi0la5+YB6w\njNDJeD0hcJ1Xef9awu0klgKPA1uARQwNXONSuQwDAyF4VVnpkiSpuNKGLoAllS3p2prXl1W2wujr\nC4GrlMjBEye6kF6SpKLyjvSR1HYuQnhtpUuSpGIydEVS27kITi9KklRkhq5IajsXwdAlSVKRGboi\nsdIlSZKSDF2RGLokSVKSoSuS4UKX3YuSJBWToSsSuxclSVKSoSsSpxclSVKSoSsSuxclSVKSoSsS\nK12SJCnJ0BWJoUuSJCUZuiKxe1GSJCUZuiKxe1GSJCUZuiJxIb0kSUoydEXimi5JkpRk6IrE0CVJ\nkpIMXZEYuiRJUpKhK5J6C+m7uqBchoGB1oxJkiS1jqErknqVrlLJapckSUVl6IqkXvciGLokSSoq\nQ1ck9SpdYOiSJKmoDF2RGLokSVKSoSuSkUKXjwKSJKl4DF2R1OteBB8FJElSURm6InF6UZIkJRm6\nIrF7UZIkJRm6IrHSJUmSkgxdkRi6JElSkqErErsXJUlSkqErErsXJUlSkqErEqcXJUlSkqErErsX\nJUlSkqErEitdkiQpydAViaFLkiQlGboisXtRkiQlGboisXtRkiQlGboicXpRkiQlGboisXtRkiQl\nGboisdIlSZKSDF2RGLokSVKSoSsSuxclSVKSoSsSuxclSVKSoSsSpxclSVKSoSsSuxclSVKSoSsS\nK12SJCnJ0BWJC+klSVKSoSsSK12SJCkpi9A1C1gFPAtcOMJxHwb6gU9ncM62N1zosntRkqRiShu6\nuoCrCcHrIOB04MBhjrsUWAqUUp6z7Q0MQKkEXV3vfs9KlyRJxZQ2dB0GrAbWAn3AbcDsOsf9PfBD\n4OWU5+sIw3UugqFLkqSiShu6dgfWJV6vr3yv9pjZwDWV1+WU52x7w00tgqFLkqSi6k75840EqMuB\nr1WOLTHM9OL8+fO37vf09NDT05NyaK0zWuiye1GSpM7R29tLb29v6s9Ju75qJjCfsKYL4CJgC2H9\nVtVzifNMAd4GvgDcnjimXC6PnwLYiy/C9Onha60lS+CKK2Dp0vzHJUmS0iuVSjCGDJW20vULYF9g\nGrABOJWwmD7p/Yn9G4E7GBq4xh2nFyVJUq20oasfmAcsI3QoXg88BZxXef/alJ/fkQxdkiSpVtrQ\nBbCksiUNF7bOyeB8bc/uRUmSVMs70kdgpUuSJNUydEVg96IkSapl6Iqgry887qceHwMkSVIxGboi\ncHpRkiTVMnRF4EJ6SZJUy9AVgZUuSZJUy9AVgaFLkiTVMnRFYPeiJEmqZeiKwO5FSZJUy9AVgdOL\nkiSplqErgtG6F/v7oVzOd0ySJKm1DF0RjFTpKpWgqysEL0mSVByGrghGCl3gYnpJkorI0BVBI6HL\ndV2SJBWLoSuCkboXwQ5GSZKKyNAVgZUuSZJUy9AVwUjdi2DokiSpiAxdEVjpkiRJtQxdEdi9KEmS\nahm6IrDSJUmSahm6IrB7UZIk1TJ0RWClS5Ik1TJ0RWD3oiRJqmXoisBKlyRJqmXoisDuRUmSVMvQ\nFYEL6SVJUi1DVwROL0qSpFqGrggMXZIkqZahKwK7FyVJUi1DVwQupJckSbUMXRE4vShJkmoZuiKw\ne1GSJNUydEVgpUuSJNUydEVg6JIkSbUMXRHYvShJkmoZuiKwe1GSJNUydEXg9KIkSapl6IrA7kVJ\nklTL0BWBlS5JklTL0BWBoUuSJNUydEVg96IkSapl6MpYuQwDA9DdPfwxdi9KklQ8hq6M9fWFwFUq\nDX+MlS5JkorH0JWx0ToXwe5FSZKKyNCVsdEW0YOVLkmSisjQlbHRFtGDoUuSpCLKInTNAlYBzwIX\n1nn/DGAl8Djwc+ADGZyzbVnpkiRJ9YzQY9eQLuBq4FjgeeBR4HbgqcQxzwFHAZsIAe06YGbK87at\nRkOX3YuSJBVL2krXYcBqYC3QB9wGzK455iFC4AJ4BNgj5TnbWltUup57Dl56KeIJJElSs9KGrt2B\ndYnX6yvfG865wF0pz9nWWta9+OabcMMNcOSRMH06zJ4dbhgmSZLaQtrQVW7i2KOBudRf9zVu5F7p\nuv9+OOssmDoVbr8dvvxlePll2GYbuPzyjE4iSZLSSrum63lgauL1VEK1q9YHgEWENV2v1fug+fPn\nb93v6emhp6cn5dBaI9fuxfvvh9NOg69+FS67DP7kTwbfu/56OPxwOOkk2G+/DE4mSVIx9fb20tvb\nm/pzRrhvekO6gaeBY4ANwHLgdIYupN8TuBf4LPDwMJ9TLpebKZq1rwcfDMWmhx4a/phVq+Dkk+GZ\nZ1KebM4cOPpo+NKX6r9/5ZXw7/8ewllXV8qTSZIkgFJ47EzTGSrt9GI/MA9YBvwG+AEhcJ1X2QC+\nCewEXAM8Rghm41Zu04vr18O994apxeHMmwcTJsBVV6U8mSRJSivt9CLAksqWdG1i//OVrRByC13X\nXQef+Qxsv/3wx0yYEKYZjzgCTjwR9tkn5UklSdJYeUf6jOXSvbh5MyxaBH/3d6Mfu+++8PWvw9y5\nsGVLipNKkqQ0DF0Zy6XS9aMfwYEHwkEHNXb8+eeH20csXJjipJIkKQ1DV8Zy6V5cuHD4xfP1dHWF\ne3h961vwwgspTixJksbK0JWx6I8BWrkS1q4NNz9txv77w6c/DTfdNMYTS5KkNAxdGYs+vbhwIZx3\nHnSPoQfi3HPDwnrXdkmSlDtDV8YaWUhfvWVW00/pef11+I//gC98YUxj47DD4I/+CH72s7H9vCRJ\nGjNDV8YaqXSVSmOsdt10E8yaBbvtNrbBlUrw+c/Dd787tp+XJEljZujKWCOhC8YQurZsgX/7t+YW\n0Nfz2c/CnXfCa3WfxiRJkiIxdGWske5FGEPouueeMDX40Y+OeWwA7LJLqJbdcku6z5EkSU0xdGWs\nmUpXUx2M11wTqlyltI/LJEwxLloE4+R5l5IkdQJDV8aiTC++9VaodJ1ySqqxbfXxj8OmTfCrX2Xz\neZIkaVSGrow10r0ITT4K6N574cMfhsmTU41tqwkTwmOBrr8+m8+TJEmjMnRlLEql68474ZOfTDWu\ndzn7bLjtNnj77Ww/V5Ik1WXoyljmoatcDqHrxBNTj22IqVNh5szwHEdJkhSdoStjmXcvrlgB224L\n++2Xemzv4j27JEnKjaErY5l3L/7kJ6HKlUXXYq0TT4Snn4Znnsn+syVJ0hCGroxlPr0YY2qxatIk\nOOssuOGGOJ8vSZK2MnRlLNPuxY0bYdUqOPLITMZW19y58L3vQX9/vHNIkiRDV9YyrXQtWQLHHttY\nihurAw4Ii+rvuSfeOSRJkqEra5mGrphTi0lnnRWqXZIkKRpDV8aa6V4ccSH95s2h+nT88ZmNbVin\nnQZ33RXuUi9JkqIwdGUss0rXAw+Eqb9dd81sbMPaZZfwaKAf/jD+uSRJKihDV8YyC10x7kI/ks99\nDhYvzu98kiQVjKErY5l1L+a1nqvq+ONDp+Rzz+V3TkmSCsTQlbFMKl3PPBOeiTh9eqZjG9GkSWFt\n180353dOSZIKxNCVsUxCV7XKFeMu9COpdjGWy/meV5KkAjB0ZSyT7sW8pxarZsyAbbaBn/88/3NL\nkjTOGboylrrStWkTPPpo6CbMW6kUFtR7zy5JkjJn6MpY6tC1bFl47M9222U+toaccUa4dcQ777Tm\n/JIkjVOGroyl7l78/vfhb/4m83E1bI894EMfgttvb90YJEkahwxdGUtV6Xr5Zbj/fpgzJ8rYGuZj\ngSRJypyhK2PNLKR/V+i65RY46STYfvsoY2vYX/0VPPggvPhia8chSdI4YujKWDOVrnd1L950E5x9\ndoRRNWm77ULwuvHGVo9EkqRxw9CVsTFPLz7+OLz6Khx9dLSxNeWCC+Dqq0d5KrckSWqUoStjY15I\nv3gxnHkmTGiTS3LIIXDwwXDrra0eiSRJ40Kb/A0/foyp0tXfH7oWzzor6tia9pWvwGWXeYd6SZIy\nYOjK0MBAuL9oI8WqIaFr2TLYay/Yf/+o42vaJz4RfkE//WmrRyJJUsczdGWo0c5FqAldixeHO8G3\nm1JpsNolSZJSMXRlqNGpRUh0L772WqgknXpq1LGN2WmnwVNPwYoVrR6JJEkdzdCVoWZDV18f8IMf\nwHHHwU47RR3bmE2aFDoZv/OdVo9EkqSOZujKUKOdi5DoXly8uD3uzTWSL34R7roL1q1r9UgkSepY\nhq4MNV3p2vQWrF0bKl3tbMcdQzC84opWj0SSpI5l6MpQ06Hr+ZfgjDOguzvuwLJwwQXhDvWbNrV6\nJJIkdSRDV4aa6l6kj80bX2vPrsV69twTZs2CRYtaPRJJkjqSoStDDVe6ymUmLriUvm13hD//8+jj\nysyXvxymGN95p9UjkSSp4xi6MtRw6FqwgImrnqBv92mxh5StD34QjjwSLryw1SORJKnjGLoy1FD3\n4h13wIIFTFr4r/QNdOUyrkwtXAg//rF3qZckqUlZhK5ZwCrgWWC4EsiVlfdXAodmcM62NGql6/HH\n4dxz4T//k4l/9r6hD7zuFDvtFBbUz50Lr77a6tFIktQx0oauLuBqQvA6CDgdOLDmmBOAfYB9gS8C\n16Q8Z9saMXRt3AgnnwxXXQWHHz70MUCd5phj4JRT4G//1odhS5LUoLSh6zBgNbAW6ANuA2bXHHMy\nsLiy/wgwGdg15Xnb0rDdi7//PXzqU3DOOVsf97P1MUCd6tvfDo8HuvnmVo9EkqSOkPYGUbsDyduU\nrwcOb+CYPYCNKc/ddoZUut58E5YvhwcfhDvvhPe/H775za3HdnSlC+A974Hvfx+OPRaOOgqmTWv1\niCRJamtpQ1ejc0ul0X5ufk9vyqG03qpXpjDxrX6YfjasXg2HHgpHHAEXXQQnnAClwf8MEyfC734H\n8+e3bLgZOASm/xj+4r5wv7EJLe7LKJfh7bfhjTfCTVzfeGNw+8MfQmkxufX3D06PJqdJS6Xwa5kw\nAbq6hu7X2+odN2HC0M+p7pdKQ7fq+ZJfm/n1NvK9kb4/lnOMpNlfgyQVSNrQ9TwwNfF6KqGSNdIx\ne1S+N0Tv2pu27k+bPJ1pk6enHFr+DpjyCh876rfwuWtD4BqhlXHyZLjkklAQ62hHHAHPPAP33Qcf\n/3h+f+mWy/Db38ILL4RtwwZ48cXw3g47hEcXVb/uvXeozE2aFLaJE8PX7u6h4QfC/pYtQ7eBgcGv\n9bbaY6o/Vy6/e7869urr6n7y1zWaev+Nh/vv3uz3x3qcJI1ja19fwdrXV6T+nLR/onYDTwPHABuA\n5YTF9E8ljjkBmFf5OhO4vPI1qVx2QXbnWr8e5syB7beHG26AqVNH/5lmlcvwm9/APffA3XfDAw/A\nzjuHe4fNmBG+Hnoo7DoulwtKktpIKfyDtOkMlcU/Y48nBKku4Hrgn4DzKu9dW/la7XB8CzgH+FXN\nZxi6Ol1/P1x6abhj/YIF4ZmSaaskL70U7gd2990hbE2aBJ/4RFhHdvTR8N73ZjN2SZKa0MrQlQVD\n13jx2GNw5plwwAFwzTXNBaP+fnj4YVi6NGyrV4dwddxxIWztvbfTXZKkljN0qX38/vfwjW+E7sY5\nc2CffQa3adNCxeqNN8J04a9/PbgtXw577RUerD1rVlgv1ugTxCVJyomhS+3nl78Ma6+efTZUrVav\nDuu/dtwR3noLDjwQDj54cJsxA3bbrdWjliRpRIYudYbNm+GVV0K4avUtJiRJGgNDlyRJUg7GGros\nNUiSJOXA0CVJkpQDQ5ckSVIODF2SJEk5MHRJkiTlwNAlSZKUA0OXJElSDgxdkiRJOTB0SZIk5cDQ\nJUmSlANDlyRJUg4MXZIkSTkwdEmSJOXA0CVJkpQDQ5ckSVIODF2SJEk5MHRJkiTlwNAlSZKUA0OX\nJElSDgxdkiRJOTB0SZIk5cDQJUmSlANDlyRJUg4MXZIkSTkwdEmSJOXA0CVJkpQDQ5ckSVIODF2S\nJEk5MHRJkiTlwNAlSZKUA0OXJElSDgxdkiRJOTB0SZIk5cDQJUmSlANDlyRJUg4MXZIkSTkwdEmS\nJOXA0CVJkpQDQ5ckSVIODF2SJEk5MHRJkiTlwNAlSZKUA0OXJElSDtKErp2Bu4FngJ8Ck+scMxW4\nD/g18CRwforzSZIkdaw0oetrhNC1H/D/Kq9r9QH/ABwMzAS+BByY4pxqQ729va0eglLw+nUur11n\n8/oVT5rQdTKwuLK/GPhUnWNeBFZU9n8HPAW8L8U51Yb8g6Ozef06l9eus3n9iidN6NoV2FjZ31h5\nPZJpwKHAIynOKUmS1JG6R3n/bmC3Ot//es3rcmUbzh8DPwQuIFS8JEmSCqWU4mdXAT2EKcQ/JSyY\nP6DOcROBO4ElwOXDfNZqYO8UY5EkScrLGmCfPE/4L8CFlf2vAf9c55gS8D3gX/MalCRJ0nizM3AP\n775lxPuAn1T2PwZsISymf6yyzcp3mJIkSZIkSVIEswhrwZ5lcGqy1pWV91cSuh3VPka7fmcQrtvj\nwM+BD+Q3NI2ikd97AB8G+oFP5zEoNayR69dDmE14EujNZVRq1GjXbwqwlDAr9CRwdm4j00huINyd\n4YkRjmnbzNJFWDA/jbC4fgXvvlHqCcBdlf3DgYfzGpxG1cj1OwLYsbI/C69fu2jk2lWPu5fQ+PLX\neQ1Oo2rk+k0mPPljj8rrKXkNTqNq5PrNB/6psj8FeJXR7y6g+I4kBKnhQlfTmSXPZy8eRvgfby3h\nTvW3AbNrjknecPURwh8ko93/S/lo5Po9BGyq7D/C4F8Aaq1Grh3A3xNu7fJybiNTIxq5fp8BfgSs\nr7x+Ja/BaVSNXL8XgB0q+zsQQld/TuPT8B4AXhvh/aYzS56ha3dgXeL1+sr3RjvGv7jbQyPXL+lc\nBv8FoNZq9PfebOCayuuR7runfDVy/fYlNDfdB/wCODOfoakBjVy/RYTH5W0gTFNdkM/QlFLTmSXP\n8mWjf4jX3jvMP/zbQzPX4WhgLvDRSGNRcxq5dpcTbv1SJvweTHMPP2Wrkes3EfggcAywLaHq/DBh\nrYlaq5Hr94+Eaccewj0r7wYOAd6MNyxlpKnMkmfoeh6Ymng9lcFS+HDH7FH5nlqvkesHYfH8IsKa\nrpHKsspPI9duBmHaA8KakuMJUyG3Rx+dRtPI9VtHmFJ8p7L9jPCXtqGr9Rq5fh8BLqnsrwH+G9if\nULVU+2rrzNJN+J9pGjCJ0RfSz8SF2O2kkeu3J2HtwsxcR6bRNHLtkm7E7sV20sj1O4Bw38QuQqXr\nCeCg/IaoETRy/RYAF1f2dyWEsp1zGp9GNo3GFtK3ZWY5Hnia8BfzRZXvnVfZqq6uvL+SUC5X+xjt\n+n2XsAC0eiPc5XkPUMNq5PdelaGr/TRy/b5C6GB8Ajg/19FpNKNdvynAHYS/954gNEao9W4lrLPb\nTKgmz8XMIkmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJElqrV0YvNfbC4QbRT5GeCTK1ZHOOQ84e4T3\nTwa+EenckiRJLXcx8L8jn6NECHUjPRatRLiT+MTIY5E0jk1o9QAkaRTVB8r2EO7aDTAfWEx4xuBa\nwh30LwMeB5YwGKBmAL2EZ9gtBXar8/kfBVYB/ZXX5xPu7L6ScEdqCA+xfQg4Lu0vRlJxGbokdaq9\ngKMJU3//F7ib8MD1d4BPEqpSVwF/DXyI8HijS+p8zscY+mDhC4HphAdGJx/3sRw4KtNfgaRCGamc\nLkntqkyoaA0ATxL+Abms8t4ThIfU7gccTHgQNISHQW+o81l7Av+VeP04cAvw48pWtQGYlcnoJRWS\noUtSp9pc+boF6Et8fwvhz7YSYZrwIw18Vimx/0lCResk4OvA/6p85gRC2JOkMXF6UVInKo1+CE8D\n7wVmVl5PBA6qc9z/MLjWq0SofPUCXwN2BP648t6fVo6VpDExdElqd+XE13r78O4KVJlQ/ZoDXEro\nPHwMOKLO5/8XYc0XhArZzYQpxl8BVwBvVN47jLBwX5IkSWNQvWXEpBGOmUAIbi7JkDRmXa0egCS1\ngS7gQEKwquck4FXCtKMkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkqRj+P8oIv3GKzttMAAAAAElF\nTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 69 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Linear Bode plot" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "freq = pylab.linspace(0, nyquist_freq, num_points)\n", + "resp = respfunc(freq)\n", + "\n", + "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 10), dpi=100)\n", + "\n", + "axes[0].plot(freq, abs(resp), 'r')\n", + "axes[0].set_xlabel('Frequency (Hz)')\n", + "axes[0].set_ylabel('Magnitude')\n", + "\n", + "axes[1].plot(freq, pylab.angle(resp), 'r')\n", + "axes[1].set_xlabel('Frequency (Hz)')\n", + "axes[1].set_ylabel('Phase (radians)')\n", + "axes[1].set_ylim([-pylab.pi, pylab.pi])\n", + "axes[1].set_yticks(pylab.linspace(-pylab.pi, pylab.pi, 9))\n", + "axes[1].set_yticklabels([\"-\u03c0\", \"-\u00be\u03c0\", \"-\u00bd\u03c0\", \"-\u00bc\u03c0\", \"0\", \"\u00bc\u03c0\", \"\u00bd\u03c0\", \"\u00be\u03c0\", \"\u03c0\"])\n", + "\n", + "None" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJeCAYAAADxx6J0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XecVNX9//HX0kERVIoNBalWxAKoGNeCgkbR2GLvLZqi\niVGjxtU0TfJVEzX2KBoNit2fBevaFRDBSlsBBQXBgigd5vfHmc0Oyy47u+zMvXfm9Xw87mNm7k75\n4Ij79nPuOQckSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkNcC/gTnA+2t4zj+BKcAEoF8+ipIkSSpW\nexACV23h7ADgqfT9AcBb+ShKkiSpmHWl9nB2M3BUxuOJQOdcFyRJkhQXTaIuoJpNgc8yHs8ENouo\nFkmSpLyLWzgDKKn2OBVJFZIkSRFoFnUB1cwCumQ83ix9bhXdO3dOVcyZk7eiJEmS1kIF0CPbJ8ct\nnD0OnAuMAAYC3xJmd66iYs4cUqk8NNSWLIEZM2DqVKiogClTYMIEePdd2GADGDgQ9tknHFtumft6\nCkRZWRllZWVRl6EG8LtLNr+/5PK7S7aSkpLu9Xl+vsPZf4E9gQ6Ea8suB5qnf3YLYabmAcBU4Afg\n5DzXt6qWLaFXr3BkWrkyBLbXXoMXXoDLLoN114XDDoMjj4Qdd4SS6qOzkiRJdct3ODs6i+ecm/Mq\n1laTJlWh7ZRTIJWC8eNh5Ej46U/Dz888E046KXTYJEmSshTHCQHJU1IC/frBn/8MkyfDv/8dhj67\nd4ezzgpDo/qf0tLSqEtQA/ndJZvfX3L53RWXpI69pfJyzdnamjsXrr0WbrklDHlefjlsumnUVUmS\npDwqCZc6ZZ257JzlUseOoZs2ZQp06AB9+8JVV4WJBpIkSTWwc5ZPn3wC550HH30Ed94JgwZFXZEk\nScqx+nbODGdReOyxcC3aSSfBFVdAixZRVyRJknLEYc0kGDYsrJf28cdhrbTp06OuSJIkxYThLCqd\nOsEjj8AJJ4SA9uKLUVckSZJiwGHNOHjhBTj22LCY7TnnRF2NJElqRF5zllTTpsH++8Mxx4QlN9xh\nQJKkgmA4S7I5c2DIkDCL8x//CDsNSJKkRDOcJd38+XDQQWFrqFtvNaBJkpRwztZMunbt4Kmnwlpo\n558f9u2UJElFw3AWR+uuGwLayy/D738fdTWSJCmPDGdx1b49PPssPPgg3HBD1NVIkqQ8aRZ1AVqD\njh3hySdh992hR48wWUCSJBU0JwQkwWuvwU9+AuXlsPXWUVcjSZLqwQkBhWjQIPi//4Mf/xi+/jrq\naiRJUg7ZOUuS886DioqwcbqL1EqSlAh2zgrZ1VeHhWr/7/+irkSSJOVIUtsvxdk5A5gxA/r3h4cf\nDhMFJElSrNk5K3RbbAF33BH24Jw/P+pqJElSI7NzllRnnQWLF8Ndd0VdiSRJWgP31iwW338PffvC\nNdfAsGFRVyNJkmphOCsmr70GRxwB770XFqyVJEmxYzgrNr/9bZgkcP/9UVciSZJqYDgrNosWwXbb\nwfXXw9ChUVcjSZKqMZwVo2efDRMEPvgA2rSJuhpJkpTBcFasjjkGNt8crroq6kokSVIGw1mxmj0b\ntt8enn8+3EqSpFgwnBWzf/0LRo6EF190701JkmLCHQKK2RlnwLx58MgjUVciSZIaKKntFTtntXnh\nBTj9dPjoI2jVKupqJEkqenbOit0++4Rrzq67LupKJElSA9g5K0RTp8KAAfDhh7DRRlFXI0lSUXNC\ngILzz4elS+GGG6KuRJKkomY4UzB3Lmy1FYweDVtuGXU1kiQVLa85U9CxI/ziF/D730ddiSRJqgc7\nZ4VswQLo2RNGjYK+faOuRpKkomTnTFXatoXf/Q4uuSTqSiRJUpbsnBW6JUugd2+45x7YY4+oq5Ek\nqeg4IUCru/POEM5efDHqSiRJKjoOa2p1xx0H06fDq69GXYkkSapDvsPZEGAiMAW4sIafdwCeAcYD\nHwAn5a2yQta8ebj27A9/iLoSSZJUh3wOazYFJgH7ArOAMcDRwMcZzykDWgIXE4LaJKAzsLzaezms\nWV9Ll4aZmyNGwK67Rl2NJElFI87Dmv2BqcB0YBkwAhhW7TlfAOul768HfMXqwUwN0aIFXHyx3TNJ\nkmIun+FsU+CzjMcz0+cy3QZsA3wOTAB+mZ/SisTJJ8P778OYMVFXIkmSapHPcJbNOOTvCNebbQLs\nANwItM1lUUWlZUu48EK7Z5IkxVizPH7WLKBLxuMuhO5Zpt2AP6XvVwDTgN7A2OpvVlZW9r/7paWl\nlJaWNl6lhezUU+GPf4QPP4Rttom6GkmSCk55eTnl5eUNfn0+JwQ0I1zgvw9h2HI0q08IuAaYD1xB\nmAjwDrA98HW193JCwNr4059g6tSw/pkkScqpuC9COxS4jjBz8w7gL8CZ6Z/dQpiheSewOWHI9S/A\nfTW8j+FsbXz9NfToAe+9B5ttFnU1kiQVtLiHs8ZiOFtbv/pVWP/sb3+LuhJJkgqa4UzZmTEDdtwR\nPvkE2rWLuhpJkgpWnNc5U5xssQUMGQK33BJ1JZIkKYOds2I2fjwceGDonrVsGXU1kiQVJDtnyt4O\nO8C228J9Nc25kCRJUbBzVuyefz5MDnj/fShJ6r8OkiTFl50z1c8++4TbF1+Mtg5JkgQYzlRSAr/4\nBfzjH1FXIkmScFhTAAsXhtmbb70F3btHXY0kSQXFYU3VX5s2cMopcMMNUVciSVLRs3Om4NNPoV8/\nmD4d2raNuhpJkgqGnTM1zOabw157wd13R12JJElFzc6Zqrz6Kpx2Gnz8MTQxt0uS1BjsnKnhBg0K\n1589+2zUlUiSVLQMZ6pSUgK//KXLakiSFCGHNbWqxYvDshqvvAK9e0ddjSRJieewptZOq1Zw6qlw\n881RVyJJUlGyc6bVTZ8OO+8cltdo0ybqaiRJSjQ7Z1p7XbvCwIFw//1RVyJJUtExnKlmZ58NN90U\ndRWSJBUdw5lqNmQIfPklvPNO1JVIklRUDGeqWdOmcOaZds8kScozJwSodl9+GZbTmDYN2rePuhpJ\nkhLJCQFqPJ06heFN99uUJClvDGdas7PPDmue2amUJCkvDGdasz32CJugv/xy1JVIklQUDGdas5IS\nOOssJwZIkpQnTghQ3ebPDwvTfvwxbLRR1NVIkpQoTghQ42vXDo44Au68M+pKJEkqeHbOlJ0xY+Cn\nP4UpU8I1aJIkKSt2zpQbO+8M660HL74YdSWSJBU0w5myU1ICp58Ot90WdSWSJBU0hzWVvcqJAZMm\nhQVqJUlSnRzWVO60aweHHALDh0ddiSRJBcvOmernzTfhxBND96wkqf/6SJKUP3bOlFsDB0KLFu4Y\nIElSjhjOVD8lJXDGGU4MkCQpR5I6LuWwZpS+/hq23BIqKmDDDaOuRpKkWHNYU7m3wQZw0EFwzz1R\nVyJJUsExnKlhzjgDbr0V7GBKktSoDGdqmEGDQjB7442oK5EkqaAYztQwlTsG3Hpr1JVIklRQnBCg\nhps3D3r2hE8+gfXXj7oaSZJiKe4TAoYAE4EpwIW1PKcUeBf4ACjPS1VqmA4dYMgQuPfeqCuRJKlg\n5LNz1hSYBOwLzALGAEcDH2c8pz3wOrA/MBPoAMyr4b3snMXFiy/Cr34FEya4Y4AkSTWIc+esPzAV\nmA4sA0YAw6o95xjgIUIwg5qDmeKktBQWLYLRo6OuRJKkgpDPcLYp8FnG45npc5l6AhsALwFjgePz\nU5oarEkTOO00JwZIktRImuXxs7IZh2wO7AjsA7QB3gTeIlyjtoqysrL/3S8tLaW0tLQxalRDnHQS\n9OkD114L660XdTWSJEWqvLyc8vLyBr++IRcJNQGOBboBVwKbAxsBdY1rDQTKCJMCAC4GVgJXZzzn\nQqB1+nkAtwPPAA9Wey+vOYubww+HffeFs86KuhJJkmIlH9ec/QvYlXB9GMD36XN1GUsYtuwKtACO\nAh6v9pzHgEGEyQNtgAHARw2oUflWuWOAJElaKw0JZwOAnwGL0o+/JgxH1mU5cC4wihC47ifM1Dwz\nfUBYZuMZ4D3gbeA2DGfJsO++8O23MHZs1JVIkpRoDRnWfBvYjdAJ6wd0BJ5N388XhzXj6C9/gWnT\n7KBJkpShvsOaDQlnxwFHAjsBw4HDgUuBBxrwXg1lOIuj2bNhq61gxgwnBkiSlJaPcAawFWFGJcAL\nrLqQbD4YzuLqsMNgv/3gzDPrfq4kSUUgl+Fsg1peW5mSvq7He60tw1lcPfssXHwxvPNO1JVIkhQL\nuQxn0wlBrISwfMY36fPrAzMIS2vki+EsrlauhO7dYeRI2HnnqKuRJClyuVxKoyshgD0H/BjYMH0c\nmD4nhR0DTj/dSQGSJDVQQ645+wDYNotzuWTnLM4qJwZ8+im0bRt1NZIkRSofi9B+Tpid2ZXQSbsE\nmNWA91Gh2mgj2HtvuO++qCuRJClxGhLOjgY6AY8AD6fvH92YRakAuGOAJEkN0tClNKLmsGbcOTFA\nkiSg/sOazRrwGS/VcC4F7N2A91KhypwYYDiTJClrDemcZf6mbQUcRtg384JGqSg7ds6S4IsvYOut\nnRggSSpq+dohoLoxwC6N9F7ZMJwlhTsGSJKKXD5ma26QcXQAhgBupKiaOTFAkqR6aUjnbDpVWzYt\nTz++AnitcUrKip2zpKicGPDgg7DTTlFXI0lS3uVjWLMVsDiLc7lkOEuSP/8Zpk+3gyZJKkr5CGfj\ngB2zOJdLhrMkcWKAJKmI5fKas42BnYA2hCC2U/q2NH1OqtnGG8Nee7ljgCRJWahP5+wk4ETCUhpj\nM84vAO4i7BaQL3bOkmbUKLj4Yhg3LupKJEnKq3wMax4GPNSA1zUmw1nSrFwJPXrAiBHQv3/U1UiS\nlDe5DGfHA/cAv6Zqtmble6SAa+rxXmvLcJZEf/sbfPABDB8edSWSJOVNLrdvqryurC01hzNpzU45\nJXTP5s2DDh2irkaSpFhy43Pl18knQ58+cOGFUVciSVJe5OOas07A6UBXqjpvKeCUBrxXQxnOkmrs\nWDj8cKiogKZNo65GkqScy+WwZqXHgFeA54CV6XMmJWVn552hc2d46ik46KCoq5EkKXYa0jkbD+zQ\n2IXUk52zJLv77rDm2TPPRF2JJEk5l4+Nz/8fcGADXicFRx4Z1jubMiXqSiRJip2GdM6+J8zcXAos\nS59LAes1VlFZsHOWdBddBEuXwjX5XIFFkqT8y8eEgDgwnCXd9Onh+rNPP4U27v4lSSpc+QhnNW1w\nPh+YASxvwPs1hOGsEBx8cDhOOy3qSiRJypl8hLO3CJuev5d+vB3wIdAOOBsY1YD3rC/DWSEYNSoM\nb44bByVJbeJKkrRm+ZgQ8DlhtuZO6WMH4BNgMPDXBryfitXgwfD99/Dmm1FXIklSbDQknPUmdMoq\nfQT0ASpwvTPVR5MmcM45cP31UVciSVJsNGQs6QHgK2BE+vVHAh2B44DXgF0arbraOaxZKObPh27d\n4L33YLPNoq5GkqRGl49rztoAPwN2Tz9+HfgXsBhYB1jQgPesL8NZIfnlL6F1a7jqqqgrkSSp0bmU\nhpKnogIGDgzLa6yzTtTVSJLUqPIxIaAX8CDhWrNp6eOTBryPFHTvDoMGhW2dJEkqcg0JZ3cCNxPW\nNNsLGA7c25hFqQiddx5cdx2sXBl1JZIkRaoh4aw18DyhPTcdKMO9NrW29tgjDGk+/XTUlUiSFKmG\nhLPFQFNgKnAu8BPCRACp4UpKqrpnkiQVsYZMCOgPfAy0B/5A2PD8r4SdA/LFCQGFaOlS6No17Byw\n3XZRVyNJUqNwtqaS7Y9/hGnT4I47oq5EkqRGkctw9gRhB4CaXpMCDs7iPYYA1xGGRW8Hrq7lebsA\nbxIWuH24ps8znBWoefOgZ0+YOBE6d466GkmS1louw9lcYCbwX+Dtaq9PAS/X8fqmwCRgX2AWMAY4\nmjBEWv15zwELCTNDH6rhvQxnhexnP4P114c//SnqSiRJWmu5XOdsY+B3wLaE7tdgQmArp+5gBuFa\ntamEGZ7LCNs/DavheT8nrKM2tx61qZD85jdwyy3w3XdRVyJJUt7VJ5wtB54GTgAGEoLWy4QZm9nY\nFPgs4/HM9LnqzxkG3JR+bHusGG25JQweDLfeGnUlkiTlXX2X0mgFHAb8BzgH+AfwSJavzSZoXQdc\nRNW1bUmdsKC19dvfhmU1liyJuhJJkvKqWT2eew+wDfAUcCXwfj0/axbQJeNxF0L3LNNOhOFOgA7A\nUMIQ6OPV36ysrOx/90tLSyktLa1nOYq1fv1gm23g3nvhlFOirkaSpKyVl5dTXl7e4NfXpzO1Evih\nlp+lCOudrUkzwoSAfYDPgdHUPCGg0p2EGaLO1ixWL70EZ58NH30ETRqyXrIkSdHL5YSAJkDbWo66\nghmEa9bOBUYRNk2/nxDMzkwf0qpKS2G99eCxx6KuRJKkvEnqNV12zorFQw/BX/8Kb70VtniSJClh\nctk5k/LvkEPg22/h5WxWa5EkKfkMZ4q3pk3DzE0XpJUkFQnDmeLv+ONhyhR4442oK5EkKecMZ4q/\nFi3gd7+DK66IuhJJknIuqVdYOyGg2CxdGjZEv/9+GDgw6mokScqaEwJUmFq0gIsvtnsmSSp4ds6U\nHEuWhO7Zgw9C//5RVyNJUlbsnKlwtWwJF11k90ySVNDsnClZliyBHj3g4Ydhl12irkaSpDrZOVNh\na9kyXHt22WVRVyJJUk4YzpQ8p50W1j0rL4+6EkmSGp3hTMnTogVceWXooDm8LUkqMIYzJdPRR8PC\nhfD441FXIklSo3JCgJLr//0/uPBCeO+9sAenJEkx5IQAFY8DD4T114d77426EkmSGo2dMyXba6/B\nccfBpElhJqckSTFj50zFZdAg2HZbuOmmqCuRJKlR2DlT8n30Eey5J0ycCBtuGHU1kiStor6dM8OZ\nCsO550JJCVx/fdSVSJK0CsOZitO8ebDVVvDKK+FWkqSY8JozFacOHeB3v4Nf/zrqSiRJWiuGMxWO\nc86BqVPhmWeirkSSpAYznKlwtGgBf/976J4tXx51NZIkNYjhTIXloINgk03gxhujrkSSpAZxQoAK\nz8SJsMceMGFCCGqSJEXI2ZoSwCWXQEUFjBgRdSWSpCJnOJMAFi4MOwfccgsMHhx1NZKkIuZSGhJA\nmzZhQdqf/QwWL466GkmSsmY4U+E68EDYbju4+uqoK5EkKWsOa6qwffYZ9OsHb7wBvXpFXY0kqQg5\nrCll6tIFLrsMTj0VVq6MuhpJkupkOFPh+/nPIZWCG26IuhJJkurksKaKw+TJsNtu8Pbb0L171NVI\nkoqIw5pSTXr1gosugtNPd3hTkhRrhjMVj/POgx9+gJtvjroSSZJq5bCmisvEiTBoELz2GvTpE3U1\nkqQi4LCmtCZ9+sAf/gDHHgtLl0ZdjSRJq7FzpuKTSsHBB4ftnf7yl6irkSQVOPfWlLLx5Zewww7w\n3//CnntGXY0kqYA5rCllo1MnuP12OOEE+OqrqKuRJOl/7JypuF1wAXz0ETzxBDTx/1UkSY3PzplU\nH3/+M8yf7+bokqTYiCKcDQEmAlOAC2v4+bHABOA94HVg+/yVpqLTvDmMGAH//CeUl0ddjSRJeR/W\nbApMAvYFZgFjgKOBjzOesyvwETCfEOTKgIHV3sdhTTWuZ5+Fk0+GsWNh442jrkaSVEDiPqzZH5gK\nTAeWASOAYdWe8yYhmAG8DWyWr+JUxPbbD846Cw49FBYvjroaSVIRy3c42xT4LOPxzPS52pwKPJXT\niqRKl14Km28eQpqdWUlSRPIdzurzG28v4BRqvi5NanwlJXDnnTBhAlx7bdTVSJKKVLM8f94soEvG\n4y6E7ll12wO3Ea45+6amNyorK/vf/dLSUkpLSxurRhWzddaBxx6DgQNhq61g6NCoK5IkJUx5eTnl\nazHJLN8TApoRJgTsA3wOjGb1CQGbAy8CxwFv1fI+TghQbr3+erj+bNQo6Ncv6mokSQkW9wkBy4Fz\ngVGEGZn3E4LZmekD4PfA+sBNwLuEACfl1+67w003wY9/DNOnR12NJKmIuEOAtCb//GcIaa+/Dhts\nEHU1kqQEcuNzqbFdcAG8+WZYC61Nm6irkSQljOFMamwrV4YFaj//POzB2apV1BVJkhLEcCblwvLl\ncOyx8MMP8PDD0KJF1BVJkhLCcCblyrJlcPjhVftxNsv3SjSSpCSK+2xNKbmaN4cHHoAFC+DEE0M3\nTZKkRmY4k+qjZUt45BGYNw+OOgqWLIm6IklSgTGcSfXVpg08/njYf/OQQ2DhwqgrkiQVEMOZ1BAt\nW4Yhzo4dwxZP330XdUWSpAJhOJMaqlkzuOsu2GYb2HtvmD076ookSQXAcCatjSZN4MYbw/DmwIHw\nwQdRVyRJSjiX0pAay733wnnnhdvBg6OuRpIUEy6lIUXl2GPhwQfh+OPh1lujrkaSlFB2zqTGNnky\nDBsGgwbB9de73ZMkFTk7Z1LUevWC0aPhm29gjz3g00+jrkiSlCCGMykX2raFkSPDQrUDBsBzz0Vd\nkSQpIRzWlHLtpZfCdWhHHw1//GNYI02SVDQc1pTiZq+9YPx4mDIFdt0VJk6MuiJJUowZzqR86NAh\n7Ml51lnhOrQbb4SVK6OuSpIUQw5rSvk2aRKcckpYwPb226F376grkiTlkMOaUtz17g2vvhomCwwa\nBFddBcuWRV2VJCkmDGdSFJo0gXPPhbFjobwc+vWDF1+MuipJUgw4rClFLZWCRx+F88+HnXaCv/8d\nunaNuipJUiNxWFNKmpISOPRQ+Ogj6Ns3BLTLL4cFC6KuTJIUAcOZFBetW8Nll8G778Inn0CPHnDN\nNbBoUdSVSZLyyHAmxc3mm8M998ALL4SJA716wW23OWlAkoqE15xJcff223DppWFD9d/8Bk49Fdq0\niboqSVKWvOZMKjSVe3M+8EDYCqpbt7AN1DffRF2ZJCkHDGdSUgwYAA8/HJbeqKiALbcMOw68/37U\nlUmSGpHhTEqarbaCO+8Mszs32QSGDIEf/Sh01rwuTZISz2vOpKRbtiysk3bjjWFrqGOPhRNPhO22\ni7oySRL1v+bMcCYVkkmTYPjwMNuzQ4cQ0o45Bjp1iroySSpahjNJsGJFuDZt+HB47DHYeWc47LCw\n2O3GG0ddnSQVFcOZpFUtXAijRsGDD8JTT8G224agdsghbhMlSXlgOJNUuyVL4Pnnq4La+uvD0KHh\n+NGPoFWrqCuUpIJjOJOUnZUrYdw4ePrpcHzwQQhoe+0FpaWwww7QtGnUVUpS4hnOJDXM11+Hrlp5\nObz8MsyaBYMGwZ57hqNfP2jePOoqJSlxDGeSGseXX8Irr4Sg9vLLYeHbvn3DYrj9+4fbbt2gJKn/\nGZGk/DCcScqNBQvgnXdg9Oiw3+fo0bBoEeyySwhtffvC9tuHjdrtsEnS/xjOJOXP55/DmDHw3ntV\nx2efQe/eIahtvz306RMed+0KzZpFXbEk5Z3hTFK0Fi6EDz+ECRNCWJs0KRyzZ4f9QHv1CmGtd+9w\nv1u3sPZaE3eTk1SYDGeS4mnRIpgyBSZPrgpskyfD9Onw7bfQpUvorlU/ttgCNtrIrpukxDKcSUqe\nhQvh009DUKs8Zsyouj9vXtiOapNNaj86dw7PadEi0j+KJFUX93A2BLgOaArcDlxdw3P+CQwFFgIn\nAe/W8BzDWYKVl5dTWloadRlqgMi+u+XLYc4c+OKLcJ1b9WPWrDC7dN48aNMmhLSOHcNR0/0NN4T2\n7aFdu3DbunVRzDr1715y+d0lW33DWT7HCZoCNwD7ArOAMcDjwMcZzzkA6AH0BAYANwED81ij8sD/\nyCRXZN9ds2aw6abhWJNUCubPh7lzQ1CbO7fq/pw54Vq4uXPhq6/C8779NhwrVoSQVv2oDG/t2kHb\ntrDOOrDuujUflT+L8UxV/+4ll99dcclnOOsPTAWmpx+PAIaxajg7GBievv820B7oDMzJT4mSEq2k\npCpY9eyZ/euWLFk1rGUelefnzYPvv6/7aNp01dDWunU4WrUKR+X9ms7V9vNWrcJwbfPmq97WdK55\n86LoAkqFLJ/hbFPgs4zHMwndsbqesxmGM0m51LIldOoUjrWRSoWg98MPIagtWBAmQixevPpt9fvf\nfLP6+cyfL1sGS5euflv93LJlodNYU4j77ruwr2rl42bNQphs2nTV+3U9bsjPSkrCjNwmTRrn/tq+\nR2WArel2TT+L6rnffQczZzbsfRuL75c3+fzkwwjXnJ2efnwcIZz9POM5TwBXAa+nHz8P/BYYV+29\npgLdc1apJElS46kgXLaVlXx2zmYBXTIedyF0xtb0nM3S56rL+g8oSZKkmjUjJMeuQAtgPLBVtecc\nADyVvj8QeCtfxUmSJBWjocAkwrDkxelzZ6aPSjekfz4B2DGv1UmSJEmSJElJNQSYCEwBLoy4FtXt\n34SZtu9nnNsAeA6YDDxLWC5F8dQFeAn4EPgA+EX6vN9h/LUiLEc0HvgI+Ev6vN9dsjQlLMT+RPqx\n318yTAfeI3x3o9PnCva7a0oY7uwKNKfma9YUL3sA/Vg1nP2VMAMXQsC+Kt9FKWsbATuk769LuCRh\nK/wOk6JN+rYZ4frdQfjdJc35wL2EBdvB7y8pphHCWKaC/e52BZ7JeHxR+lC8dWXVcDaRsLAwhF/+\nE/NdkBrsUcIOH36HydKGsCPLNvjdJclmhOWk9qKqc+b3lwzTgA2rnavXd9ckB0XlSk0L1Naxl4ti\nKHPHhzlU/cuqeOtK6IK+jd9hUjQhjDDMoWp42u8uOa4FLgBWZpzz+0uGFCFYj6Vqbdd6fXf5XOds\nbbnTeeFJ4feaBOsCDwG/BBZU+5nfYXytJAxLtwNGETowmfzu4uvHwJeEa5ZKa3mO31987Q58AXQk\nXGdWvUtW53eXpM5ZNovYKv7mEFq6ABsT/gOk+GpOCGb3EIY1we8waeYDTwI74XeXFLsR9pqeBvwX\n2Jvwd9DvLxm+SN/OBR4h7C1er+8uSeFsLNCTqkVsj6LqIkklx+PAien7J1L1C1/xUwLcQZjtd13G\neb/D+OthVROPAAAgAElEQVRA1Wyw1sBgQhfG7y4ZfkdoQHQDfgq8CByP318StAHapu+vA+xHuO66\noL+7mhaxVXz9F/gcWEq4XvBkwgyW5ynA6cQFaBBhaGw84Rf7u4TlbPwO4287wp7E4wlT+i9In/e7\nS549qWpE+P3FXzfC37vxhCWIKrOK350kSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSfm1gqo11N4F\nNo+2nEa1HfDv9P2TgOur/bycsGp/bR4grJkkqcgkaW9NSYVnIWFD9ZqUpG+Tun/gBVQFspr+DHXt\nr3cbcB7wi0auS1LMJWn7JkmFrythF5DhhC1PuhBCzmhgAlCW8dxL0s99FbgP+HX6fDlVHakOhP0J\nAZoCf8t4rzPS50vTrxkJfAz8J+MzdgFeJ6z2/RZhE/iXgb4Zz3mN0CXL1BIYCIyp489bAhxEVedw\nEvBJxp/jgDpeL0mS1KiWUxVMHgK2IAx19k//fD/glvT9JsATwB6E8PUe0Iqwj90U4Pz0814Cdkzf\nzwxnZxACHYTwNIYQBkuBb4FNCGHpDcLG0y2ACqqC3rqEgHcCcG36XC9qDmAD07VWOomw0XHmEO6C\njDor3Q+cnfH4ZWCrGt5fUgFzWFNSlBax6rBmV2AGobsFIZztRwgzEDYS7kkIZA8Di9PH49RtP0KH\n6/D04/WAHsCy9Od9nj4/nnCt1wLgC+Cd9Pnv07cPApcROnqnAHfW8FlbpF9bKQWMYNUhypeqvea3\nhGHemzLOfU74Z/JxHX82SQXEcCYpbn6o9vgvwK3Vzv2SqmvSqHZ/OVWXbLSq9rpzgeeqnSsFlmQ8\nXkH4b2Nt14MtTL/HIcARrN79Iv3akmrnqj/OtC9wGPCjGl6zcg2vk1SAvOZMUpyNInSn1kk/3hTo\nCLxCCEeVw5o/znjNdGDn9P3DM86PAn5G1f+U9gLa1PK5KcL1XxtnvFdbwrAmwO3APwkdt/k1vH4G\nsFHG4zUFsy2AG4EjWTUkkv78GWt4raQCZOdMUpRqm8VY6TnCNVdvph8vAI4jDHPeT7iw/0vCdV+V\nAejvhGUozgCezHi/2wlDhOPSz/0SOJTaZ00uA44izLhsTeiYDSZ09sYRQllNQ5qk6+pd7c9U02eU\nACcCGwCPps/NIoTN5sBmwMRaPkOSJCm2LqdqtmY+bELorK3JXcCAtfiM/YB/rMXrJSWUw5qSCkW+\n1kM7gbCsxu/qeN7fgbPW4nNOo2pWqCRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkpLrQWAOUAF8ACwAbou0IkmSpDxpGnUB\n1bQG1gNaAmcBTwPtgJOjLEqSJClfSqIuoBaTgG2BU4F1gb9n/rB79+6pioqKKOqSJEmqrwqgR7ZP\nbpbDQhqqHbAYWEYIaI9Uf0JFRQWpVCrfdakx9OhB2ezZlF16KVx0UdTVqJ7KysooKyuLugw1kN9f\ncvndJVtJSUn3+jy/Sa4KWQvbA++n73cHZkVYixpbKgVHHAG33QY33RR1NZIkxU4cO2fbAxPS98cA\njwLbETppSrpUCtq2heefhx/9CNZbD449NuqqJEmKjTiGsxsz7v8+faiAlA4YAN26wahRsPfeIaAd\ndFDUZSkLpaWlUZegteD3l1x+d8UlrhMC6pLymrOE6tYtdM26p4ffx46FAw6A+++HvfaKtjZJknKg\npKQE6pG54njNmQpdSca/nzvvDCNHwlFHwVtvRVeTJEkxYThTfqVSq4YzgD33hOHDYdgwGD06mrok\nSYoJw5niYehQ+Pe/w7VnY8dGXY0kSZExnCm/auqcVTrwQLj99nD7zjv5rUuSpJgwnCm/6prIcdBB\ncOutYZLAuHH5qUmSpBiJ41IaKnS1dc4qDRsGK1eGgPbMM7DDDvmpS5KkGDCcKb+yXQLl0ENDQBsy\nJKyH1rdvbuuSJCkmDGfKv7o6Z5UOOyzc7r8/PP449O+fu5okSYoJw5nya00TAmpy2GHQqhX8+Mfw\n4INhyydJkgqYEwIUfwceCCNGwOGHhyFOSZIKmOFM+VXfzlmlvfeGRx+F448Pt5IkFSiHNZUcu+0W\nZm8ecAAsWgRHHx11RZIkNTrDmfKroZ2zSjvuGDZOHzIEvvoKzj238WqTJCkGDGfKr2yX0liTbbeF\nV18Nszhnz4Y//GHtAp8kSTHiNWfKv8YIUt26weuvw7PPwmmnwfLla/+ekiTFgOFM+bW2w5qZOnaE\nF1+Ezz8Pi9YuXNg47ytJUoQMZ0q2ddcNC9RusAHsu2+4Dk2SpAQznCm/GrNzVql5c7jrLthjjzCj\nc+rUxn1/SZLyyHCmwlBSAldfDeefD4MGhQkDkiQlkOFM+ZWLzlmmM8+Ee+4JuwncfXfuPkeSpBxx\nKQ3lV2MspVGXwYPhpZfgoINg8mS48kpo4v+HSJKSIZ+/sQ4ARgPvA10zzl8LDMhjHYpaPtYk23pr\neOutENKOPjrsKCBJUgLkM5w9RQhhDwDbps/1BtoDb+exDkUp18OamTp2hBdegBYtYPfdYcaM/Hyu\nJElrId9jPbcARxA6aABXApem73cAxgA/ABOA6cAXwJn5LVEFpVWrcO3Z8cfDgAFhXTRJkmIsn+Gs\nP3A28DfgeGAw8DEwK31+EnAxMAfoC9wMXEYIdCoU+eycVSopgfPOg/vug2OOgWuvzc+1b5IkNUA+\nw1lvYBxwMvAY8GvgakLH7FRCt2wy8Gn6+X2Bd/NYnwrd3nvD22+H2ZzHH++OApKkWMrnbM170gfA\nWcC9wCLgj8AJwI2EQPZe+jlbEyYP1KisrOx/90tLSyktLW3sepULUXTOMm2xBbz2GpxxRliwduRI\n6NkzunokSQWnvLyc8vLyBr8+it+S7YD7gAMJnbsHgRnA4cBthGHOOwjXnHWt5T1SKYelkmn99cMK\n/htuGG0dqRTcdBNcfjnccAMcdVS09UiSClZJaEpknbmiCGd/JQxrvl7t/EvAXMI1aWPSjz8jdNWq\nM5wl1frrQ0VF2AszDsaNgyOPhP32g2uuCRMIJElqREkIZ43BcJZU7dvDtGkhpMXF/Plw2mmhozdy\nJPToEXVFkqQCUt9w5rLpUrt28MADIaDtumuY1SlJUkTsnCm/2rULi8G2bx91JTUbNw6OPRb69YN/\n/Su+dUqSEsPOmbQ2dtwR3nknXBPXty+sxWwbSZIawnCm/Ip6KY1stGkTZnDefHPoov32t7BkSdRV\nSZKKhOFMqs3QoTB+PEyeHLZ++uCDqCuSJBUBw5nyKwmds0wdO8Ijj8DPfw577QV/+AMsWxZ1VZKk\nAmY4U34lcSJHSQmcemqYLPDmm7DLLvCuO4tJknLDcKb8S1LnLFOXLvDkk3D++bD//nDppV6LJklq\ndIYz5VfShjWrKymBE06A996Djz4KS268+WbUVUmSCojhTGqIjTaChx6CsjI47DA480z4+uuoq5Ik\nFQDDmfIr6Z2zTCUlYV/Ojz6C5s1h663hrruSeV2dJCk2kvpb0h0CkqpNG5g7F9ZZJ+pKGt/YsXDW\nWeHPeNNNsM02UVckSYoBdwhQvBVS56y6nXeGt9+Go46C0lK44AL47ruoq5IkJYzhTPlV6B3Ppk3h\nnHPg/fdDh7B3b7j9dlixIurKJEkJkdQWhsOaSdWqFXzzDbRuHXUl+TF2LPzqV7BwIVx3HfzoR1FX\nJEnKs/oOaxrOlF8tW8L8+SGkFYtUCh54IOzR2b8//PWv0K1b1FVJkvLEa86kuCkpCdehTZwIffuG\na9POOw/mzYu6MklSDBnOlF+FPCGgLq1bh10FPvww7M/Zp0/Yq/P776OuTJIUI4YzKd822ghuuCHM\n7Jw4EXr2DI+XLo26MklSDBjOlF/F3Dmrrnt3uPdeePppeOqp0Em75x5YvjzqyiRJEUrqb0knBCRV\n8+bwww/QokXUlcTPyy/DZZfB7Nlh+POYY6BZs6irkiStJWdrKt6aNYNFi0JI0+pSKSgvhyuugFmz\n4JJL4LjjDGmSlGCGM8Vb06awZIlhIxuVIe3TT0Mn7bjjDLWSlEBxXkqjHzASuLza+X8CO+exDikZ\nSkvhpZfg3/+G//wHevSAf/zD2Z2SVODyGc7eBW6odq4PsA4wNo91KEpOCKi/PfeEF16AkSPh1VfD\nAra//z18+WXUlUmSciDq2ZpXAL9P3+8DTAAWE4Lc58AM4KBoSpNipn9/ePBBeOONEMz69An7eFZU\nRF2ZJKkR5TucZbZM9gEmAn2Bm4A9gPOBDwlDoI8DpwJP5LlG5ZKds7XXsyfcfDN8/DG0bw8DBsAR\nR8BrrxX+xvKSVATyGc66A8cD2wD7A78Brga+A84GTgKWAe+nn9+X0EGTVJPOneFPf4Jp08LQ5ymn\nwE47wV13weLFUVcnSWqgqFoYpwErgX+nHx8A7AbMAZoB1wKTgN61vD51+eVV8wpKS0spLS3NVa1q\nTCUlsHKl3bNcWLkSRo2Cf/4Txo2D00+Hs8+GTTeNujJJKirl5eWUl5f/7/EVV1wBMV9KY13CrM2h\n6ce7EbpmZwG3AfcRhjYfBQbW8h4upZFEqRQ0aeLQWz5MmhS2hLr3Xhg8GM46K8z+NBRLUt4lYZ2z\nK4Fy4MX040mESQArgRWE0LaAMIPzMeCSGt7DcJZEhrP8mz8f7r4bbrkl7N15xhlw4onQsWPUlUlS\n0UhCOGsMhrMkWrkyLELrd5d/qRS8+Sbceis8+igMGQJnnmk3TZLywHCm+DKcxcM334ThzltuCbs1\nnHZa2H1gk02irkySCpLhTPG1YkXYfmjlyqgrEYSQ/NZbcMcd8PDDYR21E0+EYcOgTZuoq5OkgmE4\nU3wZzuJr4cIw3Dl8OIwZAz/5SQhqgwY57ClJa8lwpvhavhxatgwhTfE1a1bYy3P48DDsecIJcPTR\n0KtX1JVJUiIZzhRfy5ZB69YhpCn+UikYOzbM9hw5MqyX9tOfwpFHwhZbRF2dJCWG4UzxZThLrhUr\n4OWXYcSIcH1a794hqB1xBGy0UdTVSVKsGc4UX0uXwjrrhJCm5Fq6FJ5/PgS1J56Afv1CUDvkEOjU\nKerqJCl2DGeKL8NZ4Vm0CJ5+Gu6/P2wdtf32cOih4ejaNerqJCkWDGeKryVLoG3bENJUeBYvhhde\ngEcegccfD9eoHXpomPm5zTbO+pRUtAxnii/DWfFYvhxefz0EtUcegRYtQlAbNgwGDgyLEUtSkTCc\nKb4WL4Z27UJIU/FIpeDdd8NEgieeCEt17L8/HHhg2EZqgw2irlCScspwpvgynAngs8/gqafgySeh\nvDxcp3bggeHYbjuHPyUVHMOZ4mvRIlh//RDSJAj/LpSXh6D25JNhOPSAA8JRWgrrrRd1hZK01gxn\niq+FC2HDDUNIk6pLpWDixBDSRo0K+37usAMMHhyOXXaBZs2irlKS6s1wpvgynKk+Fi6E116DZ5+F\n556DTz+FvfYKQW2//aB796grlKSsGM4UXz/8AB07hl+6Un3Nnh0Wv33uuRDYWrcOQW3vvWHPPd2p\nQFJsGc4UX4YzNZZUCj78MAS18nJ45RXYeONwnVppaQhrnTtHXKQkBYYzxdf334dfmD/8EHUlKjQr\nVsCECSGolZfDq68a1iTFhuFM8WU4U77UFtZ+9CPYffdwdOvmsh2S8sJwpvhasCD8gvz++6grUbFZ\nsQLGjw8TDF5/PdymUiGkDRoUbnfYAZo3j7pSSQXIcKb4+u67sN/iggVRV6Jil0rB9OkhqFWGtenT\nYeedq8LarruGRZMlaS3lKpytA3QBUsBMIOpxKcNZEhnOFGfffgtvvlkV2MaMCUOfAwZA//7hdptt\nXGtNUr01ZjhrC5wO/BToAMxJP78z8BVwL3AbEMUYleEsiebPhy5dQkiT4m7ZsnDd2ujR8Pbb4Zg1\nC3bcsSqsDRgAm23mtWuS1qgxw9kLwAjgcUIwy7QRcDBwFLBPlp+1HXABUAFckXH+b8BDwFtZvg8Y\nzpLJcKak++YbGDu2Kqy9/XbopGWGtZ13dtspSauI+zVnewKlVIWz7sBlwEn1fB/DWRJ9+y1ssUUI\naVIhSKVgxoxVw9qECbDJJqHDttNO4XbHHcO+spKKUn3DWTYXTwwCxhOGL48HdgSuA2Y0oL7qrgQu\nTt/vADwNbA1MBdoBLYEy4JZG+CxJalwlJdC1aziOOiqcW74cJk2Cd94Jx+OPh5minTqtGth22ils\nZyZJ1WST4t4Htk8fdwG3A0cSumD1VZp+3RXp2/2AS4BzCEHtKOBWYEvgImBe+vOqs3OWRN98Ey6w\n/vbbqCuR8mvFCpgypSqwjRsH774L7duHkLbTTmEpj759w6QZr2GTCkouOmfLCbM0DwFuJISlUxtQ\nW3dC560tMBT4BXAEsDEhnE0AJgOfpp/fF/h7bW9WVlb2v/ulpaWUlpY2oCTlnb90VIyaNoU+fcJx\n7LHh3MqVUFFRFdiuuy4Mia5YEUJa376w/fbhduutoVWraP8MkrJWXl5OeXl5g1+fzW/KV4BngJOB\nPYC5hGHO7Rr8qeG9mhO6ZNcTQt9NwDXAYEJwmwDsAiyt4fV2zpLo66+he/fQQZO0ulQK5swJIS3z\nqKgIf3eqh7aNNvJ/eKQEyMWEgI2BY4DRwKvA5sBewPAG1AdhzbSHgSGEjtwUwvVsPYCrgS+AO4Dp\nQNda3sNwlkRffQU9e4aQJil7S5bARx+tHtqaNg0hbbvtYNttwzpsW28NbdtGXbGkDHGfrQnh2rLX\ngGernX+J0JX7GzAm/fgz4IQa3sNwlkSGM6nxpFLw+echpL3/Pnz4YTg+/jhMPthmm6pj221hq62g\nTZuoq5aKUi7C2WHAVYTFZyufnwKiXMjHcJZE8+ZB794hpEnKjRUrYNq0ENQ++KAqtE2eHJb4qOyw\nVR59+ng9m5RjuQhnFcCPgY8bWFMuGM6SyHAmRWf5cpg6tSqsVQa3Tz6BzTcPnbWttqqauNC7d5hN\nKmmt5SKcvQ7s3tCCcsRwlkRz54b/+M+bF3UlkiotXRq6ahMnrn60bVsV1jKPLl2gSZOoK5cSIxfh\n7B+E7ZoepWrmZIpwUX9UDGdJNHduuFh57tyoK5FUl1Qq7CVaU2j75hvo1Wv10NarF7RuHXXlUuzk\nIpzdlb6tnoZOzvZDcsBwlkRffhmucTGcScn23Xc1d9sqKqBz5zDxp1evcFt5dOsGzZtHXbkUiSTM\n1mwMhrMkmjMnTPn/8suoK5GUC8uXw/TpYTeE6sesWbDZZqsGtsoQt/nmYQN5qUDlIpy1JuwIsHX6\nfmUqOqW+xTUiw1kSGc6k4rV0aZhFOmVK6LplBrc5c8L+pNWDW8+eXt+mgpCL7ZvuIczUHELYE/M4\n4jVzU0mRSrmauVSsWrQIM0B79179Z4sXhyHRyrD27rvwwAPh/ldfheC25ZZhl4TM2y239Bo3FaRs\nflOOB3YA3iNsft6csIjsgBzWVRc7Z0k0e3ZYzXzOnKgrkZQUCxeGjltFRVj2I/N2+nTYYIPVQ1vl\nbadO/g+hYiEXnbPKGZrzCftpzgY61rsyyc6ZpPpq06ZqwdzqVq4M17JlhrYnn6x6vHhxVYetenDr\n2jV086QYyuY35enAQ4RgdhewLnAZcHPuyqqTnbMk+uIL2HHHcCtJuTZ/fghq1Ttun3wCM2eGmaVd\nu656dOsWbjfbzNmlajTO1lR8ff457LST4UxS9JYtC1236dOrjmnTqu7Png0bbbR6aKs8NtvMGabK\nWmMOax5PmAzw64xzqfSbp4BrGlCfipnDmpLionnzqqBVk2XLQnctM7C9+GLV/Tlzwl6ltXXeNt0U\nmjbNwx9EhWhN4axN+rYtqy5AW8LqC9JKklQ4mjcPQatbt5p/vnQpfPbZqp23556ruj93bghoW2wR\n1nGrPCofd+kC66yTtz+OkiWpbQyHNZNo5kwYMCAMJUhSIVuyBD79tOZjxowQ7NZZZ9XAVv3o3Nk1\n3gpEYw5rXp9xv3I4s/I+wC/qVZkkScWiZcuqhXRrkkqF7lr10PbGG1WP588P17ZVD22Z3bc2bWp+\nfyXamsLZO+nb3Qi7A9xPCGhHAB/muC4VIq85k6SgpCSsw9apE+y8c83PWbQojDjMmFEV2N54A0aM\nCPc/+wzWXXf14NalSwh1XbqESQ1OXEicbH5Tvg0MApalH7sIrRrms89gt93CrSRp7axcWXP3bebM\n8N/ZmTPDzzt3DmGtMrBVvzXA5VwuFqFtD6wHfJV+3DZ9TqofA7UkNZ4mTULw6twZdtml5ucsWxaW\nL6oMa599FkLc669XPZ43zwAXM9n8k74KGAeUpx/vCZTlqB4VOoc1JSl/mjevGvKszbJlYR3KmTMN\ncDGR7W/KjQnDmCnCMOfsnFWUHYc1k2jGDNhjj/CXXpKUHJkBLrMLl3lbW4DbdNOqY5NNwmSJIpOL\nYU2AxcAXQCugV/p4pb7Fqcg5IUCSkql58zDZYIstan/O0qVhCLUysFV24N54IyyhNGtW2HmhfftV\nA9umm64e4tq3L+rfF9mEs9MJy2ZsBowHBgJvAnvnsC5JkpQkLVrUHeBWroQvvwxBbebMqtD2yiur\nnluxYvUAVz3EFfAwajax9ANgF0Ig2wHoA/wFODSHdQ0BrgOaArcDV1f7ucOaSTRtGuy1V1g9W5Kk\n2ixYUBXcqge5ymPuXOjYcfWuW/UQt+66Uf9pcjKsuRhYlL7fCpgI9K53ZdlrCtwA7AvMAsYAjwMf\n5/AzJUlSXLRtC336hKM2y5aFYdLqoe2DD1Z93Lx5zUOnmec6dozVbgzZhLOZwPrAo8BzwDfA9BzW\n1B+YmvEZI4BhGM6Sz2vOJEmNpXnzMOmgS5fan5NKwbffrt59Gz8ennyy6tx334Vh0tpCXOVkhtat\n8/JHyyacHZK+LSMsp7Ee8EyO6gHYFMhcpXQm0S54q8ZkOJMk5UtJCay/fji23bb25y1ZEmajVh9G\nHTOm6twXX4T9UDfZZNXAVj3Adeq01l24usJZM8I1Z5V9xfK1+rTseDFZofI6QUlSHLVsCd26haM2\nqRR89dWqQ6affw7jxsETT1Sdmz8/dOEyg1s91RXOlgOTgC2AGfV+94aZBWT2KLsQumerKDvttNB6\nBEpLSyktLc1LcVoLDmtKkpKqpAQ6dAhH3761P2/JEsofeYTyZ58NExsmTar/R2XxnFeBfsBo4If0\nuRRwcL0/LTvNCIFwH+Dz9OcezarXnKVSHTvCY4/BrrvmqAw1uqlTYf/9oaIi6kokScqbXMzWvKyG\nc7kcn1oOnAuMIszcvIOaJgMMHw7DhsGjj4bNtBV/ds4kSarTmn5TllB3CMvmObkQ1jl75hk44QQD\nWlJMmQJDh4YOmiRJRaK+nbM1TScoBy4gbNVUXW/gQuDletTW+IYMgbvvhkMOCdtDKN7snEmSVKc1\nhbP9gK+AGwn7ak4GpqTv3wDMISwUG60hQ+CeewxoSWE4kyRpjbL9TdkU6JC+Pw9YkZtysrb69k2j\nRsHxx8Mjj8Duu0dTldZs0iQ46CCYPDnqSiRJypvGHNbMtILQKZtD9MGsZvvvHzpohx4Kr78edTWq\njZ0zSZLWKD4bSTWG/feH//wnBLTXXou6GlXnIrSSJNWpsMIZwH77hYD2k58Y0OLGCQGSJNUp23DW\nlaqL/9sQ9teMLwOaJElKqGzC2RnASOCW9OPNgEdyVlFj2W8/uPdeA1qc2DmTJKlO2YSzc4BBwHfp\nx5OBTjmrqDENHgz33RcC2quvRl2NwHAmSVIdsglnS9JHpWZEsytAw+y7b1VAe+WVqKspbk4IkCSp\nTtmEs5eBSwjXmg0mDHE+kcuiGt2++8J//wuHH25Ai5qdM0mS1iibcHYRMBd4HzgTeAq4NJdF5URl\nQDvsMANaVOycSZJUp2zC2QrgVuBwwuSA0SRpWDPTPvvAiBF20KLihABJkuqU7bDmesAGwDvAbcC1\nuSwqp/bZxw6aJEmKrWzCWTvCTM2fAHcD/YnDhudrI7OD9tZbUVdTPOycSZJUp2zCWVNgY+BI4Mn0\nuWQOa2baZx8YPhwOPhjeeSfqaoqH4UySpDXKJpxdCYwCKgjXm3UHpuSyqLwZOhRuuw0OPBAmTIi6\nmsLnhABJkurULIvnjEwflSqAw3JTTgSGDYOlS2HIEHjhBdh666grKmx2ziRJWqNswllr4FRg6/R9\nCMOap+SqqLw74ghYsiRs+fTSS9CzZ9QVFSY7Z5Ik1SmbYc17gM7AEKCcsLfm9zmsKRrHHQdXXBGu\nRZs2LepqCpMTAiRJqlM2nbMehDXOhgHDgfuAwtxJ/NRTQwdtn33g5ZehS5eoK5IkSUUmm3C2NH07\nH9gOmA10zFlFUfvZz0JA23vvsA7axhtHXVHhsHMmSVKdsglntxEWoL0UeBxYF7gsl0VF7rzzYPHi\n0EErL4dOnaKuqHAYziRJWqNsrjm7DfiasFNAN0LX7OZ6fs5JwHhgi4xzrQn7dGZTQ/5dfHGYKLDv\nvvDVV1FXUxicECBJUp2y6Zy1Iiyd0ZWwIG0JYbbmlfX4nLuAPaud+w1wHbCyHu+TX2VloYM2dGhY\nZqNt26grSj47Z5IkrVE2XavHgIOBZcAPGUdDVLZONga2B57N+NlPgU+Ab4B3gQXAB+nnRqOkBK66\nCnbcMayHtnhxZKUUBDtnkiTVKZvO2abA/o30eZVtkz8AvyeEwzeBf6VrKQN2AM4nBLOdgCWN9NkN\nU1ICN94Ylto46ih48EFo3jzSkhLNzpkkSWuUTefsDUKXa20cA2wI/Jawafoi4GPgLMK6aSngDsI1\nbROANoThzmiDWaWmTeHuu2H5cjjlFFgZ35HYWLNzJklSndbUOXs/fdsUOBmYRlVYSlF3YPsZcHr6\nuQcQ1keD/9/enUdJVd0JHP+2CLIJiguIyjIjInpsIE7ABbF1BKKACrK4Kzo4Rj3quJsx2jO4xnjI\nKLzcrsQAABIGSURBVHFwRkFH44qIuCCaCaWDmmBEYTRCjEZ0RBGPUYkbgj1/3NdW0XbXa7qr69Wr\n/n7O6VPvvX7Lr3OP1C+/++69oZv0NKAH8ENCBa22nFJJGBG6F7A8382rq6u/266qqqKqqiomnGZq\n2xYefDC8f3bOOXDzzVaBNpdTaUiSWoFMJkMmk2ny9fm+KfuQfUesvvPebsLzxhMmtf0ZMJGQmHWJ\nntOPkBBWRucNAi5t4D41NUlVYT77DA4+OCRpV12VTAxptWQJnHUWvPhi0pFIklQ0FaEw0ejqRL7K\n2RpCt+NuhCrW7cCGZsTWDvgxMDrar11Q/WRCcrYVYdDBemAFIYF7gVBpKx1dusCTT8Lw4bDNNnDh\nhUlHlB5WziRJipUvObuTkCgtJnRL7gmc24xnrQdGNPCcWntHn68AvZrxrJa1ww7w9NNw4IHQtStM\nnZp0ROlhciZJUl75krMBZJOl2wD7onLtsgs89RQcdFBI0CZNSjqi0ueAAEmSYuVLzjY0sK1a/frB\nggUwYgR06xZWE1B+Vs4kScor31QalYSJYGt/9s7Z/qzlQ0uJgQPD3GfHHQdLlyYdTWmzciZJUqx8\nyVkbYOucny1ztru0fGgpMnw4zJwJY8bAm28mHU3pckCAJEmxGrNCgBpj/Hj48EMYNQqeew66d086\nIkmSlEImZ4V0xhnw/vswejQsWuRC6XVZOZMkKVZjlm/S5qiuhn32gaOPhvXrk46m9JicSZKUl8lZ\nodUulN6xI0yZ4jqcuRwQIElSLJOzlrDllnDvvfDOO3DRRUlHU1qsnEmSlJfJWUvp0AHmz4eFC+HG\nG5OOpjRYOZMkKZYDAlrSttuGdTj32w969YKJE5OOKHlWziRJysvkrKXtsgs8+iiMHAk9e8IBByQd\nUXKsnEmSFMtuzWIYNAjuuiuM4HzjjaSjSY5TaUiSFMvkrFhGjYJp0+Cww2Dt2qSjSY7JmSRJeZmc\nFdPUqTBpEhxxBHz5ZdLRFJ/dmpIkxTI5K7arroK+feHEE1vnHGhWziRJysvkrNi22AJmzw5dm61t\nDjQrZ5IkxTI5S8JWW8HDD8Pjj8OMGUlHU1xWziRJysupNJLSrRssWBCm1ujdG8aOTTqilmflTJKk\nWFbOktS3b6ignXoqLFuWdDQtz6k0JEmKZXKWtKFD4eabwwjODz5IOpqWZ3ImSVJeJmel4JhjYMoU\nOOqo8p5iw25NSZJimZyViiuvDN2cp51W3kmMlTNJkvIqVnJ2BvAq0DvnWGdgQZGeX/oqKmDWLHjr\nrbCSQDkq56RTkqQCKdZozZnA0DrHLgFuKNLz06FDB5g3L7yH1r8/TJ6cdESFZ+VMkqS8it2tWVs6\n2RXYHfhNzu8uAN4B1gAvA18DS2ht03306AHz58PZZ8OSJUlHU1hWziRJilXs5Ky2bPIvwFXAxcB/\nAQOAjcB04GrgB8AqYAiwocgxJm/gwNDFOW4cvPtu0tEUjlNpSJIUq1hVqeOA7QjJ2DxgXfTzK6AH\n8I/AecBs4A6gP/DHfDesrq7+bruqqoqqqqqCB52osWNh5crwuXgxdO6cdESFYXImSSpzmUyGTCbT\n5Otb8pvyTGAqoSvzcKB2Eq/5wCnAx9H+LOAXwHLgJeDQ6KcS+GkD966paQ1dZDU1MHUqfPQRzJ0b\n1uVMsyefhOnTYeHCpCORJKloKkJhotE5V0t+298CDCZ0UdYmZhOADNnE7DpCJW05oYrXFfgL8DfA\ney0YWzpUVMAtt8DHH4epNsqBlTNJkvIq5sv27QiVtNHR/v6EKTZGAfsBdwMro98tB24FlhIGBbRe\n7drBnDkwZAhUVsLEiUlH1HStodopSVIzFTM5W09IxGo9D2xT55zaxG0B0KsYQaXCjjuGNThHjoR+\n/WDQoKQjajorZ5Ik5ZXyl5hakcGDYcaMMIJz7dqko2kaK2eSJMUyOUuTyZPhuONgwgRYvz7paJrG\nypkkSXmZnKXNtGnQtSucd17SkWw+5zmTJCmWyVnabLEF3H03ZDJw661JR7N57NaUJClW61oaqVx0\n6QKPPALDhsGAATB8eNIRNZ6VM0mS8rJyllb9+sFdd4X30FatSjqaxrFyJklSLJOzNBs5Ei6+GI46\nCj7/POloGsfKmSRJeZmcpd1554WF0qdMKf3KVKnHJ0lSCTA5S7uKCpg5E955B665Julo4lk5kyQp\nLwcElIP27cPC6EOHhiramDFJR1Q/K2eSJMWyclYuevYMa3CeeiqsWJF0NPVznjNJkmKZnJWToUPh\n+uvhyCPhk0+SjqZ+JmeSJOVlclZupkyBUaPCMk8bNyYdzabs1pQkKZbJWTm68Ub46iu4/PKkI/k+\nK2eSJOVlclaO2raFBx6A++6D++9POposK2eSJMVytGa52n57mDcPRoyA/v1h0KCkIwqsnEmSlJeV\ns3I2cCDMmAHjxsHatUlHY+VMkqRGMDkrd5MmwbHHhs9vvkk6GitnkiTFMDlrDaZNg44d4YILko3D\nec4kSYplctYatGkD99wDCxfCrFnJxWG3piRJsRwQ0Fp07QqPPALDh8Oee8K++yYTh5UzSZLysnLW\nmuyxR6icTZgAq1cX//lWziRJilWs5OyXwKvAHTnHOgFPApZSimnMGDjzTBg/PkxUW2xWziRJyqtY\nydlZwL5An5xjlwA3AJZTiu2yy6BXr5CkFbOaZeVMkqRYxUrOdgSWAK9H+7sC/YH/zjnnGOAt4C/A\ny8A6QrVtpyLF2HpUVMDs2fDSS2EetGI/W5IkNahYydluwJ5AP2A74F+Bywldmi8AJxO6OauB2cBg\nYBWwD/B+kWJsXTp1CisIXH01LFpUnGc6lYYkSbGKlZydDywDXiMkaB8DbwCnAt0JXZu3A32j8zoC\n3wJfFym+1qlv3zDFxrHHwttvt/zz7NaUJClWS06lcSYwlZB4HQ58EB1/DDgB2BY4CLiK7KCASmA+\nsBewPN/Nq6urv9uuqqqiqqqqYIG3KoccEt5BO/JIeP75UFFrSVbOJEllLpPJkMlkmnx9sb8pjwO6\nATMICds4oBchgRtNqKxVAuOBQcClDdynpsYqTOHU1MCUKfDFF3D//S2XQN1zDzz6KNx7b8vcX5Kk\nElQRvlcb/eVazHnOOgAnAbdE+08QKmv3RD8dgG+A9cAKQiJ3ZBHja70qKmDmTFi1Cq67ruWfJUmS\nGlTMFQK+BH5Uz/E7c7b3jj5fIVTUVCzt28PcuTBkCFRWwujRhX+G1U5JkmK5QoCydt4Z5swJXZwr\nV7bMM6ycSZKUl8mZNrXffnDttWGAwKefFvbeTqUhSVIskzN932mnwaGHwgknwLffFu6+dmtKkhTL\n5Ez1mz4d1q2DK64o7H2tnEmSlJfJmerXti08+CDcfXf4LAQrZ5IkxTI5U8N22AEefjgskL5sWWHu\naeVMkqS8TM6U3+DBcNNNMG4cfPRR8+5l5UySpFgmZ4p37LEwcSJMngwbNjTvXlbOJEnKy+RMjXPN\nNdCuHVx4YdPv4VQakiTFMjlT47RpE9bGfOIJuPPO+PMlSVKTFHP5JqXdttvCI4/AQQfBHnvA0KGb\nd72VM0mSYlk50+YZMABuuw2OPhpWr968ax0QIElSLCtn2nxHHAGvvQZjx8Kzz0KnTo2/1sqZJEl5\nWTlT01x6KVRWwvHHw8aNjbvGypkkSbFMztQ0FRVw661hcfRLLtm86yRJUoNMztR07drBQw/BY4+F\nRC2OAwIkSYrlO2dqnm7dQnI2bBj06QOjRiUdkSRJqWblTM23224wZw6ceCK8+mrD51k5kyQplsmZ\nCmPYMJg+PYzgXLOm/nMcECBJUiyTMxXO8cfDKaeEqTa++KL+c6ycSZKUl8mZCuuKK2D33eGEE74/\nxYaVM0mSYpmcqbAqKsIKAp9+Cuec8/2EzMqZJEl5mZyp8LbaCubOhcWL4frrs8etnEmSFMupNNQy\nunaFBQtg//2hZ0846aRw3MqZJEl5WTlTy+nZMyRoF10ETz313VQamUwm6cjURLZdutl+6WXbtS6l\nlJztASwDvgJeBlYDq4CxSQalZhowIKwicPzxsHQp4D8yaWbbpZvtl162XetSSslZd+B84DVgMDAf\nOA14NMmgVADDhoXlnW6/3W5NSZJilFJy9gywAVge7Q8kVNBUDsaPD8nZkCFJRyJJUkkrhTLGmcDU\naHsu8FdgOrAS6N/ANX8C/rblQ5MkSWq2N4Hdkg6iqW4H/h7oAfw24VgkSZKKrpS6NQEqCd2anwCd\ngauTDUeSJEmSJElKiR8BK4A3gEsSjkXxZgFrgP/NOdYNeBr4I/AUsE0CcalxdgUWEUZQvwqcEx23\nDUtfe+B3wCvAH4Bro+O2Xbq0IQyMq521wPZLh7cJvYAvA0uiY2Xbdm0IAwH6AG0J/+gMSDIgxTqQ\nMC1KbnL2M+DiaPsS4LpiB6VG6wEMirY7EwbpDMA2TIuO0eeWhHd4h2Hbpc35wK8IU0uB7ZcWfyYk\nY7nKtu32A57M2b80+lFp68OmydkKwpx2EL78VxQ7IDXZPOBQbMO06Qi8COyFbZcmuwC/Bg4mWzmz\n/dLhz8B2dY5tVtuV2oCAfHYG3s3Z/7/omNKlO6Grk+ize55zVTr6EKqgv8M2TIstCD0Ma8h2T9t2\n6TEduAj4NueY7ZcONYTE+vdkpwrbrLZL08LnNUkHoIKrwXZNg87AQ8C5wLo6v7MNS9e3hG7prsBC\nQgUml21XusYAHxLeWapq4Bzbr3QdALwP7EB4z6xulSy27dJUOXuP8IJyrV0J1TOlyxpCSRdgJ8I/\nQCpdbQmJ2V2Ebk2wDdPmU+BxYB9su7TYHziC0D12L3AI4b9B2y8d3o8+1wIPA0PYzLZLU3L2e6Af\noXulHTCZ7EuSSo/5wMnR9slkv/BVeioIE0P/AfhFznHbsPRtT3Y0WAdgBKEKY9ulw08IBYi+wDHA\nb4ATsf3SoCOwdbTdCRhJeO+6rNvuMMKIsT8BlyUci+LdC6wG1hPeF5xCGMHya8pwOHEZGkboGnuF\n8MX+MmE6G9uw9O0NLCW03XLCu0tg26XRQWQLEbZf6etL+O/uFcIURLW5im0nSZIkSZIkSZIkSZIk\nSZIkSZIkSZIkSZIkSVJxbSQ7h9rLQK9kwymovYFZ0fYpwM11fp8hzNrfkAcIcyZJamXStLampPLz\nBWFB9fpURJ9pXT/wIrIJWX1/Q9z6ev8J/BNwToHjklTi0rR8k6Ty14ewCsidhCVPdiUkOUuAZUB1\nzrn/HJ37P8A9wAXR8QzZitT2hPUJAdoAN+Tc6/ToeFV0zYPA68DdOc/4IfAcYbbv3xIWgX8GGJhz\nzmJClSzXVsC+wIsxf28FMJZs5XAl8FbO33F4zPWSJEkFtYFsYvIQ0JvQ1Tkk+v1I4NZoewvgUeBA\nQvK1HGhPWMfuDeD86LxFwA+i7dzk7HRCQgcheXqRkAxWAZ8APQnJ0vOEhafbAW+STfQ6ExK8k4Dp\n0bHdqT8B2zeKtdYphIWOc7tw1+XEWet+4Mc5+88AA+q5v6QyZrempCR9yabdmn2AVYTqFoTkbCQh\nmYGwkHA/QkI2F/gq+plPvJGECteEaL8LsBvwTfS81dHxVwjveq0D3gdeio7/NfqcA/yUUNE7FZhd\nz7N6R9fWqgHuY9MuykV1rrmY0M377znHVhP+N3k95m+TVEZMziSVms/r7F8L/EedY+eSfSeNOtsb\nyL6y0b7OdWcDT9c5VgV8nbO/kfBvY0Pvg30R3eMoYCLfr34RXVtR51jd/VyHAkcDw+u55ts810kq\nQ75zJqmULSRUpzpF+zsDOwDPEpKj2m7NMTnXvA38XbQ9Ief4QuBMsv+ndHegYwPPrSG8/7VTzr22\nJnRrAtwG3ESouH1az/WrgB45+/kSs97AL4FJbJokEj1/VZ5rJZUhK2eSktTQKMZaTxPeuXoh2l8H\nnEDo5ryf8GL/h4T3vmoToJ8TpqE4HXg85363EboIl0bnfgiMo+FRk98AkwkjLjsQKmYjCJW9pYSk\nrL4uTaK4+tf5m+p7RgVwMtANmBcde4+QbLYFdgFWNPAMSZKkknUl2dGaxdCTUFnL5w5gaDOeMRL4\nt2ZcLyml7NaUVC6KNR/aSYRpNX4Sc97PgTOa8Zx/IDsqVJIkSZIkSZIkSZIkSZIkSZIkSZIkSZIk\nSZKktPl/h3fA+OvbGmoAAAAASUVORK5CYII=\n", + "text": [ + "" + ] + } + ], + "prompt_number": 70 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Logarithmic Bode plot" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "freq = pylab.logspace(-3, pylab.log10(nyquist_freq), num_points)\n", + "resp = respfunc(freq)\n", + "\n", + "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 10), dpi=100)\n", + "\n", + "axes[0].plot(freq, 10 * pylab.log10(abs(resp)), 'r')\n", + "axes[0].set_xlabel('Frequency (Hz)')\n", + "axes[0].set_xscale('log')\n", + "axes[0].set_ylabel('Magnitude (dB)')\n", + "axes[0].axhline(y=-3)\n", + "for f in halfgain:\n", + " axes[0].axvline(f)\n", + "\n", + "axes[1].plot(freq, pylab.angle(resp), 'r')\n", + "axes[1].set_xlabel('Frequency (Hz)')\n", + "axes[1].set_xscale('log')\n", + "axes[1].set_ylabel('Phase (radians)')\n", + "axes[1].set_ylim([-pylab.pi, pylab.pi])\n", + "axes[1].set_yticks(pylab.linspace(-pylab.pi, pylab.pi, 9))\n", + "axes[1].set_yticklabels([\"-\u03c0\", \"-\u00be\u03c0\", \"-\u00bd\u03c0\", \"-\u00bc\u03c0\", \"0\", \"\u00bc\u03c0\", \"\u00bd\u03c0\", \"\u00be\u03c0\", \"\u03c0\"])\n", + "\n", + "None" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAJjCAYAAACIpPyrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm8XPP5wPHPJJF9s8SSRa4lQRB70FJTKmKrfW3VVvtO\nS1U1V1dF+aFoFaVqL7FUa2nr2ok1iSWyEJVQxJoSss3vj++MO7m592bubGfOnM/79TqvOXNm5pwn\nOXdynzzfDSRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiQpNsYAk4GpwBkRxyJJkqQCdQamAQ3AMsCL\nwDpRBiRJkhSVTlEH0EGjCIncDGA+cDOwW5QBSZIkRSVuidwg4K285zOzxyRJkhInbolcJuoAJEmS\nakWXqAPooFnAkLznQwhVua+ssdxymekffljVoCRJkoo0HViz2A/HrSL3LDCMMNihK7AfcHf+G6Z/\n+CGZTCbSbezYsTVxvo58rpD3tveejr5W6DHvX+XvXXuvd+R4Pd6/evzu1fL9i813b9EiMnPnkpk9\nm7Enn0zmlVfIPPYYmbvuInPNNWTOP5/Mj35E5ogjGDtsGJmRI8n070+mRw8ya61FZscdyZx2Gpmr\nr2bs4YeT+eSTku4dRH/vyn3/4vDd6+h9avv+sUYpiVHcKnILgOOB+wkjWK8GXo00olak0+maOF9H\nPlfIe9t7T0dfK/ffUTnF7f6Veu/ae72jx2tBOWOrx+9eodeMQmy+e6kUdO8O3buT3m03WKftyRPS\nTU2QO9+cOfDWWzB1Krz6Kjz8MOmnnoKBA2HoUNhii7BtuSXpbbYpLrYIJe27197r1fzupcp+xuhl\nshmuYqaxsZHGxsaow1CRvH/x5v2L0Pz58NJL8NRTYXvsMZg7F3bYAcaMgdGjYdll2/x4KtVIJtNY\nvXhVVqlUCkrIx+LWtKo6Vsv/09TSef/izfsXoWWWgY02gmOOgeuug+nT4ZFHYJNN4PrroaEBdtsN\nbrkFPv+8lROkqxywaokVOUmSatmnn8K4cXDTTaFit88+cPzxsMEGQGjp9ddefJVakTORkyQpLv77\nX7jqKvj972H11eHUU0ntsbuJXIyZyC3JRE6SVN/mz4c774Rf/5rUC8+Tue2vsOee0MkeU3FjIrek\njPMGS5KkeCgtkYvb9CMFsSAnSUqKVAoyizJw++1w+umw/vpw/vkwfHjUoakAqRJLatZgJUmKu1QK\n9t4bXnkFttoKvv51+NnPYN68qCNThdVl06p95CRJSdHqqNW33oJjj4XXXw+DI7bcMpLYtHT2kVuS\niZwkKTHanH4kk4Fbb4WTT4bDD4exY8OcdaopJnJLMpGTJCXGUueR++9/4dBD4aOP4IYbYI2SlvZU\nmbmygyRJatvKK8O998IBB4S1XG+4IeqIVEZW5CRJirEOrewwcSLstRfstBNccIFNrTXAipwkSSrM\nyJHwzDMwbRp861vw7rtRR6QSmchJkpQk/fvDPfdAOg2bbQbPPht1RCqBTauSJMVYh5pWWxo3Do48\nEq65BnbdtaxxqTClNq3W5coOkiSpAHvsAYMGwe67w4wZcMIJUUekDrIiJ0lSjJVUkct54w3YeWfY\nYQf47W+hkz2vqsV55JZkIidJSoyyJHIQ5pnbbTcYOjQ0tTqitSpM5JZkIidJSoyyJXIAn38e1mxd\nZhm45Rbo3r1MJ1ZbnH5EkiSVR8+ecOedIYHbeWeYMyfqiLQUJnKSJKlZ165w442w+uqw/fahyVU1\ny0ROkiQtrnNnuPLKsKTXDjvAxx9HHZHaYB85SZJirKx95FrKZOCkk2D8eHjgAejbt0IXSi4HOyzJ\nRE6SlBgVTeQgnPz44+GFF+D++6FPnwpeLHkc7CBJkionlYJLLw3rtO64owMgaowVOUmSYqziFbmc\nRYvCcl5vvAF//zt061aFi9Y/m1aXZCInSUqMqiVyAAsXwn77hf1bbgmDIlQSm1YlSVJ1dO4MN9wQ\nRrEec0wVM0i1xUROkiQVrls3GDcOXnwRzjor6mgSz0ROkiR1TJ8+oZ/cuHFw0UVRR5NoXaIOQJIk\nxdAKK4TpSL72NWhogD32iDqiRHKwgyRJMVbVwQ6tee45GDMG7r0XRo2KMJB4crCDJEmKziabwNVX\nw+67w4wZUUeTODatSpKk0nz72yGJ23lnePxx6N8/6ogSw6ZVSZJiLPKm1XwnngiTJ4eBEF2sFRXC\nplVJklQbciNYzzwz2jgSxEROkiSVR+fOcPPNcMcdcOONUUeTCDatSpIUYzXVtJozcSJstx08+CBs\nuGHU0dQ0m1YlSVJtGTkSLrsszC03e3bU0dQ1K3KSJMVYTVbkcs44A559Fh54IDS7aglW5CRJUm36\n1a/C489/Hm0cdcyKnCRJMVbTFTmAd94Jkwb/5S+w7bZRR1Nz6rEidz7wKjABuAPol/famcBUYDIw\nuvqhSZKkDlllFbjuOjjoIHj33aijqTu1mMg9AKwLbABMISRvACOA/bKPY4DLqc34JUlSvu23h0MP\nDcncokVRR1NXajERehDI3eWngcHZ/d2Am4D5wAxgGuDqvJIkxUFjI3zxBfzmN1FHUldqMZHLdxjw\n9+z+QGBm3mszgUFVj0iSJHVcly5www1h9YcXXog6mroRVSL3IDCplW3XvPecBcwD2psaupa7d0qS\npHxDhsCFF8L3vheqcypZVCvabr+U1w8BdgK2yzs2CxiS93xw9tgSGhsbv9pPp9Ok0+kiQpQkSWX3\nne/AuHHw05/CeedFHU3VNTU10dTUVLbz1eL0I2OA3wLbAPnTQY8gVOdGEZpU/wmsyZJVOacfkSQl\nRs1PP9Ka99+HDTaAW26BrbeOOppIlTr9SC0mclOBrsCH2edPAsdm939M6De3ADgJuL+Vz5vISZIS\nI5aJHMDdd8PJJ8OECdCnT9TRRKYeE7lSmchJkhIjtokchClJ+vSBSy6JOpLImMgtyUROkpQYsU7k\nPvgA1lsP7roLRiVzRrF6XNlBkiQlwfLLwwUXwBFHwPz5UUcTSyZykiQpOgceCCuvHOaXU4fZtCpJ\nUozFumk1Z/p02HxzGD8eVl896miqyj5ySzKRkyQlRl0kcgC//jU89VToL5cg9pGTJEnxd+qp8PLL\n8MADUUcSKyZykiQpet26wW9/G+aWc+BDwUzkJElSbfj2t2HQILjiiqgjiQ37yEmSFGN100cu5+WX\nIZ2GyZPD9CR1zsEOSzKRkyQlRt0lcgBHHQX9+sF550UdScWZyC3JRE6SlBh1mcjNmgUjR8LEiaGp\ntY6ZyC3JRE6SlBh1mcgBnH46fPIJ/OEPUUdSUSZySzKRkyQlRt0mch9+CMOHw5NPwrBhUUdTMc4j\nJ0mS6s9yy8Epp8DYsVFHUtOsyEmSFGN1W5EDmDMnLNn1xBN1W5WzIidJkupTnz5w7LGJGL1aLCty\nkiTFWF1X5AA++CBU4yZOhMGDo46m7KzISZKk+rX88nDYYWH5Li3BipwkSTFW9xU5gLffhvXWg9de\ngwEDoo6mrKzISZKk+jZwIOy9t2uwtsKKnCRJMZaIihzASy/B6NEwYwZ07Rp1NGVjRU6SJNW/9daD\ntdeG22+POpKaYiInSZLi4cQT4dJLo46ippjISZKkeNh11zDw4Zlnoo6kZpjISZKkeOjcGY4+Gq68\nMupIaoaDHSRJirHEDHbIeecdGDECZs6EXr2ijqZkDnaQJEnJscoq8PWvO+ghy0ROkiTFy2GHwTXX\nRB1FTbBpVZKkGEtc0yrAvHlh3dUnnoA114w6mpLYtCpJkpKla1f4znfguuuijiRyJnKSJCl+DjwQ\nbr01geXIxZnISZKk+Nl009DEOmlS1JFEykROkiTFTyoF++wTqnIJZiInSZLiad99E9+82qWA9/QH\ntgQagAwwA3gS+KRiUUmSJC3NJpvA/PkwYQJsuGHU0USivYrc1sDdwCPA/sCqhGTuAODR7GtbVTg+\nSZKk1qVSoSp3221RRxKZ9uYtuRC4ApjaxuvDgaOBU8sdVImcR06SlBiJnEcu35NPwpFHxnbQQ6nz\nyDkhsCRJMZb4RG7hQlh5ZXjuOVh11aij6bBKTgicAtLAyOzz/YDLgFOAbsVeUJIkqWw6d4YxY+De\ne6OOJBLtZYCXA+sD3YHXgN7AfYR+cSngOxWPrjhW5CRJiZH4ihzAzTfDX/4Cf/tb1JF0WCWbVl8F\nRhASuVnAisCC7GcmAesVe9EKM5GTJCWGiRzw0UcwdCi8+y706BF1NB1SyabVLwjTjcwF3iQkcWSP\nzS/2gpIkSWW17LKwwQbwyCNRR1J17SVyAwgjUk9rsZ97XmmnAYuA5fKOnUkYRTsZGF2FGCRJUhxs\nuy089FDUUVRde4ncVUAfQt+4/P0+wB8rHNcQYHtCJTBnBGHAxQhgDKEPnytTSJIk+OY3E5nI1er0\nI7cBPwfuAjYBPiRU4xYBv8m+5z6gEXiqxWftIydJSgz7yGV98QUMGAAzZ0K/flFHU7BS+8i1t0TX\npXn7mbyL5H5cTiz2okuxGzATmNji+EAWT9pmAoMqFIMkSYqT7t1h1Ch49FHYZZeoo6ma9hK557KP\nXyM0Z95CSOb2AV4u8boPAiu3cvwsQuUtv/9be1lqq/8HaWxs/Go/nU6TTqc7HKAkSYqZXPNqDSdy\nTU1NNDU1le18hZTynibMHZcbqboM8BiwedmiaLYe8C/g8+zzwYSpTzYHDs0eOzf7eB8wNhtfPptW\nJUmJYdNqnscfhxNPDKs8xEQ1luh6jVCV+yD7fDngSWCtYi/aAW/Q3EduBHAjMIrQpPpPYE2WrMqZ\nyEmSEsNELs/cubDCCvD++9CzZ9TRFKSSfeRyzgWeBx7KXmgbwiCDasj/0XwFuDX7uAA4ljaaViVJ\nUgL16AHrrgvPPw9bbRV1NFVRaAa4CqF5MwOMB96pWESlsyInSUoMK3ItnHgirLoq/OAHUUdSkEqu\n7LBG3v47wJ2E6UDeaeM9kiRJ0dp8c3iq5cxk9au9DPAWoBdwN/AsIYFLEapzmwLfBuYA+1c4xo6y\nIidJSgwrci1Mnw7pNLz1VtSRFKTSgx3WJCRqXweGZo+9SRi1ehPwerEXriATOUlSYpjItZDJwIor\nwgsvwODBUUezVJUe7DAN+EWxJ5ckSaqqVCo0rz79dCwSuVK5VqkkSaovG28ML74YdRRVYSInSZLq\ny4YbmshJkiTFkoncEu85CPhp9vmqhNUVJEmSak9DA3z6KcyeHXUkFVdIInc5sCVwYPb5/7LHJEmS\nak+nTjByJEyYEHUkFVdIIrc5YTmsudnnHwLLVCwiSZKkUm2wAUycGHUUFVdIIjcP6Jz3fACwqDLh\nSJIklcE668Crr0YdRcUVkshdCowDVgR+BTwO/LqSQUmSJJVkxIhEJHKFziS8DrBddv9fQC3/zbiy\ngyQpMVzZoQ3//S+su24Y8JAqeuGEiqvkEl3LtfHe3I/Lh8VetMJM5CRJiWEi14ZMBpZfHiZPDkt2\n1ahKLtH1PCFpSxGmHPkoe3xZwnqrqxV7UUmSpIpKpZr7ydVwIleq9vrINRCStQeBXYDls9vO2WOS\nJEm1KwEDHgoZ7LAl8Pe85/8AvlaZcCRJkspk+HCYNi3qKCqqkETubeAnNFfozgJmVTAmSZKk0q2+\nOkyfHnUUFVVIIncAYeqRccAd2f0DKhmUJElSydZYo+4Tudodj1s8R61KkhLDUavt+PRTGDgQ5syp\n2SlIKjlqNeehVo5lgG2LvagkSVLF9e0LPXrAe+/BSitFHU1FFJLI/TBvvzuwF7CgMuFIkiSVUa6f\nXIITuWdbPH8MeKYCsUiSJJVXrp/c1+pzwo1CErn8FR46AZsCfSsTjiRJUhmtvjq8/nrUUVRMIYlc\nboUHCE2qM4DDKxWQJElS2ay6KjxTvw2JhSRyawNftDjWvQKxSJIkldfgwXDHHVFHUTGFzCP3RIHH\nJEmSasvgwTBzZtRRVEx7FblVgIFAT2BjwhwnGUL/uJ6VD02SJKlEQ4YkNpHbATgYGAT8Nu/4HODH\nlQxKkiSpLPr3h/nzw6TAffpEHU3ZFTKT8F7A7ZUOpIxc2UGSlBiu7FCAtdaCO++EddaJOpIlVHJl\nh4OA64EG4NT8axKaWC8s9qKSJElVk+snV4OJXKnaS+Ry/eD60Dz9CDQncpIkSbVvyBB4662oo6iI\n9hK5P2QfG6sQhyRJUmUMHpzIRC5nReAIQhNr7v0Z4LAKxSRJklQ+q6wCL70UdRQVUUgidxfwCPAg\nsCh7zKZVSZIUDyuuCO+/H3UUFVFIItcDOKPSgUiSJFXEiivCe+9FHUVFFLKyw9+AnSsdiCRJUkXU\ncSJXyLwl/yOMYJ0HzM8ey63wUIucR06SlBjOI1eADz6AYcPgww+jjmQJpc4jV/QHa5iJnCQpMUzk\nCrBoEXTrBp99Bl27Rh3NYio5IXDOxq0c+wR4E1hQ7IUlSZKqolMnWH55mD0bBg6MOpqyKiSRuxzY\nBJiYfb4+8DLQDzgGuL8yoUmSJJVJbuRqnSVyhQx2eBvYkJDMbZLdfx3YHjivQnGdALwKvAT8Ju/4\nmcBUYDIwukLXliRJ9aZOBzwUUpFbi1CBy3kFWBuYTmXmk/sm8G1gJGFwxYDs8RHAftnHQcA/geE0\nz20nSZLUugQnci8DVwA3Ezrj7UtI5rrRPIq1nI4Bfp137twMfrsBN2WPzwCmAaOApyoQgyRJqid1\nmsgV0rR6CKH6djJwEqFZ9WBCQrVtBWIaBnyDkKA1AZtmjw8EZua9byahMidJktS+Ok3kCqnIfQ5c\nkN1amlPkdR8EVm7l+FnZmJYFtgA2A24FVm/jPA64liRJSzdgALz+etRRlF0hidxw4FeEvmk9sscy\ntJ1cFWL7dl47Brgju/8MoQ/cCsAsYEje+wZnjy2hsbHxq/10Ok06nS4+UkmSFH/LLgsffRR1FDQ1\nNdHU1FS28xUyAd3jwFjgQsIghEOAzsDZZYticUcRmlHHEpLIfwKrEhLJGwn94nKDHdZkyaqcEwJL\nkhLDCYEL9K9/wS9/Cf/+d9SRLKbUCYEL6SPXg5A0pQiDDBqp7Nqr1xCqfZMIgxu+lz3+CqGZ9RXg\nH8Cx2LQqSZIKUSMVuXIrpGn1C0IFbhpwPGFeuV4VjGk+cFAbr/0qu0mSJBVu2WVrcq3VUhWSyJ0M\n9AROBH4O9CWMWpUkSYqHOq3IFd0mW8PsIydJSgz7yBVo0SLo2hW++AK6FFLHqo5S+8i19ye5h9AH\nrbWTZwgDHyRJkmpfp07Qty98/DGssELU0ZRNe4ncFoRJd28Cns4eyyV15v6SJClecs2rCUnkViHM\n93ZAdruXkNS93M5nJEmSalPfvjCn2LUMalN7048sIEzz8T1CdW4a8DBh5KokSVK89O0Ln34adRRl\ntbTeft0Jc8btDzQAFwPjKhyTJElS+SUskbseWBf4O/AzwgS9kiRJ8VSHTavtJXLfAT4DTspu+TKE\n+eQkSZLiIWEVuUKW75IkSYqHOkzkTNYkSVIymMhJkiTFlImcJElSTPXtC598EnUUZWUiJ0mSkqF3\nb/jss6ijKCsTOUmSlAw9e5rISZIkxVKvXiZykiRJsdSrF3z+edRRlJWJnCRJSgYrcpIkSTFlHzlJ\nkqSYsiInSZIUUyZykiRJMdWtGyxYELY6YSInSZKSIZUK/eTqaOSqiZwkSUqOOmteNZGTJEnJYSIn\nSZIUUzatSpIkxVT37vDll1FHUTYmcpIkKTm6dYMvvog6irIxkZMkScnRvbuJnCRJUix162bTqiRJ\nUixZkZMkSYopBztIkiTFlIMdJEmSYsqKnCRJUkxZkZMkSYopBztIkiTFlE2rkiRJMWXTqiRJUkxZ\nkZMkSYopK3KSJEkx5WAHSZKkmLJpteJGAeOBF4BngM3yXjsTmApMBkZXPzRJkhRrdda02iXqAFpx\nHnA2cD+wY/b5N4ERwH7Zx0HAP4HhwKJowpQkSbHTtSvMnx91FGVTixW5d4B+2f3+wKzs/m7ATcB8\nYAYwjVC9kyRJKkyXLnWVyNViRe5HwGPABYREc8vs8YHAU3nvm0mozEmSJBVmmWVgwYKooyibqBK5\nB4GVWzl+FnBidhsH7ANcA2zfxnkyFYlOkiTVJytyZdFWYgbwF+Bb2f2/Aldl92cBQ/LeN5jmZtfF\nNDY2frWfTqdJp9NFhilJkupKxBW5pqYmmpqayna+VNnOVD7PA6cADwPbAecSRq6OAG4k9IvLDXZY\nkyWrcplMxkKdJCkZUinw114HjB8Pxx0HzzwTdSQApFIpKCEfq8U+ckcClwHdgLnZ5wCvALdmHxcA\nx2LTqiRJ6og66yNXixW5UlmRkyQlhhW5Dpo0CQ44AF56KepIgNIrcrU4/YgkSVJldOlSVxU5EzlJ\nkpQcdda0aiInSZKSo86mHzGRkyRJyWFFTpIkKaasyEmSJMWUFTlJkqSYsiInSZIUU1bkJEmSYsqK\nnCRJUkx16QILF9bNchgmcpIkKTlSKejcOSRzdcBETpIkJUsdNa+ayEmSpGSpowEPJnKSJClZrMhJ\nkiTFlBU5SZKkmLIiJ0mSFFNW5CRJkmLKipwkSVJMWZGTJEmKqS5dTOQkSZJiyaZVSZKkmLJpVZIk\nKaa22gp69446irJIRR1ABWQymUzUMUiSVBWpFPhrL75SqRSUkI9ZkZMkSYopEzlJkqSYMpGTJEmK\nKRM5SZKkmDKRkyRJiikTOUmSpJgykZMkSYopEzlJkqSYMpGTJEmKKRM5SZKkmDKRkyRJiikTOUmS\npJgykZMkSYopEzlJkqSYMpGTJEmKKRM5SZKkmDKRkyRJiikTOUmSpJiKKpHbB3gZWAhs3OK1M4Gp\nwGRgdN7xTYBJ2dcurkKMkiRJNS2qRG4SsAfwSIvjI4D9so9jgMuBVPa1K4DDgWHZbUxVIlXVNDU1\nRR2CSuD9izfvX5w1RR2AIhRVIjcZmNLK8d2Am4D5wAxgGrA5sArQBxiffd+fgd0rHqWqyl8k8eb9\nizfvX5w1RR2AIlRrfeQGAjPzns8EBrVyfFb2eE0q9z+IxZ6vI58r5L3tvaejr9XyL4243b9S7117\nr3f0eC0oZ2z1+N0r9JpRiNt3r5D3+t2r7rmq+d1r7/Vq3r9KJnIPEppQW267VvCaNaEe/zFa2ntM\n5Mp/PhO54iTtl4mJXPnPZyJXnKR999p7vZr3L7X0t1TUQ8BpwPPZ5z/KPp6bfbwPGAu8mX3vOtnj\nBwDbAEe3cs5pwBqVCFaSJKnMpgNrRh1EsR4ijEbNGQG8CHQFViP84XLJ5tOE/nIp4O842EGSJCkS\newBvAXOB/wL/yHvtx4Sq2mRgh7zjuelHpgGXVCdMSZIkSZIkSZIkSZIkSZJq3tqElSFuJawOoXjZ\nDbgSuBnYPuJY1HGrAVcBt0UdiArWC7iO8L07MOJY1HF+5+LN33nt6ERI5hRP/Qn/OCme/KUSHwcB\nO2f3b44yEJXE71y8FfQ7r9ZWdliaa4B3CaNX840hjHKdCpzRxmd3Be7Ff5SiVMr9A/gJ8LvKhKYC\nlHr/FK2O3L9BhJkFABZWJTotjd+/eCvm/tXl77ytgY1Y/C+iM2FKkgZgGcI8dOsQ/kd5EWF5r3x3\nVTxKtaXY+5cCfgNsV8VYtaRSv39WB6LVkfv3XZorcjdVL0S1oyP3L8fvXO3oyP2r+995DSz+F7El\nYQWInB/RvEJEzjbAxcAfgJMrGZyWqoGO378TgWcJ/RyPqmRwWqoGOn7/lgN+jxWDWtBAYfevJ6GC\ncDlhJR3VhgYKu39+52pTA4XdvxPowO+8LuWLLzL5TQAAMwkrQOR7OLup9hRy/y7BSaBrVSH370Na\nX05P0Wvr/n0OHBZJROqItu6f37l4aOv+nQBcWuhJ4tZHrjWZqANQSbx/8eb9izfvX7x5/+KtLPev\nHhK5WcCQvOdDCFmt4sH7F2/ev3jz/sWb9y/eEnv/Gli8jbkLMD17vCtLdvZUbWnA+xdnDXj/4qwB\n71+cNeD9i7MGvH/cBLwNfEloVz40e3xH4DXC6I8zowlNBfD+xZv3L968f/Hm/Ys3758kSZIkSZIk\nSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZJq2kLghbxt1WjDKav1gWuy+4ew5ALZTcAm7Xz+VmC1skcl\nKXa6RB2AJLXhc2CjNl5LZR/jumj4D2lO3lr7M2TaOJ7zR+AU4MQyxyUpZjpFHYAkFaiBsJTNdYT1\nCocQEqLxwASgMe+9Z2Xf+yhwI3Ba9ngTzZWuFYA3svudgfPzznVk9ng6+5nbgFeBv+RdYzPgccL6\niE8BvYGHgQ3y3vMYofqWrxuwBfDMUv68KWBXmiuSrwGv5/05dlrK5yVJkiKzgOYk5nZgKKG5dVT2\n9dHAH7L7nYB7gK0JidpEoDvQB5gKnJp930PAxtn9/ETuSELyByHReoaQOKaBj4GBhMTqCeBrhAWu\np9OcFPYmJIPfAy7KHhtO68naFtlYcw4B3mPxZuQ5eXHm3AIck/f8YRKwwLak9tm0KqlWzWXxptUG\n4E1C1QxCIjeakPgA9AKGEZK3O4AvstvdBVxrNKFytnf2eV9gTWB+9npvZ4+/SOibNgd4B3gue/x/\n2ce/AmcTKoWHAX9q5VpDs5/NyQA3s3gz6UMtPnM6oan5irxjbxP+Tl5dyp9NUh0zkZMUJ5+1eP5r\n4MoWx06iuQ8dLfYX0NylpHuLzx0PPNjiWBr4Mu/5QsK/m231X/s8e47dgX1YsqpG9rOpFsdaPs/3\nLWAv4ButfGZRO5+TlAD2kZMUV/cTql69ss8HAQOARwiJVK5pdZe8z8wANs3u7513/H7gWJr/czsc\n6NnGdTOE/mqr5J2rD6FpFeAq4BJCJe+TVj7/JrBy3vP2krihwGXAviyeUJK9/pvtfFZSAliRk1Sr\n2hrNmfMgoY/Yk9nnc4DvEppabyEMWniP0E8tlyxdQJi640jg3rzzXUVopnw++973gD1oe/TofGA/\nwsjTHoRK3PaEiuHzhASutWZVsnGt1eLP1No1UsDBwHLAndljswiJ6TLAYGByG9eQJEmqC2NpHrVa\nDQMJFbt5/WXDAAAgAElEQVT2XAtsXsI1RgMXl/B5SXXCplVJSVCt+ea+R5iK5MdLed8FwNElXOf7\nNI+OlSRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkqQa\n9VfgXWA68BIwB/hjpBFJkiRFpHPUAXRAD6Av0A04GvgH0A84NMqgJEmSopKKOoAivAasBxwO9AYu\nyH9xjTXWyEyfPj2KuCRJkjpqOrBmsR/uVMZAqqEf8AUwn5DMvdDyDdOnTyeTyUS6jR07tibO15HP\nFfLe9t7T0dcKPeb9q/y9a+/1jhyvx/tXj9+9Wr5/cfvuFfLeqnz3unYl89BDdXX/auHeVev+AWuU\nkhjFqWkVYBQwBLgDOAG4DZjd4j2NjY2NVQ5rSQ0NDTVxvo58rpD3tveejr7W8lhTUxPpdHqpMVRD\n3O5fqfeuvdcLPV6v968ev3utHa+V+xe3714h7634d+/nPyf9+eew997tXqcakvbda+/1Qo+fc845\nAOcsNZA2xK1p9TigJ3A+8DNgX2B9QoUuJ5PNcBUzjY2N1EISruJ4/+LN+xdfjV270ti1K7zxBgwY\nEHU46qBUKgUl5GNxa1q9jJDEAfwUWJvFkzjFWC1UA1Q871+8ef/iK92pE4wZA9ddF3UoikDcKnKF\nsCInSUqOXr3gzjvhuOPgtdcgVY+/2utX0ipykiQpXyYDW24J3bpBU1PU0ajKTOQkSYq7VAqOPBL+\n8IeoI1GV1WP91aZVSVJy9OgBH3wA8+ZBQwNMneqghxixaVWSJEH//rDHHnDttVFHoioykZMkKc4y\nmeYBDkceCVdeCYsWRRuTqsZETpKkerHFFqGp1UEPiWEiJ0lSnOVX5Bz0kDgOdpAkKc66dYNPPoHu\n3cPzjz8Ogx6mTIEVV4w0NC2dgx0kSUqy/IochEEPe+7poIeEMJGTJCnuWq7m4KCHxDCRkyQpzlrr\nTrT55mHproceqn48qioTOUmS4q5lRc5BD4nhYAdJkuKsSxeYOxeWWWbx4598AkOHwmuvwUorRROb\nlsrBDpIkJVnLwQ45/fo56CEBrMhJkhRnnTvDl1+GylxLTz8NBx4Y1l/tZO2mFlmRkyQpydqqyAGM\nGgV9+sC//13dmFQ1JnKSJNUrBz3UPZtWJUmKs1QKFi5su+k0N+hh8mRYeeXqxqalsmlVkiS1rV8/\n2GsvBz3UKStykiTFWSoVVnBoq58cwPjxsP/+MG2agx5qjBU5SZLUvs02g2WXhfvuizoSlVmtJHI7\nAeOBSUBD3vGLgM2jCEiSpJqXa4FqrxqXe/344+HSSysfk6qqVhK5vxMStluB9bLH1gL6A09HFZQk\nSXVj//3huedgypSoI1EZ1UoiB/AHYB9CZQ7gZ8BPsvsrAM8AnwETgBnAO8BR1Q1RkqQa0pE+4T16\nwOGHw+WXVy4eVV2tDHYYBTwHHAisCEwEvg40AscAvwD2A64EVgd+BMwGrmrlXA52kCQlw6JFYUWH\nRYsKe/9//gMbbQRvvgm9e1c2NhWkXgY7rAU8DxwK3AWcBvyGUIk7nFCFmwL8J/v+DYAXqh+mJEk1\npKOFi1VXhXQa/vznioSj6mtlYbZIXJ/dAI4GbgDmEipx3wMuIyRvE7PvGUEYGNGqxsbGr/bT6TTp\ndLrc8UqSVBuWNtChpRNOgGOPhWOO6fhnVbKmpiaamprKdr5au4P9gBuBnQnVwr8CbwJ7A38EZgFX\nE/rINbRxDptWJUnJsGABdOsWVnYoVCYDI0fCRRfBt75VudhUkHppWs05C/hVdn8RsCdwCjCNMJo1\nV5F7A7AuLElKtkym41W1VCpU5X73u8rEpKqqtYpcOViRkyQlw/z5YTTqggUd+9xnn4X1V599Fhoa\nKhKaClNvFTlJklSoYipyAL16wcEHOxVJHbAiJ0lSXM2bF5Ky+fM7/tk33ghLd73xBvTpU/7YVBAr\ncpIkJVWxFTmA1VaDbbeFa64pb0yqKitykiTF1ZdfhmravHnFff7pp2G//WDatDCxsKrOipwkSUlV\nSkUOYPPNwyTBt99evphUVSZykiQl2Q9+AOef3/FVIlQTTOQkSYqrUityALvsAnPmwCOPlCcmVZWJ\nnCRJSdapE5x2WqjKKXYc7CBJUlx9/jksvzzMnVvaeebODRMDP/QQjBhRltBUGAc7SJKUZOVY+L5H\nDzjuOLjwwtLPpaqyIidJUlx99hkMGBAqc6WaPRuGD4dJk2DQoNLPp4JYkZMkKcnKUZEDWGEFOOQQ\nuOCC8pxPVWFFTpKkuPrf/2CllUJlrhzefhvWWw8mT4YVVyzPOdUuK3KSJCVVOaYfyTdwIBxwgH3l\nYsSKnCRJcTVnDqyySqjMlcubb8LGG8PUqbDccuU7r1plRU6SpKQqd0UOYOhQ2H13uOSS8p5XFWFF\nTpKkuPr00zDCdM6c8p532jTYckuYPh369i3vubUYK3KSJCVVJSpyAGuuCTvsAL/7XfnPrbKyIidJ\nUlx98gkMGRIqc+U2eTJsvXXoK9e/f/nPL8CKnCRJyVWpihzA2mvDzjs7grXGWZGTJCmuPv44DE74\n5JPKnP+NN2DTTUN1bsCAylwj4azISZKUVJWsyAGsthrsvz+ce27lrqGSWJGTJCmuPvooJFsff1y5\na7zzTljtYeJE12CtgHqpyG0E3AaMbXH8EmDT6ocjSVIMVLoiB2HC4cMPh1/8orLXUVFqJZF7AWg5\nxnltoBfwbPXDkSRJXznjDLj11jCvnGpKrSRyrTkH+Gl2f21gAvAFIel7G3gT2DWa0CRJqgHVqMgB\nLL88nHIK/PjHlb+WOqSWErn8n8TtgMnABsAVwNbAqcDLhGbYu4HDgXuqHKMkSbWlGokcwKmnwhNP\nhE01o1YSuTWAg4B1gR2AHwC/AT4FjgEOAeYDk7Lv34BQmZMkKbmqObivZ0/45S9DQuegwprRJeoA\nsqYTKmwA3ycMfPgceAzYCXiIkLxNyL5nOeCDtk7W2Nj41X46nSadTpc7XkmSoletptWc734XLr4Y\nbrklTEuiDmtqaqKpqals56u16Ud6E5K4HbPPv0aoxh0N/BG4kdC8eiewRRvncPoRSVIyvP8+rLMO\nzJ5dvWs+/DAcfHCYJLh79+pdt07Vy/QjOacD5+c9/xOwOfAcoSI3EfiYkPD9surRSZJUS6pdkQPY\nZhvYaKNQmVPkaq0iVw5W5CRJyfDee7DuuqEyV01Tp8KWW8KkSWGeORWt3ipykiSpUFFU5ACGDQuT\nBJ9+evWvrcWYyEmSpI47++zQX+7hh6OOJNFM5CRJiquoKnIAvXvDRRfBccfB/PnRxCATOUmSVKQ9\n94RBg+CSS6KOJLFM5CRJiqsoK3IQrv2738Gvfw2zZkUXR4KZyEmSpOINGwZHHw0nnxx1JIlkIidJ\nUlxFXZHLOessmDAB7rwz6kgSx0ROkiSVpkcPuOqqMPDho4+ijiZRTOQkSYqrWqnIAXzjG7D77vDD\nH0YdSaKYyEmSFGe1kshBGPTw4IPwz39GHUlimMhJkhRXtbYkZd++cMUVcOSR8NlnUUeTCCZykiTF\nVS01rebstBNstZVNrFViIidJksrrkkvg3nvhH/+IOpK6ZyInSVJc1WJFDqB/f7j2Wvj+92H27Kij\nqWsmcpIkqfy++U3Yf3846qja68tXR0zkJEmKq1qtyOX88pcwZQpcf33UkdQtEzlJklQZ3bvDX/4C\np50G06dHHU1dMpGTJCmuar0iB7DBBvDTn8K++8KXX0YdTd0xkZMkSZV1/PHQ0AA/+EHUkdQdEzlJ\nkuIqDhU5CDFefXWYkuT226OOpq6YyEmSpMrr3x9uvhmOOQZefz3qaOpGORK5XsDawFrZfUmSVA1x\nqcjljBoFP/4x7Lef/eXKpNhErg9wKjAemAT8CbgOeAl4FjgF6F2OACVJUh056SQYMiQ8qmTFJnJ3\nAnOAXYHVgS2BLYDVgF2Az4C7OnC+9YE/A2NbHD8/e15JktRS3CpyEOK99lp4+GH44x+jjib2uhT5\nue3aee2/wJXZrVCTgKuBdN6xNYABwFMdDU6SpMSIWyIH0Lcv3HknbL01rL8+bGHNplil9pHbiuYm\n1IOAi4ChJZ4z52fAT7P7KwDPECp9E4AZwDvAUWW6liRJ8RPnpa/WWiuMZN17b3jnnaijia1SE7kr\nCMnVBoQ+c9MITaTFyP8vxTaEZO0/wHHAa8CZwLvZa/0eOBv4Q5HXkiQp/uLYtJpv113hiCNgn31g\n3ryoo4mlYptWcxYAGWB34DLgKuDwIs6zBqGi1wfYETgR2AdYhZDITQCmEBI7CMncBW2drLGx8av9\ndDpNOp0uIiRJklRxZ58Nzz8Pxx4b+szFOTEtQFNTE01NTWU7X6l/W48A9wGHAlsD7wMvEgYvFOtQ\nYBlCH7tLCQniFcCFwPaEJG8CsBnQWvqeycS51CxJUqGmTIGdd4apU6OOpDRz5oT+ct/9buJWf0iF\nxLXofKzUitx+wIHAYYRBDqvSTqWsAL2A/YEx2edjCP3w1iRU4SZkj/ej9SROkiTFTZ8+cM89sOWW\nsOaasPvuUUcUG6X2kXsH+C3waPb5fwjzyRXrjOz5ciW1YcBGhLnpRgITs8ffoPi+eJIk1Ye495HL\nN2RIGMl6xBGhqVUFKfXu7wWcC6yUd64M0LfE85bCplVJUjK89loYMDBlStSRlM8dd4TJgp98EgYP\njjqaiou6afU8wgTAr5Z4HkmS1FH1VJHL2XPP0Odv113h0UehtwtFtafUptX/YhInSZLK6fTTYbPN\nwhxz8+dHHU1NKzWRexa4BTiA0My6F7BnqUFJkqQC1GNFDsKf6fLLoWtXOPzweE98XGGlJnL9gLnA\naEIT6y6E9VclSZKK16UL3HxzaGY988yoo6lZpfaRO6QcQUiSpCLUa0Uup2fPMC3JVlvBwIFw4olR\nR1RzSk3kehBWchiR3c/VPg8r8bySJEmwwgpw330hmVt5Zdh336gjqimlNq1eT5h6ZAzQBAwB/lfi\nOSVJUiHqvSKX09AAf/sbHH88lHF5q3pQaiK3JmHx+v8RJgLeCdi81KAkSVKBkpDIAWy4Yegzt+++\nMGHC0t+fEKUmcrllsj4hrK/aHxhQ4jklSVIhkjaac9tt4bLLYKed4PXXo46mJpTaR+6PwHLAT4C7\ngd6ECp0kSaqGpFTkcvbZB95/H0aPhscfh5VWijqiSJUjkQN4GFitxHNJkqSOSFpFLufYY+G992DH\nHUOfub5RrgwarWITuYMIAx1OyzuWIawVlgEuLDEuSZK0NEkZ7NCasWPh3Xdhjz3g73+Hbt2ijigS\nxfaR65l97ENoTu2d3c9tkiRJlZNKwe9+B8stB9/9LixcGHVEkajHND6TSWqpWZKULBMmwEEHwcSJ\nUUcSnS+/DIMfhg8Py3rFrEKZCvEWHXSxTauX5u3nmlRz+wBOvSxJkiqvWzcYNw6++U045xxobIw6\noqoqtmn1uezWDdgYmAJMBTYCupYnNEmS1K4k95HL17dv6Cd3ww1wxRVRR1NVxVbkrs0+HgNsBczP\nPr8CeKzEmCRJkjpmpZXg/vth663Dsl777BN1RFVR6vQj/YG+wAfZ532yxyRJUqVZkVvc6quHytzo\n0WEQxHbbRR1RxZW6ssO5wPOE5bmuy+7/utSgJEmSirLBBnDbbXDAAfDcc1FHU3GlJnJ/ArYAxgF3\nZPevLfGckiSpEFbkWveNb8CVV8Kuu8LUqVFHU1GlNq0CfAG8A3QHhme3R8pwXkmSpOLsvjvMng07\n7BCW8lpllagjqohSE7kjCFONDAZeJFTkngS2LfG8kiRpaazIte/73w9LeY0ZAw8/DP3rrxt/qU2r\nJwGjgDeBbxKmH/mk1KCWYgwwmTDdyRkVvpYkSbXNRK59Z54J6TTsthvMnRt1NGVXaiL3BZD7W+lO\nSLDWKvGc7ekM/I6QzI0ADgDWqeD1JEmqXa5ktHSpFFx0EQwaFAZALFgQdURlVWoiNxNYFrgTeBC4\nG5hR4jnbMwqYlr3GfOBmYLcKXk+SpNpmRW7pOnWCa6+FL76Ao4+uqwS41ERud+AjoBE4G7gqe6xS\nBgFv5T2fmT0mSVLy1FFCUnFdu8Jf/wqTJsHZZ0cdTdmUMtihC/ASsHb2eVPJ0SxdYT+xN9zQvN/y\nfyr5z9t7rSPvrdfzRHHNUs7TqVM4lnssZL+ar+c2SSoXBzt0TO/ecO+9sPnmMGwYHHxw1BGVrJRE\nbgHwGjCUMNihGmYBQ/KeDyFU5RbTeO65X+2nV1yR9EorhSf5/3Np+b+Y9p4X+1qczxPFNUs9TyYD\nixYt/ljIfrVez2kvEezceelbly6Fva8jW5cusMwyzVvXros/trVf6Otdu0L37tCjR7iWJEVlhRXg\nnnvgppsiuXxTUxNNTU1lO1+pafyjhJGq44HPsscywLdLPG9buhCSx+2At7PXPQB4Ne89mYylZtWi\nXILXVqK3aBEsXFjYtmBB4e8t9HwLFsC8eTB/fthy+4Uea+/1efNC35S5c0PSmkvqundv3vKft7ff\no0f4X3Xv3tCrV/N+y2M9e4bkWKpnTz8NJ5wA48dHHYmKlAoV1aLzsVL/a9xaI3Mls6gFwPHA/YQR\nrFezeBIn1S6bVoMFC0JC98UXzVtHns+eDW++Cf/73+LbZ58t/nzu3OakLz/h69s3zCVV6NatW9R/\nY5LUpmITuRQhYWsq4D3l9o/sJimOunSBPn3CVkmLFsHnny+e4M2ZE7aPP27ePvwQXn897H/00eKv\nffxxaHru3z80x6ywAgwY0Ly19nyFFUJzslQN9pFLvGITuSbgb8BdwJQWr61FGLm6M/CNoiOTpFJ0\n6tRchcv1k+2oTCZU9j76KFQCZ8+G999v3l55pXk/99qHH4aq38CBbW+DBoWYTPgklajYRG408B3g\nMmA9YA6hAtebMJL1BuBb5QhQkiKTSoW+dj17huSrEIsWwQcfwDvvwNtvw6xZ4fGll+CBB8L+22+H\npG+55WDoUGhoaH7MbUOHhiZhqT1W5BKvHHe/M7BCdn82sLAM5yyFgx0k1b6FC+Hdd+E//4EZM5bc\n3nwzND83NMDqq8Naa8Hw4c2PfftGGb1qxZNPwimnwFNPRR2JihT1YAcIidu7ZTiPJCVH587NTa1b\nbLHk64sWhcW+Z8yAadNgyhS4++7wOGVKSORySd3aa8N668HIkaHJ1gpNcliRSzwndJKkWtSpE6y8\ncthaJnqLFoXm2ddeC0ndq6/C3/4GEyaEBHH99UNSN3Jk2F933TCCV1LdMZGTpLjp1AkGDw7bdts1\nH89kQt+8SZNg4kT497/h4otDsjd8OGy2GYwaFR7XW8/JmeuBFbnEK8e3uAFYE/gn0DN7zk/LcF5J\nUkekUs3NtTvs0Hz8yy9DYjd+PDz+OFx0Ueibt8EGIbHbaivYemtYccXoYpdUlFITuSOBI4DlgDWA\nwcAVhJUXJEm1oFu3UIXbbLPmY59+Cs89F5K7a6+F738/NON+4xvN25AhbZ5SNcKKXOKVevcnAKOA\npwhLdQFMAtYv8bylcNSqJHXUwoWhSfaRR5q3Xr1g++1hzJjQhNuvX9RRqqXHHoMf/Sg8KpZKHbVa\n6kKEX2a3nC5UdokuSVIldO4MG24IJ54If/1rmBrlb3+DddaBK68M/fG23hp+8Qt45pkw4ELRs3CR\neKU2rT4MnEXoG7c9cCxwT6lBSZIilkqF0a7rrhvmKZs7N1Tp7r8fDj44LF+2226wxx6QTkPXrlFH\nnEw2rSZeqXe/M3A4YaUHCIvZX0W0VTmbViWp0qZOhXHj4I47wqjYnXaCPfeEHXd0qpNqeuQROOss\nePTRqCNRkUptWi1nGr8cMITQby5KJnKSVE2zZsFdd8Htt8MLL4Qq3Xe/C9tsE6ZKUeU8/DCcfXZI\n6BRLUfeRexjoS0jingP+CFxU4jklSXEyaBAceyz8619hwMSIEXDqqWG92NNPD8ckVUSpiVw/wpxx\newJ/Joxg/VapQUmSYmrQIDjttFCZu+++MOnwTjvBllvCddeFvnYqH/vIJV6piVxnYBVgX+De7DHb\nNSVJYaDEr34Fb7wBP/4x3HprmJvupJPglVeijk6qC6Umcj8jDHCYDownTAo8tdSgJEl1pEsX2HVX\nuPdeePZZ6NMnzEu3ww7w4INOoVEKK3KJV49338EOklTrvvwSbroJLrggzGF32mmw//5OY9JRDz0E\n55wDTU1RR6IiRT3YoQdwPHA58Kfsdk2J55Qk1btu3eCQQ8JAiPPOg+uvhzXWgMsugy++iDq6+LAi\nl3ilJnLXAysBY4Amwlqr/yvxnJKkpEilmptYx40LAySGDTOhkwpUaiK3JnA2IXm7DtgJ2LzUoCRJ\nCbTppnDPPYsndL//PSxYEHVktcuKXOKVmsjNyz5+AqwP9AcGlHhOSVKS5RK6O+4II11HjgwDJez/\nLC2h1ETuj4TJgH8C3A28ApxXalCSJLHZZmGS4fPOgx/8ALbfHl58MeqoaosVucQrRyL3IWGFh9UI\n1bjfd/AchwAvAkPzjvUA/l6G+CRJcZZKwS67hEERe+0FY8bA0UfDRx9FHVntMJFLtFITpe7Ad4Cz\ngJ8CY7OPHXEt8EKLYz8A/g9YVGJ8kqR60KULHHMMTJ4cpisZMQL+8hebW5P+51fJidxdwLeB+cBn\neVsxcj+NqwAjgQfyXtsfeB34iJD0zQFeyr5XkpQU/fuHEa133QUXXhgmFp4yJeqoomPTauJ1KfHz\ng4AdyhEIzZPh/ZxQ1esEPEmYo64L0AhsCJxKSOI2Ab4s07UlSXEyahSMHx+Suq9/HX7yEzjhBOhk\njxwlS6k/8U8QqmelOBBYHjgd+BYwF3gVOJowL10GuJrQB28C0JPQ5GoSJ0lJ1qVLWLf1iSfC6Nbt\ntgvruiaJFbnEK7YiNyn72Bk4FHiD5sQqw9KTu2OBI7Lv3Qm4MXv8LuBwYGVgM0JlLvcTOpIwMnZd\nYGJ7J29sbPxqP51Ok06nlxKOJCm2hg2DRx4JTa2jRsG558Jhh5ngqCY1NTXRVMYl1Yr9KW+guU9b\na+eYUcQ59yRMMHwesA8hieubvc4wQvI4Mvu+DYEftXEe11qVpKR6+eWwZuvIkWEy4T59oo6osu67\nDy66CO6/P+pIVKSo1lp9l5BQnU7oIzeTkLzlto7qChxDGKkKcBthguGfZrduhAEV84DJhObY3YqM\nXZJUr9ZdF55+Gnr0CBMLT2y3AUeKvWKbVq8jJFWPEZpGRwAnlRDHPGD7Nq6Ts3728UVg1RKuJUmq\nZz17wlVXwfXXh35zv/oVHHFE1FFVhn3kEq/YRG4dmhOrq4BnyhOOJEllctBBoSq3117wwgtw8cWw\nzDJRRyWVVbFNqwva2JckqXassw48+SS8+SbssAN88EHUEZWXFbnEKzaRG0mYlDe3rZ+3/2l5QpMk\nqQz69YO77w5rt44aFQZESHWi2KbVzmWNQpKkSurcGX7zG1hvPdh2W7jtNvjGN6KOqnRW5BLPKbAl\nSclx0EFhjda994Zx46KORipZqUt0SZIUL9tvH+Zf22UXePddOProqCMqnhW5xDORkyQlz8Ybw6OP\nhgEQH38MP2prjvkYMJFLNBM5SVIyrbFGWNpr221hwQL4yU+ijqjjXMko8UzkJEnJNXAgPPRQmDh4\n4UIYOzbqiDrGptXEM5GTJCXbKqssnsz97GdRRyQVzEROkqSVVgrJXDoNffrAD38YdUSFsSKXeCZy\nkiQBDBgADzwAW20Fyy0Hhx8edUTSUpnISZKUM2hQSOa22QaWXRb23DPqiNpnRS7xTOQkSco3bBjc\ne2+YmqRfv9B3TqpRruwgSVJLG20UlvE64IDaXpvVilzimchJktSabbaB3/4Wdt0V3nsv6mikVpnI\nSZLUloMOgu98B3bbDebOjTqaJVmRSzwTOUmS2nPOOTB0aBjF6koKqjEmcpIktadTJ/jTn+C11+D/\n/i/qaBZnRS7xHLUqSdLS9OgBt98Om28OG28c+s9JNcCKnCRJhWhogOuvDyNZZ82KOprAilzimchJ\nklSo0aPhuONg331h/vyoowlM5BLNRE6SpI4480zo2xcaG6OOxMEXqolE7mjgJWBo3rHewD+iCUeS\npHZ06gTXXhsGQDz0ULSx2LSaeLWQyP0eeKbFsTOA8yOIRZKkpVtppZDMfe97MHt21NEowWohkcvJ\n1YeHAMOBf+e9dhrwH+Bd4AXgS2A8jrqVJEVl9GjYf3847LDomjityCVeLSVyuZ/Ec4BfAKcDfwbW\nARYCFwG/BDYG3gRGAQuqH6YkSVm//CXMnBmaWaUI1EIidyCwPCFx2x6Yk91uAC4GjgL+DxgJTADW\nAqZEEqkkSfm6dg1NrGecAW+9Vf3rW5FLvKiaJo8FjiA0p+4E3Jg9fjdwCPBh9vnPCUkchERuIvAt\nQvNqmxrzRhKl02nS6XRZgpYkaQkjR8LJJ8P3vw/33WdipXY1NTXR1NRUtvPV0k/b3sCqwIXZ5+cC\nTxCSuy7AZGBNwkCITwiDJFqTyTgcW5JUTQsWwBZbwFFHwRFHVO+6N90Ed90FN99cvWuqrFIh8S86\nH6uVwQJdCRW6nbPPv0aYlmQHYEv4//buNEiK8o7j+HeXFXXxTrwDjOVtjMSYROIRFw9INJYajDEU\n8UAgXkA0RcB4oWUqpaJ4IWg8kFLxSNRoNCgVXXKAV0Q0IqtBRSMYX2wMIAE8yIunNzssO7uz7Mx0\n9/T3UzW10z19PLsPPfPj/3RPczfQFL32CnAL8BLhggdJkuJVVxeGWAcMgEGDoE+fuFukjEhKkFtD\nCG0t5gBbtVmmJeT9gVC5kyQpOfbdF0aNgtGj4ZFHKrNPz5HLvCRc7CBJUnUYNw4WLgzDnVIFGOQk\nSSqVjTeGKVNCZW7FivLvz4pc5hnkJEkqpQEDoKEhGfdiVdUzyEmSVGoTJ8L06TB/fnn3Y0Uu8wxy\nkgcyxKkAAAssSURBVCSV2nbbweWXw5gx5b99l0Eu0wxykiSVw/Dh0NwMDz1Uvn34vamZZ5CTJKkc\n6urguutg7FhYtap8+7Eil2kGOUmSyuXww6FfP5g0qTzbtyKXeQY5SZLKaeJEuOYaWLq09Nv2YofM\nM8hJklROu+4Kw4bBJZfE3RJVIYOcJEnldsEF4bZdCxeWdrtW5DLPICdJUrltvXW46OHCC+NuiaqM\nQU6SpEo491x47rnwKBUrcplnkJMkqRLq6+HSS2H8eK82VckY5CRJqpTTT4clS+Cpp0qzPStymWeQ\nkySpUurq4IorwrlyVuVUAgY5SZIqafBgWL0annii+9uyIpd5BjlJkiqptjZ8p9xll1mVU7cZ5CRJ\nqrTBg2HlSpg5s3vbsSKXeQY5SZIqrbYWLr64NFU5g1ymGeQkSYrDiSfCsmUwa9aGb8Oh2cwzyEmS\nFIcePUpTlbMil2lJCHKTgb8D0/Lm9QJmAv7rlCRVr5NOguZmePrpDVvfilzmJSHInQP0B3J588YB\nVwP+C5UkVa8ePWDcOLjyyg1b34sdMi8JQW474Hng9Wi6N7An8Me8ZU4G3gL+DcwDlhOqeDtWrpmS\nJJXBkCGwYAHMmxd3S5RCSQhyuwH7ALsDXwAuBy4iDKvOBU4lDLVOAO4E9gcWAwcASyvfXEmSSqhn\nTzjvPLjqqq6va0Uu85IQ5M4H5gOvEcJcM/AmMAzYnjC8ejuwS7RcPfA5sDqOxkqSVHIjRoSrV99+\nO+6WKGXqYtrv2cAIQkg7Gvggmv97YCiwNXAYcAWtFzzsBzwKfBl4paONT5gw4f/PGxoaaGhoKFnD\nJUkquS22CGHu2mvhxhuLX8+KXOo0NjbS2NhYsu0lqfeHANsANxHC3QlAH0LYO4ZQsdsP+D7wVWB8\nge2sXetVPJKktPngA9hnH2hqgm23LW6dKVNg/nyYOrW8bVPZ1IQgvsF5LAlDqwCbAqcAN0fTTxAq\ndvdGj02BT4A1wEJC6Duu8s2UJKlMdtghfEnw5MnFr2NFLvOqsfetyEmS0umNN+CQQ+Cdd6C+vvPl\nb74ZXn01VOaUStVSkZMkSXvsAf37wz33FLe8FbnMM8hJkpQko0fD9dd71wYVxSAnSVKSHHFE+FnM\nbbusyGWeQU6SpCSpqWmtykmdMMhJkpQ0Q4fC3LmwaFHHy1mRyzyDnCRJSVNfD2ecUdyXAxvkMs0g\nJ0lSEp1zDkyfDsuWFV7GCyIyzyAnSVIS9e4NRx4J06YVXsah1cwzyEmSlFRjxoTh1c8/j7slSiiD\nnCRJSXXQQbDlljBzZvuvW5HLPIOcJElJVVMDZ53lLbhUkEFOkqQkO/lkmDMHFi9e/zUrcplnkJMk\nKcl69QrfK3frrXG3RAlkkJMkKenOPBNuvx3WrFl3vhW5zDPISZKUdHvvHR4PPxx3S5QwBjlJktKg\nvYserMhlnkFOkqQ0OP54aGqCBQvibokSxCAnSVIa9OwZ7r86dWrrPCtymWeQkyQpLUaOhHvugY8/\njrslSgiDnCRJadGnDxx8MMyYEaatyGWeQU6SpDRpe9GDQS7TDHKSJKXJoEHQ3Awvvhgqcso0g5wk\nSWlSWwvDh4c7PTi0mnkGOUmS0mbYMHjwQVi+PO6WKGYGOSVGY2Nj3E1QN9h/6Wb/pcyOO8KAATBj\nBo3vvRd3axSjtAS5vYD5wCpgHrAEWAwcG2ejVFp+kKSb/Zdu9l8KjRwJb71F47vvxt0SxSgtQW57\n4HzgNWB/4FHgDOCxOBtVSKnfEDd0e11Zr5hlO1qmq68l+UMjbf3X3b7r6PWuzk+CUratGo+9YvcZ\nh7Qde8UsW9Zjb+BA6Ns3MefIZe3Y6+j1Sh57aQlys4FPgVei6X6EylwiVeObUWfLGORKvz2D3IbJ\n2oeJQa7020tNkKuthbFjYcstO9xHpWTt2Ovo9Uoee8mI8YWdDYyInj8ErAAmAU3AngXW+Qewa/mb\nJkmS1G2LgN3ibkQl3A4cAewAPBtzWyRJkmKXlqFVgP0IQ6sfAZsBv4y3OZIkSZIkSZIkSZIkSZIk\nJcJewBTgAcL3zyldjgNuBe4Djoq5Leq6XYDbgAfjboiK1gu4i3DcDYm5Leo6j7l08zOvA7WEMKd0\n2orw5qR08kMlPX4MHBM9vy/OhqhbPObSrajPvDRdtQpwB/Av4NU2878DLATeBMYVWPdY4HF8U4pT\nd/oP4CLgpvI0TUXobv8pXl3pv52Blht4flaR1qkzHn/ptiH9V5WfeYcSbtGV/4foQfgS4BywEfAy\nsDfhf5STgJ3abON3ZW+lCtnQ/qsBriR8j6Di093jz+pAvLrSf0NprcjNqFwT1YGu9F8Lj7nk6Er/\nVf1nXo51/xDfAmbmTY+PHvkOA64HbgF+Ws7GqVM5ut5/o4EXCec5/qScjVOncnS9/7YBpmLFIAly\nFNd/9YQKws3AjyrVOHUqR3H95zGXTDmK679RdOEzr6507YtN/hAAwD+BA9ssMzt6KHmK6b8booeS\np5j+awbOrFiL1BWF+m8lMCyWFqkrCvWfx1w6FOq/UcCNxW4kbefItWdt3A1Qt9h/6Wb/pZv9l272\nX7qVpP+qIci9D/TOm+5NSLVKB/sv3ey/dLP/0s3+S7fM9l+OdceY64BF0fyerH+yp5Ilh/2XZjns\nvzTLYf+lWQ77L81y2H/MAJYAqwnjyqdH878LNBGu/rggnqapCPZfutl/6Wb/pZv9l272nyRJkiRJ\nkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkhLtM2Be3qNPvM0pqa8Ad0TPT2P9G2Q3Agd0sP4DwC4l\nb5Wk1KmLuwGSVMBKYP8Cr9VEP9N60/CxtIa39n6HtQXmt/g1cB4wusTtkpQytXE3QJKKlCPcyuYu\nwv0KexMC0fPAfGBC3rIXRsv+GbgX+Fk0v5HWStcXgbej5z2Aq/O2NTKa3xCt8yDwOnB33j6+AfyV\ncH/EZ4HNgNlAv7xl/kKovuXbGOgPvNDJ71sDHEtrRbIJeCvv9zi6k/UlSZJi8ymtIea3QF/CcOs3\no9cHArdEz2uBx4BDCUHtFWATYHPgTeD8aLlngK9Fz/OD3EhC+IMQtF4gBMcG4CNgJ0KwmgMcRLjB\n9SJaQ+FmhDB4CjApmrcH7Ye1/lFbW5wGfMi6w8jL89rZ4n7grLzp2WTgBtuSOubQqqSk+i/rDq3m\ngMWEqhmEIDeQEHwAegG7E8LbQ8Cq6PFoEfsaSKicnRhNbwHsBnwS7W9JNP9lwrlpy4GlwN+i+Sui\nn78BLiZUCocBd7azr77Rui3WAvex7jDpM23W+TlhqHlK3rwlhL/J6538bpKqmEFOUpp83Gb6V8Ct\nbeaNofUcOto8/5TWU0o2abPeucCsNvMagNV5058R3jcLnb+2MtrG8cAPWL+qRrRuTZt5bafzHQkM\nBr7dzjqfd7CepAzwHDlJafUkoerVK5reGdgW+BMhSLUMrX4vb513gK9Hz0/Mm/8kcDat/7ndA6gv\nsN+1hPPVdszb1uaEoVWA24AbCJW8/7Sz/mJgh7zpjkJcX2AycBLrBkqi/S/uYF1JGWBFTlJSFbqa\ns8Uswjlic6Pp5cBQwlDr/YSLFj4knKfWEpYmEr66YyTweN72biMMU74ULfshcAKFrx79BPgh4crT\nTQmVuKMIFcOXCAGuvWFVonbt2eZ3am8fNcCpwDbAI9G89wnBdCPgS8DCAvuQJEmqCpfSetVqJexE\nqNh1ZBpwYDf2MRC4vhvrS6oSDq1KyoJKfd/cKYSvIvlFJ8tNBM7sxn6G03p1rCRJkiRJkiRJkiRJ\nkiRJkiRJkiRJkiRJkiRJyp7/Af4lzyISHNnnAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 71 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 71 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/python/iir++.h b/python/iir++.h new file mode 120000 index 0000000..85f0bef --- /dev/null +++ b/python/iir++.h @@ -0,0 +1 @@ +../obj/iir++.h \ No newline at end of file diff --git a/python/pyiir.i b/python/pyiir.i new file mode 100644 index 0000000..42123cd --- /dev/null +++ b/python/pyiir.i @@ -0,0 +1,142 @@ +%module iir +%include + +/* + * List conversion + * + * SWIG knows how to turn std::vector func() into a function that returns a + * Python list. We do need to tell it about all types with %template, + * however. + */ +%include +%template(vectord) std::vector; + + + +%{ +#include "iir++.h" +%} + + + +/* + * Add bounds checking + * + * The Python library is targetted at experimentation. It is therefore + * reasonable to protect its object query and creation interfaces with code + * that performs bounds checking on arguments. + */ +%exception IIR::Coeff::c { + if(arg2 < 0 || arg2 >= arg1->nc()) { + char buf[80]; + snprintf(buf, sizeof(buf), "Index of c coefficient (%d) out of bounds " + "(0,%d).", arg2, arg1->nc() - 1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + $action +} + +%exception IIR::Coeff::d { + if(arg2 < 0 || arg2 >= arg1->nd()) { + char buf[80]; + snprintf(buf, sizeof(buf), "Index of d coefficient (%d) out of bounds " + "(0,%d).", arg2, arg1->nd() - 1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + $action +} + +%exception IIR::Filter::getCoeffSet { + if(arg2 < 0 || arg2 >= arg1->numCoeffSet()) { + char buf[80]; + snprintf(buf, sizeof(buf), "Index of coefficient set (%d) out of bounds " + "(0,%d).", arg2, arg1->numCoeffSet() - 1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + $action +} + +%exception IIR::ButterworthLowPass { + if(arg1 < 1) { + char buf[80]; + snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.", + arg1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + $action +} + +%exception IIR::ButterworthHighPass { + if(arg1 < 1) { + char buf[80]; + snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.", + arg1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + $action +} + +%exception IIR::ButterworthBandPass { + if(arg1 < 1) { + char buf[80]; + snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.", + arg1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + if(arg3 > arg4) { + PyErr_SetString(PyExc_RuntimeError, "Low corner frequency exceeds " + "high corner frequency."); + return NULL; + } + $action +} + +%exception IIR::ButterworthBandStop { + if(arg1 < 1) { + char buf[80]; + snprintf(buf, sizeof(buf), "Invalid filter order (%d). Must be > 0.", + arg1); + PyErr_SetString(PyExc_RuntimeError, buf); + return NULL; + } + if(arg3 > arg4) { + PyErr_SetString(PyExc_RuntimeError, "Low corner frequency exceeds " + "high corner frequency."); + return NULL; + } + $action +} + + + +/* + * Invalid filter string handling + * + * The C++ library uses a "valid()" member in IIR::Filter to avoid throwing an + * exception if the constructor is passed an invalid filter string. Python's + * exceptions are a much nicer fit, however, so use those instead. + * + * Notes: + * 1. Doesn't match if we prepend it with the namespace. + * 2. Applies to all constructors. + */ +%exception Filter { + $action + if(!result->valid()) { + PyErr_SetString(PyExc_RuntimeError, "Invalid filter description string."); + delete result; + return NULL; + } +} + + + +%include "iir++.h" + +// vim: ts=4:sw=4 diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 0000000..1b36f11 --- /dev/null +++ b/python/setup.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python3 + +from distutils.core import setup, Extension + +iir_module = Extension('_iir', + sources=['pyiir_wrap.cxx'], + libraries=['iir++'], + ) + +setup(name = 'iir', + version = '1.0.3', + author = 'Laurence Withers', + description = 'IIR filter', + ext_modules = [iir_module], + py_modules = ['iir'], + ) diff --git a/scripts/.gitignore b/scripts/.gitignore index edcbbeb..b9fcf56 100644 --- a/scripts/.gitignore +++ b/scripts/.gitignore @@ -7,6 +7,7 @@ build.docs.none build.files.none build.firmware.gpasm build.firmware.sdcc +build.header.c build.lib.c build.lib.c++ build.make.none diff --git a/src/docs/MainPage.dox b/src/docs/MainPage.dox index f915aae..ee9a49e 100644 --- a/src/docs/MainPage.dox +++ b/src/docs/MainPage.dox @@ -1,16 +1,16 @@ /* libiir/src/docs/MainPage.dox * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ /*! \mainpage -This library allows the implementation of arbitrary IIR filters in C. It has +This library allows the implementation of arbitrary %IIR filters in C. It has functions for generating and manipulating filters in terms of coefficients, for chaining arbitrary filters together, and for generating coefficients for some -common types of filter. See \ref iir_structure for a definition of the IIR +common types of filter. See \ref iir_structure for a definition of the %IIR filter equation. \section creation Filter creation @@ -23,7 +23,7 @@ Otherwise, the library user must first create a set of coefficients using \ref iir_coeff_new(). Any number of filters can then be instantiated using that set of coefficients with \ref iir_filter_new(), or the coefficients can be chained on to the end of an existing filter instance with -\ref iir_filter_chain(). See \ref common_filters for functions to generate +\ref iir_filter_chain(). See \ref libiir_common_filters for functions to generate coefficients. \section operation Filter operation @@ -47,6 +47,11 @@ installed) which will generate a Bode plot for a given filter chain. Note the phase response can be a little rough due to the simplistic time-domain analysis of the output signal's phase. +\section cpp C++ + +Two C++ wrapper objects, \ref IIR::Coeff and \ref IIR::Filter, are supplied which provide +a true object-oriented interface on top of the C implementation. + */ /* options for text editors diff --git a/src/docs/iir_structure.dox b/src/docs/iir_structure.dox index 9831e5e..b258b96 100644 --- a/src/docs/iir_structure.dox +++ b/src/docs/iir_structure.dox @@ -1,6 +1,6 @@ /* libiir/src/docs/iir_structure.dox * - * Copyright: ©2010–2011, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ @@ -14,7 +14,7 @@ For the purposes of this library, the following notation is used: \li c[n]: array of \c x(t) coefficients \li d[n]: array of \c y(t) coefficients -This leads to a general IIR filter equation: +This leads to a general %IIR filter equation: y(t) = x(t).c[0] + x(t-1).c[1] + … + x(t-N).c[N] - y(t-1).d[0] - y(t-2).d[1] - … - y(t-1-M).d[M] diff --git a/src/docs/string_desc.dox b/src/docs/string_desc.dox index 3161afb..1c8ca01 100644 --- a/src/docs/string_desc.dox +++ b/src/docs/string_desc.dox @@ -1,17 +1,17 @@ /* libiir/src/docs/string_desc.dox * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ /*! \page string_desc Describing IIR filters as strings -This library allows the user to describe an IIR filter chain as a string, which +This library allows the user to describe an %IIR filter chain as a string, which is useful to allow configurable filtering using e.g. a configuration file. This page describes the format of such strings. -The string is first split into individual IIR filters. Each filter is written as +The string is first split into individual %IIR filters. Each filter is written as \c type(params) and separated by whitespace. The string must contain at least one filter but may contain an arbitrary number. @@ -23,7 +23,7 @@ one filter but may contain an arbitrary number. \subsection string_desc_coeff Raw coefficients -An IIR filter may be specified as raw coefficients, in which case the \c type +An %IIR filter may be specified as raw coefficients, in which case the \c type is \c raw and the \c params consists of a string: c[0],c[1],…,c[n]/d[0],d[1]…d[n] diff --git a/src/libiir++/.params b/src/libiir++/.params new file mode 100644 index 0000000..7397660 --- /dev/null +++ b/src/libiir++/.params @@ -0,0 +1 @@ +lib c++ libiir++ iir++.h diff --git a/src/libiir++/000_TopHeader.h b/src/libiir++/000_TopHeader.h new file mode 100644 index 0000000..7f4e712 --- /dev/null +++ b/src/libiir++/000_TopHeader.h @@ -0,0 +1,32 @@ +/* libiir/src/libiir++/000_TopHeader.h + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + +#ifndef HEADER_libiirpp +#define HEADER_libiirpp + +// standard includes, or includes needed for type declarations +#include +#include + +// forward declares of opaque C structures (to avoid pulling in header) +extern "C" { +struct iir_coeff_t; +struct iir_filter_t; +} + +/*! \defgroup libiirpp C++ library */ + +/*! \brief All %IIR library functions and classes. + +\ingroup libiirpp + +*/ +namespace IIR { + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/000_TopSource.cc b/src/libiir++/000_TopSource.cc new file mode 100644 index 0000000..bd47320 --- /dev/null +++ b/src/libiir++/000_TopSource.cc @@ -0,0 +1,26 @@ +/* libiir/src/libiir++/000_TopSource.cc + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + +// we need to include "iir.h" first, as it pulls in , which we need +// to take effect before "iir++.h" pulls in +#include "iir.h" + +// now remove the preprocessor definition of complex to _Complex, which is fine +// for the C header but not good for the C++ header +#undef complex + +#include "iir++.h" + +// Below are all the includes used throughout the library. +#include +#include + +namespace IIR { + +/* options for text editors +vim: expandtab:ts=4:sw=4 +*/ diff --git a/src/libiir++/100_Coeff.cc b/src/libiir++/100_Coeff.cc new file mode 100644 index 0000000..7036c2c --- /dev/null +++ b/src/libiir++/100_Coeff.cc @@ -0,0 +1,98 @@ +/* libiir/src/libiir++/100_Coeff.cc + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +Coeff::Coeff(int nc, const double* c, int nd, const double* d) +{ + cc_ = iir_coeff_new(nc, c, nd, d); +} + + + +Coeff::Coeff(const std::vector& c, const std::vector& d) +{ + int nc, nd; + double* ac, * ad; + + nc = c.size(); + ac = static_cast(alloca(nc * sizeof(double))); + for(int i = 0; i < nc; ++i) ac[i] = c[i]; + + nd = d.size(); + ad = static_cast(alloca(nd * sizeof(double))); + for(int i = 0; i < nd; ++i) ad[i] = d[i]; + + cc_ = iir_coeff_new(nc, ac, nd, ad); +} + + + +Coeff::Coeff(const struct iir_coeff_t* cc) +{ + cc_ = iir_coeff_copy(cc); +} + + + +Coeff::Coeff(const Coeff& other) +{ + cc_ = iir_coeff_copy(other.cc_); +} + + + +Coeff::~Coeff() +{ + iir_coeff_free(cc_); +} + + + +int +Coeff::nc() const +{ + return iir_coeff_get_nc(cc_); +} + + + +double +Coeff::c(int idx) const +{ + return iir_coeff_get_c(cc_, idx); +} + + + +int +Coeff::nd() const +{ + return iir_coeff_get_nd(cc_); +} + + + +double +Coeff::d(int idx) const +{ + return iir_coeff_get_d(cc_, idx); +} + + + +std::complex +Coeff::response(double freq) const +{ + return std::complex(iir_response_c(cc_, freq)); +} + + + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/100_Coeff.h b/src/libiir++/100_Coeff.h new file mode 100644 index 0000000..8a3dbf0 --- /dev/null +++ b/src/libiir++/100_Coeff.h @@ -0,0 +1,130 @@ +/* libiir/src/libiir++/100_Coeff.h + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +/*! \defgroup libiirpp_coeff Coefficient handling + +\ingroup libiirpp + +A set of coefficients for a single stage %IIR filter is represented by an +instance of the \ref Coeff class. This class is often instantiated through +one of the common filter generator functions (see \ref libiirpp_common_filters). + +*/ + + + +/*! \brief %IIR filter coefficient set + +\ingroup libiirpp_coeff + +An instance of this class represents the set of coefficients used to implement +a single-phase %IIR filter. See \ref iir_structure for more information. + +*/ +class Coeff { +public: + /*! \brief Constructor with raw coefficients. + + \param nc Number of \a c coefficients, ≥ 1. + \param c Array of \a c coefficients. + \param nd Number of \a d coefficients, ≥ 1. + \param d Array of \a d coefficients. + + */ + Coeff(int nc, const double* c, int nd, const double* d); + + /*! \brief Constructor with raw coefficients. + + \param c Array of \a c coefficients, size ≥ 1. + \param d Array of \a d coefficients, size ≥ 1. + + */ + Coeff(const std::vector& c, const std::vector& d); + + /*! \brief Constructor from C library object. + + \param cc Pointer to existing coefficient object. + + */ + Coeff(const struct iir_coeff_t* cc); + + /*! \brief Copy constructor. + + \param other Object to copy. + + Performs a full deep copy of the set of coefficients from \a other. + + */ + Coeff(const Coeff& other); + + /*! \brief Destructor. */ + virtual ~Coeff(); + + + + /*! \brief Count number of \a c coefficients. + + \returns Number of \a c coefficients, ≥ 1. + + */ + int nc() const; + + /*! \brief Retrieve \a c coefficient by index. + + \param idx Index of coefficient (0 ≤ \a idx < \ref nc()). + \returns Value of \a c coefficient at index \a idx. + + */ + double c(int idx) const; + + /*! \brief Count number of \a d coefficients. + + \returns Number of \a d coefficients, ≥ 1. + + */ + int nd() const; + + /*! \brief Retrieve \a d coefficient by index. + + \param idx Index of coefficient (0 ≤ \a idx < \ref nd()). + \returns Value of \a d coefficient at index \a idx. + + */ + double d(int idx) const; + + + + /*! \brief Evaluate response of filter at given frequency. + + \param freq Frequency (0 ≤ \a freq ≤ 1). + \returns Complex response of filter at given frequency. + + Returns the steady-state response of the filter to a pure input signal of + the given frequency. \a freq is expressed as a fraction of the intended + sampling rate (with 0.5 being the Nyquist frequency). + + */ + std::complex response(double freq) const; + +private: + // allow the Filter class to access cc_ directly + friend class Filter; + + // our internal representation is the C object + struct iir_coeff_t* cc_; + + // disallow assignment; object is immutable + Coeff& operator=(const Coeff& other); +}; + + + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/200_Filter.cc b/src/libiir++/200_Filter.cc new file mode 100644 index 0000000..d0fd159 --- /dev/null +++ b/src/libiir++/200_Filter.cc @@ -0,0 +1,119 @@ +/* libiir/src/libiir++/200_Filter.cc + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +Filter::Filter(const Coeff& coeff) +{ + coeff_.push_back(new Coeff(coeff)); + fi_ = iir_filter_new(coeff.cc_); +} + + + +Filter::Filter(const Filter& other) +{ + for(int pos = 0, max = other.coeff_.size(); pos < max; ++pos) { + Coeff* cc = new Coeff(*other.coeff_[pos]); + coeff_.push_back(cc); + + if(!pos) fi_ = iir_filter_new(cc->cc_); + else iir_filter_chain(fi_, cc->cc_); + } +} + + + +Filter::Filter(const char* desc) +{ + fi_ = iir_parse(desc); + if(!fi_) return; + + for(int pos = 0, max = iir_filter_coeff_sets(fi_); pos < max; ++pos) { + struct iir_coeff_t* c = iir_filter_coeff_set(fi_, pos); + coeff_.push_back(new Coeff(c)); + iir_coeff_free(c); + } +} + + + +Filter::~Filter() +{ + for(int pos = 0, max = coeff_.size(); pos < max; ++pos) delete coeff_[pos]; + iir_filter_free(fi_); +} + + + +bool +Filter::valid() const +{ + return fi_; +} + + + +double +Filter::f(double x) +{ + return iir_filter(fi_, x); +} + + + +int +Filter::numCoeffSet() const +{ + return coeff_.size(); +} + + + +const Coeff& +Filter::getCoeffSet(int idx) const +{ + return *(coeff_[idx]); +} + + + +void +Filter::chainCoeffSet(const Coeff& coeff) +{ + coeff_.push_back(new Coeff(coeff)); + iir_filter_chain(fi_, coeff.cc_); +} + + + +std::complex +Filter::response(double freq) const +{ + return iir_response(fi_, freq); +} + + + +std::vector +Filter::responseMinus3dB(double gain) const +{ + double* cf; + int num_cf, i; + + iir_response_minus_3dB(fi_, gain, &cf, &num_cf); + std::vector ret(num_cf); + for(i = 0; i < num_cf; ++i) ret[i] = cf[i]; + free(cf); + return ret; +} + + + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/200_Filter.h b/src/libiir++/200_Filter.h new file mode 100644 index 0000000..5a4d79f --- /dev/null +++ b/src/libiir++/200_Filter.h @@ -0,0 +1,170 @@ +/* libiir/src/libiir++/200_Filter.h + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +/*! \defgroup libiirpp_filter Basic IIR filtering + +\ingroup libiirpp + +An instance of the \ref Filter class represents the structure and state of an +%IIR filter. It may be built up with \ref Coeff instances or parsed from a +user-supplied string (see \ref string_desc). + +*/ + + + +/*! \brief %IIR filter coefficient set + +\ingroup libiirpp_filter + +An instance of this class contains both a copy of the set of coefficients used +in each filter stage (of which there may be multiple), and all of the state +(i.e. previous inputs and outputs) associated with each stage. + +*/ +class Filter { +public: + /*! \brief Construct from coefficient set. + + \param coeff Coefficient set. + + This constructor builds a filter chain with one stage. Additional stages may + be added via \ref chainCoeffSet(). This constructor will always result in a + valid object. + + */ + Filter(const Coeff& coeff); + + /*! \brief Copy constructor. + + \param other Filter to copy. + + Copies the structure (but not the state) of the filter \a other. This + constructor will always result in a valid object. + + */ + Filter(const Filter& other); + + /*! \brief Construct from filter description string. + + \param desc Description string. + + This constructor will parse the description string \a desc (see + \ref string_desc for details on what is accepted) and build up a filter + as described. + + If the parsing step fails, then \ref valid() will return false and the + object must not be used (no other functions except the destructor will be + safe to call). + + */ + Filter(const char* desc); + + /*! \brief Destructor. */ + virtual ~Filter(); + + + + /*! \brief Determine whether object is valid. + + \retval true if object is valid. + \retval false if object is not valid. + + A Filter object might not be valid if it was constructed from a description + string that contained an error. In that case this function will return + \a false and the object should not be used (no other functions except the + destructor will be safe to call). + + */ + bool valid() const; + + + + /*! \brief %Filter a sample. + + \param x Sample. + \returns %Filter output. + + Injects a single sample into the filter and returns the output value. + + */ + double f(double x); + + + + /*! \brief Retrieve number of coefficient sets. + + \returns Number of sets in use (≥1). + + */ + int numCoeffSet() const; + + /*! \brief Retrieve coefficient set. + + \param idx Index of set (0 ≤ \a idx < \ref numCoeffSet()). + \returns Reference to coefficient set. + + */ + const Coeff& getCoeffSet(int idx) const; + + /*! \brief Add another set of coefficients to the filter chain. + + \param coeff Coefficient set to add. + + */ + void chainCoeffSet(const Coeff& coeff); + + + + /*! \brief Evaluate response of filter at given frequency. + + \param freq Frequency (0 ≤ \a freq ≤ 1). + \returns Complex response of filter at given frequency. + + Returns the steady-state response of the filter to a pure input signal of + the given frequency. \a freq is expressed as a fraction of the intended + sampling rate (with 0.5 being the Nyquist frequency). + + */ + std::complex response(double freq) const; + + + + /*! \brief Search for -3dB points. + + \param gain Nominal gain of filter. + \returns List of frequency ratios at which filter magnitude is -3dB. + + This function will perform a linear scan of the frequency response of + the filter between DC and the Nyquist frequency. It will then search + for any pair of points between which the magnitude transitions the -3dB + level. It then performs a few levels of binary search in order to provide + a more accurate answer. + + Returned frequencies are expressed as a fraction of the sampling rate. + + */ + std::vector responseMinus3dB(double gain) const; + +private: + // one entry per set of coefficients, in order of execution + std::vector coeff_; + + // the C object has both state and its own copy of the coefficients + struct iir_filter_t* fi_; + + // disallow assignment + Filter& operator=(const Filter& other); +}; + + + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/300_common_filters.cc b/src/libiir++/300_common_filters.cc new file mode 100644 index 0000000..a2788ac --- /dev/null +++ b/src/libiir++/300_common_filters.cc @@ -0,0 +1,58 @@ +/* libiir/src/libiir++/300_common_filters.cc + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +Coeff +ButterworthLowPass(int order, double gain, double corner) +{ + struct iir_coeff_t* cc = iir_butterworth_lowpass(order, gain, corner); + Coeff ret(cc); + iir_coeff_free(cc); + return ret; +} + + + +Coeff +ButterworthHighPass(int order, double gain, double corner) +{ + struct iir_coeff_t* cc = iir_butterworth_highpass(order, gain, corner); + Coeff ret(cc); + iir_coeff_free(cc); + return ret; +} + + + +Coeff +ButterworthBandPass(int order, double gain, double corner1, double corner2) +{ + struct iir_coeff_t* cc = iir_butterworth_bandpass(order, gain, corner1, + corner2); + Coeff ret(cc); + iir_coeff_free(cc); + return ret; +} + + + +Coeff +ButterworthBandStop(int order, double gain, double corner1, double corner2) +{ + struct iir_coeff_t* cc = iir_butterworth_bandstop(order, gain, corner1, + corner2); + Coeff ret(cc); + iir_coeff_free(cc); + return ret; +} + + + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/300_common_filters.h b/src/libiir++/300_common_filters.h new file mode 100644 index 0000000..ddc7c20 --- /dev/null +++ b/src/libiir++/300_common_filters.h @@ -0,0 +1,144 @@ +/* libiir/src/libiir++/300_common_filters.h + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + + + +/*! \defgroup libiirpp_common_filters Common types of IIR filter + +\ingroup libiirpp + +Routines in this section can be used to build a coefficient set that matches a +common %IIR filter type. These can be passed to the \ref Filter class to build +up an actual instance of a filter. + +\note It is common to build higher-order filters by chaining together several + lower-order filters. The functions in this section do not do that directly, + but are intended to build the individual components of a chain. See + \ref Filter::Filter(const char*) for an interface which performs this + automatically. + +*/ +/*!@{*/ + + + +/*! \brief nth-order Butterworth low-pass + +\param order Order of filter (≥1). +\param gain Linear gain of filter. +\param corner Corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). + +\returns %IIR filter coefficient set. + +Uses the Exstrom labs code to compute the coefficients of an nth-order (param +\a order) Butterworth-type low pass filter with gain \a gain and corner +frequency \a corner. + +Note it is recommended to chain multiple filters together to build anything +greater than a 4th-order filter. This function won't do that directly for you. +\a gain will usually be set to be 1.0. + +The corner frequency corner is expressed as a fraction of the sampling +frequency (which is of course not known by the %IIR code). It should lie between +0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). + +*/ +Coeff ButterworthLowPass(int order, double gain, double corner); + + + +/*! \brief nth-order Butterworth high-pass + +\param order Order of filter (≥1). +\param gain Linear gain of filter. +\param corner Corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). + +\returns %IIR filter coefficient set. + +Uses the Exstrom labs code to compute the coefficients of an nth-order (param +\a order) Butterworth-type high pass filter with gain \a gain and corner +frequency \a corner. + +Note it is recommended to chain multiple filters together to build anything +greater than a 4th-order filter. This function won't do that directly for you. +\a gain will usually be set to be 1.0. + +The corner frequency corner is expressed as a fraction of the sampling +frequency (which is of course not known by the %IIR code). It should lie between +0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). + +*/ +Coeff ButterworthHighPass(int order, double gain, double corner); + + + +/*! \brief nth-order Butterworth band-pass + +\param order Order of filter (≥1). +\param gain Linear gain of filter. +\param corner1 Low corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). +\param corner2 High corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). + +\returns %IIR filter coefficient set. + +Uses the Exstrom labs code to compute the coefficients of an nth-order (param +\a order) Butterworth-type band pass filter with gain \a gain and corner +frequencies \a corner1 and \a corner2. + +Note it is recommended to chain multiple filters together to build anything +greater than a 4th-order filter. This function won't do that directly for you. +\a gain will usually be set to be 1.0. + +The corner frequencies \a corner1 and \a corner2 are expressed as a fraction of +the sampling frequency (which is of course not known by the %IIR code). They +should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling +frequency), and \a corner2 should be greater than \a corner1. + +*/ +Coeff ButterworthBandPass(int order, double gain, double corner1, + double corner2); + + + +/*! \brief nth-order Butterworth band-stop + +\param order Order of filter (≥1). +\param gain Linear gain of filter. +\param corner1 Low corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). +\param corner2 High corner frequency expressed as a fraction of Nyquist + (0 ≤ \a corner ≤ 1). + +\returns %IIR filter coefficient set. + +Uses the Exstrom labs code to compute the coefficients of an nth-order (param +\a order) Butterworth-type band stop filter with gain \a gain and corner +frequencies \a corner1 and \a corner2. + +Note it is recommended to chain multiple filters together to build anything +greater than a 4th-order filter. This function won't do that directly for you. +\a gain will usually be set to be 1.0. + +The corner frequencies \a corner1 and \a corner2 are expressed as a fraction of +the sampling frequency (which is of course not known by the %IIR code). They +should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling +frequency), and \a corner2 should be greater than \a corner1. + +*/ +Coeff ButterworthBandStop(int order, double gain, double corner1, + double corner2); + + + +/*!@}*/ +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/999_BottomHeader.h b/src/libiir++/999_BottomHeader.h new file mode 100644 index 0000000..883d3b7 --- /dev/null +++ b/src/libiir++/999_BottomHeader.h @@ -0,0 +1,14 @@ +/* libiir/src/libiir++/999_BottomHeader.h + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + +} // end namespace IIR + +#endif + +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=cpp.doxygen +*/ diff --git a/src/libiir++/999_BottomSource.cc b/src/libiir++/999_BottomSource.cc new file mode 100644 index 0000000..435f4d6 --- /dev/null +++ b/src/libiir++/999_BottomSource.cc @@ -0,0 +1,12 @@ +/* libiir/src/libiir++/999_BottomSource.cpp + * + * Copyright: ©2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 + */ + +} // end namspace IIR + +/* options for text editors +vim: expandtab:ts=4:sw=4 +*/ diff --git a/src/libiir++/build.default b/src/libiir++/build.default new file mode 100644 index 0000000..902b4ed --- /dev/null +++ b/src/libiir++/build.default @@ -0,0 +1 @@ +source src/libiir++/build.lib diff --git a/src/libiir++/build.install b/src/libiir++/build.install new file mode 100644 index 0000000..e796001 --- /dev/null +++ b/src/libiir++/build.install @@ -0,0 +1 @@ +source src/libiir++/build.install-lib diff --git a/src/libiir++/build.install-lib b/src/libiir++/build.install-lib new file mode 100644 index 0000000..54af69b --- /dev/null +++ b/src/libiir++/build.install-lib @@ -0,0 +1,35 @@ +build_target libiir++ + +# make paths (this is for Gentoo in particular) +build_dir_tree "${LIBDIR}" || return 1 +build_dir_tree "${BINDIR}" || return 1 +build_dir_tree "${INCLUDEDIR}" || return 1 + +# install library +echo "Installing libraries into '${LIBDIR}'" +install_file ${libiirpp} ${LIBDIR} 0755 || return 1 +BASE="${libiirpp_BASE}.so" +MAJOR="${BASE}.${SOMAJOR}" +MICRO="${MAJOR}.${SOMICRO}" +install_symlink "${BASE}" "${MICRO}" "${LIBDIR}" + +# install header +echo "Installing header file '${libiirpp_HEADER}' into ${INCLUDEDIR}" +install_header ${libiirpp_HEADER} ${INCLUDEDIR} 0644 || return 1 + +# install config script +echo "Installing config script into ${PKGCONFDIR}" +CONFFILE="${INSTALL_PREFIX}${BINDIR}/libiir++-config" + +do_cmd rm -f "${CONFFILE}" +do_cmd_redir "${CONFFILE}" sed \ + -e "s,@VERSION@,${VERSION}," \ + -e "s,@LIB_DIR@,${LIBDIR}," \ + -e "s,@INCLUDE_DIR@,${INCLUDEDIR}," \ + -e "s,@DEP_CFLAGS@,${libiirpp_DEP_CFLAGS}," \ + -e "s,@DEP_LIBS@,${libiirpp_DEP_LIBS}," \ + src/libiir++/config-script +do_cmd chmod 0644 "${CONFFILE}" +print_success "Done" + +# vim: syntax=sh:expandtab:ts=4:sw=4 diff --git a/src/libiir++/build.lib b/src/libiir++/build.lib new file mode 100644 index 0000000..7ece41f --- /dev/null +++ b/src/libiir++/build.lib @@ -0,0 +1,58 @@ +# These are external variables, and shouldn't clash with anything else +# libiirpp +# libiirpp_BUILT +# libiirpp_HEADER +# libiirpp_BASE +# libiirpp_DEP_CFLAGS +# libiirpp_DEP_LIBS + +build_target libiir + +if [ -z ${libiirpp_BUILT} ] +then + libiirpp_BASE=libiir++ + source src/libiir++/soversion + + libiirpp="obj/${libiirpp_BASE}.so.${SOMAJOR}.${SOMICRO}" + libiirpp_DEP_CFLAGS="" + libiirpp_DEP_LIBS="${libiir}" + SO_EXTRA="${libiirpp_DEP_CFLAGS} ${libiirpp_DEP_LIBS} -lstdc++ -lc \ + -D_GNU_SOURCE" + + echo "Building library ${libiirpp}..." + + do_cmd source src/libiir++/build.monolithic || return 1 + + MODIFIED=0 + for test in ${MONOLITHIC_TESTS} ${HDR} ${SRC} + do + if [ ${test} -nt ${libiirpp} ] + then + MODIFIED=1 + break + fi + done + + if [ ${MODIFIED} -ne 0 ] + then + echo " Compiling" + + SONAME="${libiirpp_BASE}.so.${SOMAJOR}" + do_cmd ${CXX} ${CFLAGS} -Iobj -shared -fpic -o "${libiirpp}" \ + -Wl,-soname,${SONAME} \ + ${SRC} ${SO_EXTRA} || return 1 + + # make tests work + do_cmd ln -sf $(basename ${libiirpp}) obj/${SONAME} || return 1 + + print_success "Library built" + else + print_success "Library up to date" + fi + + libiirpp_BUILT=1 + libiirpp_HEADER=${HDR} + +fi + +# vim: syntax=sh:expandtab:ts=4:sw=4 diff --git a/src/libiir++/build.monolithic b/src/libiir++/build.monolithic new file mode 100644 index 0000000..f28385e --- /dev/null +++ b/src/libiir++/build.monolithic @@ -0,0 +1,21 @@ +# These are external variables, and shouldn't clash with anything else +# libiirpp_MONOLITHIC + +SRC="obj/libiir++.cpp" +HDR="obj/iir++.h" + +MONOLITHIC_TESTS="src/libiir++/build.lib src/libiir++/build.monolithic" + +if [ -z "${libiirpp_MONOLITHIC}" ] +then + MONOLITHIC_SOURCE="$(find src/libiir++/ -name '*.h' | sort)" + make_monolithic ${HDR} Ch || return 1 + + MONOLITHIC_SOURCE="$(find src/libiir++/ -name '*.cc' | sort)" + make_monolithic ${SRC} C || return 1 + + libiirpp_MONOLITHIC=1 + MONOLITHIC_DOC="${MONOLITHIC_DOC} ${HDR}" +fi + +# vim: syntax=sh:expandtab:ts=4:sw=4 diff --git a/src/libiir++/config-script b/src/libiir++/config-script new file mode 100644 index 0000000..be6f4e4 --- /dev/null +++ b/src/libiir++/config-script @@ -0,0 +1,97 @@ +#!/bin/bash +# libiir/src/libiir++/config-script +# +# libiir++-config template. Variables are finalised at install time. +# +dep_cflags="@DEP_CFLAGS@" +dep_libs="@DEP_LIBS@" +include_dir="@INCLUDE_DIR@" +include_dir_set="no" +lib_dir="@LIB_DIR@" +lib_dir_set="no" + + + +usage() { + cat <&2 + + + +while [ $# -gt 0 ] +do + case "$1" in + -*=*) + optarg="$(echo "$1" | sed 's/[-_a-zA-Z0-9]*=//')" + ;; + + *) + optarg="" + ;; + esac + + case "$1" in + --libdir=*) + lib_dir="${optarg}" + lib_dir_set="yes" + ;; + + --libdir) + echo_lib_dir="yes" + ;; + + --includedir=*) + include_dir="${optarg}" + include_dir_set="yes" + ;; + + --includedir) + echo_include_dir="yes" + ;; + + --version) + echo "@VERSION@" + exit 0 + ;; + + --cflags) + [ "${include_dir}" != "/usr/include" ] && includes="-I${include_dir}" + echo_cflags="yes" + ;; + + --libs) + echo_libs="yes" + ;; + + *) + usage 1 1>&2 + ;; + esac + + shift +done + + + +[ "${echo_prefix}" == "yes" ] && echo "${prefix}" +[ "${echo_exec_prefix}" == "yes" ] && echo "${exec_prefix}" +[ "${echo_cflags}" == "yes" ] && echo "${dep_cflags} ${includes}" +[ "${echo_libs}" == "yes" ] && echo "${dep_libs} -L${lib_dir} -liir++" +true + + + +# kate: diff --git a/src/libiir++/soversion b/src/libiir++/soversion new file mode 100644 index 0000000..216f81c --- /dev/null +++ b/src/libiir++/soversion @@ -0,0 +1,15 @@ +# libiir/src/libiir++/soversion +# +# Copyright: ©2014, Laurence Withers. +# Author: Laurence Withers +# License: GPLv3 +# + + + +# SOMAJOR is included in the library's soname, and needs to be bumped +# after a binary-incompatible release. It is a single integer. +SOMAJOR=0 + +# SOMICRO is bumped every time there is a binary-compatible release. +SOMICRO=0 diff --git a/src/libiir/000_TopHeader.h b/src/libiir/000_TopHeader.h index 027d784..14f8421 100644 --- a/src/libiir/000_TopHeader.h +++ b/src/libiir/000_TopHeader.h @@ -11,6 +11,12 @@ /* standard includes, or includes needed for type declarations */ #include +#ifdef __cplusplus +extern "C" { +#endif + +/*! \defgroup libiir C library */ + /* options for text editors vim: expandtab:ts=4:sw=4:syntax=c.doxygen */ diff --git a/src/libiir/000_TopSource.c b/src/libiir/000_TopSource.c index ce5ca2d..0cbed21 100644 --- a/src/libiir/000_TopSource.c +++ b/src/libiir/000_TopSource.c @@ -1,6 +1,6 @@ /* libiir/src/libiir/000_TopSource.c * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ @@ -15,6 +15,5 @@ #include /* options for text editors -kate: replace-trailing-space-save true; space-indent true; tab-width 4; vim: expandtab:ts=4:sw=4 */ diff --git a/src/libiir/100_coeff.c b/src/libiir/100_coeff.c new file mode 100644 index 0000000..bbd22bb --- /dev/null +++ b/src/libiir/100_coeff.c @@ -0,0 +1,113 @@ +/* libiir/src/libiir/100_coeff.c + * + * Copyright: ©2010–2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 +*/ + + + +/* struct iir_coeff_t + * Holds a general IIR filter (i.e. the set of coefficients that define it). + * nc >= 1 and nd >= 1. + */ +struct iir_coeff_t { + int nc, nd; + double* c, * d; +}; + + + +/* iir_coeff_new() + * Allocates a new set of coefficient objects. + */ +struct iir_coeff_t* +iir_coeff_new(int nc, const double* c, int nd, const double* d) +{ + struct iir_coeff_t* coeff; + + if(nc < 1 || nd < 1) { + errno = EINVAL; + return 0; + } + + coeff = malloc(sizeof(struct iir_coeff_t)); + coeff->nc = nc; + coeff->nd = nd; + coeff->c = malloc(sizeof(double) * nc); + coeff->d = malloc(sizeof(double) * nd); + memcpy(coeff->c, c, sizeof(double) * nc); + memcpy(coeff->d, d, sizeof(double) * nd); + + return coeff; +} + + + +/* iir_coeff_copy() + * Deep copy of coefficient object. + */ +struct iir_coeff_t* +iir_coeff_copy(const struct iir_coeff_t* other) +{ + struct iir_coeff_t* coeff; + + coeff = malloc(sizeof(struct iir_coeff_t)); + coeff->nc = other->nc; + coeff->nd = other->nd; + coeff->c = malloc(sizeof(double) * coeff->nc); + coeff->d = malloc(sizeof(double) * coeff->nd); + memcpy(coeff->c, other->c, sizeof(double) * coeff->nc); + memcpy(coeff->d, other->d, sizeof(double) * coeff->nd); + + return coeff; +} + + + +/* iir_coeff_free() + * Frees memory associated with ‘coeff’. + */ +void +iir_coeff_free(struct iir_coeff_t* coeff) +{ + if(!coeff) return; + free(coeff->c); + free(coeff->d); + free(coeff); +} + + + +/* iir_coeff_get_*() + * Functions for querying the coefficients that make up an IIR filter. + */ +int +iir_coeff_get_nc(const struct iir_coeff_t* coeff) +{ + return coeff->nc; +} + +int +iir_coeff_get_nd(const struct iir_coeff_t* coeff) +{ + return coeff->nd; +} + +double +iir_coeff_get_c(const struct iir_coeff_t* coeff, int idx) +{ + return coeff->c[idx]; +} + +double +iir_coeff_get_d(const struct iir_coeff_t* coeff, int idx) +{ + return coeff->d[idx]; +} + + + +/* options for text editors +vim: expandtab:ts=4:sw=4 +*/ diff --git a/src/libiir/100_coeff.h b/src/libiir/100_coeff.h new file mode 100644 index 0000000..3e1b493 --- /dev/null +++ b/src/libiir/100_coeff.h @@ -0,0 +1,126 @@ +/* libiir/src/libiir/100_coeff.h + * + * Copyright: ©2010–2014, Laurence Withers. + * Author: Laurence Withers + * License: GPLv3 +*/ + + + +/*! \defgroup libiir_coeff Coefficient handling + +\ingroup libiir + +The functions in this module present an interface for managing the set of +coefficients of a single %IIR filter instance. See \ref iir_structure for +the definition of a coefficient set. The coefficient sets form the +basic building block of an %IIR filter instance (see \ref libiir_filter), +which may chain several %IIR filters in sequence. + +*/ +/*!@{*/ + + + +/* opaque structure */ +struct iir_coeff_t; + + + +/*! \brief Create %IIR coefficient set + +\param nc Number of \a c coefficients. +\param c Array of \a c coefficients. +\param nd Number of \a d coefficients. +\param d Array of \a d coefficients. +\returns Pointer to new general %IIR filter object. + +This function creates a new set of %IIR filter coefficients which may be used to +create filter instances through \ref iir_filter_new() or chained on to existing +instances through \ref iir_filter(). + +See \ref iir_structure for a full explanation of the parameters. + +*/ +struct iir_coeff_t* iir_coeff_new(int nc, const double* c, int nd, + const double* d) +#ifndef DOXYGEN + __attribute__((malloc,nonnull)) +#endif +; + + + +/*! \brief Create copy of a set of %IIR coefficients + +\param other Set of %IIR filter coefficients to copy. +\returns Pointer to new general %IIR filter object. + +Performs a deep copy of the set of %IIR coefficients contained within \a other. + +*/ +struct iir_coeff_t* iir_coeff_copy(const struct iir_coeff_t* other); + + + +/*! \brief Free %IIR coefficient set + +\param coeff Pointer to %IIR coefficient object. May be 0. + +Frees a set of %IIR filter coefficients previously allocated through +\ref iir_coeff_new(). Can be called on a null pointer without consequences. +Note that \ref iir_filter_new() and \ref iir_filter_chain() actually store a +copy of the coefficients, so it is possible to free \a coeff even if existing +filters are still using its coefficient values. + +*/ +void iir_coeff_free(struct iir_coeff_t* coeff); + + + +/*! \brief Query number of C coefficients. + +\param coeff Pointer to %IIR coefficient object. +\returns Number of C coefficients (≥1). + +*/ +int iir_coeff_get_nc(const struct iir_coeff_t* coeff); + + + +/*! \brief Query number of D coefficients. + +\param coeff Pointer to %IIR coefficient object. +\returns Number of D coefficients (≥1). + +*/ +int iir_coeff_get_nd(const struct iir_coeff_t* coeff); + + + +/*! \brief Get value of C coefficient. + +\param coeff Pointer to %IIR coefficient object. +\param idx Index of coefficient (≥0, < \ref iir_coeff_get_nc()). +\returns C coefficient value. + +*/ +double iir_coeff_get_c(const struct iir_coeff_t* coeff, int idx); + + + +/*! \brief Get value of D coefficient. + +\param coeff Pointer to %IIR coefficient object. +\param idx Index of coefficient (≥0, < \ref iir_coeff_get_nd()). +\returns D coefficient value. + +*/ +double iir_coeff_get_d(const struct iir_coeff_t* coeff, int idx); + + + +/*!@}*/ +/* options for text editors +vim: expandtab:ts=4:sw=4:syntax=c.doxygen +*/ diff --git a/src/libiir/200_iir.c b/src/libiir/200_filter.c similarity index 83% rename from src/libiir/200_iir.c rename to src/libiir/200_filter.c index 7ab4625..3fb7f2d 100644 --- a/src/libiir/200_iir.c +++ b/src/libiir/200_filter.c @@ -1,63 +1,12 @@ -/* libiir/src/libiir/200_iir.c +/* libiir/src/libiir/200_filter.c * - * Copyright: ©2010–2011, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ -/* struct iir_coeff_t - * Holds a general IIR filter (i.e. the set of coefficients that define it). - * nc >= 1 and nd >= 1. - */ -struct iir_coeff_t { - int nc, nd; - double* c, * d; -}; - - - -/* iir_coeff_new() - * Allocates a new set of coefficient objects. - */ -struct iir_coeff_t* -iir_coeff_new(int nc, double* c, int nd, double* d) -{ - struct iir_coeff_t* coeff; - - if(nc < 1 || nd < 1) { - errno = EINVAL; - return 0; - } - - coeff = malloc(sizeof(struct iir_coeff_t)); - coeff->nc = nc; - coeff->nd = nd; - coeff->c = malloc(sizeof(double) * nc); - coeff->d = malloc(sizeof(double) * nd); - memcpy(coeff->c, c, sizeof(double) * nc); - memcpy(coeff->d, d, sizeof(double) * nd); - - return coeff; -} - - - -/* iir_coeff_free() - * Frees memory associated with ‘coeff’. - */ -void -iir_coeff_free(struct iir_coeff_t* coeff) -{ - if(!coeff) return; - free(coeff->c); - free(coeff->d); - free(coeff); -} - - - /* struct iir_filter_t * An instantiated IIR filter. This is actually a linked list node, so that we * can create chains of filters. It also has a copy of the coefficients so that @@ -268,6 +217,25 @@ iir_filter(struct iir_filter_t* fi, double samp) +int +iir_filter_coeff_sets(const struct iir_filter_t* fi) +{ + int count; + for(count = 0; fi; ++count) fi = fi->next; + return count; +} + + + +struct iir_coeff_t* +iir_filter_coeff_set(const struct iir_filter_t* fi, int idx) +{ + while(idx--) fi = fi->next; + return iir_coeff_new(fi->nc, fi->c, fi->nd, fi->d); +} + + + /* options for text editors vim: expandtab:ts=4:sw=4 */ diff --git a/src/libiir/200_iir.h b/src/libiir/200_filter.h similarity index 55% rename from src/libiir/200_iir.h rename to src/libiir/200_filter.h index 309b5fd..eea7df9 100644 --- a/src/libiir/200_iir.h +++ b/src/libiir/200_filter.h @@ -1,23 +1,25 @@ -/* libiir/src/libiir/200_iir.h +/* libiir/src/libiir/200_filter.h * - * Copyright: ©2010–2011, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ -/*! \defgroup iir Basic IIR filtering +/*! \defgroup libiir_filter Basic IIR filtering -The functions in this module present a basic interface for representing IIR -filters, creating instances of (possibly chained) IIR filters, and filtering an +\ingroup libiir + +The functions in this module present a basic interface for representing %IIR +filters, creating instances of (possibly chained) %IIR filters, and filtering an input sample. -A general IIR filter consists of a set of coefficients, and may be created +A general %IIR filter consists of a set of coefficients, and may be created through \ref iir_coeff_new(). The filter object (the opaque struct iir_coeff_t) is then used to instantiate specific filters (the opaque struct iir_filter_t) through \ref iir_filter_new(). Each filter -instance may have an arbitrary further number of IIR filters chained on to it +instance may have an arbitrary further number of %IIR filters chained on to it through \ref iir_filter_chain(). The filter processes one sample at a time through \ref iir_filter(). @@ -26,60 +28,17 @@ through \ref iir_filter(). -/* opaque structure */ -struct iir_coeff_t; - - - -/*! \brief Create general IIR filter - -\param nc Number of \a c coefficients. -\param c Array of \a c coefficients. -\param nd Number of \a d coefficients. -\param d Array of \a d coefficients. -\returns Pointer to new general IIR filter object. - -This function creates a new general IIR filter object which may be used to -create filter instances through \ref iir_filter_new() or chained on to existing -instances through \ref iir_filter(). - -See \ref iir_structure for a full explanation of the parameters. - -*/ -struct iir_coeff_t* iir_coeff_new(int nc, double* c, int nd, double* d) -#ifndef DOXYGEN - __attribute__((malloc,nonnull)) -#endif -; - - - -/*! \brief Free general IIR filter - -\param coeff Pointer to IIR filter object. May be 0. - -Frees a set of IIR filter coefficients previously allocated through -\ref iir_coeff_new(). Can be called on a null pointer without consequences. -Note that \ref iir_filter_new() and \ref iir_filter_chain() actually store a -copy of the coefficients, so it is possible to free \a coeff even if existing -filters are still using its coefficient values. - -*/ -void iir_coeff_free(struct iir_coeff_t* coeff); - - - /* opaque structure */ struct iir_filter_t; -/*! \brief Create IIR filter instance +/*! \brief Create %IIR filter instance \param coeff Filter coefficients to use. -\returns Pointer to new instance of IIR filter. +\returns Pointer to new instance of %IIR filter. -Creates a new instance of a general IIR filter. The set of coefficients \a coeff +Creates a new instance of a general %IIR filter. The set of coefficients \a coeff is copied into the returned structure, meaning the coefficients can be freed after this function returns if they will not be needed again. @@ -98,11 +57,11 @@ struct iir_filter_t* iir_filter_new(const struct iir_coeff_t* coeff) -/*! \brief Free IIR filter instance +/*! \brief Free %IIR filter instance \param fi Filter object to free. May be 0. -Frees a previously-allocated IIR filter instance. Can be called on a null +Frees a previously-allocated %IIR filter instance. Can be called on a null pointer without consequences. */ @@ -110,12 +69,12 @@ void iir_filter_free(struct iir_filter_t* fi); -/*! \brief Add a further IIR filter to a filter instance +/*! \brief Add a further %IIR filter to a filter instance \param fi Filter instance to chain onto. -\param coeff New IIR filter coefficients to add to chain. +\param coeff New %IIR filter coefficients to add to chain. -Extends an existing IIR filter by chaining a new set of coefficients onto the +Extends an existing %IIR filter by chaining a new set of coefficients onto the end. This can be used for >4th order Butterworth filters, for example. This copies the set of coefficients from \a coeff so the coefficients can be freed after this function returns if they are no longer required. @@ -129,7 +88,7 @@ void iir_filter_chain(struct iir_filter_t* fi, const struct iir_coeff_t* coeff) -/*! \brief Create a deep copy of an IIR filter instance +/*! \brief Create a deep copy of an %IIR filter instance \param fi Filter instance to copy. \param state Non-zero to copy state as well. @@ -162,6 +121,35 @@ double iir_filter(struct iir_filter_t* fi, double samp); +/*! \brief Count number of coefficient sets in %IIR filter chain + +\param fi Filter object. +\returns Number of coefficient sets in chain (≥1). + +Returns the number of discrete %IIR filter coefficient sets used in the filter +object \a fi. There will always be at least one set. + +*/ +int iir_filter_coeff_sets(const struct iir_filter_t* fi); + + + +/*! \brief Get %IIR coefficient set from filter chain + +\param fi Filter object. +\param idx Index of coefficient set (0 ≤ \a idx < \ref iir_filter_coeff_sets()). +\returns Newly-allocated coefficient set object. + +Extracts the %IIR filter coefficient set for the given step (\a idx) in the chain +of filters in \a fi. Returns a newly-allocated filter coefficient set object +which can be freed with \ref iir_coeff_free(). + +*/ +struct iir_coeff_t* iir_filter_coeff_set(const struct iir_filter_t* fi, + int idx); + + + /*!@}*/ /* options for text editors vim: expandtab:ts=4:sw=4:syntax=c.doxygen diff --git a/src/libiir/300_common_filters.h b/src/libiir/300_common_filters.h index 897827b..dc0db48 100644 --- a/src/libiir/300_common_filters.h +++ b/src/libiir/300_common_filters.h @@ -1,16 +1,18 @@ /* libiir/src/libiir/300_common_filters.h * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ -/*! \defgroup common_filters Common types of IIR filter +/*! \defgroup libiir_common_filters Common types of IIR filter -Functions to create coefficients for various common types of IIR filter. The -coefficient structures which are returned may be used to instantiate IIR +\ingroup libiir + +Functions to create coefficients for various common types of %IIR filter. The +coefficient structures which are returned may be used to instantiate %IIR filters using \ref iir_filter_new(). The Butterworth filter code comes from the Exstrom Labs LLC code available under @@ -28,7 +30,7 @@ is a copy of the original code available in the top level of this project. \param gain Linear gain of filter. \param corner Corner frequency expressed as a fraction of Nyquist (0 ≤ \a corner ≤ 1) -\returns Newly-allocated IIR filter coefficients. +\returns Newly-allocated %IIR filter coefficients. Uses the Exstrom labs code to compute the coefficients of an nth-order (param \a order) Butterworth-type low pass filter with gain \a gain and corner @@ -39,7 +41,7 @@ greater than a 4th-order filter. This function won't do that directly for you. \a gain will usually be set to be 1.0. The corner frequency \a corner is expressed as a fraction of the sampling -frequency (which is of course not known by the IIR code). It should lie between +frequency (which is of course not known by the %IIR code). It should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). */ @@ -59,7 +61,7 @@ struct iir_coeff_t* iir_butterworth_lowpass(int order, \param gain Linear gain of filter. \param corner Corner frequency expressed as a fraction of Nyquist (0 ≤ \a corner ≤ 1) -\returns Newly-allocated IIR filter coefficients. +\returns Newly-allocated %IIR filter coefficients. Uses the Exstrom labs code to compute the coefficients of an nth-order (param \a order) Butterworth-type high pass filter with gain \a gain and corner @@ -70,7 +72,7 @@ greater than a 4th-order filter. This function won't do that directly for you. \a gain will usually be set to be 1.0. The corner frequency \a corner is expressed as a fraction of the sampling -frequency (which is of course not known by the IIR code). It should lie between +frequency (which is of course not known by the %IIR code). It should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency). */ @@ -92,7 +94,7 @@ struct iir_coeff_t* iir_butterworth_highpass(int order, (0 ≤ \a c1 ≤ 1) \param c2 High corner frequency expressed as a fraction of Nyquist (0 ≤ \a c2 ≤ 1, and \a c1 < \a c2) -\returns Newly-allocated IIR filter coefficients. +\returns Newly-allocated %IIR filter coefficients. Uses the Exstrom labs code to compute the coefficients of an nth-order (param \a order) Butterworth-type band pass filter with gain \a gain and corner @@ -103,7 +105,7 @@ greater than a 4th-order filter. This function won't do that directly for you. \a gain will usually be set to be 1.0. The corner frequencies \a c1 and \a c2 are expressed as a fraction of the -sampling frequency (which is of course not known by the IIR code). They should +sampling frequency (which is of course not known by the %IIR code). They should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency), and \a c2 should be greater than \a c1. @@ -127,7 +129,7 @@ struct iir_coeff_t* iir_butterworth_bandpass(int order, (0 ≤ \a c1 ≤ 1) \param c2 High corner frequency expressed as a fraction of Nyquist (0 ≤ \a c2 ≤ 1, and \a c1 < \a c2) -\returns Newly-allocated IIR filter coefficients. +\returns Newly-allocated %IIR filter coefficients. Uses the Exstrom labs code to compute the coefficients of an nth-order (param \a order) Butterworth-type band stop filter with gain \a gain and corner @@ -138,7 +140,7 @@ greater than a 4th-order filter. This function won't do that directly for you. \a gain will usually be set to be 1.0. The corner frequencies \a c1 and \a c2 are expressed as a fraction of the -sampling frequency (which is of course not known by the IIR code). They should +sampling frequency (which is of course not known by the %IIR code). They should lie between 0 (0Hz) and 1 (the Nyquist frequency, or ½ the sampling frequency), and \a c2 should be greater than \a c1. @@ -156,6 +158,5 @@ struct iir_coeff_t* iir_butterworth_bandstop(int order, /*!@}*/ /* options for text editors -kate: replace-trailing-space-save true; space-indent true; tab-width 4; vim: expandtab:ts=4:sw=4:syntax=c.doxygen */ diff --git a/src/libiir/400_parser.c b/src/libiir/400_parser.c index 107f51e..604aed6 100644 --- a/src/libiir/400_parser.c +++ b/src/libiir/400_parser.c @@ -103,8 +103,13 @@ iir_parser_raw(struct iir_filter_t** fi, const char* desc) * The maximum order in any single Butterworth-type filter. If the given order * exceeds this, we split the resulting filter up into multiple sets of * coefficients of this order or less. + * + * This is required since IIR filters, even with double-precision floating + * point numbers, are particularly susceptible to quantization errors. In + * practice, anything with a steeper roll-off than 2nd order can go unstable + * when presented with real input data. */ -#define IIR_PARSER_BW_MAX_ORDER (4) +#define IIR_PARSER_BW_MAX_ORDER (2) @@ -176,7 +181,7 @@ iir_parser_bw_aux2(struct iir_filter_t** fi, } iir_coeff_free(coeff); if(!order) return 0; - /* add a <4th order segment */ + /* add a <2nd order segment */ } coeff = bw(order, gain, corner); @@ -214,7 +219,7 @@ iir_parser_bw_aux3(struct iir_filter_t** fi, } iir_coeff_free(coeff); if(!order) return 0; - /* add a <4th order segment */ + /* add a <2nd order segment */ } coeff = bw(order, gain, c1, c2); diff --git a/src/libiir/400_parser.h b/src/libiir/400_parser.h index 0c0c099..ea43539 100644 --- a/src/libiir/400_parser.h +++ b/src/libiir/400_parser.h @@ -1,17 +1,19 @@ /* libiir/src/libiir/400_parser.h * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ -/*! \defgroup parser Parser for user-specified IIR filters +/*! \defgroup iir_parser Parser for user-specified IIR filters -This is a high-level interface that can instantiate a set of IIR filters based +\ingroup libiir + +This is a high-level interface that can instantiate a set of %IIR filters based on a user-specified, human-readable string. The intention of this interface is -to allow IIR filters to be specified in configuration files so that they can be +to allow %IIR filters to be specified in configuration files so that they can be easily modified by the user and easily understood/parsed by the system. See \ref string_desc for details on the string description format. @@ -21,14 +23,14 @@ See \ref string_desc for details on the string description format. -/*! \brief Instantiate an IIR filter based on a string description +/*! \brief Instantiate an %IIR filter based on a string description -\param desc IIR filter description. -\returns Pointer to newly-allocated IIR filter instance. +\param desc %IIR filter description. +\returns Pointer to newly-allocated %IIR filter instance. \retval 0 on error. -Parses the human-readable description of an IIR filter chain in \a desc, -instantiating an IIR filter object to match. Returns the new filter. If \a desc +Parses the human-readable description of an %IIR filter chain in \a desc, +instantiating an %IIR filter object to match. Returns the new filter. If \a desc cannot be parsed correctly, returns 0 and sets \a errno to \c EINVAL. */ @@ -40,15 +42,15 @@ struct iir_filter_t* iir_parse(const char* desc) -/*! \brief Instantiate a set of IIR filters based on a string description +/*! \brief Instantiate a set of %IIR filters based on a string description -\param desc IIR filter description. +\param desc %IIR filter description. \param n Number of instances to allocate. -\returns Pointer to array of \a n newly-allocated IIR filter instances. +\returns Pointer to array of \a n newly-allocated %IIR filter instances. \retval 0 on error. -Parses the human-readable description of an IIR filter chain in \a desc, -instantiating a set of \a n identical IIR filter objects to match. Returns a +Parses the human-readable description of an %IIR filter chain in \a desc, +instantiating a set of \a n identical %IIR filter objects to match. Returns a pointer to an array of new filters. If \a desc cannot be parsed correctly, returns 0 and sets \a errno to \c EINVAL. @@ -66,6 +68,5 @@ struct iir_filter_t** iir_parse_n(const char* desc, int n) /*!@}*/ /* options for text editors -kate: replace-trailing-space-save true; space-indent true; tab-width 4; vim: expandtab:ts=4:sw=4:syntax=c.doxygen */ diff --git a/src/libiir/500_response.c b/src/libiir/500_response.c index b19ee9b..8c6ad0d 100644 --- a/src/libiir/500_response.c +++ b/src/libiir/500_response.c @@ -66,6 +66,62 @@ iir_response(const struct iir_filter_t* fi, double freq) +void +iir_response_minus_3dB_aux(const struct iir_filter_t* fi, double gain, + double** freqs, int* num_freqs, int* sz_freqs, double f1, double f2) +{ + double f3; + double mag; + int i; + + for(i = 0; i < 20; ++i) { + f3 = (f1 + f2) / 2; + mag = cabs(iir_response(fi, f3)) / gain; + if(mag < 0.5) { + f1 = f3; + } else { + f2 = f3; + } + } + + if(*num_freqs == *sz_freqs) { + *sz_freqs = (*sz_freqs) ? (*sz_freqs << 1) : 4; + *freqs = realloc(*freqs, sizeof(double) * (*sz_freqs)); + } + (*freqs)[*num_freqs] = f3; + ++(*num_freqs); +} + +void +iir_response_minus_3dB(const struct iir_filter_t* fi, double gain, + double** freqs, int* num_freqs) +{ + int sz_freqs = 0, pos; + double mag, last_mag; + + *freqs = 0; + *num_freqs = 0; + + last_mag = cabs(iir_response(fi, 0)) / gain; + for(pos = 1; pos < 1000; ++pos) { + mag = cabs(iir_response(fi, 0.0005 * pos)) / gain; + if(mag < 0.5) { + if(last_mag >= 0.5) { + iir_response_minus_3dB_aux(fi, gain, freqs, num_freqs, + &sz_freqs, 0.0005 * pos, 0.0005 * (pos - 1)); + } + } else { + if(last_mag < 0.5) { + iir_response_minus_3dB_aux(fi, gain, freqs, num_freqs, + &sz_freqs, 0.0005 * (pos - 1), 0.0005 * pos); + } + } + last_mag = mag; + } +} + + + /* options for text editors vim: expandtab:ts=4:sw=4 */ diff --git a/src/libiir/500_response.h b/src/libiir/500_response.h index 5c4e62a..885e1bf 100644 --- a/src/libiir/500_response.h +++ b/src/libiir/500_response.h @@ -1,15 +1,17 @@ /* libiir/src/libiir/500_response.h * - * Copyright: ©2011, Laurence Withers. + * Copyright: ©2011–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ -/*! \defgroup response Measuring IIR response +/*! \defgroup iir_response Measuring IIR response -These functions allow measuring the IIR response (i.e. frequency and phase of +\ingroup libiir + +These functions allow measuring the %IIR response (i.e. frequency and phase of the transfer function) at a given frequency. The response is a complex number; see \c complex(5) for details. The magnitude \c cabs(z) is the gain of the filter and the phase \c carg(z) its phase delay. @@ -46,6 +48,29 @@ double complex iir_response(const struct iir_filter_t* fi, double freq); +/*! \brief Search for -3dB points. + +\param fi Filter to evaluate. +\param gain Nominal gain of the filter. +\param[out] freqs Set to point to array of frequencies +\param[out] num_freqs Number of frequencies in \a freq. + +This function will perform a linear scan of the frequency response of the +filter between DC and the Nyquist frequency. It will then search for any pair +of points between which the magnitude transitions the -3dB level. It then +performs a few levels of binary search in order to provide a more accurate +answer. + +Returned frequencies are expressed as a fraction of the sampling rate. The +function will allocate an array of \c double with \c malloc(3), so the result +should later be passed to \c free(3). + +*/ +void iir_response_minus_3dB(const struct iir_filter_t* fi, double gain, + double** freqs, int* num_freqs); + + + /*!@}*/ /* options for text editors vim: expandtab:ts=4:sw=4:syntax=c.doxygen diff --git a/src/libiir/999_BottomHeader.h b/src/libiir/999_BottomHeader.h index 731ec5a..813fb2f 100644 --- a/src/libiir/999_BottomHeader.h +++ b/src/libiir/999_BottomHeader.h @@ -1,13 +1,16 @@ /* libiir/src/libiir/999_BottomHeader.h * - * Copyright: ©2010, Laurence Withers. + * Copyright: ©2010–2014, Laurence Withers. * Author: Laurence Withers * License: GPLv3 */ +#ifdef __cplusplus +} +#endif + #endif /* options for text editors -kate: replace-trailing-space-save true; space-indent true; tab-width 4; vim: expandtab:ts=4:sw=4:syntax=c.doxygen */ diff --git a/src/libiir/soversion b/src/libiir/soversion index c779093..09c75fa 100644 --- a/src/libiir/soversion +++ b/src/libiir/soversion @@ -12,4 +12,4 @@ SOMAJOR=0 # SOMICRO is bumped every time there is a binary-compatible release. -SOMICRO=1 +SOMICRO=2 diff --git a/version b/version index ac1edc9..6e4039f 100644 --- a/version +++ b/version @@ -10,8 +10,8 @@ # VERSION contains the full version number of the library, which is # expected to be in 'major.minor.micro' format. VERMAJOR=1 -VERMINOR=0 -VERMICRO=1 +VERMINOR=1 +VERMICRO=0 # kate: replace-trailing-space-save true; space-indent true; tab-width 4; # vim: expandtab:ts=4:sw=4:syntax=sh