FROLSExampleNotebook_Matlab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Data\n",
    "\n",
    "In the cell below, we obtain a signal of the system $y(k) = 0.1y(k-1) - 0.5u(k-1)y(k-1) + 0.1u(k-2)$. You can change it as you wish, or use your own data collected elsewhere. If you want to use your own data, it is better to use [this software](https://github.com/rnwatanabe/FROLSIdentification) in the Matlab or Octave environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AIWDhsP5jURAAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMi1GZWItMjAxNiAxMToyNzox\nNSnnhUAAAAAkdEVYdFNvZnR3YXJlAE1BVExBQiwgVGhlIE1hdGhXb3JrcywgSW5jLjxY3RgAABO4\nSURBVHic7d3blqLIFgVQOaP+/5c5D3bZtiAictmXOZ+6s8w0CHbEIgBxGMfxBgBX+9/VDQCA200g\nARCEQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEg5TAMwzAMV7cC4EACiQ9kIZEdV58q/3x/\nrm4Aq3jkIFCeFRIAIQwOvVO4nzp47Kzn/30+q/C8N799zezbvZyyUC3s7l1x3rbWp8rPywopt5dh\nMz3lPT0P/tVpceOQ48wW5171qfIzcg0psWF4XeDeB+HLD9e85iPjk4NMFy4//pF3P1T58VkhJfbL\nUHH7ENeaPWk2Pc920PsSk0AqbhpajvjoQOVnJJAACEEgARCCQAIgBIFU3OyN4Je0BM6k8jMSSE0t\nf4DjzJbQ0+wNde8+rzp9zY/ve8Rf5nc+h1TfdJg9xuT9E+lrxuGamQI2WCi/H+tT5adjhVTcOI7L\nn5Od/u/yC2Avs8W2V32q/Iw8y64sR3b0pPLzskICIASBBEAIAgmAEFxDAiAEKyQAQhBIAIQgkAAI\nQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAISQOJBWfr8WZKS8aShrIBmrAMVkDSQozPEWPf25\nugFbGK5UpbbpLF8gPb6f+OPQNbZZKc6XsDxasqZ6VTjrxSnyBckC6ZFGK1+fYh+caRh8A9ar1NO6\nvflChc/KUuSZriF9m0YAJJIpkAAoLE0gWR7tQgdSmwpPLeU1pOlPVCFAdmlWSADUlmaFNF0DWRsB\nVGKFBEAIAgmAEAQSACGkuYY05eoRhSlvGrJCAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSAB\nEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBI\ndDEMV7cAWCSQrmSKBHgQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSI34HC4QmUAC\nIASBBEAIAgmAEAQSACEIJABCEEgAhCCQduOmaoBfCCSycgQAxQgkAEIQSL1YVQBhCSQAQhBI/7B0\nALiWQAIgBIEEQAgCCYAQBBL/yHUVLVdrgTUE0oFMmtfS/5CLQNqBiQ/gdwLpDBIL4COBlJicAyoR\nSACEIJAACEEgARCCQAIgBIHEzrbdauEGDUAgNdUqAFptLOQlkAAIoXsgHXHsPAwzf9ZBOsCyP1c3\n4GvDf6f2cRyvaskaw3CL3cBA6vXVc62uKdThzWFL8CKHvWRaIQ3DMB2x78ZwZAmbfKX13RWqY18q\nM2OhwsnSrJAe4/n5aPH+w2EYLj+ENNss+GXpk7Rj75X5UqsrC/XyYoarZFoh3SZj1dA9ze/BcGa0\nRIixbbUauaQj9Cq1pVkhbTO7rmJHj9XPCVeA7jtzl3c59ATawh9fWCQ9fuurK08qnHcyniVOE0jb\nxptR2sS3cXgvjPtvhRq3s1eeFspYhfPO9OpGfGkCadbH4XrkW/96qH5QhdS7V62bjPMI7CLZNaRn\nR4xVw39f+nO98a+XH97EEm2kDKTH/d/TAcwG3053h06Pe/1xc/iZ9Da7SHbK7ttPGpKLeQ06y7RC\ner6haMc0CjUJntOYUJsMcJcmkK66vfXyufvyBqyRopGn2XZT3OyDSIjJjjpImkC6c5qOLHZ5dNCF\n95HC+RJfQ3qWfcTGvFc7ZqtSGMdxuuJ5qdKXsJn9lelvQWHJVkiHsgxfpn++suHRQdPXSCNaSbNC\nMjIjsGb6ynLRzv6rOqczK6R5VgPX0v/QUMdAajjZNdzkF/ce0A8QWcdA+si0ddMJ7EQhsV7xQIow\nGPZtQ4QtWilRU4EIigfSjk6YXsPO4Mc9mBzgoV0g1buWUGlbzqHHIKZ2gRTEY048bXJMPQunbjyw\nkkC63VbPd/GnxdkWXttsiQusJJDYwuwP7E4grWL+XUMvAb8QSOfpPF/H2fY4LQFeCKRX509Y03eM\nM2ke/e3mcbYUuJxAiuKSqVkeAHEIpMrW583j41npIipdg4F3BNKSmJNdtGcRxewlIB2BVFC9p1F8\npe2GQ3YCqZ2P8/XLC5b/991fOCcVZA93KuEuez8IpKNEO7EGEFzfQNo8xf+SDe6KBg6VetJoFEiV\nTiIFfGbd0WpvHbPs9G490CiQctmrELsVtDVoTGf2s32al0DaIlrFR2vPRyvvjABaEUjb+R5VgB0J\npNyk1+1TJzTvouabTy4C6VfPA96t3gCbtQik7DP77PedZ98ogBctAimUSkGSeltSN55tftzpauZo\nvQJJPe1IZ6a2y+4LXgPBm8dUr0B6sfJi+MKJMhUPF4rz5ZZHvG/DpyS3DiSCazUUt4nQRRHacBen\nJQ8BmxSZQOrLUGGl2RMD3379I3d6Y4FAqil+0W++YzDj19rS2ezJ/9kXIJCYYYRkdOaS5dtv1cor\n8ocL6x2c1Q+kYjssO7vjEt26/ZIvl/mFpzve1Q8koukwrgJK2u2HNnvhptkLu2vlWyfdocsEEtR0\n/oRV47s/Np+N/PYe9KM3P2b3LhNI0MjCpBlq/grVmBQWricl6kyBRDiJxk83x3388+h32dCMlS/Y\ncKIv3SWu0wikNMrXIjvaVi1ZamzHRxiku2ATpyVHEEhQ3Jn3Bez44tnfirZ42r09tfPmI4EEXZwQ\nHgcJ1ZgUkvaYQKKjpMN1F0ffC3f+pfUUezNFIy8nkIAdrL834fyHWC+f9Pv9jNzRV7O+/ZKBvOEn\nkKCaX6bRLHNZ5Cf6rP+bWXr7NAIJ+or8pJwjfreMqp0gkIDdRPuEzePjol+d9TrztrqVqibQiy6B\n1GR3wi5+HC9XXR/atw3nX6oxTXUJJGCl5YcJXXW5Zcd3Of+uigVxPiUWgUDiPPEf99lW/DNR3Xbl\nXtubq98EEtSRa/aZZRM6E0iQ3gmfOTXJcgKBBLlFjorIbYusxjdLbSCQoLsa01mNrWhOIAEQgkCC\nClI/B+guV2s5gkACIkqaT0mbHYRA4mIGMA3tWPaVRpBAglIqTU90I5CKMA2RwuanmkaTqKmJCCQA\nQvhzdQO+NjwdmYzjeGFL4AgqnLaSrZCG/66TB8tmalHhC3RGeZlWSPfB+XzMOAzDMAyOIqlBhVcl\nSldKtkJ6GZkGKsWo8Bok0DZpAmnh3IXTGhSgwiFNIAFQm0ACIASBBEAIAgmAEAQSACEIJABCSBNI\nCx/I8FkNClDhkCaQ7jxYhdpU+LWSRn/SZk9lenTQOI73J6m8/PCq9sC+0lX4OHokwf4692qyFZIH\nq1CbCj/acT1qX/0u0wrpzhClNhXOvhItuZKtkIDNJF1eTfadQALm7TUJNplM+Z1AAl6dFiHFsmrD\n5uzYAwU6UyABV4o5jcZs1W1dw8I2/iOBBH29zFx7HeD/MiF+9btWcsUIJCjoPoGeOY1mf6+SkZNu\nowQStHPQPPX8Z799i4XXf/xT69d56SbobgQSEF2QJZE8O5pAgtxOmCU3v0WQIAn+RtuM40wLg7f5\nI4G0VvY9TWHHHeyvmfLCDo2wDdum2ObMEkjpOWBk2bf3rZXf0e+28fHDHe/023Z9q/wueEcgQSOR\nz+/9/qZt5/EyBFJQF45qsvv2rNqOs3nSEir5adN0Db4JJMrLOCw3OGczd7lYFaepaz7V+3trPRxo\nPYEELDkhUc68KeP3t97x13khkKCXC+dQ0/deqvakQCKTquOQ3538VbDv3s7Ncr8QSMAFOk/W316m\n6tNXAonLRBhmEdpwrfI98PGTVb8/+24aMBc+2yL1DhVI/CN1HbPBJdfzl2f/bR9NvUqERkZow44E\nElsUGwYdfLtKiLmLp1dodr9LmwsJpO/4vOqZ2m54fIfumnr73Z2NKwkkSCzXdLOBm9amCveDQEqm\ncC0CzQmkoyRNjq8+ckEfauBkPR8hIZAOkbQa4DgGxYI1H7PtoH4gdduj50vaw0mbzbXczXGo+oF0\nvvhVFb+Fu+h2PbzY6da8LV+j9tZtJpAu0LAWG24yfLT5ERJVCaQZe9VBw3qKz065RIRuj9CG49TY\nOoH02dF7ev3fL/m9lr9ruMnEpyw3EEhXUrJXKdzzVTet6nbxTCDVFHD0BmwSVbW6Fy5ae34hkP6V\ncb/6JhXeifBF40zpxgUC6Xa7okS2fc/K+Q04tA2/iNkq2lKQuxBIS9oWWfa4gjJaDTGBFE7G+vv9\nOzdPEKcl2BfMEkjbxR9Uu3+nJ8HZd+/omRQE0hc8APEuxfamaCTwrGkgxXkWQ7p5M0WDUzSSj+zH\nbpoGEpdwChFY0DGQzpzgLpxMT37rr97ul7Y5ccrlFNtBOgbSj9Ti+fQ5dCCQ/nHClHfcCau883Xe\nlgO7E0h1HBpdEQLbc5KgNoG0s7BXOIp9nSihKCR20SKQlp8jYCxFY49Qj6peo3ggVSqCStuyoMlm\ndmBX8q3igbRZpbEU9lnmG14MZaj8KYH0K5NvIvofIhNIAITQJZAu/0I8DmXHQQFdAqmVyM8rkhzA\nOwLpVfkZ88wnnZfvTGBHrQPJdFmSVRok1TqQll0ybZkr7/QDNCSQUrpqvl7/vqe1UHRBGQJpI/Pg\nbevz8TxVD5glkM5m5gWYJZD29OOx/3FZtdCwsAEZtmHAQf5c3YCvDcPw/L+jeeuJzgjluVbXFOpL\nbX/1u1BAphXSMAzTEftuDM9qPq43PHaveY/94qUyvypU6CnNCukxnp+PFu8/HIbBISSh3CvzpVZX\nFqpipq1MK6TbZKwauoS1rVaVNJ2lWSFtM7uuOoJpJJe/hTEecSZt4W8uLJIev/XVlafTKpx0Mp4l\nThNI28bbwm8Zvw9Hd8Xuf//3P3gvjGG4jeMhmbTZ7JWnxTJWx8ybXt2IL00gzfo4XCGdjPMIa5io\nPkp2DelZ27H6KOvU9Z268UcY/3r54a1xqdNNxBXSx08anXnePMi8GaQZvPCpONhRxEBa8O0nDcnF\nLoXOIgbSx9uQ0kXRCe3N1iVFpCtFiCzNNaS8aURYB5XStpviZh9EAq2kCaS7JmnUYyuL2+XRQe4j\npZWIp+wWFH76ZP4t4F/3jzct3/LwEjazvzL9LSgs2QoJvnLhZL7h0UHT10gjWkmzQjIySWe5aGf/\nVZ3TmRUSoZmfC7NzeSGQOIrpBviKQOJ6ogu4CSQAghBIlGKxBXkJJHKYTZqV8SOlIAWBtIUJDmB3\njQKpWIoU25wN9AAU0yiQOJqEAH4hkFYx1QIcTSABEIJA4i3rQuBMAgmAEAQSACEIJABCaBdIrosA\nxNQukACISSABEIJACsGJRACBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmA\nEAQSACEIpLg84A5oRSABEIJAAiAEgQRACAIJgBAEEgAhCKTLuIkO4JlAAiAEgQRACAIJgBAEEgAh\nCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRA\nCAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEELiQBqGYRiGq1uRjB7LQnlvo9NSyxpIyg6gmKyBBIU5\n3qKnP1c3YAvDlarUNp3lC6T7iB3Hcc3QNbyn9Elk4zje/2PlbrI3p/RJXskC6ZFGa1688mWQlAqn\nmEzXkL5KIwByyRRIABQW8ZTdyyng+5LI8ogyZisciBhIC6aXKwUVQA0RA0m6UJsKh1kRA2nWdAxb\nGwFU4qYGAEIQSACEIJAACGFwDQaACKyQAAhBIAEQgkACIASBBEAIAgmAENI8qeErz4+863wb4fLD\nLLr10ppHmmbpkyztPJoKf5G9yAuukF52ia+PnNWql4ZhePdY3oWfhO2TLO28VrdeqlHk1T6HND1i\navvIu0eprXkMYOFemu2H6Q+z9EmWdp5AhT+UKfKCK6SXzo3T16eZPVZ60a2X1mxvlj7J0s7jqPBZ\nBYq8VCAt1Gi0lemhxr9m/1UvTWXpkyztPJoK3yBFt9S8qQEeoh0Dwu7KFHmpFRKsFO3UOewuY5EL\nJNqJc4ICDpK0yJ2yo5GF+7KghtRFLpBoIfKHAWEXBYpcIFFf6mNGWKNGkZe6hrSwJ1LvpH1166U1\nAzVLn2Rp57Ua9lKZIi8VSHfBn40RRLdeWjPksvRJlnZeq2EvFSjyao8Ous11cb1tXGnhvs8+vbQ8\n5GYftTL7r3FkaecJVPhdmSIvuEIK/myMIPTSVJY+ydLOa+mlWcG7peAKCYCMCq6QAMhIIAEQgkAC\nIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIpMqGYVj5jZDr\nXwlxqPBiBFJuH4fZyu+78rVYxKTCWxFIZW04HnQISSIqvB6BVNlXR4UOIUlHhRfz5+oGsN3jcO/+\nH8/jbfZI8OWHs+NzGHyrPVGo8G6skBL7alxNB7DTFwSnwruxQqpgdtxODydfXna/XPz8w3EcjWEC\nUuFNWCEVtHLIjeP47pzG3i2CPanwqgRSI8Yhtanw7ARSC4/jxOGva9sD+1LhNQikLl5OXxi0FKPC\nC3BTQy/TEesWWCpR4alZIbXgaJHaVHgNAgmAEKxnc3s+Knw5WTH9TMb016d73ykOQlHhrVgh5bYw\ntD4+RmV2rO7VMNiFCm/FTQ3prT/cc2BIRiq8DyukmjaPTEOaFFR4SQKpsq9OUDibQToqvBiBVNaG\nI0EHjySiwusRSMWtPCp06xFJqfBK7CQAQrBCAiCE/wPBStX4nnO0sAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "x =\n",
       "\n",
       "     4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "Fs = 100; % Sampling frequency of the data acquisition, in Hz.\n",
    "u = randn(2000, 1);    \n",
    "y = zeros(size(u));\n",
    "t = (1:length(y))/Fs;\n",
    "for k = 5:length(y)\n",
    "   y(k) = 0.1*y(k-1) - 0.5*u(k-1)*y(k-1) + 0.1*u(k-2); \n",
    "end\n",
    "subplot(1,2,1);plot(t, u); title('input'); xlabel('t(s)')\n",
    "subplot(1,2,2);plot(t, y); title('output');  xlabel('t(s)')\n",
    "x=2+2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Prepare the data\n",
    "\n",
    "Here we throw away the first 100 samples to avoid transient effects. This prepare could also involve some kind of normalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input = u(100:end);   \n",
    "output = y(100:end);  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Identification Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$\\require{dsfont}$\n",
    "The model structure that we use to identify the nonlinear  system is called as NARMAX and has the following structure:\n",
    "\n",
    "$\n",
    "\\begin{equation}\n",
    "y(k) = F\\left[y(k-1), y(k-2),..., y(k-m_y), u(k-d), ..., u(k-d-m_u), \\xi(k-1), ..., \\xi(k-m_\\xi) \\right]  \n",
    "\\end{equation}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "where $u$ is the input signal, $y$ is the output signal and $\\xi$ is the residue signal. The format of the NARMAX model used is the polynomial form:\n",
    "\n",
    "$\n",
    "\\begin{equation}\n",
    "  y(k) = \\beta_0 + \\sum\\limits_{i_1=1}^n\\beta_{i_1}x_{i_1}(k) + \\sum\\limits_{i_1=1}^n\\sum\\limits_{i_2=i_1}^n\\beta_{i_1i_2}x_{i_1}(k)x_{i_2}(k)+...+\\sum\\limits_{i_1=1}^n...\\sum\\limits_{i_l=i_{l-1}}^n\\beta_{i_1i_2... i_l}x_{i_1}x_{i_2}... x_{i_l} + \\xi(k)\n",
    "\\end{equation}\n",
    "$\n",
    "with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$\n",
    "\\begin{equation}\n",
    "x_m(k) = \\left\\{ \\begin{array}{l l}\n",
    "u(k-m) & 1\\leqslant m \\leqslant m_u\\\\\n",
    "y(k-(m - m_u)) & m_u+1\\leqslant m \\leqslant m_u + m_y\\\\\n",
    "\\xi(k - (m - m_u - m_y)) & m_u+m_y+1 \\leqslant m \\leqslant m_u+m_y+m_\\xi\n",
    "\\end{array}\\right.  \n",
    "\\end{equation}\n",
    "$\n",
    "\n",
    "where $m_u$ is the number of input  past values to be considered in the model search, $m_y$ is the number of output  past values to be considered in the model search, $m_\\xi$ is the number of the residue past values to be considered in the model search and $l$ is the maximal polynomial degree to be considered to form the combinations between the signals $u$, $y$ e $\\xi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Calling each combination $x_{i_1}x_{i_2}... x_{i_l}$ as $p_i$, the Equation above can be rewritten as:\n",
    "\n",
    "$\n",
    "\\begin{equation}\n",
    "  y(k) = \\sum\\limits_{i=1}^M \\beta_ip_i(k) + \\xi(k)\n",
    "\\end{equation}\n",
    "$\n",
    "\n",
    "where $M$ is the number of combinations used in the model search. In matricial form, the equation above is:\n",
    "\n",
    "$\n",
    "\\begin{equation}\n",
    "  Y = PB+\\xi \n",
    "\\end{equation}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "where:\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "  Y&=\\left[y(1),y(2), ..., y(N)\\right]^T\\\\\n",
    "  B &= \\left[\\beta_1, \\beta_2, ..., \\beta_M \\right]^T\\\\\n",
    "  \\xi&=\\left[\\xi(1), \\xi(2), ..., \\xi(N) \\right]^T\\\\\n",
    "  P = \\left[p_1,p_2,..., p_M\\right] &= \\left[\\begin{array}{cccc}\n",
    "  p_1(1)&p_2(1)&...&p_M(1)\\\\\n",
    "  \\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
    "  p_1(k)&p_2(k)&...&p_M(k)\\\\\n",
    "  \\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
    "  p_1(N)&p_2(N)&...&p_M(N)\n",
    "  \\end{array}\\right]  \n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $^T$ stand for the matrix or vector transponse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Parameters of the identification process\n",
    "\n",
    "For the identification of the system described at the beginning of this notebook we use the parameters in the cell bellow. As in this demonstration of the FROLS method we know our system, the used parameters are exactly what we need,  but normally we do not know what parameters to use! Normally, to find these parameters, we need to try some values before finding the right ones. If you use $m_u, m_y$ and $m_\\xi$ values higher than necessary, the algorithm will work but the time of processing will be longer. If you use them lower than necessary, the algorithm will not find the correct model. Try to modify the parameters below to see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 2; \n",
    "my = 1; \n",
    "me = max(mu,my); \n",
    "degree = 3; \n",
    "delay = 1; "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dataLength = 500; % number of samplings to be used during the identification process. Normally a number between 400 and\n",
    "% 600 is good. Do not use large numbers.\n",
    "pho = 1e-2; % a lower value will give you more identified terms. A higher value will give you less.\n",
    "delta = 1e-1; "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the P matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "numberOfCandidates = round(findPMatrixSize(mu - delay + 1, my, degree));\n",
    "begin = randi([1 length(input) - dataLength - 1], 1);\n",
    "u = input(begin+1:begin + dataLength);\n",
    "y = output(begin+1:begin + dataLength);\n",
    "[p, D]= buildPMatrix(u, y, degree, mu, my, delay);\n",
    "y = y(max(mu,my) + 1:end);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The terms that the FROLS algorithm will use to try to find the model are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Columns 1 through 6\n",
       "\n",
       "    '1'    'u(k-1)'    'u(k-2)'    'y(k-1)'    'u(k-1)u(k-1)'    'u(k-1)u(k-2)'\n",
       "\n",
       "  Columns 7 through 10\n",
       "\n",
       "    'u(k-1)y(k-1)'    'u(k-2)u(k-2)'    'u(k-2)y(k-1)'    'y(k-1)y(k-1)'\n",
       "\n",
       "  Columns 11 through 13\n",
       "\n",
       "    'u(k-1)u(k-1)u(k-1)'    'u(k-1)u(k-1)u(k-2)'    'u(k-1)u(k-1)y(k-1)'\n",
       "\n",
       "  Columns 14 through 16\n",
       "\n",
       "    'u(k-1)u(k-2)u(k-2)'    'u(k-1)u(k-2)y(k-1)'    'u(k-1)y(k-1)y(k-1)'\n",
       "\n",
       "  Columns 17 through 19\n",
       "\n",
       "    'u(k-2)u(k-2)u(k-2)'    'u(k-2)u(k-2)y(k-1)'    'u(k-2)y(k-1)y(k-1)'\n",
       "\n",
       "  Column 20\n",
       "\n",
       "    'y(k-1)y(k-1)y(k-1)'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disp(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The first 15 lines of the $P$ matrix are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Columns 1 through 7\n",
       "\n",
       "    1.0000    0.4590    1.9524    0.0381    0.2107    0.8962    0.0175\n",
       "    1.0000   -0.3993    0.4590    0.1903    0.1595   -0.1833   -0.0760\n",
       "    1.0000   -0.2254   -0.3993    0.1029    0.0508    0.0900   -0.0232\n",
       "    1.0000   -0.2803   -0.2254   -0.0180    0.0786    0.0632    0.0051\n",
       "    1.0000    1.3067   -0.2803   -0.0269    1.7075   -0.3663   -0.0351\n",
       "    1.0000   -0.2844    1.3067   -0.0132    0.0809   -0.3716    0.0037\n",
       "    1.0000   -0.7917   -0.2844    0.1275    0.6268    0.2251   -0.1009\n",
       "    1.0000    0.2406   -0.7917    0.0348    0.0579   -0.1905    0.0084\n",
       "    1.0000   -0.0946    0.2406   -0.0799    0.0090   -0.0228    0.0076\n",
       "    1.0000    0.7932   -0.0946    0.0123    0.6292   -0.0751    0.0098\n",
       "    1.0000    0.9594    0.7932   -0.0131    0.9205    0.7611   -0.0126\n",
       "    1.0000    1.0871    0.9594    0.0843    1.1818    1.0430    0.0916\n",
       "    1.0000    0.6472    1.0871    0.0586    0.4188    0.7035    0.0379\n",
       "    1.0000    0.6326    0.6472    0.0956    0.4001    0.4094    0.0605\n",
       "    1.0000    1.9033    0.6326    0.0440    3.6225    1.2039    0.0838\n",
       "\n",
       "  Columns 8 through 14\n",
       "\n",
       "    3.8118    0.0744    0.0015    0.0967    0.4114    0.0080    1.7497\n",
       "    0.2107    0.0874    0.0362   -0.0637    0.0732    0.0303   -0.0841\n",
       "    0.1595   -0.0411    0.0106   -0.0114   -0.0203    0.0052   -0.0359\n",
       "    0.0508    0.0041    0.0003   -0.0220   -0.0177   -0.0014   -0.0142\n",
       "    0.0786    0.0075    0.0007    2.2313   -0.4786   -0.0459    0.1027\n",
       "    1.7075   -0.0172    0.0002   -0.0230    0.1057   -0.0011   -0.4856\n",
       "    0.0809   -0.0363    0.0163   -0.4962   -0.1782    0.0799   -0.0640\n",
       "    0.6268   -0.0275    0.0012    0.0139   -0.0458    0.0020    0.1508\n",
       "    0.0579   -0.0192    0.0064   -0.0008    0.0022   -0.0007   -0.0055\n",
       "    0.0090   -0.0012    0.0002    0.4991   -0.0596    0.0077    0.0071\n",
       "    0.6292   -0.0104    0.0002    0.8832    0.7302   -0.0121    0.6037\n",
       "    0.9205    0.0809    0.0071    1.2847    1.1339    0.0996    1.0007\n",
       "    1.1818    0.0637    0.0034    0.2711    0.4553    0.0245    0.7648\n",
       "    0.4188    0.0619    0.0091    0.2531    0.2589    0.0383    0.2649\n",
       "    0.4001    0.0279    0.0019    6.8945    2.2914    0.1595    0.7615\n",
       "\n",
       "  Columns 15 through 20\n",
       "\n",
       "    0.0341    0.0007    7.4420    0.1452    0.0028    0.0001\n",
       "   -0.0349   -0.0145    0.0967    0.0401    0.0166    0.0069\n",
       "    0.0093   -0.0024   -0.0637    0.0164   -0.0042    0.0011\n",
       "   -0.0011   -0.0001   -0.0114   -0.0009   -0.0001   -0.0000\n",
       "    0.0098    0.0009   -0.0220   -0.0021   -0.0002   -0.0000\n",
       "    0.0049   -0.0000    2.2313   -0.0225    0.0002   -0.0000\n",
       "    0.0287   -0.0129   -0.0230    0.0103   -0.0046    0.0021\n",
       "   -0.0066    0.0003   -0.4962    0.0218   -0.0010    0.0000\n",
       "    0.0018   -0.0006    0.0139   -0.0046    0.0015   -0.0005\n",
       "   -0.0009    0.0001   -0.0008    0.0001   -0.0000    0.0000\n",
       "   -0.0100    0.0002    0.4991   -0.0082    0.0001   -0.0000\n",
       "    0.0879    0.0077    0.8832    0.0776    0.0068    0.0006\n",
       "    0.0412    0.0022    1.2847    0.0692    0.0037    0.0002\n",
       "    0.0391    0.0058    0.2711    0.0400    0.0059    0.0009\n",
       "    0.0530    0.0037    0.2531    0.0176    0.0012    0.0001"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disp(p(1:15,:))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The system identification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "global l q g err A ESR M0;\n",
    "s = 1;\n",
    "ESR = 1;\n",
    "M = size(p, 3);\n",
    "l = zeros(1, M);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The function that implements the FROLS algorithm is written below, with brief explanations about some parts inserted in the code. Actually this function can handle more than one data collection of the same system (that is the reason of the m in mfrols), but in this notebook we will not use multiple trials. Just ignore the outside loop and the variable j.\n",
    "\n",
    "function beta = mfrols(p, y, pho, s)\n",
    "    \n",
    "$\\text{The global variables are used due to the lack of pointers in Matlab/Octave}$\n",
    "\n",
    "    global l;\n",
    "    global err ESR;\n",
    "    global A;\n",
    "    global q g M0;\n",
    "    beta = [];\n",
    "    M = size(p,2);\n",
    "    L = size(p,3);\n",
    "    gs=zeros(L,M);\n",
    "    ERR=zeros(L,M);\n",
    "    qs=zeros(size(p));\n",
    "    \n",
    "    for j=1:L\n",
    "        sigma = y(:,j)'*y(:,j);\n",
    "        \n",
    "        for m=1:M\n",
    "            if (max(m*ones(size(l))==l)==0)\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\text{The Gram-Schmidt method was implemented in a modified way, as shown in Rice, JR(1966)}$\n",
    "\n",
    "$\\text{For each } m, \\text{ we compute: } \\\\\n",
    "qs_m = p_m - \\sum\\limits_{r=1}^{s-1}\\frac{q_r^Tqs_m}{q_r^Tq_r}q_r$\n",
    "                 \n",
    "                 qs(:,m,j) = p(:,m,j);\n",
    "                 for r=1:s-1\n",
    "                     qs(:,m,j) = qs(:,m,j) - (squeeze(q(:,r,j))'*qs(:,m,j))/...\n",
    "                         (squeeze(q(:,r,j))'*squeeze(q(:,r,j)))*squeeze(q(:,r,j));\n",
    "                 end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    "$\\text{We compute for each candidate that was not chosen:}$\n",
    "\n",
    "$gs_m = \\frac{y^Tqs_m}{qs_m^Tqs_m}\\\\\n",
    "\\text{ERR}_m = gs_m^2\\frac{qs_m^Tqs}{y^Ty}$\n",
    "\n",
    "                gs(j,m) = (y(:,j)'*squeeze(qs(:,m,j)))/(squeeze(qs(:,m,j))'*squeeze(qs(:,m,j)));\n",
    "                ERR(j,m) = (gs(j,m)^2)*(squeeze(qs(:,m,j))'*squeeze(qs(:,m,j)))/sigma;\n",
    "            else\n",
    "                ERR(j,m)=0;\n",
    "            end\n",
    "        end \n",
    "    end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    \n",
    "$\\text{And here we chose to be part of the model the candidate term with the highest ERR.}$\n",
    "\n",
    "    ERR_m = mean(ERR, 1);\n",
    "    l(s) = find(ERR_m == max(ERR_m), 1);\n",
    "    err(s) = ERR_m(l(s));\n",
    "    for j=1:L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\text{After the model term has been chosen we compute the values of the matrix } A$\n",
    "\n",
    "\n",
    "$a_{r,s} = \\frac{q_r^Tp_{\\text{chosen term}}}{q_r^Tq_r}, 1 \\leqslant r \\leqslant s - 1\\\\\n",
    "a_{s,s} = 1$\n",
    "\n",
    "        for r = 1:s-1\n",
    "            A(r, s, j) = (q(:,r,j)'*p(:,l(s),j))/(q(:,r,j)'*q(:,r,j));    \n",
    "        end\n",
    "        A(s, s, j) = 1;\n",
    "        q(:, s,j) = qs(:,l(s),j);\n",
    "        g(j,s) = gs(j,l(s));\n",
    "    end    \n",
    "\n",
    "    \n",
    "$\\text{After this process the matrix } A \\text{ will have the format: }\\\\\n",
    "A=\\left[\\begin{array}{cccc}\n",
    "  1&a_{12}&...&a_{1s}\\\\\n",
    "  0&1&...&a_{2s}\\\\\n",
    "  \\vdots&\\vdots&\\ddots&\\vdots\\\\\n",
    "  0&...&1&a_{s-1,s}\\\\\n",
    "  0&0&...&1  \\end{array} \\right]$\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  ESR = ESR - err(s);  \n",
    "  \n",
    "$\\text{The stop criterium is based on the ERR value of the chosen term. If it is lower than } \\rho,\\text{ the process of search for terms stops.}$\n",
    "    \n",
    "    if (err(s) >= pho && s < M)\n",
    "       s = s + 1; \n",
    "       clear qs \n",
    "       clear gs\n",
    "       \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "      \n",
    "$\\text{Here is the recurrent call of the frols algorithm.}$\n",
    "\n",
    "       beta = mfrols(p, y, pho, s);\n",
    "    else\n",
    "       M0 = s;\n",
    "       s = s + 1;\n",
    "       for j=1:L\n",
    "\n",
    "$\\text{ At the end, we compute the coefficients } \\beta \\text{ of each term by doing:}\\\\\n",
    "B = A^{-1}g$\n",
    "\n",
    "            beta(:,j) = A(:,:,j)\\g(j,:)';\n",
    "       end       \n",
    "    end   \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "q = []; err=[]; A=[]; g=[]; \n",
    "beta = mfrols(p, y, pho, s);\n",
    "la = l(1:M0)';\n",
    "Da = D(la)';"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[Dn, an, ln] = NARXNoiseModelIdentification(input, output,  degree, mu, my, me, delay, dataLength, 1, pho,  ...\n",
    "        beta, la);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[a, an, xi] = NARMAXCompleteIdentification(input, output, Da, Dn, dataLength, ...\n",
    "        1,  delta, degree, degree);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Identified Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'u(k-2)'    '0.1'\n",
       "\n",
       "    'u(k-1)y(k-1)'    '-0.5'\n",
       "\n",
       "    'y(k-1)'    '0.1'"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i = 1:length(Da)\n",
    "    disp({Da{i}, num2str(a(i))}) \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that the found terms should be the same of the system difference equation on the top of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Generalized Frequency Response Function (GFRF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If your identified system has terms with inputs and outputs, the GFRF will be non-null for \n",
    "degrees higher than the maximal polynomial degree. In this case, a good number\n",
    "is to add one to your maximal polynomial degree. If you have only\n",
    "inputs or outputs terms, the GFRFdegree will be the maximal\n",
    "polynomial degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "GFRFdegree = 3;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Hn = computeSignalsGFRF(Da, Fs, a, GFRFdegree);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GFRF of order 1: \n",
       "\n",
       "            /   pi f1 i \\\n",
       "         exp| - ------- |\n",
       "            \\     25    /\n",
       "  - ---------------------------\n",
       "       /    /   pi f1 i \\     \\\n",
       "       | exp| - ------- |     |\n",
       "       |    \\     50    /     |\n",
       "    10 | ---------------- - 1 |\n",
       "       \\        10            /\n",
       "GFRF of order 2: \n",
       "\n",
       "             /   pi f1 i \\    /   pi f1 i \\    /   pi f2 i \\\n",
       "          exp| - ------- | exp| - ------- | exp| - ------- |\n",
       "             \\     25    /    \\     50    /    \\     50    /\n",
       "  - --------------------------------------------------------------\n",
       "       /    /   pi f1 i \\     \\ /    /   pi f1 i   pi f2 i \\     \\\n",
       "       | exp| - ------- |     | | exp| - ------- - ------- |     |\n",
       "       |    \\     50    /     | |    \\     50        50    /     |\n",
       "    20 | ---------------- - 1 | | -------------------------- - 1 |\n",
       "       \\        10            / \\             10                 /\n",
       "GFRF of order 3: \n",
       "\n",
       "    /    /   pi f1 i \\                      /   pi f3 i \\                \\\n",
       "  - | exp| - ------- | exp(-#2) exp(-#1) exp| - ------- | exp(- #2 - #1) | /\n",
       "    \\    \\     25    /                      \\     50    /                /\n",
       "  \n",
       "     /\n",
       "     |\n",
       "     |    / exp(-#2)     \\ / exp(- #2 - #1)     \\\n",
       "     | 40 | -------- - 1 | | -------------- - 1 |\n",
       "     \\    \\    10        / \\       10           /\n",
       "  \n",
       "     /    /             pi f3 i \\     \\ \\\n",
       "     | exp| - #2 - #1 - ------- |     | |\n",
       "     |    \\               50    /     | |\n",
       "     | -------------------------- - 1 | |\n",
       "     \\             10                 / /\n",
       "  \n",
       "  where\n",
       "  \n",
       "           pi f2 i\n",
       "     #1 == -------\n",
       "             50\n",
       "  \n",
       "           pi f1 i\n",
       "     #2 == -------\n",
       "             50"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for j = 1:GFRFdegree\n",
    "        disp(['GFRF of order ' num2str(j) ': '])\n",
    "        pretty(Hn{j})\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Undefined variable Hn.\n"
     ]
    }
   ],
   "source": [
    "Hfun = matlabFunction((abs(Hn{1})));\n",
    "freq = [linspace(-50, 0, 50) linspace(0, 50, 50)];\n",
    "plot(freq, Hfun(freq))\n",
    "xlabel('f (Hz)'); ylabel('H_1')\n",
    "x = 2+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AETFRooOilTWQAAACJ0RVh0Q3JlYXRpb24gVGltZQAxOS1KYW4tMjAxNiAxOToyNjo0\nMK4r91wAAAAkdEVYdFNvZnR3YXJlAE1BVExBQiwgVGhlIE1hdGhXb3JrcywgSW5jLjxY3RgAACAA\nSURBVHic7Z3tUeQ800ZbT70BLBksGSwZLBncZAAZLBFsEQFkABlABksGkAFkABno/dF2jyx1t+T5\n9Md1amtrxiN7zIzHx91qySHGSAAAAMCp+d+pdwAAAAAggpAAAABMBAgJAADAJICQAAAATAIICQAA\nwCSAkAAAAEwCCAkAAMAkgJAAAABMAggJAADAJICQAAAATAIICQAAwCSAkAAAAEwCCAkAAMAkgJAA\nAABMAggJAADAJICQAAAATAIICQAAwCSAkAAAAEwCCAkAAMAkgJAAAABMAggJAADAJICQAAAATAII\nCQAAwCSAkAAAAEwCCAkAAMAkgJAAAABMAggJAADAJICQAAAATAIICQAAwCSAkAAAAEwCCAkAAMAk\ngJAAAABMAggJgKkTQjj1LgBwDCAkACZNCCHGeOq9AOAYQEgATAU1Emq3UQgBsRSYNbj4AmASIBIC\nAL8BAAAAkwApOwCOyl4Sa6O2gDwemAsQEgBHJcZ45LREjBFOArMAKTsADgI7AL8vANqBkAAAAEwC\npOwA2In5FlvPd8/BUkGEBAAAYBIgQgKgCcQTABwaREgAAAAmASIkAHKmHwkddA8RC4JT8X+n3gEA\npgWm8Fn5nw9OyNp/ewAAACYCUnZgpSAxNRZ8XODQIEICALSCfCY4KDi8wMLBFD4AzAUICQAAwCRA\nHxJYCOgTOhX45MG+QIQEwPxAXw5YJIiQwMzA9TgASwXXWQAAACYBIiQwXRAJzRrEsmAsiJDAREE3\nCQBrAxESOD3qdTRsBMDagJDAiUEktDaQxwMWOBcAAI4NrkKACiIkcCTQxb1H5v5JwkZABfdDAkcC\n5yAAgA8iJLBnEAkdgUXaHYcNQCYXADAV0Le0cvD1gy3BbR0AAPsFQgIAADAJ0IcEKqBPCJwWHIHr\nARESAACASYAICWzAdehcwDcFFgmEBDp2KXDC+REAsDsQ0hrBZKZzB18WoW9piaAPaXVgqAcAYJrg\n3AQAAGASIGW3WJDQAOsEh/18QYQEAFgayEvPFHxtswdT+AAAlgGEBMD8QAQAFgn6kGbDlPuEJrtj\nAAhT/gUBBtdZYA/ggh0AsDuIkCbHHK/jYCMwEWb32wEpuIX55MDJHQCwThAhnRJczQGwX7a4nptj\nTmKpIPV/MtDvAgAAKTgnAgAAmARI2R0cJATA3sERdQTwIR8fCOngxBgXH4bipwuWR4wRB/aRgZD2\nBiIhcDQWf4kzEfA5Hxn0IQEAlgNqhWbNfsYhpZFB+9FQXasMOKZwqGEyUwDWDJx3OPbwyW6njepa\nVvoLhwIA4LTASQdi14+1DBdaAojqWmKjsk114/sCkRAAAByTPRQ1ZKfsxjN4y1rbbXlfrKE6DgAA\npsNOQnKKyvb+0kFBdRyYFzhiJwjqbHdnomXfHJpk3+6BcmhIB+8OfocAIKeyO9Od7ZtHpWVnOuf7\nbuy7KhvgGAKzAwctWCQTjZDIuOh2rsT58sTPB+JnfCDwwYKJMMFgfYK7NFkmGiGp4Q4v9L0iTkIk\nBACYAnxSwvmnhSlGSDsaBZlcAFbLNH/709yrCTJFIQEAAFghOwnJT53t9yUAAFgS6Fsq2UOEpBZn\n77KWWvPdvmUAFg9+CwsAt7coOcZcdk6FgrWWX0233X6Cw4E+2yODDxwskolOHaQWJrRXK+C6Aywb\n2AgskmVeZ+H6EYB1soDf/pqndZ7oOKQdWed3CQBYAGs+fc3+agIAAMAywDgkAAAAkwBCAgCAGbCG\n21sgZQfA/FhA1z0AJYiQwB5Y/IUbAOAILFNIOD+CZYPwCAhLOt0tU0jgyOD8CCbCks7OjSxpCqJl\nCgnnRwDAeljMGQ9dowAAACbBMiMkAABYOXMsE0eEBAAAYBIgQgJgfszuyheAFiAksAdwfgQA7M4y\nhYTzI1g2yLSDRbJMIYEjg/MjmAi4GJ01yxQSzo8AADA7UGUHAABgEiwzQgIAAJAyi2QmhAQAAGAS\nQEgAzI9ZXO2CSTGL3hkICewBnB8BALuzTCHh/AiWRHk8z+JqF4CxLFNI4Mjg/Hg4cLfyUeBidNYs\nU0j4AYM5op5M2w/mOc7uDEAKLr4AmASIhABYZoQEwOyAjcBBmUX0DCEBcFSOn1ibxZkIAELKDoA5\nMja/h3wgmAWIkMAewDV4yaRKDGAjMAuWKaTpnAjAaokxHk4DEAxYJMsUEjgyaz4/TioSGsV899xh\neX/RqkBmGQCwHNBbNmvw5QHQBF964/cCwOGAkAAAYPnMInZEHxIAAxbZszIKfALgVMzAmQCAjFlc\n7QIwFkRIYA/M9IJ6prsNwFJZppBwogFVZh1knGTP8bMCh2aZQgJHZuJn9h1v6wCYGOP0nTT9PQQO\nyxQSzjVAmHUkNDXwSYKDgt8qAACASbDMCAmsEBQrnwp88rNgFt/R/516BwDYD4j1TwU+ebAvECGB\nmYHrcZrJ1S6YFLO4boCQwB445vnxoLd1AACckGUKCdePCwCRkMMslIxvEIxlmVV2qPQFYJ3gtz9r\nlhkh4YicEbiOBgAwuJoAABwVBDHAYpkREpggiIQAM4spiJbHLD5zjEMCRwIXxUDAwQBUECGB/TOL\na7FZg08YjGUWFwEQEtgD6fkRPQRgO2BZsEwh4cg+Dritw6lY5IeMviWwTCGBI5BGQos8P4Ljs/uB\nBKXNmmUKCefHI4APGQCwX5DuBxX4khPHCTgVOALXA4QEAADLZxbVRstM2YEtwMBVAMBpmYEzAQAZ\ns7jaBWAsiJBWxyEiIYRW4Pggpl8ey7zOwvXjkcEHDgDYHURIi+WY14+wEZgIiJlmDS5sAZgfqIS2\nwCczazDb9+zBL3A9yOX/7/iPkCnV4A8En8xMwdcGwNTJPJTyGi4JlyOggVlIGhHSnJjFIQX2Reoh\nFk+JhEoELYH5gxPcbICNVoIfD6ULxVK/4hsRvYcLHCEO0Pb0wTluiszOPbPb4QmSlYeVNmLYSawi\n9lDKe7ggnHPBbFnmeWTW58c57vwc93kiiIcyu7yHC9VJWVSkglAJzJRlnkdwfgQTJ42HLLWkThIP\n/Ywf/OAznKsrcpBEa42T8NufNfjyTgYy2uuEv3fxCtlqocQuaXtrRW6cGovWd4BtIST8EqcDhATA\nMZCQqEUtjG+jdMVMRVkD/MbBXICQDg6uv1YOHwA/4hcRfYczyy6pk1gwvEp1LXKNRXASIKKZJDMx\nDungTP8gAIdAQiLxis/P+PEeztVVfsSvz8JJrCKKkYg+Q/CdNIuT0dTAh3Z88InvDURCgPFVpIY7\nn72KrLWyFT/DOWVHmuakVFqE0+t44KQjg48b7AH8bqnw0Hc4a1UL0cAuITgrdo/UTztxkrJZgpPA\n1NnPAZrWsLZvsGWtbLRg48YPen5EJFSyZiE58ZDvpO6RoRZ1a94q/Yo/44cSPyUNVvtNgemzh6Oz\nvAFJyzZb1lJvbdK4cfzqwKHpjs8Y/ZhmG7UMN9i1/+obnwVnRW+zfZsF/zqO8NvH6eVw7PrJluFC\nSwDRstZ2W94jiISAyuY6qTnPltbLEbWq5Uf8ytsL5Yq8Vx+RiOi8ZqzlHtjHsQWcdCD2ICQ1sqkK\nqboWfAAmxSBeb86zMZuQqEUt/db09uqKIXQqElQnhUBv/cILnFLB5Nip7Nu5W7DjpJa1jh8M4ccJ\nVJR4yDqAY/wunNSp6CPSedDt8hUphCzYIqLOHGfGWrIi81G0+dA2+5Y8fYs47FfFLL7unXbR0oav\nk5a1ysfOBndkFt8TODJ6Xi55udr9s1GRYDmJ+nAnVZFw4a5Fmo3Sd7Q2228cB/++mHhSZxYnukkP\njM1iqd2/b/Urmf6XBI6GHg+VR0gs4o8ERUVt76074y0qcRKr6F/fY2S9F8dJ6mb7jc/iPDULJv4x\nTnz3mP+degcqxAS/ZQjBSQbSTC4QZor/yU+ckJAu7f6lj9N/1nJ2g2qIj9iJJOUs0Fno1GLxlqwo\n7WWVf5HOtc//PNB5oOdIF5VvZ9ZfH1gSM0jZtW+ZDp/fAyqzkH3LafeeiIhuiyXZQuFv8viuZSfS\nzqQ0ymEu3Wjmomifchny9OBz8vTK2PJF3+xqBt/gTJnFr2MiTDpltwVZtxM4DlP7vVkHwHP/4Kp/\n8Ng/uCGi3jr3lCM2Uv2UquhvsXCwhEOZMk5i/kW6sM3RCL/F83Ajz8WWL4bNnheSu5vgXxHjQj7b\nI7C0CAmskMbrD7YRq+ixePWmaENE97aiGH71r/EqGa4qX9rwL+ZxUmqOKze5d6mpKEXipIugNFtE\nkIRT/6yZ6DgkCAn4lBJK7ybE9xH6N2xw2T/IbJSqKCXVUsltsinegpq+yyKqJidlZOawnHTZfyCO\nkHh1p80inATmC2ZqAPPACYMaVfQ2bCOUKiKiq2QjvHqmFivGomHjjFuiv5quqGqsrEMoc9JloMd+\nyY0W/ciKVDcWfmLH4cjB3Cxix2PMZaeK5KBz2YEFUH776YjTs/7p2bBNpiIiuhwaSxAtqbGRuh2m\nVJHATrJsRFpflNrgb0vwxE7iwOhx+KNQnXQV6D4SEd3WjIVf2bE4piRmIaQ9FDVwl136dF9rbbdl\nMF8yCUUifp6pKF0iWkr9xDrJAqOUC6KP/vE5EQ17j0obVbkp3jdTzm3yLvxGjrfu7Dzhpg5QcnSP\nxe/iMdJVYh3WzH3/9H74qnAV6O8mZ46f2xE45oc8iy8Uhx3YA7ucv/yShGxuuLNiibpQwqbSRhwY\nfRTLWUuWii6HGb+yF0qNqNIicjUOK2v20qo/p+9KKfYrtcRxkgRG+YaKuvC/w2Z3ODmAY7PMYw7X\nd0dm1AduGSj+TNp8UvxB4Ttvo4pHVRTvTSh6mEoVMef9g5akH1vtcRgYlXCopPZRMWX1eUqZAFTH\nQimwnG6GgVHJ7WYQUm4jBk46BWvuLF/aOCRwElpqWPJVfvWvvucqIqL4Y/M/EYVvij8pfA6CITVa\nosRGlCT93mo2SjebCkbtguIlF1sl98pNXRgvpSV8TOa2K6NisFNRFc7dEek2micLuBid+/7vwuy/\nPDBB/CxcqiLSAqN8a9+KsajBRoONEJEhpHMt8JLeILULioguiL4KewmXydud2zXlThqQERtZYVbp\nJEoK2SUCM7N2VLPR3IKkBQhpzeDLA/shldB/8ZmIXsKV3XwDyyYNjDbb/N40GCz/pPir81lW72DZ\niLccvnMnlTaSTVHNRtKyrDIv3yU1ypXRuZWqJUsGXrpOKvmXlAUOEDPdBvoTiYgejHwdEd2h6A4c\nFQgJ7IQTDMVrIqLwREQU/0tWeaH4O3n62rcvbFSqiHobbZ6+b+rrfBvJZtP6Ov2Gery1HxS+dXNY\nERVpNpL3EhzPOf1SpZNERW9G/vCy2NTGUn+ST0t10l2g637hE04Up2fHvqVZxI4z2EUwNVIJXcdH\nInoKaU9HpyIiCk+5iogGNiKi8Nq14VepN0GLjbqFnPpTd7XwHNGmVsK3kTTOyiKsiOqfbSOGnWTZ\nSLbv9EtdDsvTy1ILK39YxnAdoqXMSamNGDhp5kBIYFGIh67jY26gPxQeiIYqIjcw6ha+DtpIS7Iz\ndcqOvVP8TeE1d5Jqo+6lbyKnF6qI1aSKwXIYr+jbSC38E7i0wemXYpzBVeTGSRX+xI2TShsxcBI4\nMDjCQIXUQ1QEQ1Uk+mm30SDX97Nf7tpINkjJed+xEW9K0n1CaSNZhWo24sDOytfJz0wV0kXRGWbF\nOmwsJ8y6GI7PpaS2whqVlaPaiCCkiTKL0KeR5fwlKUv6hk6I1z/Uj44JtxT/JqvcUfwz3MhDv8rv\n4XLXRpslT/3qNRulW6YGG8kW0vKELSIqKlJ8ZUVDGbpV04Bl9TkNKymcIIlRB/9WRwp3fV1w0qxY\nzBlvIX9GxmK+npOgpubiI/HDVEVETTbiNmE4O1uLjbrlT0SakFQbkSukzEayna+ajayIiowUX3p+\ntzq3qmnAtPpcfV91uK6UeFixWnUuJbK09BTod6TXSf+48NufNfjyQIc5gUKDisiwUdaGEi1lQnJs\nxJvt3JYOp/1t/CGvFK8pPGlle1qY1W3NsFG2YuYkJ8X3Yduoa0NEbhqQ7Kko5FWxWtnMcRK5A3IV\nrmNnI2bCToKQZg2+PNCpqCxVaCH+7RzTaKPcatLDdK3t2JO22V91G21WT+d6cGykCcxacVBo7mvM\nfJGoQUi+jVraqE6SAvRRswJ2zMFJgJnjFEQ4qtaLGRL1o114YKs85SVxODGAKCw1h2ojGgqpW8Ja\nutYaFzbabPy3uuMDG8lGSOrIXSERjQiqukJz30Za4d+mAXlpQGLb/SQiCp+eb5iW2j/xkBR9WO+r\nOkkmYSLqtQQnzYpZxI4z2EWwX5Q7Ub1RuCCqqYhIsVG2SrfBNhsR10Q8UrhRRKLaiCR397tYXtho\nsLWajTYt24Iq8vuc+t1TnSQ2kq0pxX4/kwaak7JKimqkpe5Go5NoOO5YDiCcPWbELIT0v1PvADge\nIYT8hkNvRKTYKD4rgZFjI159s4Xi7nKOjYi7qZ6S5a6N4j3F+838Dt1y30Z/O4XkL73na8XrrqDO\nsVHWciyZjYgo/hpkxjIbEVH8mafOsmxh/FHk1vpmlo14obrWV79i+o/3PKRxkh1kgwkyfRsRZvte\nA85ZI1x0TqJeS93jPtyJz0qajgobyVp5jNXX15U2ytiUTlybbdhGXft7Crdm7m6zypMerpFmo27L\n115E5azYvfo92KsscVfaKKO0Ubedn12c1AnGzhaqmwqfppPCsJqc+u2Hb2Oi2yQso+TomsX5Dkyc\nGQRxWzCL4PQIKNm5vo87nBPRxkY0lJMs2ayYOKlqo3QhkW4jCY/y5Td215EWY3W7d62sQkMhhbth\nkZ6xCpFZaF6umHU7ZTbaLH/tTv2+58iw0abNZ6XjSqotyk1ZTqI0BVd2oRkTL2V/tYDfHdiFZZ64\nVy6k4X3fiZ85KiLDRpmuNttssxGlkVZRy2Da6B+FS624zoixKmURZW16X0pQdVgqsO4lbS1xkmWj\nrtlrJTaqCkmdED1vY8yPLlsoh+gS9YWLVpFhw/LufX9tUqOn+gGu/Lc/d5bZh7TOIzL0DBf2D867\nf/EtN03VRtQLTKw2aGzbiN8rvm1CGarZiIjiv838Dj7htuu7Srugupe0ZB13JrXYqB3uTKrb6I/e\nj9U1eKd4TfF6c7envMEnxf8o/ud1XPn7QJz3k8Z9xMaajL/0Lav37ZXl8m/zV/zq/pWHIgBVcDUx\ne8qf/Vf8QURn4Tv9brs4KVLWPBMPaTbqFmYx1jNRzUb5Fu7rNtosudwMiTWTflnlxXX/2Ok64lFT\n19pLdkRFtaCKtMK/7tXX5A9pTAP+HDb4LKapVXt3+h0Irw2RlrUn9lR+GV1cNSwXVMFJBjSCooa5\nIh76iD/Pwyf1HiKis/BNRJmN5Kk84IWdqJLSBt9G1MdJqZby9oWNKCkub4Rzd2R0QW1HuEs8et2w\nD387J/l9To3wnOgtacBBg898Vov4X+6kamwkm+q2YCQP1dyd2KhcK3PSJh3a7wwXQ85xhObymEUy\ncwa7CIQsGGpUUbZEXSgbrtpos/y881nmJNVG3Xa+iIjCmVZBrs1u7QgpC482b31dCY8G1e3XyUv+\nWtf6S92K9/0u/S5efdUnm+geGzaSIKm00abNi9d3VQZJ6abCS30WpY2H+o1b3WBZnERDJ8m6As45\nJwFCAvthi1w8f6tpYJRsTVnIy7/ij7PwnZfe2TaSFSnN4BnTUIuQaOgky0bEWbsPCudaTYQWmXU7\n0GCjTeNroqqNns2ISmy02bHfydPCRt3yh3rUVd5QKm/w0j3wUoU/iSRHV04e6A/71TY7yknZRmSs\nWPen4cwDCpCymyilhN76X/xFeE+fypKP5JI4GJ3jvo2I6Cv+CBxmvTXZiFL5GZKgoY2IKH51Tqra\niEh3kr7KVd9YnUav2D02DVG9kEF1UmYjosHoKMtGJLm7a/3VUTiZOh4CRYbVuryfapFrpUgk3Wa1\neH3z9HX49KnbPiGPBzSWKaRZBKclmYSe439X4YVGqoiIOI/HzUJ/hujSa66NGH4cRs5D0G3fTdYN\n2n9RUOcJ4FUu814rcZITHg0ajyyW03dD05gPO4nItBFRvYaQVZeWJuYNXpLa9N/j9lBInSQq6l7i\nzGfblnndbAubV5/Me1wRtASGzPLEXWVeQhqbkUv99FFUU52HzzftCpbVpWbwvorKqrPwzVs+D59Z\nkJSFR9l2sowfGTbqXjojIiVCymyUvjUZZRTUh0eD/UyHxNprkREhOSk+0sKjTbMGIakz+JVbDrfG\njaOSHbac1NWaP7Tl/aw9MbacMvhMrO1cF0v6zyfV8+4/23n99kHGMiOk6R+RpYQe4zUR3YQnInoe\nnkKuwsu/5MRwkZ0PEhwbpQGTJNmIFBulfMSfnP2Tyrr2jJ9PONuEYlbWrp3MRpTESb6NNn9XQ0Ql\niTvHRvVdfUhm8FNnlU223OUAs5rv2q52mbo/3f++k8iwkb9xda3uk2nYWrym8LDZQ+q1hIBp5eBq\n4niokdBj//O1VERE/4aXqZfhVZpxA4lmShuVWT5ZaKlIwqMUTgO2ZPxkI10vlBYeiY02jXsnWeER\nEYXzPghT6+t2C6oyJ1U0Ztso3HYxXznfBCU22iy5KYr91Flo0wI5tYfsd/9Y67tSnTTI+13nr272\n53e3WSaNbBrjIW9hsqvdPO5JrhJnp7UBIR0QJxd33/8Qb7X+BOk9+ldkTFIbCVlXkyCBUbmc38JJ\n1mVYQlJtJJtyknXZWuykqo02jbMhsfZa5AyWKrORcoJ2CjRsIYmNuqeXxW2i1Glqr/vHRuAlQnLC\nI3aSU0mROokzdbmAr8uVur0iIwM51kntQEt7ZxbJzGWm7E6FaiB2j4jnPvll88LH4rd7pf12L8Mr\nFSEUEV2FF97CRXiiYQ+TYyPe1Fl4oeEwJstGrMYQXhsP6bPw/S/+DuG1dFJpI0pyfSqpjaTxptDc\ntlE5M0X3kr2Wz6aQb0/DdSVx56QBN4USdrKuS07aHVdCS9Jv07hWXqHvzLWin24P1fCuKOsfzDWF\nPN5qmIEzp8zYeoQsMCpVxIm7rNm/+FsNjCixUbYF0gImRoSUboRTYb6NmMvESVZ4xDbatE+cpNpI\n1iIyM296VvC5LiTeT6U8YfcUXzk6Si3Q6PtI1DmTumY3FcNVhaTeSD5v8+BtJAuSusBIKhudwOt6\nuORp81jXT5uT0pFqmw3ifLVoIKRWLPf8ifdE9NBf0f1Nfm134U5dwqgiudd+96wlJzbKuBmGSiml\njWRTqo2oEBL1sRq5vVC5w/o5GpxVeIcvwnuZSbO0R4bAKLFR9zQM52gYn+IrV0ydpNqoe4nnm/Bt\n9NbN9ac3uE3GHbsFgeHKNQdHY67VJLhRDOE6STw0+EysvW1zElF+3xMBJ65Fskwh7Z4t9UOfP/H+\nQctiiHtKFcnyrI0U1zkqSpN+aTmDZaO0fTZoyQmzbsKTWstQ9mNR76RqeJS254FH1irZfqaByHYR\nlT5FhRtRyds11k0MhkbZhYK+kNhGRKaQ0nFXqpCyHi/VSWl8VneS40Vty4y58zs4iYahEtHmXlyb\nO5ss8Qy2WiCkrn21zZ/ktyI2+jP8AaWWymzkKIrswKhcLn1Rvo3S9hJ5VJN+mZMsGxHRZXhlh5W1\nCc4q1fCIaRfSW/xVRlRkCIkovzVUudbYugl2kh8eybQUeoOb4cy2tUmSMiep9ReZk7JsoSUkid4r\nNYTD8UMDz23rpPKe9y3IRBuUaGkWXffAYplFDdUjsjTQdfKTfepTA7LwKdyksvmj/fK4gbx0F26J\n6G/8a6mIiO7CHbe/DbdU1DuoliI7PlPDrPv458LI+FERZj3G6/PwRP20rb6NuP1Z4SRnFbVxZiMi\neou/Qnjn2oGqw6SxYNmIabFRO12BQ20cleWk1EZEeeJOnZOiC/IaqgG7jRR9V2riLp921u/QGqrI\npxN28XaUKjBzv6rY4YccLjc2IkLhw0JY0dVEKqH/ikP+JTm6r4ufWqooeZz1HjmW8m2UNWasmElW\n4cZpcGMJjNd6CLdlLYMaZvHWqjZKG8sdmLaIqJzii8aIipKgyrFRCGRFVKQJ6aytkK9L8WlO4vBo\n8/SsiFSsmZbuKzMkdRNM+LXpf7xKitRJ+pwUTuLO2aYdJKVhUP1mWm1OMm9qvJrT2sJYuJAyCYl1\nREjOEtYS66dUlLykeoh6YaQbyXqPHIFVbZS25xO9byNpnDrJ6YVigaldQU4ZRYvDMif5QiK7m0od\ndNXFK66QaH91E+mKpZMyG3ULeydZNqJmIdVjo6taJUUvCdMKapXBs7dl1Umyov4uh3ESs+zz2yhm\nkcycwS5uB6uojISY0kNWG1VF1DtGJFf2MDlhluOwVGBZBs9Zq8VG0pid5IdHEvxlxXLOKmQIyY+o\nHBtJuUdjUGVN1seIjaTxjnUT2YqjhOTYqGt2UbMRq7caQrlZNSkQMN8lSR5SmVirOWkzsvUfEX8g\nh3cSkXJ3Y4KZZiKkZfYhEVGMMYTwEq4y5YiKfsd/r+HSb0BET+GStB4mWYsfPPRaEqmMJV2RH9yG\nG6mX8x12WwjMWutPvL8a5voyxEbc+FKLk9pRHcadSeSOlJK1sp4nx2EOmY32gj5Wtz9Xqjai2hzn\n3boX+p0P803ZTpK7UoWLSjGFY74uEqrFavrGWWBDc3QzcVi9YjW5bp5eFg2KwEg8RP0HhU6mWfzt\nM3DmjmT1C7+LZP9ruKQk1ikbSBvGj7qciCrLCma9QdtlBTNTpiNq/Vyik6wr2/PI3LERFdWCqhYh\nSePGoOomPKkF306Kb4uhUWTEVZsJkDQbdSvyHOdOvq7ff9VJyliosvtnmMXyiykcJ/k9VWqQJEOF\nvLLDWjyUGUiptlfrKoefefc5G1N1LP68N1P+d+odOAiphGKMfPD9jv9U2fDCsmIf5wAAIABJREFU\nl3BlNUhxbMSrP4UbSc0JqY14I//FZxYDS8UJqn7Hf2pNXRmKXcdHLgr3bcTvflNM7aLaiLo4qW4j\naSkvVde6yO7mRkSajYSWFN9jvM5OQFZ49BZ/NVbWccvBNu0Vqzb6ij++4o/yop6GNtJXr814VN6x\nN75RdjxWE4bWpjLi82DL4abbMv8r45iWd+R/8WPzT3nfD8q+jm7hMPSMXxS/unuAyb/Ne/WM3ktw\nSJYZIanZUj74MuVw6PMrvhHRe7goG3CbX/1Pk9ukdlHjqtdwmfYGORpzVPQUbmSzr+Ey66ZyQjHf\nRum7p7UG/lppY2HHiIobZ107+4qoZNIgJ1k3qm4iLeRzNNY4pWx54yjShNQywYQESZZC0iBJtVEZ\nJKWbqnZoMfpmreLDIhjqMpDWH6iGRA1xEnG/3UBFRMP5G7sdWOKZcHYsM0JSjy0OlST59hou2TQi\nG36cZuekjSzhNtLPZMVVEipxOKLu5Eu4+hXf1IiKhjaiYZzk28iJqLI9UeMki/bGWZxUbazGSSrt\nLRu3tkUo5tNNKVvrJSLucxoGSWp4FGN3ynZio/jRBxZGQCNBkhUbxbeNGMpNpa9mdDm9t6aQa7PB\nf0owREboYy1vjJOIujiJ/xF1oZJ4iANWQsw0DRZb1GDBxQ7UR0Ulv+Lba3+qcNq8GOFUys/48RLO\nVSGxjeQtnsJFOghX3fLv+O+hqLBIEYdxSyviSfkvPt8UVYIppcYEJ6hycNbywyNe62IYUVlrPcbr\n0FY30T62l0fgUsMMfuo05+W0SV0dhFQfGBfoXYGDm6lrwc/UcYEDGZ1bZflDpyIZ++WslRzIcrHX\njRQu/ignTlK2bwisdBKPFctmlD8zZpcXJyFmOj7LjJB8+Dh7V7P4/fIfTj8A0Xu4+Bk/aFjsIHBc\nxQ1+xo+X4vJSbCRwqES9VKpdWSVOREWuVxyytdIgqZria++m4iCpxUZURFTOWu08xuuz4Xlq/yV5\nxiR+HCf5XUfVS/ZqTxhvP0alGi2jMdCRyKZb68PcsnQmhcuud43/jSV1Tzjf/GM4ykn/ZQspcVK2\nnIg+4k/5N/gzlxUzzeIPWaOQqE/flU56Dxc/4hfb6Ef8UhuIjX7Gj5/xI3MSq+hnco2XOuklXJU2\nYthJjopky2qKr0Sc5NiId+ZXfCuzfOparJnG2GhUlq+l2Vja6ybESVUbZfYSslkqGhN3TPVC/CP+\nVEMBSmy0daFE2swxFifuJM+Wv2o7iWijosEqXyMSdLyEJZQWKWT5N4H1U13IWuKbTzKipbf4i/9R\nYqZZnNNnzUqFxKROYtNkgVHmJFbRz+HPMXWSBEYZ7CQ5+6s7w++uhlzZlksnWSZjJ/k24seqk1T8\nSGu7FB+vpdqrXEs00xhUZVhrWaaprkjGnEniJOemG4xzlguhuz+96iRlLJTbL6UqJw3RmqKo5uRh\nNfjznZT+y2KaFhqdxAtTJxHRR/wplzKiJSJ6jv/NV06zyEAus8puFHxgOTm67/5aV5UN8xnOWxr4\nNpK3K2v2yi1/hnMp5LPiqtdw+SN+fYczvx8r3Q0ZUeRrTFqmqGtxeV7VRmljeclfy0nWHbqQr5wA\naYs5ztMVL7Vb8YqNhPPwWZ1g4qzWLxXCsIJOe99y7gPqIzmnQ0s6k9L2m23axYfZyCrB/AO1fbYa\nq8PF9F3ZFpxF98UyixpGTZIhZQ4+FdlYA/DSBkTvIZROyiKzH/HrNXGSE3U9hXOyCyvYRrzBl8JJ\naubwV3x7KOra1bW4ZVaJrq4ldRN75yY8tQzv/RPvpRSiGlSVRRMqaSlEdUpZq2giXbGrgxiewa1b\nJrZQr5J4G5fNa0GybXq5oO0kSjyUmb4xKoqRgtZYDYl4jl11qkZ1oTUrscT02QkEftqaVafshBjj\nt5by/w5n3+GMviJ9xU8jiy+yoai32TQgoqLjqswTUuekS7JtVEVsJBtMayusfqxRtGf5HHapm1Bf\n2q5uImuZLbQK+c7Ct2MjH39Fy0aSuHMKGThxVz0lVox1YTaL0ahtO/f2qlvR6FcblZT7ij+MyReI\nv5HsHyUdQpJ/e4u/ysEJPH6uXHgVXtIlz/E/7kl9jNfyj1/iqVLCkOpfBIRlCmmLK5TSSaKi7rnm\npIFsSHFS3oAGTlJtxLCT6mFZMrLKR5xUtVFmL0EtDmypm+COt73XTYwaR+VziEK+4V2j6jmif/F3\ny7mLnVQ9ce94GmQnWdIqnSRGaX+LcNb94xXVj8hZnv6B6TAjKvTD8VC5kR2dxFqSJawlucUz/+t3\nD3JqZZkpu+3g3B0borNRxlf8DIEl8WnnJriN2YDYSQ0dV8nbZWRR12sIaeIuC48Ezt05NnIE6Wis\naiN+/DN+PITzxvI81Yjl23E+sGW+CYft1vLJJqoop4hVwyNO3BF5ybrz8Pkc/zsLL9Up+EJ4t/qQ\n/AbUh2hZV3+KjI7iw1yZZ9YyWdg0S1+SsuxsFWs50cBA6fKL8F7e8rFc2C/P03TspGyeeyLKnMSU\nl0T38Y/c1lnb59xJyO+lQEgDNv1JpY2Yr/gpA77trXzyFFrbstGh5iQ16hInWTYS3sOFqpbURlm3\nk2MjLnBo+7MGVDWm2suYb6K1buJGmwCpRDqTWgr5GidEFydtneWj3kZOg2rZejqXUnez3YZiChW2\nizm+yig9kAbV7adbEydJwCSdgo2DxiwnZUh4lMZJ6RWPOoFWNtG+OKmcgJ+K+6JZYdPeRTWL20/M\nYBePTwjBFBIRndWEFAK9RbqwnRQCfUQiovNQykMJzs42TlJslG62Ieqis3phhbTn87gjpPf+Zglq\nfZ1aHMg36agGVZ9tQpK1qG2+ib0U8mVrNd44yr+NYbriTXgqlZDZ6KoIkvT594aFEvqs53YxxXn4\ntMIspjKtXxy0ly2fh8/2ojgqPJTuvHVzrJblZSovOzZU/agLSbuv5m14UBeWO7YFY0/dENKM0Z3E\nKnqLRGT6hm3EqG3ERszQSXqqkDoneTaiipAGWx46ycnUcb9axUb9u2czt1odYJzMbMnyZU7yNdZe\nASFO8pN1o+7TIU5yZpWl5nvDl0JSY6PUSc59C/15ZlMhqbFR6aR0U86d5ikRSfnnOCnHcuZTTqPt\n4p50uTz2Qx9rudM4vf9L2YDhCzIqrp+chZRMGNY4Ln6mJ3ak7HS63F3qhrPENET0FqnMy4Vam8xG\nRPQRv/2OK+arlgbs31q2lpJv+SuqBeijeC9GXb5ocdIoMo2libt9pfgakTnOy5ccjVVt1HjPw8d4\nfa4FSRnSmeRk6mQKvmoDMvquPuLP0DspdJdkv9LVzwwnsU6sbZ4ZTkoLGbLJLy40J7n9Q+9qRYNa\nM6lqRl1eLpSeJPFQdiFVDu4ul8tC0mJ91UPX8dG6n0B2C55y3WmyzAhpX8HpxkmZjYQ0Bgq1NqWN\nhHO344pxUoXZW1+0Rl3cA+SFR31qUc3ylT6WFJ8XHvWfRhm4WBUcW2f5jpzio7Zb8TbeyVDiJL/r\n6Cq8+F0jF66QpIHvPy5wcAYOZ+OHKOnmsbZcxkmqilLa4yTx0KgQp3F5Vsvg64dszaTL1UEU2UzK\nrJ9yemVVV9kdQWdxqkeE5LGpcVBNQ30MJI/9NpaNWjgL9E8Gyte289Yadb1r4VQLio2oC5LIHUS8\nHdsV8m0dHqVkhXx7IYuTrLiK4yQi8m3kv5fcSNdpwNs/Dy+OOchVWkoWrr3FX+eGkyROKj2kjlGt\nUqbj1D98VDzEqAZS3cMXTy0hUXpDGWrTTHqXNQmPBInpy8If6mOmiWtpmRHSfglW6MNcuMaSBmQL\n6TzQcyQiujKKKcRGRHRZyxP2b1oXEm/2UnfSJjzqd3Jwl0I3eVgPj/qW6UXcFn1O5VqfbULiFdvr\nJtonQGq8k2F7nxPZQroKL9xpcRserLgh7SZRKxqyQgm1D0l2tRpmkeEtK05i1anusZxU7oYTDNEO\n8ZATBgmqk/zl8jg10FO4UW8rU4ZEqYScw1teTdPO0z/bQ0hNmE66SFziN+A2pZPOkwakOSm1EXPZ\nkCekXoR+uSBvuXBSbqN+Vze31vULCM8VJ+kVGSHwZV3LzExHKOTbrm6CtsryNd5a9yHcqkISGzGl\nk9Qb6WY1ZuWkOGX1QVZA6FcKODm9cstE9Bz/uwovToLOclL6dLvaBHV5KqHy2NjCSVkijo92yz3O\nrc7SLchjNR5KX6I5SCgFQmpFcdJF4RK/ARVOOi8a0NBJpY0YdlJL6OaHR/3W8j4nv6+rWs5OVDrJ\nLBF0s3yDtcYU8rWk+OgwhXyN94an5iliVSdlQqKhk6wRVOKk0kay2bQ4W53brTSfbMrp0BIniYrS\nN210kjo8KGWLerlSQtu5J30qjy1bVJ2USuh3/PcaLp0jM+suIs1DKPteGgMnlbKhoZOsBskpW2kg\nzb6iaSPmsmYjJw2oRV3pdON7EBINnFQdPlWfkKJvXK2boG2zfHuf45waZuRTU3zlipmTShsx7CR/\nuiM/DUhJv1S1soADozLM8qMoS4RO4k48tN04oXK5Ewkxu+TisquTsU7iB+Wkyb6T+IFzPoeQlkaQ\nEa9Enku4jdPgI3o24jZEFRs9RrppyBNS4SQ76uJ7VXg26lN89fFVtBFSxUZvkS4aZkjq2x+okI8a\n5jin5A4dwt4L+cjQmDjJshFzGx78EoxGIfk1BWwIp3PLyh/ehCdnrbL4kB/sq3NIHlfr4hqXlxm5\nsn2LkzIPvYZLdRb/zEkzzcs5oMpuBF3RnSMSInp2bcQNfBtVYRsR0WOki4Y84XOks6R+3fLcv/jt\nFKb7qCt+NEyzxLy1zZBE9UI+R36OwzZzD9bI7ruxl0K+DCuo+hPvr8ItaTMCCNVZAGTw5lV48Lum\nLsODEyGRq7QUDozEEI/x+spw0nP877J3Er+F3znUUiyXSai9Ls5ZTomHyisedb5gy0lWPMQ5utJJ\nv+O/l2Qy5cV4SFimkA4XnNaddFWb0Pcq0N9IV24Idd8HIqU8xEaM6qTt4C2fG046T3bm33DAr6Mx\njgUt0hTo0EmVOSkaaR6r69jImOP8olo3wWs9NGT50sryXWZ6vQ0PvO6N0TVVzrqWkcZe9/GP6iSp\nwrC8QkTP8b+LPu+nCsNZN1ORrLJdoXZZfj3WPUIqIUswFukqEtZIiZB1u+fSSTy1//I8JCzz9hMH\nJcZoWodl87fWgMhsIzYiovtIlyPnq7eCMw6S/E4ph/NixX/eDQkHKz4bLcuKjDfzplM5Rkt12lm5\n6UYlqOpbpmw3x7mPc9ON6opWFRYlNiKiP/G+3Fpmo/v4JxvDVGYC7+OftIjgMrymNYGP8doaBXVl\n28jiKrzIDlijspxBRUTEdy2Re5fsawb3h3DL//6Lz/yP7Anp/YnqX8IVH1H8jxf+im/WHWTYSUT0\nGi75X4xxwTYi9CFtTR4nsV3+JkvutBK7v8Mld8VG7ouv43Y4COlR+75uah1XVOuXSrd8U5QCOmtR\nw/gqfveW4VPUFwf6fU59y6xezllri6FRVJtVtqXPiYbdTr7G/EK+tL+hrK9T47A05WUVQXDhNdmZ\nwNvw4JSnZ31C2aacDi1eUdSVlbCPqpTjB+onsF3nUJaRczqBqsvTeMiaaJ/sOGnxIVEGhLQ9GyeV\npmHuhmdkv41qI+a2G8Gq24i5qdmI3/2uIQ1IQycdQkj+8KkWG9E4IdFWhXzO6YMLQNT6OnWapZZC\nPr4iVk+s5RQy6Unfn7Dcz9SpM1JnDciNdcRJarWFYxdyLdg4anUvhdpUq9Xewkn8IDsYGp0kMdPa\nzs8Q0k50EwuppmHYN5aN0gaWjZjbmo3+RHownJS9e+YkJ+riHiDHRveRbu0+pzJA7Ivf6sOnrBFL\nRQVHUyHfqLG61AnMtxE/zpxUuU+HW8iXnozKc6U6rwzfd6fa8+HbiLtPHDHwg5bKvUa7pI0dHapO\n2m88lD5NY5rd4yEuk7FuPOYfJHQYD6Hse/mEYJtGuKu1uWuwEZEpJLYRUzpJdeFdUsDte84vPefd\nLp1klRFe1cbzpsOn9pXiM8bqVu/l0SIkGjppOyGlNmJSJzldR0/hpmojnvHMytel76IGN+nd5JxA\n52/8exfuqsXoqrdanLTfeOhAwRAVUXi7kw6qIgZCWgUVJ901hFDXkZ7cfB375kGTR2oj5qHol9pO\nSDeB/kY9xUeJjWQn0xSfnzys2khapoV8flB1gLG61fvqCuykaoqPioiKNBsx7CS/kIGLhvfV7ZQ5\nKbURUzrpJjylbRwnqXdQTV+1lMkP9hsP7SUYIldFjCUk6p10BA/NCwhpD+hOYhVdRyKiJ7sD6bpf\nrjrpduibzEmljWgoJD9VSHbURa6QLo1d5VO/PwMFZ/m02WD1iS2qWb6LPudZHslqPfp57VaHYfOH\nlJexzq06toioyBYSEb2GS99GMohyX91O4qTSRow4SQKjsoFz327fWOXd7fxq7PbaBGYX96TL5XEW\ncI9yEqsIp9+MZQrp+MFp7qTUNEzpJLVNFnaovknL4coG1DvJsRHVhHSTrFv2OTnBHLlzWEiKb5SQ\nWlJ8pNVNONUWzRFVOkeDf3v4RiHR8N7wjo1YnNYEM+mKmZOsuKql24lP4qqNmLtwV20gVitbVqMo\n2m1WUycdt8fOIcs9VSchJPKBkPb5pptKtmvt3Z+GZ3mrjZy1VdlQ7yTLRtKmmifkt1PTgKU709Jz\nR0h+sk5WvM3r5SpTMak4KT6qCYnsuomiP0wGMFZuD3+m38bQvLuHNl8ZwzYiIlVImY1kFQkL/LjK\ntxH3OVm+YRuRKySqScuKovyZERpn1N69UDtdKI/LzqGxTuIHizzf7pFlztRwkm99czc/1TS8/K7P\ng1ltGMdGTNVG//XvVZK++3WkmyINWK71N3bTRvg2+mPPQJGVEd4nE0z4w6e4QctUFM/J/BG+jf71\nbcq6iXKtj/geAhFVbETKjeH9oMqyUcrP+PE0dJJqI97aQ7gk7T5vsiIR/RefH4zRTmKy6/h4pznp\nLtzJinfh1srpUa+W6l/ndBFlqK7it7AGGo+axUeWZxKyBNMIPDSWZUZIpyKEmmmI6KnW5qlmo4dA\nRF789F//0ktDnpCSOEm1UbquUw2YFl+UlX7qircNE9FaWT5rrava0Kg01BsTUZmBTjaZ+llbRBUj\nhWCVM5TX43KrUEdjr7aQnsJN1nFVrSzP4qTURrJKmY7LMod+pDUqGKKiTygt2ThcbcLWwRDBQ+OB\nkPZMxUlPgX5HerXbSAPHN78jEeltUhsxL7XOLWoT0l2g/yK92GrJii/Srh1fY1UbSUv/1h6yFhlC\nKsdUXQ4L+XyNabfW1W/Iexa4gKoiJKLSSaWNGD7N+TYSEZZqUfNRac+/pTGu5yZbHtTXfKttVGlR\nEkJtMWq1xbjpn+nn4igpp3wPF2Pdkw8kIKKpeghl3+tCblmvO+kpENHGJVYH0u9+ueWb38nC10ID\n/xWrvDR0XMnu+eGRJSSr+KJlBgprPK81i1I1y3dlD9dVB/le9nUTbUFVfifD0kbMmX5jeEptxCRO\nsmxEycnOCqqyVGE1eqD+vOx0OFHbUCdyc27sJOl8ai9YSJ823u3bcVI5xZxV/LbUeAhCWi+5k56G\nIiHNSWqbzDe/iy/rNUmUlTZiXho6rnwh3Q3TgNVSQNlhoqbxVaOE1JLio2K4rj/lRIuN+sYyosi7\nG+9bpAsjotLK03nSIMdG6dwQZVClFlNweYI/dR419Dk5N3aSRJ/jrRZjWaNWnX1rzEzKY/Uj2tpJ\nM/LQvNiPkLrOfCIa8/WMWkvij/F7t2ca92TjpNI0zOtwEJLVRs7aagNp4wiJiF4aEom8KXVAVZkG\nTM/7jpBa+pxIG8y7RSFfuVY6NMqxEZcsNkZU1M/g59uIuSgiKrvcvMVGsp00qNpi+ldKKs7VnqfS\nZFm3U5noK52UtqlW98njbJvlvslLav1h+rT6KY1yEiR0BPYgpNQrTMs2R60ljRv3diLBabfblkio\nd5JlI0pk42+EqGKjX5HeG/KEVDiptJFskzNjfm1FY58THSzF1ygkKuY4d9ZyhHSmFF9sIipn8NNX\npDNlBtjSRgw7ybeR3HVQ7UQpw6yWGgHp+7Ek4fRLWfNKiI22iIeyJc7tVrdzUnaXEzmrTOQMszx2\n/VjLcKElgBi71kmEtHtMFoIrEua11qalAdlCYhsxpZNUF74mA6p8z1VLz6lI8ZEdVPkpPtq2kI/a\nbrpBI+c4bxQSdU6qC6nfQtZPvt3cENlGWuaGkBmm/biq2udEbr9UJi0aFsu1OMkKg1K2dpIMXGVg\nnSOzByGVW6j6YNRa2+UDJ4LnJBbJj0jfbpufkT5rDYjo06ho+DVc+N6WJ+Q+p5aoy1JLum5jio/c\nLN92hXyOkC6NQcEftTnOJZIrCr63vKdGsR0ZFVstHFe7RtR0k/QGOaNxyZ2PnJJbxlnpMhmX6ndN\nqVugtmCofTLAFiftxUAImPbFTp+jFUO0hDuNaw1K1w4mpIMeT7qTXgP96BdaQhLZEOlOShtQ4aTS\nRtQmJGqOulQhWcV+LTNQbFHI59tIQr2y4Lu0EXPTYCPmahDZeDbiMcVaOYMZaflDcZNq9ay+zumL\naqwap4aTfjbPntrnpKbsGo3FFRmysDr/QqOTynuz7uuHPzaLA1QmLaR0yeGEdISrm9xJqY2Y0kmZ\nbKhwUtmAEiepNmLeax1XNDINWC09l7WoeTzvvgr50pL3VDOWkG7cKWXVvivu/vFtJFvIhuI6lRHO\nUNzscA3d3OS+jaRHRD1Hl3GDOKBcRVYUbTjzkFptrBo/8dDYsa5qPJQtWbMtZhHGTVdI/tP2zWZt\nTvWVbJxU2oj5Hg5CKmVDiZOsBtyGyLQR816z0dg0YLX0nFwhqeN5Gwv5yBDSrT0jnx8eOUGVMzeE\nda+mciN8BPo26msxyvq67er0si6lxiI9qg3F5Qd+ls9pQH2opGbktp53Tpj++feYQEhHFZI0c973\ntN9HF9SrNmK++94ySzbU+6bawBHSe/A6rrZIA1LvJD88srJ8O6b4qIioSLMRcxeI2uY4p2JKWafc\nnAvHy8mNtqsaT4v9EidVyiJiVJ2kjuhkJ1VtxFiRR3a/7XIAUDZVtm8RJ8Zq0c80T7XTGawyFyYq\npO2K99KWzNQOhRCM8Ij5rsmGiD7t2IhfpUjk5uucvqut04C85S1SfGTHVe0pPioiKtVG5ArpxhiD\nxd0/vo1kC9kU5paQuOtLm3NIHw7F+DaSx0ntuDW/APWjaloqobNQSS0W4OVSN+FHNmrkZMVDNB/9\nWEBL7UxaSBbtI5YmeCiYTvoORHwtfGbm64g635gNZLlR0eD0Xe2SBpSoa1SKj7btc3KCKt9G1k03\nVBvJWi02ku3I5EaOjWT11EmqjRh2knoYB6NWgojs25hS30HFd8HIip6dG8pRrUiPdsjgpTJLJTSp\n3y84KAu5/YRV4jJZLeVsbEREX0oYNJBNrDUgokjvQxOUNiKiH7Eb5+TYqIps+UekF7cbifkV6aE5\nxfdQm/t8C7Kbbvg2su6pUdqIiB4jXbiF5ul2niOd9U7ybcSSK92j2sglKyXnx+/hjGsiWjJ4JVkG\n72XbDF766qR/rTtz8u6DyTKDcUjyKu1wmE6qKDMPkgY2Es4GNd9U7naoNaBNnKTaKN2BrdOAatSV\n1ss5ycPqOCdp2V7IVw2PhOoc53fu5EZORxS1CUk2Re79MrI0IFE3GyzZMdNWw2xpfAYvU1FKWuFW\nvRc4M4Xf5tGYwVXyKZjHTA2Nm91uh09C56SuikE9R5xtSt102VDnJK8BEbVVUnjJOjsNaHnuux9R\n5HdlqSk+2raQ76UP9aiYJba0EVOdUlYdiuvbSGoxRhXpkeakc6OafOQwW3no3LtWhuLSmAxedYyt\n1V0kj0/+SzwtCJUyjjGXnTPi1VlLfZcljWILQQ2MUs6IyJUNEfk2ooqQ/I4rRXXDwMXbbEOfE5Hi\nJKeQj9yhUdaMfJaNyBXSnTsU1yqLyIZPEW0KzS0bpV1Kw/o600ZvcRMqpVil5L2T1NnHnVlEqyqS\nMCh76tfF4RR8EmYhv2PM9u1XMVhrNW5kR04bMNWcVBVSIPogOvdCKOLTzbkijzxVWDjJShVW04C+\nkN6LjGVjio8MIb3YM/L5NpLqwTKicqr7WmwksJaqNpIlRN3cRY6N0qfUa8mx0bAGLy1ncO6zUGrJ\niYpaVDT9s+HJOfSUMdP/Cmawiwfl5DXihpNYRb1LrP4hkrOJ2iZtQLmTtuy4ki37UVfsmrVU+lFb\nis8aq6vaiOH0nZWsy6rGiTZDYh0bcS0GUT4rqz+Zntp15BTvtdgoXU72ZK9GDZ5TgFfOOM4PnKkQ\nsoKFlJWfYcYyheTNCVm7kFJOFS0VTjobioQ034Ram7IBbZyk26h/9136pb5DvhvVSj9Z0U/xpabM\nXDJWSE/bjnzKupSkJqJlMj2iQXbOsRFnBYnyYKiawcuKyEmrepBqiGLEEmmWkl6lbKisNRMr4iGw\nCxBSzkm0lDiptBFzvgk7iBraqA24DTX0XVXzhKSkAXMb9e2bUnzRK+Qrc4mSbXNsJO9b1teNFZJa\nff7Qz9HQPnfRs1FEzlwVtXwyqqk9g8dz6zUOWrJHLFkjZ7MS8GzYEM4nYBcgJJ3jB879O1oiod4l\nToOWNlUhnRG9EV005AlJSwNaa1Ur/YwsX2kjhsfqVm0kG0nnOPc7ol61uSTGTqbnz13U2KVEtYny\nnAweaSUPRsDU3a6JiJJyO0tF6tysOI0clInUBh+BhQyM3Tt8w4tD3/aifMfdtnHe0IDPMmdmrXnX\n4E2r3ysDrw/6TurlKvV+BoMVh+N5LRv5vBZW+xXpqZ/jvDpu93ccGMi3EXvutqivs2zEy9VyBits\neu4zctUM3kVfE8EMs3OqjaTuTh589mZKa75LFUnv0RrOkidnNgP8dwYJFGEQAAAPVElEQVQRUhNH\nOxRCsLJt50R8OrisNbDaiI2Yi6LjioYNaBgnVdOA1odTrfQzsnyOjT5TC1Je8L3FUNwy9ffaXxxU\nbcQ8JSOQfBvJU6JNd5EVG2UZPEG1kTq69lwLmJJyO6Gsu5NyhhKcN8AhgJBGcBwtFU7ik0KaJyl9\ncz5sULY5L2RDiZPOtFelTXRtRK6QtL8lnd3VH/Or8mnPyOfbqCvWKGrHrY6oV7dLSa2VIHtU06gM\nnhUwWU5SbXQ+rHSwR8haKipDIpwuJsWoSm6UfS+Q4/QtJU4qTcNcJg2ooY3jG7JflTa+jd76ZtUU\nX79KN0uFE1R9dUFbWc5gVve5I5+se2o4NnK6lKxScs4KEuXBkG+pxhq8suRBUG2kziYug5aGZOV2\n2a0lCCqaKqNugzD9L/E0u2h1luw4ePaYHDpa6v9w1TTMZa2BtHF84wvpon8LJ0+YpQFbUnxUG8xr\nDNe1bESukNTU32c/8sm3kTwl2sxd5NgofSqyaYyZ7vpDfWzJAyNOygIjpqwCH2btpKghY1I/OmAx\nC9m0MGkhjZ1e6PgcNGAKIewmpEuiZ6Ir1zc8jkRtc+HmCTkyU9dqSfFZhXzBHK5bsVHsVm8s0iO7\nTu/VzeC12EgW0lYZPNK6jlQbpdPlXSY/Fn9AUtGBlM6sKlXdU/utgTVwSiHtfQLWE3LImY1K5bCK\nxCVWvk4GLVq+SUc1XhWxjvq+1RwgNaf4qHCSaiPmrMFGshHajEBybCRTP2TBkD/MlgonqTaikRk8\nVlF6j0Eic9ASK6qcvJVoM26ppC+3K1U0bDXF3xcYxXwDppMJabubkU/2gz5carFw0uVQJFQ4KdVV\n2ibVQNkgbaPaKN24b6N/fbNSS6XJLjbBjWMjnrKvtUiPOi1VbSRPpbvIsZG8lIZKvo3Sp0Td3Hqq\njZx58Ji0h8lRkZDN3XCmKSphmj8rsDUznYLoBOd3CSack/h296KdCHvfyd5JqmmYq6QoXL/vWe+b\nC7sBt6GGPCEZTqpm+ay1qGYj2UJzkR4ZQvrUwqbqMFs1s0fjM3ik3RpDHVdLxXR5TGmjBhVlfUXp\nVEAT/zWBVXGygbFZ/9AsTNPI3kex9fJ2RPLs6kra+DaqktruQot1MpP9Izpvy/KRMVY3m0jpg4i6\nobgVG/XpOKI8GLLCph/9vZTKcgZLVGREQqqNpFTvLjGTZSN1QnEOki6D8hJjqyhLzb2Hi2X81kA7\nszjHnixCIrd/aNYRUsoed9h1kkyx7Dd4JLqxYyxuQEQ3RgeSnwbcLst3MdyIP8ms4BfpxXyJdBdZ\nMVM2AZKfwcssJaGSk75T58ETsq6jzEY3RVFD1ofkKKpgRr8dsDtZFmr6588JdcmkH9ZihER77V7q\nN5UVI1AvEiK6KRpwm8fkaemkrAENneSnCt8abER2sq4Ms/xJZqmPurjMr3DP2C4ltXivnABCGFuD\nZ9komzicKQOjG22SobK4Lqv5lumCeub1kwG7UD3hTPksCiEdlb3sfxIqlSJhbpIGpLW5GdbgWRth\ntu5zyir9qCjkc7qUfBulT+XDdG3EI23bh9nKtA4tNXjpKFpG6uuqNqJeSJzBy6gGRpmKhsz3lwLG\nMraKYZrnUgjpBOz+V/QHnyoS5qatQUsbJw34mDTLtKSWV0hERbXRuGQU6akLGctGxZR9ksFzbCSI\nlpzSBms5kRIJlTZSR8g+xk0hw0jm/gMBjUx23oCtwWzfJ4CTubtoqTY1+I39Uso90a27BfbNlaEW\nMRk/uBxm+VSNPduBEZNm/8pgyAqbLEJRKPFF1DumxUbU2+tV+7SdgEmmJkr7iqzAKOUuqWgQJ6lR\n0TAkWsb5CFRZnoRSTln27S+f1zikrdlFS1qcxCK575/eFg3UNmWDckmWfNsuy3fVTx5BWu8RaWGW\noNroPC8KH2TwrFtsxL5BVkSu9TOlpXqf/bnAD5jKyVtVrMCIUVVUsLCfA7CY6biisUxFSC3zMiwj\nX6eyy9EWQujL5yjRTMrtMLdWthFvOYm+9hQf2RUQ1twQ1bKIqo3ShbHBRkJ/oh87aImKTN2rdnuL\nF+1+tbuxyJ8AKMmyIGv43ic0l50/s6rVZmFsJ90QgqEigZVTbbN1n1NmshutFFDN+zGOjf71D6gt\ngyddSuqoJqfqQZtZXLVRev/AFDUwUudrEJwZgwpGHRXLyyWshGVn5Kqc7KhtlP8Kv57toqV+LSv6\nEawQivouJcs3ThqwTPFRUhChdkExjpAuteVbZPCkiJyMuVxpmOVL8G0kvIcug5dRBkaWim4V/TA7\nVr6s5Cczd1YYCVngkJ0uW0RLw1BJjYpuh0vKNplvnCyfpAGtuKpapJf1S6U1EX4Gj/ZRgxcqARPR\n4M4XVNze4l0rsVPLH8aAn+RKWEm30CggpKkzVkvJ1ZaVoEtjpmobJ8unVkwIN8OOq90zeFnMlCbx\nzpPHKedFm6ScwQqYslkbmHYVpbf+GwN+iWsAwZAPyr6nztiZ8dL2NVqM5TTg1dWSPCqKwkVCTgZP\nKCOkcgk/tVRERcDEj+Vj0Uob1EJwDpKy7iInKhrjoYmUFIGDssJ+h61BhDQbxh7WRa9SGfRY+Tpn\niZUGpN0yeGXMJEk8tTOJ4ZfeiskdzrXGH1qvEnUledbsQfsGP7f1AA9tAYQ0S9qvc6eRwUtjJjWD\n51c9WCUPcvemCundFsbyX+z27SVcKS/zV+DGo/iJrQdIaEcgpBkzXku71OBZDTJXjarBY6oxU1ny\nkMB3WxjegO5H/Crvheog1knRDdQAflOrAhLaIxDSvBlVqJNoyYp40oyc2qZsUHWYEzMJo0oeiN7i\n4C7dLfdFLUOZj5jNvsNa2sJDC/gRoW9pLKiROwQQ0kLYR7QkVEfR+hm8LOrySx7ShUwpp4S34g+8\nKCSk3kG1yr+IGwsBn/UEQ1ZVlD99we6fCarsFkJ7Md6wDK8seSA7AKqWRThT5wlWyYOQSOi5QT8t\nbbJmMg9COnGca6Nln32Aw3oktAXljb93vbcOPuLl0Z5M2HfJQ0t2jopRt5TMuKrNrNM1GU6ow83K\nGwhtdnb4t9uzIZTgRyGsdrqHlWfkWi5tDzHd6EqPtjWwVfdSSSmPrIdJbfPozn1XKiq5ox0jWnKU\nU3VSAzj+q6zHSQiGhJYv/RA3ZFjLobZmGi9b+EgqzJTWapdk1XR+G6sBEWmTkFoNtp0tG4c6KIGE\nSuSM4Xw4B7qHKoS0FvwDxbquaZvxQcjmpquR3fZ7Z3AwgxYwf4+P86tP71dHEBLYkf3dPf1I4Pic\nJnMsE195t1A76gdV3kCVICSwCyu8vxRYOQiG9kgqmwMJCWXfy8e5MJzjdS4APugWmi8Q0vJxfpNj\npxIHwOeEJXnw0AKAkEBXTgMtgd3hY+loRxEktDAgJEDU0FcJQCOHPnggoUPTOOFL+/RC7UBIYACS\neGCaoEbutKjnhCwa3r0EF1V2wASnALAvtsvjIRg6Fapa/JlVrTbj3hdfM/CBlsBeaHQSJDQRGivm\n9/t9QUigFeTxwIGAhAADIYFxQEtgXyD4BhkoagDjQNUD2IUsGDryTFRg4iBCAtsDLYEWqr0R67nD\nBfDBcQB2BYkXoLLFHblwCK0cCAnsDZxTAMoTwC5ASHNlmr98REsrZJqHIpgjENIsmf6NJBAtLR54\nCOwdCGl+lOf6yZ79J7tjYDsgIXBQUPY9S7JzwWTLZ1EjvgAgIXA0IKSZ4YhnsrWzqZYIJ7WZgO8L\nHB8ICRwJ3OFi+iAYAqcFQgLHBnm8SQEJgekAIYHTAC2dEEgITBMICZwSaOmYoFsITBwICZweVD0c\nDgRDYEZASGAqoOphXzTeWg2AqTHRQmHgoJZ3T7bmezsQLW0BPjQwdxAhzZJMP9McFbsLiJYaWVVG\nblV/7DpZ1GX1epj+XHb7BVpKWed5eW3H/DqBkObKCs9KK9fSCr9xYUbzN4JdgJDAzFhVT8maJZSy\nhn5TQOhDArNj8d1LkFDGHOdvBNsBIYG5srBBtauK/ABQgZDAvJm1lhAMAZACIYElMCMtQUIAWEBI\nYDlMVkuQEAAtQEhgaUxnZrwp7AMAMwJCAsvkVMV4CIYA2BoICSycI+TxIKGDEmO0Kr/xaS8MCAms\ngkNoCRm5Y7L4+RsBYaYGsEJ20RKCoVOBuezWAIQEVkp7fAMJTQR8EYsHQgJrxwqYkJED4MhASAAQ\nISMEwARAUQMAuOc3AJMAQgLrxUrKTXCuBwDWAFJ2AOhASwAcGQgJAADAJPjfqXcAAAAAIIKQAAAA\nTAQICQAAwCSAkAAAAEwCCAkAAMAkgJAAAABMAgyM3Ya9zHKWDnPxh7yMHRCDATQAgDmCCGk0e7wR\ny4GcARUBAOYIhLQlMcbdw6ODgjuYAQDmBYR0Mg4axyBIAgDMDvQhjUPCjn3ddXS/q5eThMJMAIC5\nACGNI8aIVNhJsD5258Z66qvgmKD8B4wCQtqSHQ/iXVa3TsHZ8nW6M/uTESOekFmU/6zwNzJlIKRj\n4/wAtvttVK8Zl3RG9v+W8qMIIRztE0AMp7LjH3ic8p/FfwtzAUUN82ZV13ctZ40yTDzY7mxJGcOd\nak/mAsp/1gMipAnhZ8a3WHExZLUkjHqbV2v1o31EU47hjgnKf8AWIEKaMWvrki1ji6mFF8uI4fbC\nUv8ucFAQIc2VtdmIyWKLE+5JxoxiuCOD8h/QDoQ0S9ZmI/Uv5fPIpE7lav/QdHZvXqD8Z4VASGBa\nZOea2Z0jJhvDrRx8F7MAQjo2u+cH1OyQbLy6BByCucRwywDlP0sFQtoPY0eWHOEkNdNLwumfMuYe\nw60Q5E7nAoQ0GrWnNEvUjD30/fYrqctyaDmhOKHnCj8xIMBGMwJCOgFHq+pZ9o9QPdGUFwd7f99l\nf6oLAzaaFxDSHtguPDpo1m6m+ToL6Ywpl1ebHXPeIMRwAOwCBsbuk3bHHOEEtLBzXGO9xqTSm1YM\n5zwFwu7fXVr+k3GItwO7gwhpb2wR8RwoSFpqWdd0ZK++6cRjuEkxkUSCvMVBtw/aWeaZ6/gs1QFg\nFI2mWdVs3xnytzcWoza23H2X1vZFTBOcRvcAbATAKKYWIeH3OxHQhwQAmDQHTakhXzcp0Ie0H9Ar\nAMAhOMIYCfxapwOEtAdwQANwUFD+sxIgJADApDmcM2CjqYE+JAAAAJMAQgIAADAJICQAAACTAH16\nAAAAJgEiJAAAAJMAQgIAADAJICQAAACTAEICAAAwCSAkAAAAk+D/Ac100ENaNnWXAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Hfun = matlabFunction((abs(Hn{2})));\n",
    "freq1 = [linspace(-50, 0, 20) linspace(0, 50, 20)];\n",
    "freq2 = [linspace(-50, 0, 20) linspace(0, 50, 20)];\n",
    "[F1,F2] = meshgrid(freq1,freq2);\n",
    "surf(F1, F2, Hfun(F1,F2));\n",
    "xlabel('f_1 (Hz)'); ylabel('f_2 (Hz)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nonlinear Output Frequency Response Function (NOFRF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The frequency response of the output signal $y$ of a system to a given input $u$ is:\n",
    "\n",
    "$Y(j\\omega) = \\sum\\limits_{n=1}^N Y_n(j\\omega)$\n",
    "\n",
    "where:\n",
    "\n",
    "$Y_n(j\\omega)=\\frac{1/\\sqrt{n}}{(2\\pi)^{n-1}}\\int\\limits_{\\omega=\\omega_1+...+\\omega_n}H_n(j\\omega,...,j\\omega_n)\\prod\\limits_{i=1}^nU(j\\omega_i)d\\sigma_{\\omega n}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example of the use of NOFRF, below the expressions above will be used to predict what will be the output of the system to the input signal $u(t) = cos(2\\pi 8t)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t = 0:1/Fs:10000/Fs-1/Fs;\n",
    "u = cos(2*pi*8*t);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we compute the FFT of $u(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AETFRsNaDa2XwAAACJ0RVh0Q3JlYXRpb24gVGltZQAxOS1KYW4tMjAxNiAxOToyNzox\nM/LpNcYAAAAkdEVYdFNvZnR3YXJlAE1BVExBQiwgVGhlIE1hdGhXb3JrcywgSW5jLjxY3RgAAA0g\nSURBVHic7d3Rcqs2GIVRq5P3f2V6wRwPAwQEtuMtaa2LzjlxaA1N8uUHGZdpmh4A8G3/ffsJAMDj\nIUgAhBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACL8fPsJrJVSnn+epunS\n5y/VbAtAjqwJaVWX32IDQH+CJqQ5P8vJppRSSqmZdcxDAK3LmpBWXanMjBoBdCBlQjo4O3cwJD23\nqr/y5DQgMLLk3+BTgvSK3StPBwc9+f/Hp1WeAu3V4Lv/GP4IDL77j/jfyHsI0mNz5emLzwSAe7Ku\nIV01/bP64GOMLA2wi8BA2g4SAN0QpLEMfgJ98N1/DH8EBt/9fIIEQISUIN1bFDe/cvYzzwiAP5US\npNlbbh10uuwbgEBBy76nadpOPKuurGKzu8l2KwDyZU1IN24dtP0cNQJoUdCENDvOye6jCgTQgawJ\nCYBhCRIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACI\nIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARB\nAiCCIAEQQZAAiPDz7SewVkp5/nmapnub39gQgO/KmpCWNdr+9ermADQkaELaDjellFKKcQdgBFkT\n0qo9l1JkPAJoWsqEdJCTmiHpOV3VZOnFy1QArWjrN/WUIL3i6kIGEQIGsfpxF96nrFN2N1hWB9CH\n5oMEQB/aDpLxCKAb/VxD2n5EqAAa0vaEBEA3UiakgxXbB4PO9iGzEUCjsiakF28dBEC7Uiakx78h\naRWh3UX0BiCA/mRNSK/cOgiApgVNSLPjCJ0mSsMAGpU1IQEwLEECIIIgARBBkACIIEgARBAkACII\nEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABE\nECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAA\niCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIP99+AmullOef\np2m6sdWlDQEIkTUhrbqy+utvm2w/rWZDAKIETUhzRZbDzRybmnFntdX8T3MSQEOyJqRVQk6Lsm1Y\nzVYABEqZkA5Osr191rl3mQqgOW1dv0gJ0j27Odkdm063AujP6sddeJ/aDtKSuQegaVnXkN4l/LcA\nALb6CdL0z/xXTQJoSz9BenK+DqBFbQdp91WxALQoJUgWxQEMLiVIsxu3Drq9FQBRgpZ9T9O0PQW3\nu4j++cF5k8dehMxVAG3JmpBu3ARoubLu0oYARAmakGbHLfntUQUCaF3WhATAsGonpBdXCphgADhm\nQgIgQu2EZMQB4KNMSABEECQAIry07Ht3pYOTewDccCdIxyvuno8qEwD1rgWp8l1Zl7fzkaWPKuXh\nAAN9uPw6pMrb+Ty3KqVoEgCnLkxIN7oyb6JJnzDPRu5sDnSjdpXdK0VRIwBOWfYNQITaIC3fqcgb\nhwPwdiYkACLcX/Z9PCS5bgTAJX+xqAEATt1Z9m0ZNwBvd2FRw/PPV2tkBQQApy5MSDduBSRFAFS6\n9gZ9ywXfp/eyO/00AHi6tsrumaVHxfQjRQDUu/P2E8vVDQePAkC9l96gT3sAeBd3agAggiABEEGQ\nAIhw+R1jj7mqBMA9b56QvBIWgHuuvTD21PzKWXMSAFe9eUJavnIWAOpZ1ABABEECIIIgARBBkACI\nIEgARPjI65As+wbgqjffqQEA7nnp7Se2zEYA3PPmOzUAwD0WNQAQQZAAiCBIAET4o1V2LkEBcMyE\nBEAEq+wAiGBCAiDCm18Y+7rlxar6sWx1ics8B9CcrAlp1ZWalRTzm6Yf/3sAyBc0IW1vzDrH5mDc\neYZntdX8T3MSQEOyJqRVQiqLcm8rAKKkTEgHJ9nePuvcu0wF0Jy2rl+kBOmeezkRIWAQqx934X3K\nOmX3Ft4kEKBFvQUpvP8A/KbtU3ZLuyvuAGhFD0GySAGgA80HyWAE0IeUa0gHObn6wlgAWpQSpNmN\nWwc91AigC0Gn7KZp2t6YbncRfeXKeqECaEjWhOQmQADDCpqQZscRUiyAXmVNSAAMS5AAiCBIAEQQ\nJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACI\nIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARB\nAiDCz7efwFop5fnnaZpubHt1KwASZE1Iyxpt/wpAx4ImpO18U0oppVROPOoF0LSsCWnVnvoUqRFA\n61ImpIOinA5Jz0crs/TKZSqAhrT1y3pKkP6SCAGDWP24C+9T1ik7AIYlSABEECQAIggSABEECYAI\nKUE6WPlmURzACFKCNHPrIIBhBb0OaZqm7T0XdhfRm5kA+pM1Id27dRAAHQiakGaVdwm68SgAybIm\nJACGJUgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIg\ngiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIE\nQARBAiCCIAEQQZDaNk2PUr79JADeQZAAiCBITSrlMU3ffhIAbyVIAEQQJAAiCBIAEQQJgAiCBEAE\nQQIggiABEEGQ2uNFSECXBAmACIIEQARBasz2fJ37qwJ9ECQAIghSS35bzjAPSTVzUhl7mBp89x/D\nH4HBdz/fz7efwHssv86mHpegzft3sGfzQ6efBhCrhyCtfusppfTRpOVuVe6QLAHtaj5Ic42WBSql\nZDbp6tmC23uwzNL2QS9jAjI1H6TH5hzdNE3HZ4q/dRr5jzOw+5+bUz3wifSR9302+BEYfPfTtR2k\ng/AcDknf+ZIM+U5wXRfI1HaQbgg8lQfAw7JvAEIIEgARBAmACIIEQARBAiBC20E6WDJnNR1AW9oO\n0mx766BvPRMAbku8xc5V2wJ1sFMAo+lhQtreOuhbzwSA23qYkADoQA8TEgAdECQAIggSABEECYAI\nggRAhIHeD2n5cqVB1hZu3999++isswOyemna7t51vPsPR+Cfg2+BXnf/tzsDbPcx8AiMMiG5m8NK\nrweklLLdl9OPdLP7D0dg4fgdpSs/s1eZR2CI1yFtf0s6Hh368PwK++03o/4OyO4ubz/Y6+7Pftu7\nxzBH4Om3b4G+d79mX2KPwCgT0lB3c9j9HXml4wNSs2u97v7uj5WhjsDT8bdA37tfszuZR6D/INWP\n7d2Y/tl9dMADsjT47j/GOAKVl47qH2rFcxfKwm+fc7D5twy0qIHuhfyW90UHP4LHOTij7e/W7vWh\nJg6IING5hr4b3ytwDdUfGPZ/98ruRcR8/Z+yY2QNfSt+lOMwiN3T9fNfm/gaMCHRp4NFhoNYrakr\npf8ltcaj1gkSvRnzVNWBaZqa+O34XX67jO+LIZ8g0ZXBByM/eWla/9eQDr45x/y+7fiA1NSo492v\n1PERmDaWH3/++WDzv3iWn3TjNYiVD/2N/oM0y7xPxhd1fEBqvqk63v1H3d71fQRODbX7u3Nz5hHo\n/zrnbHu4R9jxSy8P7OCA1L84v8vdn1XeW7PjI7D027dAx7u/+wXQyhEYZULKvE/GFw1+QDre/d2b\ndJx+pKcjUKPj3a/5v7/9YMgRGGVCAiDcKBMSAOEECYAIggRABEECIIIgARBBkACIIEgARBAkACII\nEgARBAmACIIEQARBAiCCIAEQQZAAiPDz7ScAbTh4f/Tlu8AdvCni6aOvfz40zYQE507f4PlDzZAi\nhiJIUGv3zVhPW/W6P/hPQAJBgld9dI4xJDEO15DgxHNA2V7ReXF2Odh81aFSijLRPRMSnFAC+Bsm\nJKi1W6ZXcnW8YG/5aS4jMQJBgpsOInGvH6dLxs1q9M0pO4hgBgITErzf8Qtjb2wIIzAhwfe5IwM8\nBAm+To1gJkjwTWoET4IEQASLGuCm118etLoHxOpffvoR6IwJCV7i5qrwLiYkOFc/nRx/5vYWDPef\nE3THhAT3/VlRpIsRCBK86qOn1JyvYxyCBC/5g9nFeMQgBAne4ENzjBuqMhRf7gBEMCEBEOF/Z+Cd\nzvGmi8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fres = 0.1; % the frequency resolution of the NOFRF\n",
    "fmin = 0; % lower frequency limit of the NOFRF\n",
    "fmax = 50; % upper frequency limit of the NOFRF\n",
    "\n",
    "\n",
    "U = computeSignalFFT(u', Fs, fres); f_u = -Fs/2:fres:Fs/2;\n",
    "plot(f_u, abs(U));xlim([0 50]); xlabel('f(Hz)');ylabel('|U(f)|')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AETFSog67TtmAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxOS1KYW4tMjAxNiAxOTo0Mjoz\nMg9GYkAAAAAkdEVYdFNvZnR3YXJlAE1BVExBQiwgVGhlIE1hdGhXb3JrcywgSW5jLjxY3RgAAA0e\nSURBVHic7d3RcuK4AkVRdCv//8u6D9RQtG2MIHZ8JK31MNUNMWO5EzbCwim11hsAXO1/V+8AANxu\nggRACEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIA\nEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIMLPSY9bSnn8udb6xbav\ntnp+5C8eHIBMp8yQFs1Y/PU3D7t+qKMeHIBrlcNnGOv5zf6MZ3Pz9ddv3v7qiwHozikzpEUe2lP0\ndrrz3SMDkO/gc0g7RSnlzWzsce+Bb/Ed8jgAw0h+HX/WoobDfXcQN7d6m8axGb7hX70Xl5l8+Lf4\nl+l9L/v+6OwUAMk6DlJ46gH4SDdv2T2zuA5gPJ0F6Teft/39hmMw/Kt34UqGf/UusKenIJkYAQzs\n4HNIO6n4ZUXUCGBsPV066KZGAOM6/i27Wuv6mgub1wH6oi6v2iZUAL0LunQQADM7a1FD41WC2u9V\nNYCxdfzBWABGIkgARBCkvrl8EjAMQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIg\nARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAE\nQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARfi78\nf5dSHn+utX6x7adbARDrshnSc43WfwVgNtfMkNbzm1JKKaVxxqNeAOO5bIa0aE97itQIYEgXzJB2\nivJ2kvS4tzFLvzlNBdC7vl7BX7mo4Q+IEDCzxXNgeJ8s+wYggiABEEGQAIggSABEECQAIlwQpJ2V\nbxbFAUzLpYMAiHDN55BqretrLmyulzdnAphEZ5cOAmBUV16pofEqQV/cC0B3rLIDIIIgARBBkACI\nIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARB\nAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBECEn5Met5Ty+HOt9aitnr/g0wcHINkpM6RFNtYV\n+W6rzcdpfHAAwh0/Q7oX4nniUkoppexPZd5u9QjP4mvu/zVPAujdKTOkRR4aa9Gy1XePDEC+g2dI\nO2+g7cxjvtvq0/1RL2A2fZ3UOGtRw+Fqreu3/tZv9K23+oudA4i0eA4M71M3Qbo9NWlx41X7A8CB\nevocklV2AAPrZoa0+e6cVXYAw+hjhvTqXJEOAQyjjyABMLyDg/TdgjfL5ADo49JB9yw1XkwIgB4d\nv6ihZXH2+pxQ45LuzQKZRQEMoJtLB9VaNxc1qBHAGM5a9r3fiVf3vq2L/ACMyio7ACIIEgARBAmA\nCIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggS\nABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQ\nJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSN0r5eo9ADiCIAEQQZA6Vsqt1qt3AuAgPyc9bnl6\nI6k2P2u2bFX+fYuq/cEBSHbKDGnRjNJ2lqNlq/WNjQ8OQLjjZ0j3QjxPXEoppZT9qUzLVptfc9yO\nA3ClU2ZIi/Y0vqv2xVa1Vm/ZAYzh4BnSzpRlZ5LUstV6evTp/kgXMJu+3kY6a1HDqR6H+G1jRAiY\n2eI5MLxPnQVpc+GD6gAMoLMg3SxqABhUZx+M3Vz4IEsAA+gsSACMSpAAiHBwkHbWFxx+FwAj6ebS\nQZuni6yyAxjG8avsaq33q/4sbnz+6zokLVttfg0AY+js0kHrrzE9AhjDWZ9D2u/Eq3sbr1/35T4B\nEMwqOwAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBB\nkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIg\ngiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACL8nPS4pZTHn2ut\nZ2x1/+L2Bwcg2SkzpOeurP96yFaNjwlAL46fIa0nLqWUUsr+VOa7rQAYxikzpEVFGqPSvpXpEcB4\nDp4h7aRiZ7rz0VaPuVRLlr47lQUwhr5evp+1qOEkny5kECFgZpsv6GP1tOzbsjqAgfUUJAAG1k2Q\nTI8AxtblOaT1LUIF0LtuZkgAjO3gGdLOauydSUzLVuvNzY0ARtLrpYMAGMzx55Du051FTl59uPWj\nrQAYWJeXDgJgPGetstvPyat7P4qQYgGMxCo7ACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQ\nJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEmDFI//6edAAizBgkAAIJEgARBAmACIIEQARB\nAiCCIAEQQZBIYTk+TE6QAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACII\nEgARBAmACIIEQARBAiCCIAEQQZAAiPBz0uOWp1//WWs9cKvy7y8WbX9wAJKdMkNaNKO0/W7qt1uV\nUjZv/HwHAYhz/AzpXojnics9JPtTmbdbPcKz+Jr7f82TAHp3ygxpkYfGWrRs9d0jA5Dv4BnSzhto\nO/OY77b6dH/UC5hNXyc1zlrUcLjvciJCwMwWz4Hhfep72ff6zBMAneo4SOGpB+Aj3bxl92xzxR0A\nXessSBYpAIyqpyCZGAEM7OBzSDup+OVdagQwtp4uHXRTI4BxHf+WXa11fdG5zbXwzze2bHV73Tah\nAuhdZ5cOAmBUZy1q2M/Jq3u/O8/EMEq5+XeGaXX8wVgARiJIRDA3AmYMUq03Vx0CSDNjkAAIJEgA\nRBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQ\nAIggSABEECQAIggSABEECYAIgtS9Wm+lXL0TAL8mSABEECQ4njkrfEGQelXKrdardwLgOIIEQARB\nAiCCIAEQQZAAiCBIAEQQJDiFld/wKUEagYs1RLEiH74jSARRVpiZIAEQQZC4nve4gJsgdcozODAe\nQerPZo2cfQF6N2OQSs/P3Dtzo8YmhQ//7LKGD/9shn/1LrDn5+od+Njzt1Sd5n2rx6D3R/x4Nu/o\nwKwTex9FR0PYNMYo4C91FqTFC5xSyqhNWrySax/l/Su7yNL+Tno2h9n0FKR7jZ4LVEpJbtJv3h74\n5Zies7S+89Pn+t+3YXNPdh5zvf+p/8j/WCS2xwkrXKinIN1W79HVWr9+U/gP3ky+/GlocwfuCf9o\n+L8/r/PdoXje6rh/r8/G/tlDby02uWVdRujE4f/S5T8vXC53erGwnh7t335zAhNgJfk5v7MZ0keS\njzsACzMu+wYgkCABEEGQAIggSABEECQAInQTpJ0lc1bTAQygmyDdrS8ddNWeAHCsbj4Ye7cuUF/7\nD8Arnc2Q1pcOumpPADhWZzMkAEbV2QwJgFEJEgARBAmACIIEQARBAiDCyL8Pae35Y0yTLC/c+QWG\nt6EPyOIja29/haPhn75PV2j8BZ4jDf/V5QJe/XbTV/deYqIZkqs8LIx6QEop67G8vcXwz92tK+wM\naobh78s8ArN8Dmn9Qml/6jCGxzdZy69+H+OAbA55feOow7+9HtptjuE/e/X9P/bwW8YSewQmmiFN\ndZWHzZfJCwMfkJahDTn8zaeVeYb/bP/7f+zhtwwn8whMEaT2mfsw6n82753wgDwz/C/u6kvjqaP2\nu3rxGEJ58uprdja/ylyLGhheyAu9q+w8Bc9zZGYb79rm+aEuDoggMb6OfiAPFLiG6g/M+W+9tnkS\nMd8Ub9kxs45+Gs/jIExi8736+1+7+B4wQ2JYO4sMZ7BYU1fK+EtqTY96J0gMaM53q16ptXbx6vgo\nr07j+07IJ0iMZuaJkWdeujbFOaSdn885f3QHPiAtNRp4+C0GHn5deb798eedzf9iL8/0xQcQG+/6\nG1ME6S7zUhkXGviAtPxcDTz8lqENPPwWUw1/c96ceQTGP8/5sD7iM4z9o08IDnBA2j+fP+Twb83X\n1hx1+Auvvv8HHv7mN0AvR2CiGVLmpTIuNPkBGXX4m1foeHvLMMNvNPDwW/711zeGHIGJZkgAJJto\nhgRAMkECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAET4P01tEaqF\nUDd+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "% Initialize variables for the NOFRF computation\n",
    "fv = -Fs/2:fres:Fs/2;\n",
    "f_out = fmin:fres:fmax;\n",
    "f_inputMin = 4;\n",
    "f_inputMax = 12;\n",
    "NOFRF = zeros(size(f_out));\n",
    "% NOFRF computing\n",
    "for i = 1:3\n",
    "       if  logical(Hn{i} ~= 0)\n",
    "           HnFunction = matlabFunction(Hn{i});\n",
    "           for j = 1:length(f_out)\n",
    "               if i == 1\n",
    "                   validFrequencyIndices = abs(fv-f_out(j))<=1e-3 & f_out(j)>=f_inputMin(1) & f_out(j) <= f_inputMax(1);\n",
    "                   if ~isempty(U(validFrequencyIndices))\n",
    "                        NOFRF(j) = 2 * HnFunction(f_out(j)) * U(validFrequencyIndices);\n",
    "                   end\n",
    "               else\n",
    "                   NOFRF(j) = NOFRF(j) +  2 * computeDegreeNOFRF(HnFunction, U, Fs, i, f_out(j), fres, f_inputMin, f_inputMax);\n",
    "               end\n",
    "           end \n",
    "       end\n",
    "end \n",
    "plot(f_out, abs(NOFRF));xlim([-1 50])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AETFSorfGY0EAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxOS1KYW4tMjAxNiAxOTo0Mjo0\nMzcAxBEAAAAkdEVYdFNvZnR3YXJlAE1BVExBQiwgVGhlIE1hdGhXb3JrcywgSW5jLjxY3RgAAAz3\nSURBVHic7d3rburIFoVRqpX3f2X3D7QRso1vYHtW1Rg6OkqHkE2RhI9lG1OGYXgAwN3+u/sGAMDj\nIUgAhBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACL83X0DjiilvD4ehuHA\ntXZdEYALlOoel0ddeWxIy/QqG68IwGUqC9IzLe+3efqZ2avMXmv5igBcqb59SKOEbCzKsWsBcJma\n9iF92vL2vOiHjVn4hwCqlvx0vKYgHXPs3p+91m+zVwVL7kSHq+5wyY/4Z9v1bbL7idU9TwBcrMcg\nhT9HAOhT+5vs3jm4DiBWL0E69lrakQ4zZsmd6HDVHS45XxdBMhgB5KtpH9JCTva+MBaANDUF6Wl0\nSMLGIxTUCCBcZZvshmEopSyfJnX2kG6nswMIV9+E5CRAAE2qbEJ6Wo6QYgHUqL4JCYAmCRIAEQQJ\ngAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACII\nEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABE\nECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAA\niCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIg\nARBBkACI8Hf3DRgrpbw+HobhwHWn13r/nu/2fn8AzpM1IY3K8SkkALQnaEKazjellFLKxjlmtV7m\nIYBkWRPSqBnbU6RGALVLmZAWirI6JL0unf0mr0/u2jv1za4sgBB17fhICdLZZvdOLZRGhIAGjB7K\nwvvUS5Aek71TN94SAKay9iGdYfhn9MmHLAEkaT9IAFRBkACIIEgAREgJ0nkHvG15lRIAt0sJ0tNl\npw5aPewbgIsFHfY9DMN0mpk9iH5XSGa/7d5vAsDZsiakY6cO2vttf/idAfiVoAnpaeNZgvZeqkAA\n4bImJAC6JUgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAE\nQQIggiABEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmA\nCIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggS\nABEECYAIggRABEHqSyl33wKADwQJgAiCBEAEQQIgwt/Gryvf7XwYhuGbqwPQPBMSABG2TkhGHABO\nZUICIIIgARBh6ya7WbNHOti4B8ABR4K0fMTd61JlAmC7fUF6T9FCb55f9vx/WQJgi92vQ9oSmNfX\nlFJKKZoEwKodE9KBrjyvokkArNp6lN03RVEjAFY57BuACFuD9NwbNP0YAH7ChARAhOOHfS8PSafu\nN9p49PnCde3WAkhzxUENvzUKoY2HAG04ctj3jYdxT+ebXS91Ui+AWDsOanh9fGwr2a+M/vXtKVIj\ngGQ7JqQDe19+24CF77Y6JL2Pd3v/rZzNlQC71PVEfN8b9L3PGavnslv9smSV3myAd6OHsvA+7TvK\n7pWlx4aFeUwHYLsjbz+xvPlLhwA44Ks36NMeAH7FmRoAiCBIAEQ48jqkuyxsIbTxEKB2OyakhCY9\nnDoIoFFHTq564zgyDMP0nAuzB9qbmQDqsuPkqu9He984lxw7dRAA4b56YexdMdh4lqADlwJwlyNH\n2YWMSgC05OALY1fPIWQQAWAXr0MCIMLBCamB83kDEGV3kKQIgDPs22SnRgCcZOuEJEUAnOrI65AA\n4Od2bLJTIwDOs+PUQafeDgA653VIAEQ4clDDAQYsAJaZkACIsHVCMuIAcCoTEgARBAmACIIEQARB\nAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQAIggSABEECQAIggSABEECYAI\nggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIA\nEQQJgAiCBEAEQepOKXffAoA5ggRABEHqSCmPYbj7RgB8IEgARBAkACIIEgARBAmACIIEQARBAiCC\nIAEQQZAAiCBIAEQQJAAiCBIAEf7uvgFj5e1k1MPmM68tX6t8OMH19u8PwNmyJqRROT6F5CfXAiBK\n0IT0DMn71FJKKaUszzHbr2UeAkiWNSGNmrExIVuupUYA4VImpIXtbAtD0pZrvb5m196pY7uyAKLU\ntQsjJUhnm93PtFAaEQIaMHooC+9TL0F6TPYz3XhLAJjK2od0huGf0ScfsgSQpP0gAVAFQQIggiAB\nECElSMcOeNty0fN1st/cNgAukBKkp8tOHbR62DcAFws67HsYhuk0M3sQ/fsnV681+wXT7wzAvbIm\npJNOHeS8dgD5giakp+VUfLp0NTAKBBAua0ICoFuCBEAEQQIggiABEEGQAIggSABEECQAIggSABEE\nCYAIggRABEECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZB2KOXuWwDQLkEC\nIIIgARBBkKiMDafQKkECIIIgARBBkACIIEgARBAkACIIEgARBAmACIIEQARBAiCCIAEQQZAAiCBI\nAEQQJAAiCBIAEQQJgAiCBEAEQQIggiABEEGQqI93MYcmCRI1KeUxDHffCOAcgrSVZ+UApxIkACII\nEgARBAmACIIEQARBAiCCIAEQQZAAiCBIAEQQJAAiCBIAEQQJgAiCBEAEQQIggiDB1Zw5HmYJUi+8\nkxAQTpAAiCBI3RkGm4yARIIEN/CcAKYEqQujHUjbh6SS9MB5zW6ws5ecuTMv6gd9jQ6XnO/v7htw\nnfffvyHwIeFazyZ1fzcAQXqZkEbPhrp6cvQpPM8mvf8vWWP5rH1PXtU3nlhdTEjP/LxPRaWUUkrz\nc9LzUWNhlaOL5h5lhpCHntlNjtX9AEPuTMjU/oPy4/GYbc9CkEopj0cLd8v3P9vwbFf3+P5+X1Z3\n49/VPuF1LPovuv0JaWHr3OKjbQt/bT95yOhq8+bZmrkvm1kIUdoP0gHJzyAAWtXLQQ0AhBMkACII\nEgARBAmACIIEQIT2g7RwyJyj6QBytB+kp55PHQRQhehX7f7QtECdLBygFr1MSKP8qBFAml4mJADC\n9TIhARBOkACIIEgARBAkACIIEgARvB/SJu8vY2r4uMTpe71PL31q4E4YvTRtdkWNLfnR66qfFn69\nW1ryp1f9z75r9sKltzAhrXOWh0dbd0IpZXr7Vz9T9ZIfva76Zfmdozd+ZUsyV+11SCumz6qWx4h6\nvX4jPz2TauNOmF3m9JMtLfnp04oeTa/65dOvd3tL3nL7Y1dtQlrX/FkeZp87jzR2J2xZTktLnn24\naX7VL8u/3u0tecsSMlctSEu2j/lVG/6ZvbSTO+Fdh0t+tLvqjbuOtl+U7HWzy5tPX7Nw9bs4qIG+\nhDwTvNjCw3Hbd0gPa5ya3T9UxZ0gSFDTX+z3Ao+tOklXP9aR2Z2F+Wyyo3cV/bn+XM9rb9Ls5vfn\nf1bxszYh0a+FAwsbNjqmrpQ2D7XteTyqlyDRo342Wy0YhqGKZ83f+LRLv9sfejhBojsdDkYehamC\nfUhLFv6A+/nbbuxO2FKjxpa8UWOrHibeP//6eOHqV9zKXzvwmsKNF11DkNZlnmPjYo3dCVv+8Bpb\n8mPbitpb9armlzw7H2euus39mb81/VG1eqfteglhpXfC9hftN7Pkp43n3Gxs1e8+/Xo3tuTZH3Qt\nqzYhrcs8x8bFOrwTGlvy7Mk4Vj9T+6q3aGzJW37K00+GrNqEBEAEExIAEQQJgAiCBEAEQQIggiAB\nEEGQAIggSABEECQAIggSABEECYAIggRABEECIIIgARBBkACI8Hf3DYBaLbwb+vt7wS287eHqpd9/\nPVTEhARHrL7l80nNkCIaJkhw3OzbsK626nsX/BNwPUGC3zt1jjEk0Sr7kGC314Ay3aPz5eyycPVR\nh0opykRjTEiwmxLAGUxIcNxsmb7J1fIBe+9fZjcS7REk+JmFSBzrx+oh42Y1WmKTHYQyA9EbExJc\nYfmFsQeuCO0xIUEiZ2SgQ4IEcdSIPgkSZFEjuiVIAERwUAP8zPcvDxqdA2L0zVc/A1UzIcGPObkq\nHGNCgiO2TyfLXzk9BcPx2wSVMyHBL11WFOmiPYIEv3fqJjXb62iVIMGPXTC7GI9okiDBKU6aY5xQ\nlYb55QYgggkJgAj/A758paCvz/fGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = zeros(size(u));\n",
    "t = (1:length(y))/Fs;\n",
    "\n",
    "for k = 5:length(y)\n",
    "   y(k) = 0.1*y(k-1) - 0.5*u(k-1)*y(k-1) + 0.1*u(k-2); \n",
    "end\n",
    "Y = computeSignalFFT(y', Fs, fres); f_y = -Fs/2:fres:Fs/2;\n",
    "plot(f_y, 2*abs(Y));xlim([-1 50]); xlabel('f(Hz)');ylabel('|Y(f)|')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab_kernel"
  },
  "language_info": {
   "codemirror_mode": "Octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "octave"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}