{ "cells": [ { "cell_type": "markdown", "metadata": { "nbsphinx": "hidden" }, "source": [ "# Random Signals and LTI-Systems\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": [ "## Power Spectral Densitity\n", "\n", "For a wide-sense stationary (WSS) real-valued random process $x[k]$, the [power spectral density](../random_signals/power_spectral_densities.ipynb#Power-Spectral-Density) (PSD) $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ is given as the discrete-time Fourier transformation (DTFT) of the auto-correlation function (ACF) $\\varphi_{xx}[\\kappa]$\n", "\n", "\\begin{equation}\n", "\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sum_{\\kappa = -\\infty}^{\\infty} \\varphi_{xx}[\\kappa] \\; \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,\\kappa}\n", "\\end{equation}\n", "\n", "Under the assumption of a real-valued LTI system with impulse response $h[k] \\in \\mathbb{R}$, the PSD $\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the output signal of an LTI system $y[k] = \\mathcal{H} \\{ x[k] \\}$ is derived by taking the DTFT of the [ACF of the output signal](../random_signals_LTI_systems/correlation_functions.ipynb#Auto-Correlation-Function) $\\varphi_{yy}[\\kappa]$\n", "\n", "\\begin{align}\n", "\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) &= \\sum_{\\kappa = -\\infty}^{\\infty} \\underbrace{h[\\kappa] * h[-\\kappa]}_{\\varphi_{hh}[\\kappa]} * \\varphi_{xx}[\\kappa] \\; \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,\\kappa} \\\\ \n", "&= H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\cdot H(\\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega}) \\cdot \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = | H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2 \\cdot \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\n", "\\end{align}\n", "\n", "The PSD of the output signal $\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of an LTI system is given by the PSD of the input signal $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ multiplied with the squared magnitude $| H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2$ of the transfer function of the system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example - Pink Noise\n", "\n", "It can be concluded from above findings, that filtering can be applied to a white noise random signal $x[k]$ with $\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = N_0$ in order to create a random signal $y[k] = x[k] * h[k]$ with a desired PSD\n", "\n", "\\begin{equation}\n", "\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = N_0 \\cdot | H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2\n", "\\end{equation}\n", "\n", "where $N_0$ denotes the power per frequency of the white noise. Such a random signal is commonly termed as [*colored noise*](https://en.wikipedia.org/wiki/Colors_of_noise). Different application specific types of colored noise exist. One of these is [*pink noise*](https://en.wikipedia.org/wiki/Pink_noise) whose PSD is inversely proportional to the frequency. The approximation of a pink noise signal by filtering is illustrated by the following example. The PSDs $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ and $\\Phi_{yy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ are estimated from $x[k]$ and $y[k]$ using the [Welch technique](../spectral_estimation_random_signals/welch_method.ipynb)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDM5OS42MTM3NSAzMjUuOTgwNzUgXSAvUGFyZW50IDIgMCBSIC9SZXNvdXJjZXMgOCAwIFIKL1R5cGUgL1BhZ2UgPj4KZW5kb2JqCjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMSAwIFIgPj4Kc3RyZWFtCnic7ZpNjxa5Ecfvz6foI3ugcdnllzqCSJByybKLspdc0OwsCwI2wG4IUj58/lXuLreZgeEhWSKiDEIzU+Mu2+V6+VX7oeXZiZYnS1ie4f/bhZYH+P/kFPDbi1MSWQulmvHb88NvKeZVWsBPzzHy+NvPp9NPpzt3oeINnnlwOqW8RrFn6ppZf1C1cWWahM+PwpDWknOX+uMHoU3yarmqOuW6UlxiKSvh++vL5Yfl5XLnbtTlxOVP2Bq2ubbl7SmsRQJJCbUx9nzn/uXfn15cfvfg3nLx5vpFX7eS0/enh8urfYYAQ/oMun2V3KCM17AZ4HQPJ/H2dO/RcuePtFBYHv10inVNHErMS4yr5OXRj6dbYQ3fLI+eLX94ZHP/flbAgVAiKum48oPwTDvcpO7jlmi8Um0htVjrZIz8ZYxB3FYpLWP2w/KP0jPNcbPCjxuEmNZAmTPlzAeL0Jdyjxjy2kIWOa7/IDzTHjepu8EaUteWYyvcUpDJHF/IQWKJa2FE6+TfR+m5BrlR4Q25I7OuN9ZGoR5NEr+UhySSlVtoVKbcd5Cem0pvVPhxkyDTrDnWIpQazyb5Ql6SalkTtTid6UF4rkFuUHeDOUpbU8kUuBQ6miMND3mFx8JyO2DyGAQxpupkVa9CVbp4AY0nHYaFBR0Z/Ie19AGY7/b9y2eP//Lb949fvrn95xeXTx4v9385PdR/v6+xP5U9PruM01rlfV3d0Ae71TUuKa2VS6ubyY6FbU21idX4YMYvZnsz6umKUa8x6IunL39784UNKmGt8JnZoC4806BHZbtBj7quNWgrqJMSP8Wg+SswKHECQ7c8W3RIzzTppG636aTtWqMSh5WEyqdYlb8Gq4IKIgV5z6ouPdeqR3Vu1aO2660qCTNz+BSrpq/AqjGHlVNKs1WH9EyrTup2q07arrVq5LKWVPhTrBq/AquO0v/iA73weVXqqM7L1FHb9VYVrfHtkwoV/ZesCoYzq1KkXFrRGgBnKdtXxS81YAsx5SjLd+8b/BDOBbYTqonURm0NzTwvRbQcMbYKIyEftmpSdFslwgAqhBFLF4ZVYk7JhpaV8z6UQiZSIUiJx/NodHVoBoGnXSpRGieV8oo1b2qrIGHoZDBIf+dikwVmG4uVowPYNCSmVnRsScjePhsIsIlK9e3MrqHFLDAspG2tvrNYuZJqqLBt2aVwlMYm5TVtW0trCmBt3VvFUW17i/0sdG9VVtn2BviE25sdWlxr3DUwNhRVb8trpl1DwGNB19vgaGGX4hg0ap6fhHTEbjMR4L5K09raPjYFfOneBJmh7lI8L1X1wrO5jNlwzjqWAibxzcWSkhQT5w2dVUwIGa4mRpuZXAk2XW00BfOO/Ty4VjZx2pOjGrnkJs3EZY3BdVfKkkyM1l9847mmaFNGRG9zJSFRstGR11z3Y0U/GGoXVyB992PBqaDaa7DD4KvwPjoX4WabRyzU5CtpJScxcUObtItrqZFMt+6NhjdHsnAgHDQFd1EM7hZEQDQZpsJBmm6ERPHt4PDwTcWICa77lIKWv5mSnPVMdiWM/GCGzYjWvO+yceG+QA1oD7eYqJhqhEXxcKNUorkZIcewx1tDELHNWBCyNMQx2x4rZnenZI6IDRPzXq8gRt5M1WZEaGR3S67MfSENB+VHhi40dPs1mN7zScAy+rkjOiS7F4cauY9G4I6wayiDdsCIj+wnmdDyk20ScEKHEItFuhhReswfmXTKiFgoMnTDamzijP7aHbPABcXEDYljnA2OwpSgxw6eQwL8yGISNt2aR0stNXEsJi76OmN3tUC12JTa0g/HhJ2DjUaLGj3OUMl6lo0RweqOmZEZoy0QsdDC8MBYLeJjUgPt4gxfCjZlQs3x7VCCffpofSHnKSZyzrYSTFT9dEAhyY5Bk14eCT8wzGJiROUhbSCP27qzFjGfUlJPVBGxIH5oLKgbpgSxUH07CfWRbd2IBZbdVPAHJFoVF0RA85QbSusWLMjKIzfSbkEEQxuJvyFX2UoQDMV3mUPD2au4agx4sDYRSzIoxAg63zzWXWzKimAdu4Q/kK0b0dDIj7hgieY+iIYyIg0JqbsP+r0kw31ybM3ESHLDNQlZ2kYjGsSrW9CXgrYSQbj6LqNWdtsOgiTnPQBRoqttB+UEoeFBUpCoq4mxWndNTqFsYoSruyZSkljMJwRJddeUUsUqtb2IGmdZpWU2Me9QbM6GhSQT1x1O1O0KlBYTi9MJTqcGsRkRIzueJC0jSDQmzs4nUQtatY4FhdgJBaMxZbMZgYzBHbYh0VswwIedUXDwVLNRQ0Lo7JDCmgaz+RpK6YFSUBhQdlTcfX3PsFKbuQmqpnMKRoN+jDISQucIKiCgfi8UnFSQS1EWWxcnRxXSWbLFdsrlyCqghGTrRuiI7xK42BMyjtppBaUNSLyJ2XElIYgszFKpTiv6FPKuja3BcSWqc3dPSwic5ntE3jUSx+/OK5a+MplBEDfsB1lKZgMh5HYHFk3XSI3mf40dWHSLrQNsQtzswILijMQvdpCImwEspUSkuOV4B4eNwOygBP6weOtH/h2KXwgcz4IsLcpS2j1k7s1S+xi6A8/VRYBIR3RPoLMC5rSgGegOSCTkVEu5g92ruk2WTBO7ZyvobC+eDvCOY4u1WCYamK4EhgRkJX9gOiiO0wb6A9NR/ZCF2FLFwPQgCNsNnAemB9bkxeaiA9ODYi9SJx8xXQu5Zj1zGMd0koYIT7XD+47phOyXsAkLcMd0fTlRCRVUjphOopHJwbDKMZ0aWkRUEpvNMV1BBWwSjIgc00l1SRZbr2M6wXnhpWzZzjGdsJqccs8ajulAAeQSdM1lwnQVE746NjqmK4oRGMTOfmA6kjZMgePnCdMJdi1g+pnSFQkrkl3nQ6d0FeOUI8tE6YhqcATHDlpO6YBO7fTEct2gdEA+FsJb0+GUjsNHKUfVsk2mYOem4oxyx0QzvJNe6SA4rPcZ8E5GBq0Yfw1419ENGbpzo8M7zAGIwz+Z4F3BHKuJHVUd3vVCEtXTMGawOynQAKV7z+HsroNDgD/Hid11kgb+yDSxuy0JcdUJ29ldxSy1Wm4c8K7vJhEwYSP9Hd51u9g7W0Ya8E6ajABxnemd3lWMdpuEJnpXsRbyDsdO72phxHr/ZMOgd1R/pP+a++ad3ilpz8E5zPRup4cyHNJE7zoamZ67Czq9qxjG7p3woHc0ckBw6c7t8E4RpxRBdjO8q7ji4G3vA95VjMrHG9Pv8K4uiixI1poOeAenIWuWDiAD3pEpVrVfbhO8qxKsttj7ugHvNrqkfjgD3ikGhfeWZnjHGvQ8pDcGDu8qRh5pvTFweNdARN0PliYHvKNXgBlAzDTBu4qRxoVmeNcgByREnuEdcYdCg4amTfCuYhxQb60GvCMXIYWGZpQw4F3FjDRtLw4GvGu6wbqjNckD3ruSvOl2eLfR6MsN1ga8o6zCJtprT/CuCQ7FlNoM78iWyHA4yjLBuypBO8N9gQ7vhKpUYXTDrwHvm7hYfhjwrlOqdeLE7joYbhB7z+HsrnlZcHxGpIPdSZ0R2uyF12B3CsW6CHsnMdhdxZnL1ig6u1spwCyJJnZXMRJpozyxu66kCPUiO9hdR2Nof3E22F3Fuo7Or87uKkZznzryOburGJW+2roHuxPqOvy4V+XB7joa6bRnkwHvKs4VZSdP8K5K4JAh0gTvpgS9ZogTvOtoDN4+SuXwrqOTvpOOE7yrYYk2nHR2R3VAvpaQJ3RHECBWpcPEQHct06UiUcUJ3XV0y1J7N+PormLkztBndHRXcUaApDyhuyGAAI95QncVM9CA0oTuunU0RM3a24HuOlqgvNYJ3m1K/M0oaMC7rRuG6ufu8I4GYrVkFyd6t+2gBIWZ3tWAiMjWF+j0ritBlYiGj4PebSUZCCoTveuUap1t9E7vOhpnwL2pdEzHBnRGI7LrpfNdgjL5uEf4/M/CXXO5/pHLjHOu58+7Jgm6n+1+hGxXT/wuI/YPTyGMata7DJTa7VG90/j2l7eXr5fv/3Z58evrx8+X+5cv3zz99d3+YYg7d5Oayz8wCZP1j0siT4rfg7WxwJavSPupVoQlYmeIE9h6Hypo/aP2CgdhCdvAi9NBqgjtWodUIcdi8ziVjKl8VS67OO5Axlot/7egixliRM/VqY7CctB6dQfPj9Kx2eNUeUx11a4X+gHQe6dr+8nPuRU66ae3Gu3zNG2etvMLV6XvOdi4ejvcpKlP7c/xvmr7qIxekDXzO3jbuFu79U+7SkNlSSmbn99WDgw+MC7Vxv3Dxlmrk/tHcPpd3C5HJxjI9M/P+0R/vWUj9R278jvpTV/Y59jGXNoQLQuFQWn4WrR02HTTap7ZQAQqEuD299vzyDHvNx+Yd9t60opEfLheXDY51Qa6MqP0Pzy1P9jNTtL7EH/i5bbu9N5F5Qc1/bg/wLXkMf6e32u+96mmaZHbOV695/z256f9ltMeR3HbnsdPSV8zYfMBQVMWuEZveW74lJS5eH9dspzxuuTg05r6r/q0S8/06VL+oz797gM+/e7/Pv2/69MPT/8CwuPvagplbmRzdHJlYW0KZW5kb2JqCjExIDAgb2JqCjM0MjMKZW5kb2JqCjE2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzI3ID4+CnN0cmVhbQp4nD2SS47EMAhE9zkFF2jJ/Pw5T0az6rn/th92ehZRoQBFUbhrSpPo8lKXHCrZm/zoFVPFzeTviYbENHE30MEp97WjCP4EuA5m7gzRSolBvXWJyvKRsfPX19OXB+/L22GwUXqWmFefkrEGC/J0dPEB5tza7n+V7yuTjiy9Sk9F70u9WJaoDzH7YvUp3MbmG9t6MJ5M8Ws2UauZU9baIroM1EVKN7HVC1DtUYHjXEdvFW0DEIE7QQeVCbVqbMFISLRF1i/MaUL3xP/tAbpf2sR17IPYtI1nfEUaQ9DduEdZw8wx6nxqc8Pa09km65TiTC7XsjhWCTRKy1hqS/Csel0UY4MZaKUQzdFqbGVodujrZXg/yAxvO4qxqKVKz543p+BcfezXFHrwPIsdwZuc2C3o55CsQgYtNmHqLF2v8Dnnff1+ACo9fUYKZW5kc3RyZWFtCmVuZG9iagoxNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE1MiA+PgpzdHJlYW0KeJw9T8sRQyEIvFvFNsCMIILW8zI5mf6vAU08sbK4H/GBCmIWkFSHuMOt4sWFtaOr41OkHbQKi4PmBpa/ErGCvIHr1ukYPWabeIrnxhg6Y4awcMyRDFnaxTPuti9Fjg2Cu2FWoekplj7kemEbB1J6s7RdoaZwLqYiY4Zx22mfWyXSjtMvSv2ariK9w9SStYvGn31/Abz3NmcKZW5kc3RyZWFtCmVuZG9iagoxOCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDkyID4+CnN0cmVhbQp4nD2MsQ3AMAgEe6b4BSJhjG3YJ0rl7N/mLSdp4PQP19KgOKxxdlU0HziLfHhL9YSNxJSmlUdTnN3aFg4rgxS72BYWXmERpPJqmPF5U9XAklKU5c36f3c9x6sbugplbmRzdHJlYW0KZW5kb2JqCjE5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTM5ID4+CnN0cmVhbQp4nD2PsQ3FMAhEe6a4BZAAGxvPk+hX/vu3wXGSAvF0oDvwYRCw1SzpaFLgteNUshpgF/zJpIHVBNotoRVoXUooDlo66whE2xb16Qd9rpN45FKxpGovtb4pYrk79I7RbVrAo2dO9q7Q5uByT0ZBJk7KU9ahkzR9NpkH1CLJZTza+9D8/pn0uwC7vC3bCmVuZHN0cmVhbQplbmRvYmoKMTQgMCBvYmoKPDwgL0Jhc2VGb250IC9EZWphVnVTYW5zLU9ibGlxdWUgL0NoYXJQcm9jcyAxNSAwIFIKL0VuY29kaW5nIDw8IC9EaWZmZXJlbmNlcyBbIDEwMSAvZSAxMDYgL2ogMTIwIC94IC95IF0gL1R5cGUgL0VuY29kaW5nID4+Ci9GaXJzdENoYXIgMCAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdIC9Gb250RGVzY3JpcHRvciAxMyAwIFIKL0ZvbnRNYXRyaXggWyAwLjAwMSAwIDAgMC4wMDEgMCAwIF0gL0xhc3RDaGFyIDI1NSAvTmFtZSAvRGVqYVZ1U2Fucy1PYmxpcXVlCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDEyIDAgUiA+PgplbmRvYmoKMTMgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDk2Ci9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL0l0YWxpY0FuZ2xlIDAgL01heFdpZHRoIDEzNTAgL1N0ZW1WIDAgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9YSGVpZ2h0IDAgPj4KZW5kb2JqCjEyIDAgb2JqClsgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAKNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCAzMTggNDAxIDQ2MCA4MzggNjM2Cjk1MCA3ODAgMjc1IDM5MCAzOTAgNTAwIDgzOCAzMTggMzYxIDMxOCAzMzcgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNgo2MzYgNjM2IDMzNyAzMzcgODM4IDgzOCA4MzggNTMxIDEwMDAgNjg0IDY4NiA2OTggNzcwIDYzMiA1NzUgNzc1IDc1MiAyOTUKMjk1IDY1NiA1NTcgODYzIDc0OCA3ODcgNjAzIDc4NyA2OTUgNjM1IDYxMSA3MzIgNjg0IDk4OSA2ODUgNjExIDY4NSAzOTAgMzM3CjM5MCA4MzggNTAwIDUwMCA2MTMgNjM1IDU1MCA2MzUgNjE1IDM1MiA2MzUgNjM0IDI3OCAyNzggNTc5IDI3OCA5NzQgNjM0IDYxMgo2MzUgNjM1IDQxMSA1MjEgMzkyIDYzNCA1OTIgODE4IDU5MiA1OTIgNTI1IDYzNiAzMzcgNjM2IDgzOCA2MDAgNjM2IDYwMCAzMTgKMzUyIDUxOCAxMDAwIDUwMCA1MDAgNTAwIDEzNTAgNjM1IDQwMCAxMDcwIDYwMCA2ODUgNjAwIDYwMCAzMTggMzE4IDUxOCA1MTgKNTkwIDUwMCAxMDAwIDUwMCAxMDAwIDUyMSA0MDAgMTAyOCA2MDAgNTI1IDYxMSAzMTggNDAxIDYzNiA2MzYgNjM2IDYzNiAzMzcKNTAwIDUwMCAxMDAwIDQ3MSA2MTcgODM4IDM2MSAxMDAwIDUwMCA1MDAgODM4IDQwMSA0MDEgNTAwIDYzNiA2MzYgMzE4IDUwMAo0MDEgNDcxIDYxNyA5NjkgOTY5IDk2OSA1MzEgNjg0IDY4NCA2ODQgNjg0IDY4NCA2ODQgOTc0IDY5OCA2MzIgNjMyIDYzMiA2MzIKMjk1IDI5NSAyOTUgMjk1IDc3NSA3NDggNzg3IDc4NyA3ODcgNzg3IDc4NyA4MzggNzg3IDczMiA3MzIgNzMyIDczMiA2MTEgNjA4CjYzMCA2MTMgNjEzIDYxMyA2MTMgNjEzIDYxMyA5OTUgNTUwIDYxNSA2MTUgNjE1IDYxNSAyNzggMjc4IDI3OCAyNzggNjEyIDYzNAo2MTIgNjEyIDYxMiA2MTIgNjEyIDgzOCA2MTIgNjM0IDYzNCA2MzQgNjM0IDU5MiA2MzUgNTkyIF0KZW5kb2JqCjE1IDAgb2JqCjw8IC9lIDE2IDAgUiAvaiAxNyAwIFIgL3ggMTggMCBSIC95IDE5IDAgUiA+PgplbmRvYmoKMjQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNTkgPj4Kc3RyZWFtCnicPVJJcsMwDLv7FXwCd0nvSacn5//XAvTUlxBjigAIpneLytn4aStZfuTHLjtLIrd8B22T+4qqB6RugswlxyVXiFnK50qQWLSUmVifqQ7KzzWoVfjCT8xMTIckEUvIST2KsH5eB/egfr2k81tk/KNjg9JbkkwFnRrulaOU2LBUnxrkHjTdlFafmXZlByyNN2SlJnX69dPjB0swpS8S0UGIxJ/kcocsCykH8Xau3kB4V7sg1VMrDztEmo+R3lIFqzkzAUtG0w4jYG3WmCBVArxFawi0FuLyhU4rIj45N3QwTVWehciwUuFYe89oDRe6hrzScHqScf/5H7wr3tfvH9rzZmQKZW5kc3RyZWFtCmVuZG9iagoyNSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MCA+PgpzdHJlYW0KeJw9kEsSwyAMQ/ecQkfA+H+edLpK7r+tDZ1ssBiE9MB9YiKjFieCr8SHBqXDJPBsFYR7MNkRcoTkBE2GsoMkcQ0NBqXCpmOZ78mmddJKrLzRftl3NGaddIotRYd2If/n9SLco+Aa6xk8D2AxyNpKpeyZMFplpq7yqOi1H9PhPQ9Eq8Xl9Qau8NpHN6koKkvq/kR3NNj+kbf7Ht8fmWU4JAplbmRzdHJlYW0KZW5kb2JqCjI2IDAgb2JqCjw8IC9CQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIzNwovU3VidHlwZSAvRm9ybSAvVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJw9UbtxxTAM6z0FRuBH/Gied5cq2b8NKDkpeIApEQTkpyzRhZ9niOD7We7/yAOSrVBthCc0FZEN08DnSRFYbqQm3F7c54RslSP24lwgwhDtPAlppAsWOxkL3hc/j6seZqy5Yfy+M5p9VHTVUR28ew7jZk0/TpTd682sjlub+3TvrhOHa0gmn/cfnJRKp5csgzpLuLA2mhrW47woxljMOP4nqrBNsrajCsHSJUgq0IAYShLGgMUt/iInWg4L2psbaeudyU6qNIqGF6MM3qD1RjiKdJF8mGsrg7GpmDa++eQlN+j7Z7+fr18Da1rrCmVuZHN0cmVhbQplbmRvYmoKMjcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjUgPj4Kc3RyZWFtCnicRY87EgMhDEN7TqEjgH/AeTaTir1/G8s7SRosjCU/ois69srDY2PKxmu0sSfCFu5SOg2nqYyviqdnXaDLYTJTb1zNXGCqsMhuTrH6GHyh8uzmhK9VnhjCl0wJDTCVO7mH9fpRnJZ8JLsLguqUjcrCMEfS90BMTZunhYH8jy95akFQmeaNa5aVR2sVUzRnmCpbC4L1gaA6pfoD0/9Mp70/3PQ9gAplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9CQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIzOAovU3VidHlwZSAvRm9ybSAvVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJxFkE2ywzAIg/c+BUcwmB/7PJ15q/T+20okfV1kYAxI+jKWh0QeeQ877EwsUPFq5hJa8hq6VbyOaE5x211XOifs5hRumJ++MVVMqKJ7CVU1Q+ijcTD5Ol7DffXbe/ix3oplfRX51cmZrZxW7dS1vdkxDTeYjjd3XqqQgKokog8dX/+OV6eoIjZsN9xIMcHo4IuF3Q3jyahi1eW2ReNhgmmY9XpgjF9xEFLxC6CX7g2coQ9w2Rf413F6jTB0DLiBovjWeRQLkehRS9uS9c7AziqFG8zImzs0VYhBVVC1YT3g88GeP/xr/H0AqH5ggAplbmRzdHJlYW0KZW5kb2JqCjI5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzM4ID4+CnN0cmVhbQp4nDVSS5JbQQjbv1PoAq5q/s15nJrV5P7bCOysoIEWEpAWOMjESwxRjXLFH3mC8TqBv+vlafw+3oXUgqci/cC1aRvvx5o1UbA0YinMPvb9KCHHU+PfEOi5SBNmZDJyIBmI+7U+f9abTDn8BqRpc/ooSXoQLdjdGnZ8WZBB0pMaluzkh3UtsLoITZgbayIZObUyNc/HnuEynhgjQdUsIEmfuE8VjEgzHjtnLXmQ4XiqFy9+vY3XMo+pl1UFMrYJ5mA7mQmnKCIQv6AkuYm7aOoojmbGmtuFhpIi9909nJz0ur+cRAVeCeEs1hKOGXrKMic7DUqgauUEmGG99oVxmjZKuFPT7V2xr99nJmHc5rCzUjINznFwL5vMESR73TFhEx6HmPfuEYzEvPldbBFcucy5JtOP/SjaSB8U1+dcTZmtKOEfquSJFdf4//zez88/kDd9sQplbmRzdHJlYW0KZW5kb2JqCjMwIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzA0ID4+CnN0cmVhbQp4nD2SO5LDMAxDe52CF8iM+JPk82Qnlff+7T4yyVaASYkAKC91mbKmPCBpJgn/0eHhYjvld9iezczAtUQvE8spz6ErxNxF+bKZjbqyOsWqwzCdW/SonIuGTZOa5ypLGbcLnsO1ieeWfcQPNzSoB3WNS8IN3dVoWQrNcHX/O71H2Xc1PBebVOrUF48XURXm+SFPoofpSuJ8PCghXHswRhYS5FPRQI6zXK3yXkL2DrcassJBaknnsyc82HV6Ty5uF80QD2S5VPhOUezt0DO+7EoJPRK24VjufTuasekamzjsfu9G1sqMrmghfshXJ+slYNxTJkUSZE62WG6L1Z7uoSimc4ZzGSDq2YqGUuZiV6t/DDtvLC/ZLMiUzAsyRqdNnjh4yH6NmvR5led4/QFs83M7CmVuZHN0cmVhbQplbmRvYmoKMzEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA1NCA+PgpzdHJlYW0KeJwzNjZXMFAwNDJX0DUyNlUwMjRQMDczUUgx5IIxc8EssGwOF1whhAmSz4GrzOFKAwBMkA8VCmVuZHN0cmVhbQplbmRvYmoKMzIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzAgPj4Kc3RyZWFtCnicNVFJbsMwDLzrFfOBAOIuv8dBT+3/rx3SCWBgaEuczREbGxF4icHPQeTGW9aMmvibyV3xuzwVHgm3gidRBF6Ge9kJLm8Yl/04zHzwXlo5kxpPMiAX2fTwRMhgl0DowOwa1GGbaSf6hoTPjkg1G1lOX0vQS6sQKE/ZfqcLSrSt6s/tsy607WtPONntqSeVTyCeW7ICl41XTBZjGfRE5S7F9EGqs4WehPKifA6y+aghEl2inIEnBgejQDuw57afiVeFoHV1n7aNoRopHU//NjQ1SSLkEyWc2dK4W/j+nnv9/AOmVFOfCmVuZHN0cmVhbQplbmRvYmoKMzMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMjcgPj4Kc3RyZWFtCnicNU87sgMhDOs5hS6QGYxtYM+zmVQv92+fZLINEv5I8vRERyZe5sgIrNnxthYZiBn4FlPxrz3tw4TqPbiHCOXiQphhJJw167ibp+PFv13lM9bBuw2+YpYXBLYwk/WVxZnLdsFYGidxTrIbY9dEbGNd6+kU1hFMKAMhne0wJcgcFSl9sqOMOTpO5InnYqrFLr/vYX3BpjGiwhxXBU/QZFCWPe8moB0X9N/Vjd9JNIteAjKRYGGdJObOWU741WtHx1GLIjEnpBnkMhHSnK5iCqEJxTo7CioVBZfqc8rdPv9oXVtNCmVuZHN0cmVhbQplbmRvYmoKMzQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDUgPj4Kc3RyZWFtCnicRVC7jUMxDOs9BRcIYP0se553SJXbvz1KRnCFIVo/kloSmIjASwyxlG/iR0ZBPQu/F4XiM8TPF4VBzoSkQJz1GRCZeIbaRm7odnDOvMMzjDkCF8VacKbTmfZc2OScBycQzm2U8YxCuklUFXFUn3FM8aqyz43XgaW1bLPTkewhjYRLSSUml35TKv+0KVsq6NpFE7BI5IGTTTThLD9DkmLMoJRR9zC1jvRxspFHddDJ2Zw5LZnZ7qftTHwPWCaZUeUpnecyPiep81xOfe6zHdHkoqVV+5z93pGW8iK126HV6VclUZmN1aeQuDz/jJ/x/gOOoFk+CmVuZHN0cmVhbQplbmRvYmoKMzUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDcgPj4Kc3RyZWFtCnicTVG7bUQxDOvfFFzgAOtreZ4LUl32b0PJCJDCIKEvKaclFvbGSwzhB1sPvuSRVUN/Hj8x7DMsPcnk1D/muclUFL4VqpuYUBdi4f1oBLwWdC8iK8oH349lDHPO9+CjEJdgJjRgrG9JJhfVvDNkwomhjsNBm1QYd00ULK4VzTPI7VY3sjqzIGx4JRPixgBEBNkXkM1go4yxlZDFch6oCpIFWmDX6RtRi4IrlNYJdKLWxLrM4Kvn9nY3Qy/y4Ki6eH0M60uwwuileyx8rkIfzPRMO3dJI73wphMRZg8FUpmdkZU6PWJ9t0D/n2Ur+PvJz/P9CxUoXCoKZW5kc3RyZWFtCmVuZG9iagozNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDkwID4+CnN0cmVhbQp4nE2NQRLAIAgD77wiT1BE0P90etL/X6vUDr3ATgKJFkWC9DVqSzDuuDIVa1ApmJSXwFUwXAva7qLK/jJJTJ2G03u3A4Oy8XGD0kn79nF6AKv9egbdD9IcIlgKZW5kc3RyZWFtCmVuZG9iagozNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiFtCNEGUglgQpWYmZhBJOAMilwYAybQV5QplbmRzdHJlYW0KZW5kb2JqCjM4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNDUgPj4Kc3RyZWFtCnicMzK3UDBQsDQBEoYWJgrmZgYKKYZclhBWLhdMLAfMAtGWcAoingYAn30MtQplbmRzdHJlYW0KZW5kb2JqCjM5IDAgb2JqCjw8IC9CQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDM3Ci9TdWJ0eXBlIC9Gb3JtIC9UeXBlIC9YT2JqZWN0ID4+CnN0cmVhbQp4nOMyNDBTMDY1VcjlMjc2ArNywCwjcyMgCySLYEFk0wABXwoKCmVuZHN0cmVhbQplbmRvYmoKNDAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjEgPj4Kc3RyZWFtCnicRZBLEsMgDEP3nEJH8EcGfJ50ukrvv60hTbOAp7FABncnBKm1BRPRBS9tS7oLPlsJzsZ46DZuNRLkBHWAVqTjaJRSfbnFaZV08Wg2cysLrRMdZg56lKMZoBA6Fd7touRypu7O+Udw9V/1R7HunM3EwGTlDoRm9SnufJsdUV3dZH/SY27Wa38V9qqwtKyl5YTbzl0zoATuqRzt/QWpczqECmVuZHN0cmVhbQplbmRvYmoKNDEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMTQgPj4Kc3RyZWFtCnicPVC7EUMxCOs9BQvkznztN8/Lpcv+bSScpEI2QhKUmkzJlIc6ypKsKU8dPktih7yH5W5kNiUqRS+TsCX30ArxfYnmFPfd1ZazQzSXaDl+CzMqqhsd00s2mnAqE7qg3MMz+g1tdANWhx6xWyDQpGDXtiByxw8YDMGZE4siDEpNBv+tcvdS3O89HG+iiJR08K755fTLzy28Tj2ORLq9+YprcaY6CkRwRmryinRhxbLIQ6TVBDU9A2u1AK7eevk3aEd0GYDsE4njNKUcQ//WuMfrA4eKUvQKZW5kc3RyZWFtCmVuZG9iago0MiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDgwID4+CnN0cmVhbQp4nEWMuw3AMAhEe6ZgBH4mZp8olbN/GyBK3HBPunu4OhIyU95hhocEngwshlPxBpmjYDW4RlKNneyjsG5fdYHmelOr9fcHKk92dnE9zcsZ9AplbmRzdHJlYW0KZW5kb2JqCjQzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjM2ID4+CnN0cmVhbQp4nE1QS25EIQzbc4pc4EkkIQHOQ9VV5/7bscNU7SqGGH9ID+myVR7rU2J1iezypU2XyjJ5FajlT9v/UQwCbv/QyEG0t4ydYuYS1sXCJDzlNCMbJ9csH487TxtmhcbEjeOdLhlgnxYBNVuVzYE5bTo3QLqQGreqs95kUAwi6kLNB5MunKfRl4g5nqhgSncmtZAbXD7VoQNxWr0KuWOLk2/EHFmhwGHQTHHWXwHWqMmyWcggSYYhzn2je5QKjajKeSsVwg+ToRH1htWgBpW5haKp5ZL8HdoCMAW2jHXpDEqBqgDB3yqnfb8BJI1dUwplbmRzdHJlYW0KZW5kb2JqCjQ0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTQ3ID4+CnN0cmVhbQp4nD1PuQ0DMQzrPQUXOMB6LFvzXJDqsn8bykZSCCJA8ZFlR8cKXGICk445Ei9pP/hpGoFYBjVH9ISKYVjgbpICD4MsSleeLV4MkdpCXUj41hDerUxkojyvETtwJxejBz5UG1keekA7RBVZrknDWNVWXWqdsAIcss7CdT3MqgTl0SdrKR9QVEK9dP+fe9r7CwBvL+sKZW5kc3RyZWFtCmVuZG9iago0NSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE0OSA+PgpzdHJlYW0KeJw1j0sOAyEMQ/c5hS8wUn6EcB6qrqb33zZhWgkJC9svwRaDkYxLTGDsmGPhJVRPrT4kI4+6STkQqVA3BE9oTAwzbNIl8Mp03zKeW7ycVuqCTkjk6aw2GqKMZl7D0VPOCpv+y9wkamVGmQMy61S3E7KyYAXmBbU89zPuqFzohIedyrDoTjGi3GZGGn7/2/T+AnsyMGMKZW5kc3RyZWFtCmVuZG9iago0NiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDQ5ID4+CnN0cmVhbQp4nDM2tFAwUDA0MAeSRoZAlpGJQoohF0gAxMzlggnmgFkGQBqiOAeuJocrDQDG6A0mCmVuZHN0cmVhbQplbmRvYmoKNDcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNTcgPj4Kc3RyZWFtCnicRZC5EUMxCERzVUEJErAI6rHH0Xf/qRf5SrRvAC2HryVTqh8nIqbc12j0MHkOn00lVizYJraTGnIbFkFKMZh4TjGro7ehmYfU67ioqrh1ZpXTacvKxX/zaFczkz3CNeon8E3o+J88tKnoW6CvC5R9QLU4nUlQMX2vYoGjnHZ/IpwY4D4ZR5kpI3Fibgrs9xkAZr5XuMbjBd0BN3kKZW5kc3RyZWFtCmVuZG9iago0OCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzMiA+PgpzdHJlYW0KeJwtUjmOJDEMy/0KfmAA6/Lxnh5M1Pv/dElVBQWqbMs85HLDRCV+LJDbUWvi10ZmoMLwr6vMhe9I28g6iGvIRVzJlsJnRCzkMcQ8xILv2/gZHvmszMmzB8Yv2fcZVuypCctCxosztMMqjsMqyLFg6yKqe3hTpMOpJNjji/8+xXMXgha+I2jAL/nnqyN4vqRF2j1m27RbD5ZpR5UUloPtac7L5EvrLFfH4/kg2d4VO0JqV4CiMHfGeS6OMm1lRGthZ4OkxsX25tiPpQRd6MZlpDgC+ZkqwgNKmsxsoiD+yOkhpzIQpq7pSie3URV36slcs7m8nUkyW/dFis0UzuvCmfV3mDKrzTt5lhOlTkX4GXu2BA2d4+rZa5mFRrc5wSslfDZ2enLyvZpZD8mpSEgV07oKTqPIFEvYlviaiprS1Mvw35f3GX//ATPifAEKZW5kc3RyZWFtCmVuZG9iago0OSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMxNyA+PgpzdHJlYW0KeJw1UktyQzEI279TcIHOmL99nnSyau6/rYQnK7AtQEIuL1nSS37UJdulw+RXH/clsUI+j+2azFLF9xazFM8tr0fPEbctCgRREz34MicVItTP1Og6eGGXPgOvEE4pFngHkwAGr+FfeJROg8A7GzLeEZORGhAkwZpLi01IlD1J/Cvl9aSVNHR+Jitz+XtyqRRqo8kIFSBYudgHpCspHiQTPYlIsnK9N1aI3pBXksdnJSYZEN0msU20wOPclbSEmZhCBeZYgNV0s7r6HExY47CE8SphFtWDTZ41qYRmtI5jZMN498JMiYWGwxJQm32VCaqXj9PcCSOmR0127cKyWzbvIUSj+TMslMHHKCQBh05jJArSsIARgTm9sIq95gs5FsCIZZ2aLAxtaCW7eo6FwNCcs6Vhxtee1/P+B0Vbe6MKZW5kc3RyZWFtCmVuZG9iago1MCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE3ID4+CnN0cmVhbQp4nDM2tFAwgMMUQy4AGpQC7AplbmRzdHJlYW0KZW5kb2JqCjUxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTMxID4+CnN0cmVhbQp4nEWPyw0EIQxD71ThEvIZPqmH1Z7Y/q/rMJpBQvhBIjvxMAis8/I20MXw0aLDN/421atjlSwfunpSVg/pkIe88hVQaTBRxIVZTB1DYc6YysiWMrcb4bZNg6xslVStg3Y8Bg+2p2WrCH6pbWHqLPEMwlVeuMcNP5BLrXe9Vb5/QlMwlwplbmRzdHJlYW0KZW5kb2JqCjUyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzM4ID4+CnN0cmVhbQp4nDVSOa7dQAzrfQpdIIB2zZznBal+7t+GlF8KQ7RWipqOFpVp+WUhVS2TLr/tSW2JG/L3yQqJE5JXJdqlDJFQ+TyFVL9ny7y+1pwRIEuVCpOTksclC/4Ml94uHOdjaz+PI3c9emBVjIQSAcsUE6NrWTq7w5qN/DymAT/iEXKuWLccYxVIDbpx2hXvQ/N5yBogZpiWigpdVokWfkHxoEetffdYVFgg0e0cSXCMjVCRgHaB2kgMObMWu6gv+lmUmAl07Ysi7qLAEknMnGJdOvoPPnQsqL8248uvjkr6SCtrTNp3o0lpzCKTrpdFbzdvfT24QPMuyn9ezSBBU9YoaXzQqp1jKJoZZYV3HJoMNMcch8wTPIczEpT0fSh+X0smuiiRPw4NoX9fHqOMnAZvAXPRn7aKAxfx2WGvHGCF0sWa5H1AKhN6YPr/1/h5/vwDHLaAVAplbmRzdHJlYW0KZW5kb2JqCjUzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ4ID4+CnN0cmVhbQp4nC1ROZIDQQjL5xV6QnPT77HLkff/6QrKAYOGQyA6LXFQxk8Qlive8shVtOHvmRjBd8Gh38p1GxY5EBVI0hhUTahdvB69B3YcZgLzpDUsgxnrAz9jCjd6cXhMxtntdRk1BHvXa09mUDIrF3HJxAVTddjImcNPpowL7VzPDci5EdZlGKSblcaMhCNNIVJIoeomqTNBkASjq1GjjRzFfunLI51hVSNqDPtcS9vXcxPOGjQ7Fqs8OaVHV5zLycULKwf9vM3ARVQaqzwQEnC/20P9nOzkN97SubPF9Phec7K8MBVY8ea1G5BNtfg3L+L4PePr+fwDqKVbFgplbmRzdHJlYW0KZW5kb2JqCjU0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggODggPj4Kc3RyZWFtCnicNYy7EcAwCEN7T8EIBouP98mlSvZvg+3QgKR394KDOkHyuBspnC5u2Vd6G4+TniYAsfRMQ+3fYEXVi1oULV9uY9BiKr4/+iQglnXyXjj0kBLeH8UXHXsKZW5kc3RyZWFtCmVuZG9iago1NSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzOCA+PgpzdHJlYW0KeJw9j0EOAzEIA+95hT8QKXZCWN6zVU/b/19Lmt1e0AiMMRZCQ2+oag6bgg3Hi6VLqNbwKYqJSg7ImWAOpaTSHWeRemI4GNwetBvO4rHp+hG7klZ90OZGuiVogkfsU2nclnETxAM1Beop6lyjvBC5n6lX2DSS3bSykms4pt+956nr/9NV3l9f3y6MCmVuZHN0cmVhbQplbmRvYmoKNTYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMTAgPj4Kc3RyZWFtCnicNVDLDUMxCLtnChaoFAKBZJ5WvXX/a23QO2ER/0JYyJQIeanJzinpSz46TA+2Lr+xIgutdSXsypognivvoZmysdHY4mBwGiZegBY3YOhpjRo1dOGCpi6VQoHFJfCZfHV76L5PGXhqGXJ2BBFDyWAJaroWTVi0PJ+QTgHi/37D7i3koZLzyp4b+Ruc7fA7s27hJ2p2ItFyFTLUszTHGAgTRR48eUWmcOKz1nfVNBLUZgtOlgGuTj+MDgBgIl5ZgOyuRDlL0o6ln2+8x/cPQABTtAplbmRzdHJlYW0KZW5kb2JqCjIyIDAgb2JqCjw8IC9CYXNlRm9udCAvRGVqYVZ1U2FucyAvQ2hhclByb2NzIDIzIDAgUgovRW5jb2RpbmcgPDwKL0RpZmZlcmVuY2VzIFsgMzIgL3NwYWNlIDQwIC9wYXJlbmxlZnQgL3BhcmVucmlnaHQgNDYgL3BlcmlvZCA0OCAvemVybyAvb25lIC90d28gL3RocmVlCi9mb3VyIC9maXZlIC9zaXggNjYgL0IgNjggL0QgODAgL1AgODMgL1MgOTcgL2EgOTkgL2MgL2QgL2UgMTA1IC9pIDEwOCAvbAoxMTAgL24gL28gL3AgMTE0IC9yIC9zIC90IDExOSAvdyAxMjEgL3kgMTI0IC9iYXIgXQovVHlwZSAvRW5jb2RpbmcgPj4KL0ZpcnN0Q2hhciAwIC9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnREZXNjcmlwdG9yIDIxIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDIwIDAgUiA+PgplbmRvYmoKMjEgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zIC9JdGFsaWNBbmdsZSAwCi9NYXhXaWR0aCAxMzQyIC9TdGVtViAwIC9UeXBlIC9Gb250RGVzY3JpcHRvciAvWEhlaWdodCAwID4+CmVuZG9iagoyMCAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzQyIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjMgNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjEyIDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTIgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwNQo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTgyIDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoyMyAwIG9iago8PCAvQiAyNCAwIFIgL0QgMjUgMCBSIC9QIDI3IDAgUiAvUyAyOSAwIFIgL2EgMzAgMCBSIC9iYXIgMzEgMCBSIC9jIDMyIDAgUgovZCAzMyAwIFIgL2UgMzQgMCBSIC9maXZlIDM1IDAgUiAvZm91ciAzNiAwIFIgL2kgMzcgMCBSIC9sIDM4IDAgUiAvbiA0MCAwIFIKL28gNDEgMCBSIC9vbmUgNDIgMCBSIC9wIDQzIDAgUiAvcGFyZW5sZWZ0IDQ0IDAgUiAvcGFyZW5yaWdodCA0NSAwIFIKL3BlcmlvZCA0NiAwIFIgL3IgNDcgMCBSIC9zIDQ4IDAgUiAvc2l4IDQ5IDAgUiAvc3BhY2UgNTAgMCBSIC90IDUxIDAgUgovdGhyZWUgNTIgMCBSIC90d28gNTMgMCBSIC93IDU0IDAgUiAveSA1NSAwIFIgL3plcm8gNTYgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAyMiAwIFIgL0YyIDE0IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL0NBIDAgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTIgPDwgL0NBIDEgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTMgPDwgL0NBIDAuOCAvVHlwZSAvRXh0R1N0YXRlIC9jYSAwLjggPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgL0YxLURlamFWdVNhbnMtT21lZ2EgMjYgMCBSIC9GMS1EZWphVnVTYW5zLVBoaSAyOCAwIFIKL0YxLURlamFWdVNhbnMtbWludXMgMzkgMCBSID4+CmVuZG9iagoyIDAgb2JqCjw8IC9Db3VudCAxIC9LaWRzIFsgMTAgMCBSIF0gL1R5cGUgL1BhZ2VzID4+CmVuZG9iago1NyAwIG9iago8PCAvQ3JlYXRpb25EYXRlIChEOjIwMjAxMjE2MTUzNDMyKzAyJzAwJykKL0NyZWF0b3IgKE1hdHBsb3RsaWIgdjMuMy4yLCBodHRwczovL21hdHBsb3RsaWIub3JnKQovUHJvZHVjZXIgKE1hdHBsb3RsaWIgcGRmIGJhY2tlbmQgdjMuMy4yKSA+PgplbmRvYmoKeHJlZgowIDU4CjAwMDAwMDAwMDAgNjU1MzUgZiAKMDAwMDAwMDAxNiAwMDAwMCBuIAowMDAwMDE3NTUyIDAwMDAwIG4gCjAwMDAwMTcyMjIgMDAwMDAgbiAKMDAwMDAxNzI2NSAwMDAwMCBuIAowMDAwMDE3NDA3IDAwMDAwIG4gCjAwMDAwMTc0MjggMDAwMDAgbiAKMDAwMDAxNzQ0OSAwMDAwMCBuIAowMDAwMDAwMDY1IDAwMDAwIG4gCjAwMDAwMDAzOTcgMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDAzODk1IDAwMDAwIG4gCjAwMDAwMDU0NTggMDAwMDAgbiAKMDAwMDAwNTI1MCAwMDAwMCBuIAowMDAwMDA0OTE3IDAwMDAwIG4gCjAwMDAwMDY1MTEgMDAwMDAgbiAKMDAwMDAwMzkxNiAwMDAwMCBuIAowMDAwMDA0MzE2IDAwMDAwIG4gCjAwMDAwMDQ1NDEgMDAwMDAgbiAKMDAwMDAwNDcwNSAwMDAwMCBuIAowMDAwMDE1ODAwIDAwMDAwIG4gCjAwMDAwMTU2MDAgMDAwMDAgbiAKMDAwMDAxNTExMiAwMDAwMCBuIAowMDAwMDE2ODUzIDAwMDAwIG4gCjAwMDAwMDY1NzMgMDAwMDAgbiAKMDAwMDAwNjkwNSAwMDAwMCBuIAowMDAwMDA3MTM4IDAwMDAwIG4gCjAwMDAwMDc1MDkgMDAwMDAgbiAKMDAwMDAwNzc0NyAwMDAwMCBuIAowMDAwMDA4MTE5IDAwMDAwIG4gCjAwMDAwMDg1MzAgMDAwMDAgbiAKMDAwMDAwODkwNyAwMDAwMCBuIAowMDAwMDA5MDMzIDAwMDAwIG4gCjAwMDAwMDkzMzYgMDAwMDAgbiAKMDAwMDAwOTYzNiAwMDAwMCBuIAowMDAwMDA5OTU0IDAwMDAwIG4gCjAwMDAwMTAyNzQgMDAwMDAgbiAKMDAwMDAxMDQzNiAwMDAwMCBuIAowMDAwMDEwNTc2IDAwMDAwIG4gCjAwMDAwMTA2OTMgMDAwMDAgbiAKMDAwMDAxMDg2MyAwMDAwMCBuIAowMDAwMDExMDk3IDAwMDAwIG4gCjAwMDAwMTEzODQgMDAwMDAgbiAKMDAwMDAxMTUzNiAwMDAwMCBuIAowMDAwMDExODQ1IDAwMDAwIG4gCjAwMDAwMTIwNjUgMDAwMDAgbiAKMDAwMDAxMjI4NyAwMDAwMCBuIAowMDAwMDEyNDA4IDAwMDAwIG4gCjAwMDAwMTI2MzggMDAwMDAgbiAKMDAwMDAxMzA0MyAwMDAwMCBuIAowMDAwMDEzNDMzIDAwMDAwIG4gCjAwMDAwMTM1MjIgMDAwMDAgbiAKMDAwMDAxMzcyNiAwMDAwMCBuIAowMDAwMDE0MTM3IDAwMDAwIG4gCjAwMDAwMTQ0NTggMDAwMDAgbiAKMDAwMDAxNDYxOCAwMDAwMCBuIAowMDAwMDE0ODI5IDAwMDAwIG4gCjAwMDAwMTc2MTIgMDAwMDAgbiAKdHJhaWxlcgo8PCAvSW5mbyA1NyAwIFIgL1Jvb3QgMSAwIFIgL1NpemUgNTggPj4Kc3RhcnR4cmVmCjE3NzY5CiUlRU9GCg==\n", "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2020-12-16T15:34:32.782596\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.3.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.signal as sig\n", "\n", "fs = 44100\n", "N = 5 * fs\n", "\n", "# generate uniformly distributed white noise\n", "np.random.seed(1)\n", "x = np.random.uniform(size=N) - 0.5\n", "# filter white noise to yield pink noise\n", "# see http://www.firstpr.com.au/dsp/pink-noise/#Filtering\n", "a = np.poly([0.99572754, 0.94790649, 0.53567505]) # denominator coefficients\n", "b = np.poly([0.98443604, 0.83392334, 0.07568359]) # numerator coefficients\n", "y = 1 / 3 * sig.lfilter(b, a, x)\n", "# estimate PSDs using Welch's technique\n", "f, Pxx = sig.csd(x, x, nperseg=256)\n", "f, Pyy = sig.csd(y, y, nperseg=256)\n", "\n", "# PSDs\n", "Om = f * 2 * np.pi\n", "plt.plot(\n", " Om, 20 * np.log10(np.abs(0.5 * Pxx)), label=r\"$| \\Phi_{xx}(e^{j \\Omega}) |$ in dB\"\n", ")\n", "plt.plot(\n", " Om, 20 * np.log10(np.abs(0.5 * Pyy)), label=r\"$| \\Phi_{yy}(e^{j \\Omega}) |$ in dB\"\n", ")\n", "plt.title(\"Power Spectral Density\")\n", "plt.xlabel(r\"$\\Omega$\")\n", "plt.legend()\n", "plt.axis([0, np.pi, -60, -10])\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's listen to white and pink noise" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy.io import wavfile\n", "\n", "wavfile.write(\"uniform_white_noise.wav\", fs, np.int16(x * 32768))\n", "wavfile.write(\"uniform_pink_noise.wav\", fs, np.int16(y * 32768))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**White noise**\n", "\n", "[./uniform_white_noise.wav](./uniform_white_noise.wav)\n", "\n", "**Pink noise**\n", "\n", "[./uniform_pink_noise.wav](./uniform_white_noise.wav)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cross-Power Spectral Densities\n", "\n", "The cross-power spectral densities $\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ and $\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ between the in- and output of an LTI system are given by taking the DTFT of the [cross-correlation functions](../random_signals_LTI_systems/correlation_functions.ipynb#Cross-Correlation-Function) (CCF) $\\varphi_{yx}[\\kappa]$ and $\\varphi_{xy}[\\kappa]$. Hence,\n", "\n", "\\begin{equation}\n", "\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sum_{\\kappa = -\\infty}^{\\infty} h[\\kappa] * \\varphi_{xx}[\\kappa] \\; \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,\\kappa} = \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\cdot H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\n", "\\end{equation}\n", "\n", "and\n", "\n", "\\begin{equation}\n", "\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sum_{\\kappa = -\\infty}^{\\infty} h[-\\kappa] * \\varphi_{xx}[\\kappa] \\; \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,\\kappa} = \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\cdot H(\\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega})\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## System Identification by Spectral Division\n", "\n", "Using the result above for the cross-power spectral density $\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ between out- and input, and the relation of the [CCF of finite-length signals to the convolution](../random_signals/correlation_functions.ipynb#Definition) yields\n", "\n", "\\begin{equation}\n", "H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})}{\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})} = \\frac{\\frac{1}{K} Y(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\cdot X(\\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega})}{\\frac{1}{K} X(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\cdot X(\\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega})} \n", "= \\frac{Y(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})}{X(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})}\n", "\\end{equation}\n", "\n", "holding for $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\neq 0$ and $X(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\neq 0$. Hence, the transfer function $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of an unknown system can be derived by dividing the spectrum of the output signal $Y(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ through the spectrum of the input signal $X(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. This is equal to the [definition of the transfer function](https://en.wikipedia.org/wiki/Transfer_function). However, care has to be taken that the spectrum of the input signal does not contain zeros.\n", "\n", "Above relation can be realized by the discrete Fourier transformation (DFT) by taking into account that a multiplication of two spectra $X[\\mu] \\cdot Y[\\mu]$ results in the cyclic/periodic convolution $x[k] \\circledast y[k]$. Since we aim at a linear convolution, zero-padding of the in- and output signal has to be applied." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example\n", "\n", "We consider the estimation of the impulse response $h[k] = \\mathcal{F}_*^{-1} \\{ H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\}$ of an unknown system using the spectral division method. Normal distributed white noise with variance $\\sigma_n^2 = 1$ is used as wide-sense ergodic input signal $x[k]$. In order to show the effect of sensor noise, normally distributed white noise $n[k]$ with the variance $\\sigma_n^2 = 0.01$ is added to the output signal $y[k] = x[k] * h[k] + n[k]$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDQxMi40MjYyNSAzMjUuOTgwNzUgXSAvUGFyZW50IDIgMCBSIC9SZXNvdXJjZXMgOCAwIFIKL1R5cGUgL1BhZ2UgPj4KZW5kb2JqCjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMSAwIFIgPj4Kc3RyZWFtCniclVpNbx43Dr7Pr5hj9mBZpL6PDdIaWGAPTozdw2IPgeu6ycppayfbv78P550PaiK/8zZAUOepJFIkRT7kmMbPA42Pox0/4++fI403+Ps4WPzrafDExnPkgH9V9S/HwZRsUwBqm3/9Ogy/DNc/4IgX7LkZBp8N2WlPMsE7LMKxNhjmBqwatM7EEE7oul2BsxA+CXmEwlDeZKgPkYIM0RvncklJi1WgN3YWO7zF/f8c3t6N1z/RSHa8+2WIZCgUlqWcRmZTwnj38/DG/m28+zz+eDdMwgZy1oQSrG2kaPS8GGJv2IXoqITklBzaCyrR5EDWtYIUeiCoWOMLxUjRU1GCeCeIo4O5k41RC9LoeUEcokk+xUKOY1CC3E6Q42JcdLY0gjR6XpBjZyxuzmSJSAnye0ElmEiFuBWk0ANBuRgOJQUb4SQlKKyC/sAuO17ZUfRAkEYJV+xKmfDD/RMOlOdk/BTH2Hv9E8+nv/mvHDLgkNtZ3zXegzNMhV2U9+JN3IGv2d/w6AteKdbAZX6OW7OzyiqFLP4fQiizFqPQM3IQfkwOhsGfWQy/JsYjVC0FvCYtZkPPiCEvrvYxMMvrPwnyrwkqDo+PI56TFrSh5wTB1QlOLin7MsuJr8hhb43LLlPWchR6Rg5cY0KOxTuPnDcLyq8JyslQ8iX6RtCGnhOUnXEp21ys41NiMSpqJWCvJHSpGGLKMC5FyV2IdhW2aR+zv04xG5FLE04+FYdpxSL8zb+nFc4USy6m6WTbDfwAg+MFh94h/1mexh+jROMoFWn6IZroY0lIMBHBLppi19W7h88f//ntw8cvL1ffvnyyzvL47jc8q1vs/74GuQADwkIxGsJ/nx/Gf41fRjJBaohYI8QcWfQ2Ps5/EnyTbMgZOZvL+P5mvH738L9P9w/vb96O9y+qvOi3q1BUz5RSnB7VhyFZA/+F3KzeQC9FxOU0L4bisGZ27WqNIm1SdKflGYWFQsqxWa5Q/BRtcuW0HOUhe4KWzfIGJQRbmU8vSYxSCrfLNxRuSjadgvPDQBYptiCJtklNw4hW1Iw4X5YIhqZMpbWNhpE44H+Xlw3JOO+QS3YbNhheINxhkYA8jft469sNCsZTC86nRcJW1JsNCg4oXShHYdmQwGZs3Gm0oQgsVEu/2AhJLnvsbx2sYQQEqv0cPCT7xRithzXsDeqPW9cnGDtJ9W/XbzAbm2GxxUKo+Ynxilont7CDqELzBoSu83lfujTMBpRqMWiKBoKZdz5WcIBpiXOcN2TJ4SW7nUUbOPmIBzZvQPyiZCDLthsUDNqKG/jlBhufajesMBck0GBzmEUwOCnYTmwlbCjjxUUYNc7vDOGBRMTsWy8rmHM21hFCZNkBQgUm4eNuxwpLki+4A8+OEwriELm+9bSCWSwQSl6eM4N1lVRsoN2OFT5pxZnmaJUKFqUk7y6+wXA0sgEe3Oxu3M1I2ITWtgrmjPQYwMmWm4PG5BB2KWxDxRdBiNm8HCnBQ8G4M+0Gi4AYLdPsbsVqmx0bDFnRbRIi6AFobdwZVsEe5BAilhsgll3InNonpGE22bn1UYuGiGQpl82GBi5sg1s8BytHWyRLtRsaOBCtbhDOBsqWd65uYBjLzwEL15iMC+XW0S0MHuCSWzbA73jtOe02bDDyGHm4Yd5AUoFpV3oaFHykrMtBwCHNl9bNGmbj5D3kZcPaUrQbNpjQG+HH2WvwB6wNerbTSMHQDZ1HXlTyiCskatu6WcNIxcEFn5YNaA5cwcvabdhgD9YFhrN4AYkqo/e1rZtbOCF1hjmQUIik35POqNmgYNgrIHf7ZQMoIOov7fysYBaHpLRsSKBunh3t/KxgWMkGv6RK9BIweNpxBoUiokBNVoWQcUoBz9r5WcEgA86GhfHcjpfTvAsJ3vPj6yvHPRXUK4fXVz4/zu2ikO+VGrKMGVZ6yAj2+6fx+h9WGOy62OBmSF8lZvB/XAcrg00oxAdrrwhlF6c66WuOFoN+BpCCaJFjLzgYIYF87yVxHC1GE4d2AVkbih8tDgYZy1vn0YwdrZUqgKcvUe0OlWDknmCF1NKhxsJboAYJUzhcmixMHDNy7qEOKBIODQyhHT3Uwcscyl/mDqE1iZAskwsXGCIJHz7VrUMLgzEgB2QExgUxYTOjE/HSLRwslp6FpCwgKi44OHvwAwtyfqgFC8OzUXLAoe/wLhGZ6KI8gvPQFleo3RE/pQj7XXC/IEMqd8H1ZAqFnglEMx2+aBTrFCTbu+NMcSW2SNPTs4eLUQmRKmBnouMQQudgS0CrwYehWYRRCdOx+TiCsslFOl+E0aHzrtghr2BhRFopF+QVEAuhjyC0FxgO7RNYFJM/9IhEfYIeoGnHIWRBQ4tDcknHlsOjJnkeSETHlpMaBENEcnRoOeRuaVFQEtJx2kSTi1MtOHG6JAOgHwuWPV2wGPVGGDmc5w6NjF7KS7NaIh0bWYoeXp5w2UvuZwtKAp70cW2SDOfZphTYHweR1NOM8pT52NcSGGjekeD42H9Ci13MDC2OrQwGhSKCY20+tpzEfUL6LDlqy93+JYrF49/XaRrexIkaORAjAgvywgEktmfChLR+6TStywv/MgFU2kELGT2+vxla8a/I7wztanc4V7szuNodtdXuRK32B2e1Px6r/SFY7Y+6an+gVbtjq9qfTtX+EKr2Z021P1Gq/blR7c+Han8KVPuzntqMbvTEXE9oNKznMA2upi0aVzOVFleTE43r+YjG9RSkwddhRwPrmYbGewON2p9b1P50ovZnELU/aaj9eULtjw1qfzhQuyOA2u/0a7+fr/2uvfZ789rvwGu/z679brr2e+ba74xrt/+tB+lM0pRdP223+an7Wb33rRxndT63P73yuV1WX/rBvlm7nnHmZDvd5/S5nqbk+7h+teLpO7UnrM+nzI1gOe2VD1c/vnz99PTx68PP46en37/Vlwek8pfff/vy8rB8z7r+wcmvA6y/ywCbnX6TAfEVJAvIZzOUAUc8fQDOs44LBmeUIhPFgDy+gN5atQwNtDT0GszTsvtBYXgS6y8vKDQhPovF2VoMtDitVfqs2L3WfEXlvcgILRVK+gBy69JNlAKjOnVFZ/1rg60X1YJmg3xvy3v5bYy3y29jTJVV3LurrrvoxetYFcOSOLnE7bFd0Gy/3dF+oferSsiGiMYpYL4+f1tD4+Cz3neznJ2uvOqVZLZo5cenBkbNtJPzTkNW0VixnvEc69GmaI53HXgzyGV3eZ7NdtlsajfFOnvu0NzfL/cXs2T5PJtD0bhkbXD3cgID5cR+Akuaqh5yJmqIs1N4CjNxgaR6Ao+gNMJ3gS7iUP5Q5sA6NLj64HTCBoMtMHhBc2w2CV1SDI0OKCwoCn46QemLGoyiIL8poq6W7Xa11Q4beD+Zx4YQPbWLUS/B6Kg9FxWDUQb2OgSD8iKz+EZh37lbAyodfMdm+tjNvEoH7Qulr/LbdjXt4k48zMnBTkn5/Nvdolye78OS75c3fDv8H/NOu6QKZW5kc3RyZWFtCmVuZG9iagoxMSAwIG9iagoyNjMyCmVuZG9iagoxNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIxMyA+PgpzdHJlYW0KeJw9UDGSAzEI6/0KngAIsP2evbkq+X8bYWdS7IpBSEYUQlSQ/GWVzFL5s5E5BVB5j9gg9RqhXywyPk+1BeES06hPIgRknxEzBXOzU4K1Lu48TEk4NZyLgEQqv90M2ikklPPLqb/4jN6jK2+nSvGkDiDjRhdVsR3cYIotojZjnmdbXLCFk+w1fP4q48plXYE228SZP9mFuuR5AGQyGY+LJVVhF7lu+e3sLRmccmrdyGQyCMP2NSPXRTtM9Rk4zxBY1FQc52YYuPQ4Iuj0Pf8z/j+cDk48CmVuZHN0cmVhbQplbmRvYmoKMTggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4NyA+PgpzdHJlYW0KeJw1jcENwDAIA/+ZwiPEECDZp+qr3f9bSNSPfbKMsVjoUEtxCsI7LjZO3fg2iUNPZgFlJI1lsFFUJ4fEJ2RakrEWs8W+nREQrw7FdqLH/idPuz+4ThnECmVuZHN0cmVhbQplbmRvYmoKMTUgMCBvYmoKPDwgL0Jhc2VGb250IC9EZWphVnVTYW5zLU9ibGlxdWUgL0NoYXJQcm9jcyAxNiAwIFIKL0VuY29kaW5nIDw8IC9EaWZmZXJlbmNlcyBbIDEwNCAvaCAxMDcgL2sgXSAvVHlwZSAvRW5jb2RpbmcgPj4gL0ZpcnN0Q2hhciAwCi9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnREZXNjcmlwdG9yIDE0IDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL1N1YnR5cGUgL1R5cGUzIC9UeXBlIC9Gb250IC9XaWR0aHMgMTMgMCBSID4+CmVuZG9iagoxNCAwIG9iago8PCAvQXNjZW50IDkyOSAvQ2FwSGVpZ2h0IDAgL0Rlc2NlbnQgLTIzNiAvRmxhZ3MgOTYKL0ZvbnRCQm94IFsgLTEwMTYgLTM1MSAxNjYwIDEwNjggXSAvRm9udE5hbWUgL0RlamFWdVNhbnMtT2JsaXF1ZQovSXRhbGljQW5nbGUgMCAvTWF4V2lkdGggMTM1MCAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTMgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTYgMCBvYmoKPDwgL2ggMTcgMCBSIC9rIDE4IDAgUiA+PgplbmRvYmoKMjMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA3OSA+PgpzdHJlYW0KeJxNzbsNwCAMBNCeKTwC4P8+UaqwfxsbIkJjP+lOOsEOFdzisBhod7ha8aVRmH3qmRKSUHM9RFgzJTqEpF/6yzDDmNjItu+3Vu4X3hscGQplbmRzdHJlYW0KZW5kb2JqCjI0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzA0ID4+CnN0cmVhbQp4nD2SO5LDMAxDe52CF8iM+JPk82Qnlff+7T4yyVaASYkAKC91mbKmPCBpJgn/0eHhYjvld9iezczAtUQvE8spz6ErxNxF+bKZjbqyOsWqwzCdW/SonIuGTZOa5ypLGbcLnsO1ieeWfcQPNzSoB3WNS8IN3dVoWQrNcHX/O71H2Xc1PBebVOrUF48XURXm+SFPoofpSuJ8PCghXHswRhYS5FPRQI6zXK3yXkL2DrcassJBaknnsyc82HV6Ty5uF80QD2S5VPhOUezt0DO+7EoJPRK24VjufTuasekamzjsfu9G1sqMrmghfshXJ+slYNxTJkUSZE62WG6L1Z7uoSimc4ZzGSDq2YqGUuZiV6t/DDtvLC/ZLMiUzAsyRqdNnjh4yH6NmvR5led4/QFs83M7CmVuZHN0cmVhbQplbmRvYmoKMjUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA3MSA+PgpzdHJlYW0KeJwztjRQMFCwMFPQNTQ2VDCyNFYwNzNQSDHkAgqBWLlcMLEcMMvMEsQyNDdDYumaGUJlkVgg43K4YAbnwMzL4UoDAPG0FiMKZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY3ID4+CnN0cmVhbQp4nDO2NFAwULA0V9A1NDZUMDYwUTA3M1BIMeSCMXPBLLBsDhdMHYRlBmIYGZogscyAxoEl4QyQGTlw03K40gDOgxXTCmVuZHN0cmVhbQplbmRvYmoKMjcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMjcgPj4Kc3RyZWFtCnicNU87sgMhDOs5hS6QGYxtYM+zmVQv92+fZLINEv5I8vRERyZe5sgIrNnxthYZiBn4FlPxrz3tw4TqPbiHCOXiQphhJJw167ibp+PFv13lM9bBuw2+YpYXBLYwk/WVxZnLdsFYGidxTrIbY9dEbGNd6+kU1hFMKAMhne0wJcgcFSl9sqOMOTpO5InnYqrFLr/vYX3BpjGiwhxXBU/QZFCWPe8moB0X9N/Vjd9JNIteAjKRYGGdJObOWU741WtHx1GLIjEnpBnkMhHSnK5iCqEJxTo7CioVBZfqc8rdPv9oXVtNCmVuZHN0cmVhbQplbmRvYmoKMjggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDUgPj4Kc3RyZWFtCnicRVC7jUMxDOs9BRcIYP0se553SJXbvz1KRnCFIVo/kloSmIjASwyxlG/iR0ZBPQu/F4XiM8TPF4VBzoSkQJz1GRCZeIbaRm7odnDOvMMzjDkCF8VacKbTmfZc2OScBycQzm2U8YxCuklUFXFUn3FM8aqyz43XgaW1bLPTkewhjYRLSSUml35TKv+0KVsq6NpFE7BI5IGTTTThLD9DkmLMoJRR9zC1jvRxspFHddDJ2Zw5LZnZ7qftTHwPWCaZUeUpnecyPiep81xOfe6zHdHkoqVV+5z93pGW8iK126HV6VclUZmN1aeQuDz/jJ/x/gOOoFk+CmVuZHN0cmVhbQplbmRvYmoKMjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzOTIgPj4Kc3RyZWFtCnicPVJLbgUxCNvPKbhApfBNcp6p3u7df1ubzFSqCi8DtjGUlwypJT/qkogzTH71cl3iUfK9bGpn5iHuLjam+FhyX7qG2HLRmmKxTxzJL8i0VFihVt2jQ/GFKBMPAC3ggQXhvhz/8ReowdewhXLDe2QCYErUbkDGQ9EZSFlBEWH7kRXopFCvbOHvKCBX1KyFoXRiiA2WACm+qw2JmKjZoIeElZKqHdLxjKTwW8FdiWFQW1vbBHhm0BDZ3pGNETPt0RlxWRFrPz3po1EytVEZD01nfPHdMlLz0RXopNLI3cpDZ89CJ2Ak5kmY53Aj4Z7bQQsx9HGvlk9s95gpVpHwBTvKAQO9/d6Sjc974CyMXNvsTCfw0WmnHBOtvh5i/YM/bEubXMcrh0UUqLwoCH7XQRNxfFjF92SjRHe0AdYjE9VoJRAMEsLO7TDyeMZ52d4VtOb0RGijRB7UjhE9KLLF5ZwVsKf8rM2xHJ4PJntvtI+UzMyohBXUdnqots9jHdR3nvv6/AEuAKEZCmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDcgPj4Kc3RyZWFtCnicTVG7bUQxDOvfFFzgAOtreZ4LUl32b0PJCJDCIKEvKaclFvbGSwzhB1sPvuSRVUN/Hj8x7DMsPcnk1D/muclUFL4VqpuYUBdi4f1oBLwWdC8iK8oH349lDHPO9+CjEJdgJjRgrG9JJhfVvDNkwomhjsNBm1QYd00ULK4VzTPI7VY3sjqzIGx4JRPixgBEBNkXkM1go4yxlZDFch6oCpIFWmDX6RtRi4IrlNYJdKLWxLrM4Kvn9nY3Qy/y4Ki6eH0M60uwwuileyx8rkIfzPRMO3dJI73wphMRZg8FUpmdkZU6PWJ9t0D/n2Ur+PvJz/P9CxUoXCoKZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDkwID4+CnN0cmVhbQp4nE2NQRLAIAgD77wiT1BE0P90etL/X6vUDr3ATgKJFkWC9DVqSzDuuDIVa1ApmJSXwFUwXAva7qLK/jJJTJ2G03u3A4Oy8XGD0kn79nF6AKv9egbdD9IcIlgKZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiFtCNEGUglgQpWYmZhBJOAMilwYAybQV5QplbmRzdHJlYW0KZW5kb2JqCjMzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNDUgPj4Kc3RyZWFtCnicMzK3UDBQsDQBEoYWJgrmZgYKKYZclhBWLhdMLAfMAtGWcAoingYAn30MtQplbmRzdHJlYW0KZW5kb2JqCjM0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjU1ID4+CnN0cmVhbQp4nEWRS5IDIAhE956CI4D85DyZmlVy/+00mEw2dpeo/YRKI6YSLOcUeTD9yPLNZLbptRyrnY0CiiIUzOQq9FiB1Z0p4sy1RLX1sTJy3Okdg+IN566cVLK4UcY6qjoVOKbnyvqq7vy4LMq+I4cyBWzWOQ42cOW2YYwTo81Wd4f7RJCnk6mj4naQbPiDk8a+ytUVuE42++olGAeCfqEJTPJNoHWGQOPmKXpyCfbxcbvzQLC3vAmkbAjkyBCMDkG7Tq5/cev83v86w53n2gxXjnfxO0xru+MvMcmKuYBF7hTU8z0XresMHe/JmWNy031D51ywy91Bps/8H+v3D1CKZogKZW5kc3RyZWFtCmVuZG9iagozNSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MSA+PgpzdHJlYW0KeJxFkEsSwyAMQ/ecQkfwRwZ8nnS6Su+/rSFNs4CnsUAGdycEqbUFE9EFL21Lugs+WwnOxnjoNm41EuQEdYBWpONolFJ9ucVplXTxaDZzKwutEx1mDnqUoxmgEDoV3u2i5HKm7s75R3D1X/VHse6czcTAZOUOhGb1Ke58mx1RXd1kf9JjbtZrfxX2qrC0rKXlhNvOXTOgBO6pHO39BalzOoQKZW5kc3RyZWFtCmVuZG9iagozNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIxNCA+PgpzdHJlYW0KeJw9ULsRQzEI6z0FC+TOfO03z8uly/5tJJykQjZCEpSaTMmUhzrKkqwpTx0+S2KHvIflbmQ2JSpFL5OwJffQCvF9ieYU993VlrNDNJdoOX4LMyqqGx3TSzaacCoTuqDcwzP6DW10A1aHHrFbINCkYNe2IHLHDxgMwZkTiyIMSk0G/61y91Lc7z0cb6KIlHTwrvnl9MvPLbxOPY5Eur35imtxpjoKRHBGavKKdGHFsshDpNUENT0Da7UArt56+TdoR3QZgOwTieM0pRxD/9a4x+sDh4pS9AplbmRzdHJlYW0KZW5kb2JqCjM3IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggODAgPj4Kc3RyZWFtCnicRYy7DcAwCER7pmAEfiZmnyiVs38bIErccE+6e7g6EjJT3mGGhwSeDCyGU/EGmaNgNbhGUo2d7KOwbl91geZ6U6v19wcqT3Z2cT3Nyxn0CmVuZHN0cmVhbQplbmRvYmoKMzggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzYgPj4Kc3RyZWFtCnicTVBLbkQhDNtzilzgSSQhAc5D1VXn/tuxw1TtKoYYf0gP6bJVHutTYnWJ7PKlTZfKMnkVqOVP2/9RDAJu/9DIQbS3jJ1i5hLWxcIkPOU0Ixsn1ywfjztPG2aFxsSN450uGWCfFgE1W5XNgTltOjdAupAat6qz3mRQDCLqQs0Hky6cp9GXiDmeqGBKdya1kBtcPtWhA3FavQq5Y4uTb8QcWaHAYdBMcdZfAdaoybJZyCBJhiHOfaN7lAqNqMp5KxXCD5OhEfWG1aAGlbmFoqnlkvwd2gIwBbaMdekMSoGqAMHfKqd9vwEkjV1TCmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA0OSA+PgpzdHJlYW0KeJwzNrRQMFAwNDAHkkaGQJaRiUKKIRdIAMTM5YIJ5oBZBkAaojgHriaHKw0AxugNJgplbmRzdHJlYW0KZW5kb2JqCjQwIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTU3ID4+CnN0cmVhbQp4nEWQuRFDMQhEc1VBCRKwCOqxx9F3/6kX+Uq0bwAth68lU6ofJyKm3Ndo9DB5Dp9NJVYs2Ca2kxpyGxZBSjGYeE4xq6O3oZmH1Ou4qKq4dWaV02nLysV/82hXM5M9wjXqJ/BN6PifPLSp6FugrwuUfUC1OJ1JUDF9r2KBo5x2fyKcGOA+GUeZKSNxYm4K7PcZAGa+V7jG4wXdATd5CmVuZHN0cmVhbQplbmRvYmoKNDEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzIgPj4Kc3RyZWFtCnicLVI5jiQxDMv9Cn5gAOvy8Z4eTNT7/3RJVQUFqmzLPORyw0QlfiyQ21Fr4tdGZqDC8K+rzIXvSNvIOohryEVcyZbCZ0Qs5DHEPMSC79v4GR75rMzJswfGL9n3GVbsqQnLQsaLM7TDKo7DKsixYOsiqnt4U6TDqSTY44v/PsVzF4IWviNowC/556sjeL6kRdo9Ztu0Ww+WaUeVFJaD7WnOy+RL6yxXx+P5INneFTtCaleAojB3xnkujjJtZURrYWeDpMbF9ubYj6UEXejGZaQ4AvmZKsIDSprMbKIg/sjpIacyEKau6Uont1EVd+rJXLO5vJ1JMlv3RYrNFM7rwpn1d5gyq807eZYTpU5F+Bl7tgQNnePq2WuZhUa3OcErJXw2dnpy8r2aWQ/JqUhIFdO6Ck6jyBRL2Jb4moqa0tTL8N+X9xl//wEz4nwBCmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMTcgPj4Kc3RyZWFtCnicNVJLckMxCNu/U3CBzpi/fZ50smruv62EJyuwLUBCLi9Z0kt+1CXbpcPkVx/3JbFCPo/tmsxSxfcWsxTPLa9HzxG3LQoEURM9+DInFSLUz9ToOnhhlz4DrxBOKRZ4B5MABq/hX3iUToPAOxsy3hGTkRoQJMGaS4tNSJQ9Sfwr5fWklTR0fiYrc/l7cqkUaqPJCBUgWLnYB6QrKR4kEz2JSLJyvTdWiN6QV5LHZyUmGRDdJrFNtMDj3JW0hJmYQgXmWIDVdLO6+hxMWOOwhPEqYRbVg02eNamEZrSOY2TDePfCTImFhsMSUJt9lQmql4/T3AkjpkdNdu3Csls27yFEo/kzLJTBxygkAYdOYyQK0rCAEYE5vbCKveYLORbAiGWdmiwMbWglu3qOhcDQnLOlYcbXntfz/gdFW3ujCmVuZHN0cmVhbQplbmRvYmoKNDMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNyA+PgpzdHJlYW0KeJwzNrRQMIDDFEMuABqUAuwKZW5kc3RyZWFtCmVuZG9iago0NCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMSA+PgpzdHJlYW0KeJxFj8sNBCEMQ+9U4RLyGT6ph9We2P6v6zCaQUL4QSI78TAIrPPyNtDF8NGiwzf+NtWrY5UsH7p6UlYP6ZCHvPIVUGkwUcSFWUwdQ2HOmMrIljK3G+G2TYOsbJVUrYN2PAYPtqdlqwh+qW1h6izxDMJVXrjHDT+QS613vVW+f0JTMJcKZW5kc3RyZWFtCmVuZG9iago0NSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzOCA+PgpzdHJlYW0KeJw1Ujmu3UAM630KXSCAds2c5wWpfu7fhpRfCkO0VoqajhaVafllIVUtky6/7UltiRvy98kKiROSVyXapQyRUPk8hVS/Z8u8vtacESBLlQqTk5LHJQv+DJfeLhznY2s/jyN3PXpgVYyEEgHLFBOja1k6u8Oajfw8pgE/4hFyrli3HGMVSA26cdoV70PzecgaIGaYlooKXVaJFn5B8aBHrX33WFRYINHtHElwjI1QkYB2gdpIDDmzFruoL/pZlJgJdO2LIu6iwBJJzJxiXTr6Dz50LKi/NuPLr45K+kgra0zad6NJacwik66XRW83b309uEDzLsp/Xs0gQVPWKGl80KqdYyiaGWWFdxyaDDTHHIfMEzyHMxKU9H0ofl9LJrookT8ODaF/Xx6jjJwGbwFz0Z+2igMX8dlhrxxghdLFmuR9QCoTemD6/9f4ef78Axy2gFQKZW5kc3RyZWFtCmVuZG9iago0NiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI0OCA+PgpzdHJlYW0KeJwtUTmSA0EIy+cVekJz0++xy5H3/+kKygGDhkMgOi1xUMZPEJYr3vLIVbTh75kYwXfBod/KdRsWORAVSNIYVE2oXbwevQd2HGYC86Q1LIMZ6wM/Ywo3enF4TMbZ7XUZNQR712tPZlAyKxdxycQFU3XYyJnDT6aMC+1czw3IuRHWZRikm5XGjIQjTSFSSKHqJqkzQZAEo6tRo40cxX7pyyOdYVUjagz7XEvb13MTzho0OxarPDmlR1ecy8nFCysH/bzNwEVUGqs8EBJwv9tD/Zzs5Dfe0rmzxfT4XnOyvDAVWPHmtRuQTbX4Ny/i+D3j6/n8A6ilWxYKZW5kc3RyZWFtCmVuZG9iago0NyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE3MSA+PgpzdHJlYW0KeJxNkE0OQiEQg/ecohcwofMDj/NoXOn9t3bw+eKC9EshQ6fDAx1H4kZHhs7oeLDJMQ68CzImXo3zn4zrJI4J6hVtwbq0O+7NLDEnLBMjYGuU3JtHFPjhmAtBguzywxcYRKRrmG81n3WTfn67013UpXX30yMKnMiOUAwbcAXY0z0O3BLO75omv1QpGZs4lA9UF5Gy2QmFqKVil1NVaIziVj3vi17t+QHB9jv7CmVuZHN0cmVhbQplbmRvYmoKNDggMCBvYmoKPDwgL0JCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNzIKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnic49I1sjBVsDAwUMjl0jUyNAYzc7h0LY0VzAzNQCxDM0MY08jEUsHcGMw0NjaHiZoYmMIVQM2CqjU1gxgLZeZwpQEAk4MVTgplbmRzdHJlYW0KZW5kb2JqCjQ5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjEwID4+CnN0cmVhbQp4nDVQyw1DMQi7ZwoWqBQCgWSeVr11/2tt0DthEf9CWMiUCHmpyc4p6Us+OkwPti6/sSILrXUl7MqaIJ4r76GZsrHR2OJgcBomXoAWN2DoaY0aNXThgqYulUKBxSXwmXx1e+i+Txl4ahlydgQRQ8lgCWq6Fk1YtDyfkE4B4v9+w+4t5KGS88qeG/kbnO3wO7Nu4SdqdiLRchUy1LM0xxgIE0UePHlFpnDis9Z31TQS1GYLTpYBrk4/jA4AYCJeWYDsrkQ5S9KOpZ9vvMf3D0AAU7QKZW5kc3RyZWFtCmVuZG9iagoyMSAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMgL0NoYXJQcm9jcyAyMiAwIFIKL0VuY29kaW5nIDw8Ci9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0NiAvcGVyaW9kIDQ4IC96ZXJvIC9vbmUgL3R3byAvdGhyZWUgL2ZvdXIgL2ZpdmUgL3NpeCA1NgovZWlnaHQgNjkgL0UgOTEgL2JyYWNrZXRsZWZ0IDkzIC9icmFja2V0cmlnaHQgOTcgL2EgMTAwIC9kIC9lIDEwNSAvaSAxMDggL2wKL20gL24gL28gL3AgMTE0IC9yIC9zIC90IC91IF0KL1R5cGUgL0VuY29kaW5nID4+Ci9GaXJzdENoYXIgMCAvRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250RGVzY3JpcHRvciAyMCAwIFIKL0ZvbnRNYXRyaXggWyAwLjAwMSAwIDAgMC4wMDEgMCAwIF0gL0xhc3RDaGFyIDI1NSAvTmFtZSAvRGVqYVZ1U2FucwovU3VidHlwZSAvVHlwZTMgL1R5cGUgL0ZvbnQgL1dpZHRocyAxOSAwIFIgPj4KZW5kb2JqCjIwIDAgb2JqCjw8IC9Bc2NlbnQgOTI5IC9DYXBIZWlnaHQgMCAvRGVzY2VudCAtMjM2IC9GbGFncyAzMgovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250TmFtZSAvRGVqYVZ1U2FucyAvSXRhbGljQW5nbGUgMAovTWF4V2lkdGggMTM0MiAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTkgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMjIgMCBvYmoKPDwgL0UgMjMgMCBSIC9hIDI0IDAgUiAvYnJhY2tldGxlZnQgMjUgMCBSIC9icmFja2V0cmlnaHQgMjYgMCBSIC9kIDI3IDAgUgovZSAyOCAwIFIgL2VpZ2h0IDI5IDAgUiAvZml2ZSAzMCAwIFIgL2ZvdXIgMzEgMCBSIC9pIDMyIDAgUiAvbCAzMyAwIFIKL20gMzQgMCBSIC9uIDM1IDAgUiAvbyAzNiAwIFIgL29uZSAzNyAwIFIgL3AgMzggMCBSIC9wZXJpb2QgMzkgMCBSCi9yIDQwIDAgUiAvcyA0MSAwIFIgL3NpeCA0MiAwIFIgL3NwYWNlIDQzIDAgUiAvdCA0NCAwIFIgL3RocmVlIDQ1IDAgUgovdHdvIDQ2IDAgUiAvdSA0NyAwIFIgL3plcm8gNDkgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAyMSAwIFIgL0YyIDE1IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL0NBIDAgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTIgPDwgL0NBIDEgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTMgPDwgL0NBIDAuOCAvVHlwZSAvRXh0R1N0YXRlIC9jYSAwLjggPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgL0YxLURlamFWdVNhbnMtdW5pMDMwMiA0OCAwIFIgL00wIDEyIDAgUiA+PgplbmRvYmoKMTIgMCBvYmoKPDwgL0JCb3ggWyAtOCAtOCA4IDggXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMSAvU3VidHlwZSAvRm9ybQovVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJxtkEEOhCAMRfc9RS/wSUtFZevSa7iZTOL9twNxQEzdNNC+PH5R/pLwTqXA+CQJS06z5HrTkNK6TIwY5tWyKMegUS3WznU4qM/QcGN0i7EUptTW6Hijm+k23pM/+rBZIUY/HA6vhHsWQyZcKTEGh98LL9vD/xGeXtTAH6KNfmNaQ/0KZW5kc3RyZWFtCmVuZG9iagoyIDAgb2JqCjw8IC9Db3VudCAxIC9LaWRzIFsgMTAgMCBSIF0gL1R5cGUgL1BhZ2VzID4+CmVuZG9iago1MCAwIG9iago8PCAvQ3JlYXRpb25EYXRlIChEOjIwMjAxMjE2MTUzNDMzKzAyJzAwJykKL0NyZWF0b3IgKE1hdHBsb3RsaWIgdjMuMy4yLCBodHRwczovL21hdHBsb3RsaWIub3JnKQovUHJvZHVjZXIgKE1hdHBsb3RsaWIgcGRmIGJhY2tlbmQgdjMuMy4yKSA+PgplbmRvYmoKeHJlZgowIDUxCjAwMDAwMDAwMDAgNjU1MzUgZiAKMDAwMDAwMDAxNiAwMDAwMCBuIAowMDAwMDE0NjM2IDAwMDAwIG4gCjAwMDAwMTQwOTMgMDAwMDAgbiAKMDAwMDAxNDEzNiAwMDAwMCBuIAowMDAwMDE0Mjc4IDAwMDAwIG4gCjAwMDAwMTQyOTkgMDAwMDAgbiAKMDAwMDAxNDMyMCAwMDAwMCBuIAowMDAwMDAwMDY1IDAwMDAwIG4gCjAwMDAwMDAzOTcgMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDAzMTA0IDAwMDAwIG4gCjAwMDAwMTQzODIgMDAwMDAgbiAKMDAwMDAwNDEwMSAwMDAwMCBuIAowMDAwMDAzODkzIDAwMDAwIG4gCjAwMDAwMDM1NzAgMDAwMDAgbiAKMDAwMDAwNTE1NCAwMDAwMCBuIAowMDAwMDAzMTI1IDAwMDAwIG4gCjAwMDAwMDM0MTEgMDAwMDAgbiAKMDAwMDAxMjcwNSAwMDAwMCBuIAowMDAwMDEyNTA1IDAwMDAwIG4gCjAwMDAwMTIwNDEgMDAwMDAgbiAKMDAwMDAxMzc1OCAwMDAwMCBuIAowMDAwMDA1MTk2IDAwMDAwIG4gCjAwMDAwMDUzNDcgMDAwMDAgbiAKMDAwMDAwNTcyNCAwMDAwMCBuIAowMDAwMDA1ODY3IDAwMDAwIG4gCjAwMDAwMDYwMDYgMDAwMDAgbiAKMDAwMDAwNjMwNiAwMDAwMCBuIAowMDAwMDA2NjI0IDAwMDAwIG4gCjAwMDAwMDcwODkgMDAwMDAgbiAKMDAwMDAwNzQwOSAwMDAwMCBuIAowMDAwMDA3NTcxIDAwMDAwIG4gCjAwMDAwMDc3MTEgMDAwMDAgbiAKMDAwMDAwNzgyOCAwMDAwMCBuIAowMDAwMDA4MTU2IDAwMDAwIG4gCjAwMDAwMDgzOTAgMDAwMDAgbiAKMDAwMDAwODY3NyAwMDAwMCBuIAowMDAwMDA4ODI5IDAwMDAwIG4gCjAwMDAwMDkxMzggMDAwMDAgbiAKMDAwMDAwOTI1OSAwMDAwMCBuIAowMDAwMDA5NDg5IDAwMDAwIG4gCjAwMDAwMDk4OTQgMDAwMDAgbiAKMDAwMDAxMDI4NCAwMDAwMCBuIAowMDAwMDEwMzczIDAwMDAwIG4gCjAwMDAwMTA1NzcgMDAwMDAgbiAKMDAwMDAxMDk4OCAwMDAwMCBuIAowMDAwMDExMzA5IDAwMDAwIG4gCjAwMDAwMTE1NTMgMDAwMDAgbiAKMDAwMDAxMTc1OCAwMDAwMCBuIAowMDAwMDE0Njk2IDAwMDAwIG4gCnRyYWlsZXIKPDwgL0luZm8gNTAgMCBSIC9Sb290IDEgMCBSIC9TaXplIDUxID4+CnN0YXJ0eHJlZgoxNDg1MwolJUVPRgo=\n", "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2020-12-16T15:34:33.228346\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.3.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 1000 # number of samples for input signal\n", "\n", "# generate input signal\n", "# normally distributed (zero-mean, unit-variance) white noise\n", "np.random.seed(1)\n", "x = np.random.normal(size=N, scale=1)\n", "# impulse response of the system\n", "h = np.concatenate((np.zeros(20), np.ones(10), np.zeros(20)))\n", "# output signal by convolution\n", "y = np.convolve(h, x, mode=\"full\")\n", "# add noise to the output signal\n", "y = y + np.random.normal(size=y.shape, scale=0.1)\n", "\n", "# zero-padding of input signal\n", "x = np.concatenate((x, np.zeros(len(h) - 1)))\n", "# estimate transfer function\n", "H = np.fft.rfft(y) / np.fft.rfft(x)\n", "# compute inpulse response\n", "he = np.fft.irfft(H)\n", "he = he[0 : len(h)]\n", "\n", "# plot impulse response\n", "plt.figure()\n", "plt.stem(he, label=\"estimated\")\n", "plt.plot(h, \"g-\", label=\"true\")\n", "plt.title(\"Estimated impulse response\")\n", "plt.xlabel(r\"$k$\")\n", "plt.ylabel(r\"$\\hat{h}[k]$\")\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**\n", "\n", "* Change the length `N` of the input signal. What happens?\n", "* Change the variance $\\sigma_n^2$ of the additive noise. What happens?\n", "\n", "Solution: Increasing the length `N` of the input signal lowers the uncertainty in estimating the impulse response. The higher the variance of the additive white noise, the higher the uncertainties in the estimated impulse response." ] }, { "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*." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.13" } }, "nbformat": 4, "nbformat_minor": 1 }