{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Realization of Recursive Filters\n",
    "\n",
    "*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct Form Structures\n",
    "\n",
    "The output signal $y[k] = \\mathcal{H} \\{ x[k] \\}$ of a recursive linear-time invariant (LTI) system is given by\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = \\frac{1}{a_0} \\left( \\sum_{m=0}^{M} b_m \\; x[k-m] - \\sum_{n=1}^{N} a_n \\; y[k-n] \\right)\n",
    "\\end{equation}\n",
    "\n",
    "where $a_n$ and $b_m$ denote constant coefficients and $N$ the order. Note that systems with $M > N$ are in general not stable. The computational realization of above equation requires additions, multiplications, the actual and past samples of the input signal $x[k]$, and the past samples of the output signal $y[k]$. Technically this can be realized by\n",
    "\n",
    "* adders\n",
    "* multipliers, and\n",
    "* unit delays or storage elements.\n",
    "\n",
    "These can be arranged in different topologies. A certain class of structures, which is introduced in the following, is known as *direct form structures*. Other known forms are for instance [cascaded sections](cascaded_structures.ipynb), parallel sections, lattice structures and state-space structures.\n",
    "\n",
    "For the following it is assumed that $a_0 = 1$. This can be achieved for instance by dividing the remaining coefficients by $a_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Direct Form I\n",
    "\n",
    "The [direct form I](https://en.wikipedia.org/wiki/Digital_filter#Direct_Form_I) is derived by rearranging above equation for $a_0 = 1$\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = \\sum_{m=0}^{M} b_m \\; x[k-m] + \\sum_{n=1}^{N} - a_n \\; y[k-n]\n",
    "\\end{equation}\n",
    "\n",
    "It is now evident that we can realize the recursive filter by a superposition of a non-recursive and a recursive part. With the elements given above, this results in the following block-diagram\n",
    "\n",
    "![Direct form I filter](direct_form_i.png)\n",
    "\n",
    "This representation is not canonical since $N + M$ unit delays are required to realize a system of order $N$. A benefit of the direct form I is that there is essentially only one summation point which has to be taken care of when considering quantized variables and overflow. The output signal $y[k]$ for the direct form I is computed by realizing above equation.\n",
    "\n",
    "The block diagram of the direct form I can be interpreted as the cascade of two systems. Denoting the signal in between both as $w[k]$ and discarding initial values we get\n",
    "\n",
    "\\begin{align}\n",
    "w[k] &= \\sum_{m=0}^{M} b_m \\; x[k-m] = h_1[k] * x[k] \\\\\n",
    "y[k] &= w[k] + \\sum_{n=1}^{N} - a_n \\; w[k-n] = h_2[k] * w[k] = h_2[k] * h_1[k] * x[k]\n",
    "\\end{align}\n",
    "\n",
    "where $h_1[k] = [b_0, b_1, \\dots, b_M]$ denotes the impulse response of the non-recursive part and $h_2[k] = [1, -a_1, \\dots, -a_N]$ for the recursive part. From the last equality of the second equation and the commutativity of the convolution it becomes clear that the order of the cascade can be exchanged."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Direct Form II\n",
    "\n",
    "The [direct form II](https://en.wikipedia.org/wiki/Digital_filter#Direct_Form_II) is yielded by exchanging the two systems in above block diagram and noticing that there are two parallel columns of delays which can be combined, since they are redundant. For $N=M$ it is given as\n",
    "\n",
    "![Direct form II filter](direct_form_ii.png)\n",
    "\n",
    "Other cases with $N \\neq M$ can be considered for by setting coefficients to zero. This form is a canonical structure since it only requires $N$ unit delays for a recursive filter of order $N$. The output signal $y[k]$ for the direct form II is computed by the following equations\n",
    "\n",
    "\\begin{align}\n",
    "w[k] &= x[k] + \\sum_{n=1}^{N} - a_n \\; w[k-n] \\\\\n",
    "y[k] &= \\sum_{m=0}^{M} b_m \\; w[k-m]\n",
    "\\end{align}\n",
    "\n",
    "The samples $w[k-m]$ are termed *state* (variables) of a digital filter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transposed Direct Form II\n",
    "\n",
    "The block diagrams above can be interpreted as linear [signal flow graphs](https://en.wikipedia.org/wiki/Signal-flow_graph). The theory of these graphs provides useful transformations into different forms which preserve the overall transfer function. Of special interest is the *transposition* or *reversal* of a graph which can be achieved by\n",
    "\n",
    "* exchanging in- and output,\n",
    "* exchanging signal split and summation points, and\n",
    "* reversing the directions of the signal flows.\n",
    "\n",
    "Applying this procedure to the direct form II shown above for $N=M$ yields the transposed direct form II\n",
    "\n",
    "![Transposed direct form II filter](direct_form_ii_t.png)\n",
    "\n",
    "The output signal of the transposed direct form II is given as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = b_0 x[k] + \\sum_{m=1}^{M} b_m x[k-n] - \\sum_{n=1}^{N} a_n y[k-n]\n",
    "\\end{equation}\n",
    "\n",
    "Using the signal before the $n$-th delay unit as internal state $w_n[k]$ we can reformulate this into a set of difference equations for computation of the output signal\n",
    "\n",
    "\\begin{align}\n",
    "w_n[k] &= \n",
    "\\begin{cases}\n",
    "w_{n+1}[k-1] - a_n y[k] + b_n x[k] & \\text{for } n=0,1,\\dots,N-1 \\\\\n",
    "-a_N y[k]  + b_N x[k] & \\text{for } n=N\n",
    "\\end{cases}\\\\\n",
    "y[k] &= w_1[k-1] + b_0 x[k]\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Computation of output signal\n",
    "\n",
    "The following example illustrates the computation of the impulse response $h[k]$ of a 2nd-order recursive system using the transposed direct form II as realized by `scipy.signal.lfilter`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAEYCAYAAAAj/u7rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2UXWV96PHvjyTgWCopoIEMWLHS1FiE1BRFurwDtQ2+XElpWULfsK0NvUv7tmqE1K6+2FqwtLXt1aVSq6BWopeLKbewjEo6tVVEwEEC0igCQiYIJXjExJGQmd/94+zRk8m8nD05Z87Ze76ftWbl7H2es/fv2Xsy+3ee59nPjsxEkiQtbof1OgBJktR7JgSSJMmEQJIkmRBIkiRMCCRJEiYEkiQJEwJp0YuI10XEf/Y6Dkm9ZUIg9YmIeCAiXt7rOCQtTiYEkrouIpb2OgZJszMhkPpQ0Yz/2Yh4R0Q0IuK+iHhpsf6hiHg0Ii5qKX9VRLwnIj4VEd+OiH+PiB8u3ntORGTrRTkihiPi9dPsN4p9PhoR34qIOyPix4v3joiIv46IByPikWJ/A23E/zjwp8X6X4+IeyLimxGxtSXG2fY7Y92K918aEbcWn7s1Il46pZ5/XsTy7Yj4ZEQcW7z3tIj4cETsLo7xrRGxonjvqIj4p4h4OCJGI+IvImLJ/M+o1P9MCKT+9WLgTuAY4CPAZuAngecBvwy8MyKObCn/S8CfA8cCdwD/PI99/izwMuBHgeXAa4HdxXtvL9afVsQwCPzxHPHfBzwLeFtErAf+EDgPeCbwH8A1bex3xrpFxNHADcA/0DxOfwvcEBHHtHz2F4FfK+I4HHhTsf4i4CjgxOKzvwWMFe9dDewv6rmmiO+gBEqqExMCqX/dn5kfyMxx4KM0L1xvzcwnM/OTwD6aF6xJN2TmZzLzSeAtwBkRcWLJfT4F/CDwY0Bk5j2Z+XBEBPCbwO9n5uOZ+W3gL4ELZtnWrsz835m5PzPHgIuBy4pt7i8+f1rxbX/a/bZRt1cBX83MDxX7uQb4L+B/tnz2A5n5lSKGj9FMaCbregzwvMwcz8zbM/OJopXgFcDvZebezHwUeMccdZUqz4RA6l+PtLweA8jMqetaWwgemnyRmXuAx4GVZXaYmduAdwLvAh6JiCsj4hk0v9E/Hbi9aF5vAJ8o1s/koSnLPwz8fcvnHwcCGJxlv3PVbSXw9Sn7+TrN1otJ32h5/R2+f8w+BGwFNkfEroj4q4hYVsS5DHi4Jdb30mxhkGrLhECqj++1BhRdCUcDu4C9xeqnt5Q9bqaNZOY/ZOaLgBfQbMLfCDxGMwF5QWYuL36OyswjZ9oOMPVRqg8BF7d8fnlmDmTm52bZ71x120XzAt7q2cDoLHFN1vOpzPyzzFwNvBR4NfCrRZxPAse2xPmMzHzBXNuUqsyEQKqPV0bET0XE4TT722/JzIcy879pXiB/OSKWRMSvAz8y3QYi4icj4sXFN+W9wHeB8cycAP4ReEdEPKsoOxgR60rE9x5gU0S8oPj8URFx/mz7natuwI3Aj0bEL0bE0oh4LbAa+Ne5gomIsyLilGKw4BM0uxDGi66KTwJ/ExHPiIjDIuJHIuJ/lKirVDkmBFJ9fAT4E5rN6S+iORBv0m/S/Ma9m+Y38M/NsI1n0Lzwf5Nm0/tu4K+L9y4B7gU+HxFPAJ8GVrUbXGZ+nObAxM3F5++i2Vc/135nrFtm7qb5zf4Pis+8GXh1Zj7WRkjHAdfSTAbuAf4d+HDx3q/SHID45SKma4Hj262rVEWRObVVT1LVRMRVwM7M/KNex9Jpda6b1E9sIZAkSf2VEETE+4uJSe6a4f2hYvKRO4qf2e6BliRJbeqrLoOIeBmwB/hgZv74NO8PAW/KzFcvdGySJNVZX7UQZOZnaA4akiRJC6iKDxw5IyK+RPP+4zdl5t1TC0TEBmADwMDAwItOPLG9ydomJiY47LC+ypE6pq51q2u9oL51s17VU9e61bVeAF/5ylcey8zZJg47SNUSgi8CP5yZeyLilcAW4OSphTLzSuBKgLVr1+Ztt93W1saHh4cZGhrqXLR9pK51q2u9oL51s17VU9e61bVeABExdQbPOVUqNcrMJ4ppS8nMG4Flk08ukyRJ81ephCAijiseskJEnE4z/t2zf0qSJM2lr7oMIuIaYAg4NiJ20pyZbBlAZr4H+AXgf0XEfprzql+Q/XSbhCRJFdVXCUFmXjjH+++k+UQ0SZLUQZXqMpAkSd1hQiBJkkwIJEmSCYEkScKEQJIkYUIgSZIwIZAkSZgQSJIk+mxioirYMjLKFVt3sKsxxsrlA2xct4r1awZ7HZYkSYfEhKCELSOjbLpuO2NPjQMw2hhj03XbAUwKJEmVZpdBCVds3fG9ZGDS2FPjXLF1R48ikiSpM0wIStjVGCu1XpKkqjAhKGHl8oFS6yVJqgoTghI2rlvFwLIlB6wbWLaEjetW9SgiSZI6w0GFJUwOHHzztXeyb3yCQe8ykCTVhAlBSevXDHLNFx4E4KMXn9HjaCRJ6gy7DCRJkgmBJEkyIZAkSZgQSJIkTAgkSRImBJIkCRMCSZKECYEkScKEQJIkYUIgSZIwIZAkSfRZQhAR74+IRyPirhnej4j4h4i4NyLujIifWOgYJUmqo75KCICrgHNmef8VwMnFzwbg3Z3Y6ZaRUc68fBuv+8Rezrx8G1tGRjuxWUmSKqOvEoLM/Azw+CxFzgU+mE2fB5ZHxPGHss8tI6Nsum47o40xAEYbY2y6brtJgSRpUemrhKANg8BDLcs7i3XzdsXWHYw9NX7AurGnxrli645D2awkSZWytNcBlBTTrMuDCkVsoNmlwIoVKxgeHp5xg5MtA9Otn+lzjeIzs2233+zZs6dS8barrvWC+tbNelVPXetW13rNV9USgp3AiS3LJwC7phbKzCuBKwHWrl2bQ0NDM25w8PPbpk0KBpcPMNPn3r3jZgCGhs5oN+6eGx4enrE+VVbXekF962a9qqeudatrvearal0G1wO/Wtxt8BLgW5n58KFscOO6VQwsW3LAuoFlS9i4btWhbFaSpErpqxaCiLgGGAKOjYidwJ8AywAy8z3AjcArgXuB7wC/dqj7XL+mOQThzdfeyb7xCQaXD7Bx3arvrZckaTHoq4QgMy+c4/0E3tDp/a5fM8g1X3iQRqPB1kvO7vTmJUnqe1XrMpAkSV1gQiBJkkwIJEmSCYEkScKEQJIk0Wd3GdTNlpFRrti6g12NMVZ6O6MkqY+ZEHTJ5EOTJp+TMPnQJMCkQJLUd+wy6BIfmiRJqhITgi7ZNcNDk2ZaL0lSL5kQdMnK5QOl1kuS1EsmBF3iQ5MkSVXioMIu8aFJkqQqMSHoosmHJgF89OIzehyNJEkzs8tAkiSZEEiSJBMCSZKECYEkScKEQJIkYUIgSZIwIZAkSZgQSJIkTAgkSRImBJIkCRMCSZKEzzJQH9kyMsoVW3ewqzHGSh8GJUkLyoRAfWHLyCibrtvO2FPjAIw2xth03XYAkwJJWgB2GagvXLF1x/eSgUljT41zxdYdPYpIkhYXEwL1hV2NsVLrJUmd1VcJQUScExE7IuLeiLh0mvdfFxH/HRF3FD+v70Wc6ryVywdKrZckdVbfJAQRsQR4F/AKYDVwYUSsnqboRzPztOLnfQsapLpm47pVDCxbcsC6gWVL2LhuVY8ikqTFpW8SAuB04N7MvC8z9wGbgXN7HJMWyPo1g1x23ikcvqT5Kzm4fIDLzjvFAYWStED66S6DQeChluWdwIunKffzEfEy4CvA72fmQ1MLRMQGYAPAihUrGB4ennPnjcYY4+PjbZcFOl62m/bs2dPzGOayHDjpGQCHsenFh8G3vsrw8Fdn/UwV6jVfda2b9aqeutatrvWar35KCGKadTll+f8B12TmkxHxW8DVwNkHfSjzSuBKgLVr1+bQ0NCcO3/3jptpNBq0WxZgaOiMjpbtpuHh4bbq1mtlj1dV6jUfda2b9aqeutatrvWar35KCHYCJ7YsnwDsai2QmbtbFv8RePsCxKVD4GRDklQN/TSG4Fbg5Ig4KSIOBy4Arm8tEBHHtyy+BrhnAeNTSZOTDY02xki+P9nQlpHRXocmSZqibxKCzNwPvBHYSvNC/7HMvDsi3hoRrymK/U5E3B0RXwJ+B3hdb6JVO5xsSJKqo5+6DMjMG4Ebp6z745bXm4BNCx2X5sfJhiSpOvqmhUD142RDklQdJgTqGicbkqTq6KsuA9XL5N0Eb772TvaNTzDYwbsMJu9eGG2MMfj5bd69IEmHyIRAXbV+zSDXfOFBAD56cWfmYvBRyZLUeXYZqHK8e0GSOs+EQJXj3QuS1HkmBKoc716QpM4zIVDlePeCJHWegwpVOd28e0GSFisTAlXS5N0LjUaDrZcc9MBLSVJJJgTSIuBTJyXNxYRAqjnnbZDUDgcVSjXnvA2S2mFCINWc8zZIaoddBlJFtTsuYOXyAUanufg7b4OkVrYQSBU0OS5gtDFG8v1xAVtGRg8q67wNktphC4FUQbONC5jaSlB23gbvSJAWJxMCqYLKjgto96mT3pEgLV52GUgV1K3nOXhHgrR4mRBIFdStcQHekSAtXiYEUgWtXzPIZeedwuFLmv+FB5cPcNl5pxxys75Pkixvy8goZ16+jdd9Yi9nXr5t2oGdUhXMOYYgIo5uYzsTmdnoQDyS2tTuuIAyNq5bdcAYAli8dyS0M7jSMReqk3YGFe4qfmKWMkuAZ3ckIkk945Mkm9q90Je522Nyu97BoX7VTkJwT2auma1ARIx0KB5JPdaNloeqafdCX2bMha0J6nftjCFo5y/C4vyrIamW2r3Qlxlz4R0c8zM5RuOkS29wjEaXzZkQZOZ3ASLiL6a+FxFLWstIUh20e6Evc7eHd3B8X7sX+TIzcpo4HLoydxkMRsSFkwsR8Szg050PSZJ6q90LfZm7PbyDo6nMRb7dVpUy29TMysxUeDGwNSK+BiTwAeCSTgYTEecAf09zkOL7MvPyKe8fAXwQeBGwG3htZj7QyRh6pUqDjaoUK1Q33tHGGIOf39bX8Vb12M4Vb5nBlZNjLhqNBlsvOXvGfZe9g6Oux7bMQMx2W1XmO7izCv/HYOF+F9q57fCDwBeBEeANwEeA/cD6zLy3U4EU3Q/vAn4G2AncGhHXZ+aXW4r9BvDNzHxeRFwAvB14badi6JUqDTaqUqxgvN1UpVihfLydHlxZJsmo87Et03XS7pM6qzi4s92L/ELGG5k5e4GInwZObfl5DnAr8Fngrsy8tiOBRJwB/GlmriuWNwFk5mUtZbYWZW6OiKXAN4Bn5iyVWLt2bd52221z7v8DF/42xz58Py88ce5pF7788BMArD7+GR0pO/Jggyf3jx+0/oilS1jz7OXTfuaB3XsBeM4xPzBnDA/s3suTTz7JqpVz122u7c4n1k4fr9ay+/fvn/WcLcSx7WTZsvF289jOVXa+x7ZTv4tly3bzd7ed38Uy2+y339u5zlmZeMuUfWzPk9z32F4mJr7/J/6ww4LnHvsDHHvkEV3f/2QM9z+2l/GJ5IilSzjx6IED9j2fsu3Wa2q89x01yHtfeC7Q7J767KUzt0hFxO2ZuXbGAtOYs4UgM28CbmrZyVJgNc3k4CVARxICYBB4qGV5J/Dimcpk5v6I+BZwDPBYa6GI2ABsAFixYgXDw8Nz7vzw3MfhhyWNxtzzKy3JCYCOlZ3ul3Ny/Uyfe+I7xXaXPDVnDE98Z4LM9uo213bnE2unj1dr2cPmOGcLcWw7WbZsvN08tnOVne+x7dTvYtmy3fzdbed3scw2++X39lv7kkf2TpAJtz/wOMc+PTjq8IOnpCkT79FPS76xF1q/xkU0108tuxRY8fTgkb3JRMKyw4Jjnx4s3T9Go+Xbf5ltlon1W/uSbxT1nyzztf/ew97vfOeg41Cm7NcbEwckAwATE8nXH9vL0v0HtmrMFO9oY6yta1sZc7YQLJSIOB9Yl5mvL5Z/BTg9M3+7pczdRZmdxfLXijK7Z9puuy0EAMPDwwwNDc2/EvN05uXbpm0Wmy0DfO17bwbaa8p87XtvLvo3X9FW2dm2O59Yu2muc7YQx7aTZfvt+M5mvse2U7+LZct2+9h28u9H2Vi3jIyWmkyqneM1takammMephs0OZ94O90n3u42y8TarbInXXoD0115A7j/8lfNe7sHbGseLQRz3mUQEV/sRJk27ARObFk+geYMidOWKVoqjgIe78C+e6pbD6rphirFCsbbTVWKFaoVb5lYJy/c+8ab3/rnGmG/ZWSUkQcb3HL/47Penldm3oSyx3b9mkE+e+nZ3H/5q/jspWd3pC+83W1261bRsmMjpjPd+oX8vW3nLoPnR8Sds7wfNC/Mh+pW4OSIOAkYBS4AfnFKmeuBi4CbgV8Ats02fqAqJn9xqzCiuEqxQrXjHW2M9fXUwVU+tv0eb5lYy4ywnyl5aN3npDIXuKoe27n+j7U7qLFs2TJ3nCzksW0nIfixNspM38lRQjEm4I3AVpq3Hb4/M++OiLcCt2Xm9cA/AR+KiHtptgxccKj77Rfr1wz25X+e6VQpVqhuvL3qwiqjH47t5DfefeMTnHn57LeQ9UO87Wo31jIX7jLJQ5kLXJl4+0G7/8fKXLi7eZFfqGPbzqDCrwNExNnALwEN4C7gTpp3GTzZqWAy80bgxinr/rjl9XeB8zu1P0nVVuYbb12VuXCXSR588mW5C3e/XuTLKDMx0YdpzkOwFHghsB54AfC8LsQlqWbKfJNvV9kJaeqozIW7TPJQpe6rbipz4e7Hi3wZZRKCezPz48Xr/9ONYCTVU7e+yft8gHLfTMt+669S95UOXZmZCm+OiD/IzL/pfliS6mQ+U8u205pQtp+7rtr9ZlqlwX9aeO20EFxNcxKi44B1EfE7wJeKnzsz09YCSbOaz9Sy7bQm2M9dXtWbtdU9hzpT4enYfSBpDmW+yZdpTfAbr9Q5ZcYQAM3bA2neYTDb3ASS9D1lvsmXHRfgN16pM+acqVCSDtX6NYNcdt4pDBYtAoPLB6adAhfKzeImqXNKtxBI0nx0YzIYSZ1jC4G0ANqdP14HtiYEs7cmSOocWwikLnM2vfIcFyAtPFsIpC4r89Q4SeoVEwKpy5xNT1IVmBBILbrR1++oeUlVYEIgFWbq6z/UpGDjulUMLFtywDpHzUvqNyYEUqFbff2OmpdUBd5lIBW62dfvqHlJ/c4WAqlgX7+kxcyEQCqU7et3siFJdWKXgVQo8+Q8JxuSVDcmBFKLdvv6yzyiV5KqwC4DaR6cbEhS3ZgQSPPgAERJdWNCIM2Dkw1JqhvHEEjzUGYAoiRVgQmBNE9ONiSpTuwykCRJJgSSJKlPEoKIODoiPhURXy3+/aEZyo1HxB3Fz/ULHackSXXVFwkBcClwU2aeDNxULE9nLDNPK35es3DhSZJUb/2SEJwLXF28vhpY38NYJEladCIzex0DEdHIzOUty9/MzIO6DSJiP3AHsB+4PDO3zLC9DcAGgBUrVrxo8+bNbcWxZ88ejjzyyHnUYOFddktzRrxNL557IpzLbhljfHycP3rp7HX73K6neP/2fexPOOZpwc//6DJeunJZR+Ltliqds7LqWjfrVT11rVtd6wVw1lln3Z6Za8t8ZsFuO4yITwPHTfPWW0ps5tmZuSsingtsi4jtmfm1qYUy80rgSoC1a9fm0NBQWxsfHh6m3bK99u4dNwMwNHRGW2UbjcasddsyMsqHbtrO/iI/3P3d5EP3jLP6+av7+ta6Kp2zsupaN+tVPXWtW13rNV8LlhBk5stnei8iHomI4zPz4Yg4Hnh0hm3sKv69LyKGgTXAQQmByvNhPZK0uPXLGILrgYuK1xcB/zK1QET8UEQcUbw+FjgT+PKCRVhzPqxHkha3fkkILgd+JiK+CvxMsUxErI2I9xVlng/cFhFfAv6N5hgCE4IO8WE9krS49cXUxZm5G/jpadbfBry+eP054JQFDm3R2LhuFZuu235At4EP65GkxaMvEgL1ng/rkaTFzYRA3+PDeiRp8eqXMQSSJKmHTAgqaMvIKCMPNrjl/sc58/JtbBkZ7XVIkqSKMyGomC0jo2y6bjv7xicAGG2Msem67SYFkqRDYkJQMbNNICRJ0nyZEFSMEwhJkrrBhKBinEBIktQNJgQVs3HdKgaWLTlgnRMISZIOlfMQVIwTCEmSusGEoIKcQEiS1Gl2GUiSJBMCSZJkQiBJkjAhkCRJmBBIkiRMCGpv8kFIO7454YOQJEkzMiGoMR+EJElqlwlBjfkgJElSu0wIaswHIUmS2mVCUGM+CEmS1C4TghrzQUiSpHb5LIMaa30Q0mhjjEEfhCRJmoEJQc1NPghpeHiYoaGhXocjSepTdhlIkiQTAkmSZEIgSZLok4QgIs6PiLsjYiIi1s5S7pyI2BER90bEpQsZoyRJddYXCQFwF3Ae8JmZCkTEEuBdwCuA1cCFEbF6YcKTJKne+uIug8y8ByAiZit2OnBvZt5XlN0MnAt8uesBSpJUc/3SQtCOQeChluWdxTpJknSIIjMXZkcRnwaOm+att2TmvxRlhoE3ZeZt03z+fGBdZr6+WP4V4PTM/O1pym4ANgCsWLHiRZs3b24rxj179nDkkUe2V6GKqWvd6lovqG/drFf11LVuda0XwFlnnXV7Zs44Jm86C9ZlkJkvP8RN7ARObFk+Adg1w76uBK4EWLt2bbY7IU+dJ++pa93qWi+ob92sV/XUtW51rdd8VanL4Fbg5Ig4KSIOBy4Aru9xTJIk1UJfJAQR8XMRsRM4A7ghIrYW61dGxI0AmbkfeCOwFbgH+Fhm3t2rmCVJqpN+ucvg48DHp1m/C3hly/KNwI0LGJokSYtCX7QQSJKk3jIhkCRJJgSSJMmEQJIkYUIgSZIwIZAkSZgQSJIkTAgkSRImBJIkCRMCSZKECYEkScKEQJIkYUIgSZIwIZAkSZgQSJIkTAgkSRImBJIkCRMCSZKECYEkScKEQJIkYUIgSZIwIZAkSZgQSJIkTAgkSRImBJIkCRMCSZKECYEkScKEQJIk0ScJQUScHxF3R8RERKydpdwDEbE9Iu6IiNsWMkZJkupsaa8DKNwFnAe8t42yZ2XmY12OR5KkRaUvEoLMvAcgInodiiRJi1JfJAQlJPDJiEjgvZl55XSFImIDsKFY3BMRO9rc/rFAXVsf6lq3utYL6ls361U9da1bXesFsKrsBxYsIYiITwPHTfPWWzLzX9rczJmZuSsingV8KiL+KzM/M7VQkShMmyzMEeNtmTnjGIYqq2vd6lovqG/drFf11LVuda0XNOtW9jMLlhBk5ss7sI1dxb+PRsTHgdOBgxICSZJUTl/cZdCOiPiBiPjBydfAz9IcjChJkg5RXyQEEfFzEbETOAO4ISK2FutXRsSNRbEVwH9GxJeALwA3ZOYnOhxK6W6GCqlr3epaL6hv3axX9dS1bnWtF8yn2zwzuxGIJEmqkL5oIZAkSb1lQiBJkkwIJkXEORGxIyLujYhLex1Pp9RpuueIeH9EPBoRd7WsOzoiPhURXy3+/aFexjgfM9TrTyNitDhvd0TEK3sZ43xExIkR8W8RcU8xNfnvFuvrcM5mqlulz1tEPC0ivhARXyrq9WfF+pMi4pbinH00Ig7vdaxlzVK3qyLi/pZzdlqvY52PiFgSESMR8a/FculzZkJA80AC7wJeAawGLoyI1b2NqqPOyszTanC/7VXAOVPWXQrclJknAzcVy1VzFQfXC+AdxXk7LTNvnOb9frcf+IPMfD7wEuANxf+rOpyzmeoG1T5vTwJnZ+apwGnAORHxEuDtNOt1MvBN4Dd6GON8zVQ3gI0t5+yO3oV4SH4XuKdlufQ5MyFoOh24NzPvy8x9wGbg3B7HpCmKSagen7L6XODq4vXVwPoFDaoDZqhX5WXmw5n5xeL1t2n+sRqkHudsprpVWjbtKRaXFT8JnA1cW6yv6jmbqW6VFxEnAK8C3lcsB/M4ZyYETYPAQy3LO6nBf+7C5HTPtxdTOtfNisx8GJp/pIFn9TieTnpjRNxZdClUrlm9VUQ8B1gD3ELNztmUukHFz1vR9HwH8CjwKeBrQCMz9xdFKvv3cWrdMnPynL2tOGfviIgjehjifP0d8GZgolg+hnmcMxOCpumeqlSLzJHmdM8/QbM75A0R8bJeB6S2vBv4EZpNmw8Df9PbcOYvIo4E/i/we5n5RK/j6aRp6lb585aZ45l5GnACzdbT509XbGGj6oypdYuIHwc2AT8G/CRwNHBJD0MsLSJeDTyambe3rp6m6JznzISgaSdwYsvyCcCuHsXSUa3TPQOT0z3XySMRcTxA8e+jPY6nIzLzkeKP1wTwj1T0vEXEMpoXzH/OzOuK1bU4Z9PVrS7nDSAzG8AwzTESyyNicqr7yv99bKnbOUX3T2bmk8AHqN45OxN4TUQ8QLO7+2yaLQalz5kJQdOtwMnFqMzDgQuA63sc0yFbJNM9Xw9cVLy+CGj3QVl9bfKCWfg5Knjein7MfwLuycy/bXmr8udsprpV/bxFxDMjYnnxegB4Oc3xEf8G/EJRrKrnbLq6/VdLcho0+9krdc4yc1NmnpCZz6F57dqWmb/EPM6ZMxUWituD/g5YArw/M9/W45AOWUQ8l2arADQfZPWRKtcrIq4Bhmg+svQR4E+ALcDHgGcDDwLnZ2alBujNUK8hms3OCTwAXDzZ714VEfFTwH8A2/l+3+Yf0uxrr/o5m6luF1Lh8xYRL6Q5AG0JzS+MH8vMtxZ/SzbTbFIfAX65+EZdGbPUbRvwTJrN7HcAv9Uy+LBSImIIeFNmvno+58yEQJIk2WUgSZJMCCRJEiYEkiQJEwJJkoQJgSRJwoRAkiRhQiBJkjAhkNQFEfHyiPhQr+OQ1D4TAkndcCrN2dEkVYQJgaRuOBUYiYgjIuKqiPjLYq54SX1q6dxFJKm0U2k+xXAr8L7M/HCP45E0B59lIKmjiscCPwZ8nebDfW7ucUiS2mCXgaROW03zkeL7gfEexyKpTSYEkjrtVOBzNJ/N/oGIWNHjeCS1wYRAUqedCtyVmV8BLgE+VnQjSOpjjiGQJElVwb/xAAAAKUlEQVS2EEiSJBMCSZKECYEkScKEQJIkYUIgSZIwIZAkSZgQSJIk4P8DuyOxd4tc2hAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as sig\n",
    "\n",
    "p = 0.90*np.exp(-1j*np.pi/4)\n",
    "a = np.poly([p, np.conj(p)])  # denominator coefficients\n",
    "b = [1, 0, 0]  # numerator coefficients\n",
    "N = 40  # number of computed samples\n",
    "\n",
    "# generate input signal (= Dirac impulse)\n",
    "k = np.arange(N)\n",
    "x = np.where(k==0, 1.0, 0.0)\n",
    "\n",
    "# filter signal using transposed direct form II\n",
    "h = sig.lfilter(b, a, x)\n",
    "\n",
    "# plot output signal\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.stem(h)\n",
    "plt.title('Impulse response')\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$h[k]$')\n",
    "plt.axis([-1, N, -1.5, 1.5])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Compute the step-response $h_\\epsilon[k] = \\mathcal{H} \\{ \\epsilon[k] \\}$ of the filter by modifying above example.\n",
    "\n",
    "Solution: The step response can be computed by changing the input signal to \n",
    "\n",
    "```python\n",
    "x = np.where(k>=0, 1.0, 0.0)\n",
    "``` \n",
    "\n",
    "or by cumulative summation of the impulse response due to the relation $h_\\epsilon[k] = \\sum_{\\kappa =0}^k h[k]$ of the step response to the impulse response\n",
    "\n",
    "```python\n",
    "he = np.cumsum(h)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "**Copyright**\n",
    "\n",
    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2018*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}