{ "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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEYCAYAAABfgk2GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3Xl8VNX5+PHPk3Wyh0BWQtjDvi+CyCa4oyhqtZsiKmqp1larpdalVVt/tdXWr9alFZeqVatVcccFXEBFQHaQfQkEQhayb5M5vz/OJISQkMk6M8nzfr3mlZl779z7nEzyzLnnnnuOGGNQSinV8QV4OwCllFLtQxO+Ukp1EprwlVKqk9CEr5RSnYQmfKWU6iQ04SulVCehCV+pdiQivUTEiEiQt2OpTUQmi8j33o5DtS1N+KpeIrJHREpFpEhEDovIsyIS6e24qolIiIj8VUQy3DHuEZG/tfEx54rIl218jOrfe6GIHBWRFSJyvYi06f+qMeYLY8yAOnHMbMtjqvanCV+dzPnGmEhgNDAW+J03gmigNrwQG9N4IAqYBqxpx7DqJSKBrbCb840xUUBP4AHgduDpVtiv6uQ04atGGWMOAO8DQwFEJEVEFotIrojsEJFr3csd7tppN/frO0TEKSLR7tf3VtfCRSRURP4iIvvcZxBPiEiYe900d839dhE5BDxTT1jjgDeMMQeNtccY83z1SncNdaGIbBaRPBF5RkQctdbPEpG1tWrRw2ut6yEi/xORIyKSIyKPisgg4AlgovuM4qh722dF5HEReU9EioHpInKeiHwnIgUisl9E7mnm7z3fGLMYuAy4UkSqf/+e/O5uEZEsEckUkatqle1c9++kUEQOiMittd/nfv5vIA14213W20TkXRG5sXZ8IrJeRC5qTtmUd2jCV40SkR7AucB37kUvAxlACnAJ8EcROd0YUwZ8C0x1bzcV2AtMqvX6M/fzB4B0YCTQD+gO3FXrsElAHLaWO7+esL4GfiUiPxORYSIi9WzzY+AsoK/7WL9zl2cUsAi4DugKPAksdifSQOAdd9y93HG9bIzZAlwPfGWMiTTGxNY6zo+A+7FnGl8CxcAVQCxwHnCDiFxYT3weMcasxP6+J7sXefK7i3Evvxp4TES6uNc9DVznPoMYCnxaz/F+CuzDfYZnjPkz8Bzwk+ptRGSEe//vNrdcyguMMfrQxwkPYA9QBBzFJr9/AGFAD6AKiKq17Z+AZ93P7wUeAYKAQ8AvsAnKAZRiE6xgk2LfWvuYCOx2P58GVACOk8QXCCwAlgPlwEHgyjrxX1/r9bnATvfzx4F76+zve+wX0kTgCBBUzzHnAl/WWfYs8Hwjv8u/AQ+7n/cCTH37rxX3zHqWfw3c4eHvrrT2/oEsYIL7+T7sF110nf1PAzIaisP9+eUB/d2v/wL8w9t/p/po2kNr+OpkLjTGxBpjehpjfmaMKcXW6nONMYW1ttuLre2BrcFPw7b7bwA+wibSCcAOY0wOEA+EA6vdTSpHgQ/cy6sdMfaMoV7GmCpjzGPGmEnYmvT9wCJ300u1/XViTHE/7wncUn1s9/F7uNf3APYaY5we/YZOPA4icoqILHU3CeVjzwy6NWF/9ekO5OLZ7y6nTvwlQPUF94uxX357ReQzEZnoycHdn8UrwE/cF5B/CPy7JQVS7U8Tvmqqg0CciETVWpYGHHA/XwEMAC4CPjPGbHavP5djzTnZ2FroEPcXSqwxJsbYC8TVPB7G1RhTaox5DFsDHVxrVY86MR50P98P3F/r2LHGmHBjzH/c69IauFDcUEx1l78ELAZ6GGNisG3/9TU5eURExmET/pd49rtrkDHmW2PMbCABeBN4taFN61n2HLaZbAZQYoz5qolFUV6mCV81iTFmPzap/8l9kXY4tp34Bff6EmA1trmlOsGvwNZyP3Nv4wL+CTwsIgkAItJdRM7yNA4Rudl9oTFMRIJE5EpsG/p3tTZbICKpIhKHbQ55xb38n8D17pq4iEiE+0JrFLASyAQecC93iEj1NYjDQKqIhDQSXhT2LKhMRMZj2/ibTESiRWQW9prJC8aYDS353YntyvpjEYkxxlQCBYCrgc0PA31qL3AneBfwV7R275c04avm+CG2Lfog8AZwtzHm41rrPwOCscmz+nUU8HmtbW4HdgBfi0gB8DH2zMBTJdjEcwhb610AXGyM2VVrm5eAJcAuYCdwH4AxZhVwLfAo9qxgB7Z9HmNMFXA+9mLoPuzF0svc+/sU2AQcEpHsk8T2M+APIlKIvZjaUC26IW+737sf+0X1EHBVrfUt+d39FNjjft/12Bp7ff4E/M7dbHRrreXPA8Nwf8Er/yLG6AQoquMRkT3ANXW+iFQLicgVwHxjzGnejkU1ndbwlVIeEZFw7NnLU96ORTVPixK+iFwqIptExCUiY+usWyj2ppzvm9I2q5TyPe7/4SPYtv2XvByOaqYWNem4u8C5sDeu3OpuG0VEBgP/wd72noJtY0x3t48qpZTyghbV8I0xW4wx9Y2wNxt7d2K5MWY39gLT+JYcSymlVMu01RCt3bF3BlbL4NiNOccRkfm4b513OBxjevboTmTRbgBKw5JwBvnMAI3N4nK5CAjoeJdKOmK5tEz+oyOWqyVl2rZtW7YxJr6x7RpN+CLyMXZsjrruMMa81ZzgajPGPIX7ItCAAQPM1pWfwMND7MqJc+Gs+1t6CK9atmwZ06ZN83YYra4jlkvL5D86YrlaUiYR2evJdo0mfGNMc8bEPsDxdzmmcuxOzJOrKDn2POPbZhxaKaVUfdrqnGgxcLl79MHeQH+O3YRzcpXuhB/XFw6uBWdFG4WolFKdS0u7ZV7kHkN7IvCuiHwIYIzZhL27cDN2YKcFHvfQqU74faZCVTkc2tCSEJVSSrm1tJfOG8aYVGNMqDEm0RhzVq119xtj+hpjBhhj3vd4p9UJv/cU+1ObdZRSqlX43mXuylL7M64vRKdChmctQUoppU7O9xJ+9UXbkAhIHas1fKWUaiW+l/Crm3SCwyB1HBzdB4WHvRuTUkp1AD6c8MMh2T2v9JEt3otHKaU6CN9O+HF97fOcnd6LRymlOgjfS/gVJRAQBEEhEJUMQQ7I3dX4+5RSSp2U7yX8ylJbuwcICIC4PprwlVKqFfhgwi+xF2yrxfXRJh2llGoFPprww4+9jusDebvBpUPpK6VUS/hgwi89PuF37QtVFVDg2dhrSiml6ud7Cb+iGELq1PBB2/GVUqqFfC/hV5bWacPXrplKKdUafDDhF0NwxLHX2jVTKaVahQ8m/Do1fO2aqZRSrcJHE3748cu0a6ZSSrWY7yX8uhdtQbtmKqVUK/C9hF+3SQe0a6ZSSrUC30v4VeXHX7QF7amjlFKtwKcSvhiXfVK3hq998ZVSqsV8KuGDsT/qtuFHJUNQmCZ8pZRqAZ9K+Mdq+HUSfkAAxPXWhK+UUi3gUwm/poZft0kHoFt/yNKZr5RSqrl8KuEfq+FHnLgyeYTtmll6tH2DUkqpDsKnEv5Ja/gpo+zPzHXtF45SSnUgPpXwa2r4dS/aAiSPtD8Pftd+ASmlVAfiUwn/WA2/noQfHgexaZC5tn1DUkqpDsKnEr6YkyR8sM06WsNXSqlm8amEDw10y6yWPBLy9kBpXrtFpJRSHYVPJfxjNfx6LtrCsQu3B7VZRymlmsrHEn5jNfwR9qe24yulVJP5VMIHY2e3CmggrPA46NJL2/GVUqoZfCrhi3E1XLuvljxSm3SUUqoZfCrhg2k84aeMgqN7oSS3fUJSSqkOwqcSvhhX/Tdd1ZbivgFL2/GVUqpJfCrh2xp+Az10qiWPBAT2fdMuESmlVEfRooQvIg+KyFYRWS8ib4hIbK11C0Vkh4h8LyJnebQ/T9rww2Jts86uZS0JXSmlOp2W1vA/AoYaY4YD24CFACIyGLgcGAKcDfxDRAIb350HbfgAfaZBxrdQVtDMsJVSqvNpUcI3xiwxxjjdL78GUt3PZwMvG2PKjTG7gR3A+Mb2J8aDJh2AvtPBVMHe5c0LXCmlOqGgVtzXPOAV9/Pu2C+AahnuZScQkfnAfIARycEcyiti67JlJz2QuCo5LSCEzM9fYEemB18QXlRUVMSyRsrjjzpiubRM/qMjlqs9ytRowheRj4GkelbdYYx5y73NHYATeLGpARhjngKeAhjTPdQkpfYiadq0xt94YDKp+dtJ9WRbL1q2bBnTfDzG5uiI5dIy+Y+OWK72KFOjCd8YM/Nk60VkLjALmGFM9WA4HAB61Nos1b2sEa76Z7uqT59p8NGdUHAQolM8e49SSnViLe2lczZwG3CBMaak1qrFwOUiEioivYH+wMpG9+dpGz7YdnzQ3jpKKY8VllWyYmc2xeXOxjfugFrahv8oEAp8JCIAXxtjrjfGbBKRV4HN2KaeBcaYqsZ314SEnzAEwrvZhD/yR82LvhEul0EE3GWjwunipW/2siu7mH4JkaQnRjGmZxeCA1vW2cnlMlS6XIQEBtQcqz2VVVax/XARWw8V0Cc+gjE94zx6nzGGd9ZnUljm5NKxqS3+PbSW3OIKIkODCAlq/Xgy8kpYs+8o/eIj6ZcQ2SbHaK5dR4r4aPNh4iJCSIsLp29C5Em3N8bw7Z48ln6fxcWjU+nXyPatpcplKCpzEhMe3OJ9uVyGgIDj/2eqGxqq/5eyCst4b30mSzYfZuXuXJwuw8CkKJ6bN57EaEeLY6hWVlnF81/tYV1GPtPS4zlzcFKTy7jpYD6lFVWMSutCYEDr54IWJXxjTL+TrLsfuL/JOw3xsEknIAD6TLUJ3xhoxURZVlnFM8v38PiyHUSHBTNnVHfSk6J4+KNt7DxSTHhIICUV9vsrISqUH5/Sk/NHJON0GQrLKgkJDCQ+KpTY8GDyylyszzjKgbxSdmUXsy+nhMToUEamxRIbHsK76zN5a+1BsovKAYgICWRSv25cOKo7pw9MwBF8Ym/WTQfzue+dLWTmlzIlPZ7TByYwoU/Xerc9GZfL8Min23n00x04XaZm+U8mpLHwnEFEhB7/51Hl3iYwQMgqcXHFopV8sT0bgGeW7+bu84dwWv9uTYqhWk5ROUs2H2ZvTgn780oQoHtsGKlx4Zw3LJm4iJAT3pORV8JL3+wjPiqU/glRHMwv5bXVGazcnUtIYACDkqMY2j2G4akxDO0ew4DEKIJqfSkVlzvZn1dCpdN+4e4vdJGZX0qX8JATfpcul+HFb/byp/e31nz2wYHC1PR45k/py7heXU74sq6scpFXXEFeSSW5xRXklVSQW1yBIziQlFgHCVEOyp1VFJU5ySupJKuwjOzCcvonRnH6wISa37/LZSgsc1JU4aSozEl2UTlHCssprnASGRqEIziQxesO8t6GTIw5LgS6hQkTM79jxsAEzhicSERoEPkllSzblsWzK/bw3b6jNZ/fnbMG86PxaQBkFZbz6dYs3t94iE0H8hmYHMXotC6M7x3HuF5xNb+f/JJKSiurSIrxLHFuPljA7a+vZ8OBfM4eksSNM/oxJCUGsH9fmfml7M8tJb+0kqAAIShQ6B4bRu9uEQQFBnC4oIxv9+Ty5uZyHlz/BVsPFTI8NYbzh6fQs2s4H246xIebDuNyGVLjwgkLDmDt/qO4DPRPiOTqyb3p1TWC+97ZzJx/rOD5q8fTN755X3TlzioO5ZeRU1zB9sOFPPLJDg4cLSUuwv5fLwzYwMxBiVxxak8m9ukKwNGSSgIChJiw478IisudPL+5nE8/+BKALuHBTBuQwMgesQxOiWZwcvQJ/4/NIabuX4gXjU0JNKvefhrGzAXg4NFS/t8HW+nTLZLRPWNJinZQWO4kv7SSPdnFRHz/Oj/Ydy+/jPgzn5X1oUdcOBeOTOHcYck4ggMpr6wit6TCJpHcErIKy8kuKqfC6WJkj1gm9OlKlCOIPTkl7Mst4Yh7/dKtWWTmlzFtQDxVLsOXO7IxBnp1Deeu8wczfUAChwvKWbv/KP9ZuY/Pth3xuIzdIkPILa6gOr8GBwqnD0xgeGos5ZVV5BRXsGTzYY4UlhMSFMDQlGhG9IglMdpBcGAAO48U8fLKfcSGhzAiNYavduVQVukiLNh+UZw+MIHTByY0+g9YXO7kllfX8cGmQ8wansx5w5LpnxjJK9/u519f7iYlJoxZI5IZ3j2WwADhg42ZfLIli8JyJ47gACqdLsJCgrjt7AEkRTu4790t7Mst4cKRKfz+gqHEhAfjchm+3pVDZn4ZIUEBOIID6RoZQnxkKNFhwYjYOJ5bsZfnv9pDSUUVwe5/cPv5l1FR5SIyNIj5U/pw9Wm9a/7oV+/N47p/ryK7qOK4cvXpFsEFI1MorahifUY+Gw/mU1hmT9/DQwIZlRbLgMRoNh7IZ82+vOO+6KoFBQin9InjjEGJxEc52J5VyBfbs1m9N4/J/btx88z+HDhaxvr9R3l9TQZ5JZUMTo4mJdZBaHAgxeVOdh0pJiOvhHp275HQoABGpMaSXVxORm4pFVWuk24fFRrETyb25IqJPSmvdLEvt4TvDxXy4ept7CsJIquwnLDgQPolRLLpYD4uAz3iwpg/uQ9T0uO5442NfLkjmx5xYeQUVdR8qfXsGs6Ynl34/lAhWw8VUuUyNrYesRwuKGNvjm3JHZUWy4Uju9M9Noz80koqq1xMTo+v+Swz80t5bsVe/vXFLmLDgzlvWDL/++4AhWVOukaEUFZZRWllVYO/r5CgALqEB3O4wFaMQgJhTM+uDEyO4ptduWzOtPfkRIYGccbgRGLCgtmfW0JuSQWT+8dz/vBk+idG1exvQ0Y+Vz27ktKKKs4amsR5w5Iprqjig42ZfLMrl5TYMIakRDMkJZrBKTEMSo6ipKKKPdnFbDpYwGfbjvDVzhxKK481XAxOjuaO8wZxat+urMvI5511B3ltTQZHSypJjnFQVO6ksMyJiN32lN5dCQ0OoKC0ks+2HeFAXilXTerNqLRYPt2axWfbjpBbXFHz93DmkCTmjOpOubOKVXvyyC2p4IE5wwkJCkBEVhtjxjb2d+V7Cf+Dl2H4pQD86f0tPPnZLkQ4oeYCkOSo5DOuZXn0uXzU61bW7T9a88HXJzQogG6RoQAcOFp6wnoRiAsPoX9iJDfPTGeC+1v5UH4ZGw/kMzm9G6FBJ9aidx4pYtWeXCJCg4hyBFNeWcWRonKOllRyJGM3k8YMJznGQa9uEUSGBlFc7mTDgXwOF5QxpX88XerUXqtchhU7s/ns+yOsz8hnw4H8mj+sAIGfTOjJr85IJzbc/qN8tSuHpVuz+GRLVk25BidHc87QJC4ek0pKbBjGGPbnlrJ6Xy5r9h5l2bYsDuSV8ttzB3H1ab2Pq52u2pPLn97fyoaM/JpEExsezBmDEuneJYySiioO7N/P7y6fTHKM/Ycud1bxj6U7eXTpDuIjQ/nh+DQWrzvAziPFDX4etX/vs4ancMPUvgxIiqo5lXW5DN8fLuThj7axZPNhIkICGd2zC33jI3lp5T6SYxw8feVYosOC2XG4iPDQIEakxhxXFpfLsDe3hPUZR/lu31G+3ZPL94cKGZwSzal9uzEkJRpHcCBBgcKq79aT2iedPdnFfLzlcE3sItAzLpzrp/blsnE9jtt/aUUVr63ez9vrMykqc1LmrMIRFEif+Aj6dIsgIdpBXEQIseHBxEWE0CU8hNKKKg4eLSWrsBxHcABRjmBiwoJJjHYQGx7Mmr15vL/xEOszjpIcE0aPuHDio0KJDA0kIjSIbpGhxEeFEhESRFF5JYVlTvomRBLtOLH5YNmyZUyZMpVVe/N4a+0Bth0uZEKfrjW1x9q/6+e+2sOKnTmkdgkjLS6cU3p3ZVByVE15i8udrNyTyxfbslm9L4+UGAdDu8cgAovXHmTrocITjj+uVxdCgwJZvtNWmi4encqdswYRGx5CfmklL3y9l4NHS3EEBxIeEkhKbBg9uoTTJSIYlwsqqqrYm2O/vLIKyxmSEs24XnEc2f4dM0+fXnOcHVlFHDxayvjecR6f6e7LKeGRT7ezZNMhCtyVgvioUCb368bhwjI2HSzgaEllve9Niwtn2oB4hqfG0jUihK6RIQxNiTmheamssoq31x3k061ZxEeFkhYXTklFFV/tzGH1vjyMMUQ7gkmNC2dWShnXXjSj5r3GGA4VlLEls4Bl3x9h8bqDNfFU58R3bzqNISkxfpzwP3kLBs2iymWY9MCnDE6J5m+Xj2T9/nxySyqIcgQR7Qiy/wSRoch/59obsH61FQKD2H64sKbG7QgOJCYsmJ5dw0mLCycmLLjmj/dQfhnf7M6hvNJFz67h9OwaQbfIkONO+VtDa3S1crkM5U4XFU4XgYFCZAOndsYYtmcV8enWLD7Zcphv9+QhAiNSY8nIKz2u2WhkWiw3TO130iaYCqeLbYcLKSp3nnCtoqFyrc84yq9eXceOrCJGpMbU1FgqnC5KKqrILa7gSGE5BWXVf7jClP7djqt91WfNvjz+tyaDVXvy+P5wIaf0juPxH4854cvSE8aYeq+V1C3T7uxiSiqc9I2PbHJzma9oz+6LO7IKKS6vIjY8mMoqwwcbM1m87iAVThezR3bn4tGppHX14E56D7RmuSqcLr7elUNYSCCja7WdG2PIzLeJf2tmARGhQfSOj6BffCQ94lpejrrXCBsrU7mzihU7c4gJC6bS6eKyp77mmbnjmD4wweOE35o3XrUO92iZX+/K4VBBGXecN4hoR3DDiWnYJbD5Tdj9GfSbQf/EqEaTB0BSjIPZI+u9F8znBAQIYSGBhIWcPOmICOmJUaQnRnH91L7szy3htdUZfLbtCFPSuzE6rQtjenYhPTHKowtCIUEBDO0e06RYh6fG8u5Np3E4v5wecWGtdhF6dFoXRqd1AWytKTSo+Re4PX1f724eXk9SAPRLOP7/7uen9+fnp/f3UjSeCwkKYEp6/AnLRYSU2DBSYsM4Y3Biqx+37tlAY0KDApk+IAGwzd0AhwrKmrQP30v47rF0/rfmAFHu9riT6ncGhMbAhteg34yTb9vJ9IgL55dnpPPLM9Lb9bihQYGtVpOrj7/WtpVqLfFRoYjA4SYmfN/pU1YtOIxS98WTc4YlNf7PHeyAQefDlreh8sR2eaWU6miCAwPoGhHaERJ+BEs2H6K4ooqLRqU2vj3YZp2KQti+pG1jU0opH5EYHVrTa8lTPpjww3jjuwN0jw3jlN6e3QBE7ykQkQAb/tu2sSmllI9IinZwKN/fa/gh4azak8fpAxM8v6gREAhD58C2JVCW37bxKaWUD0iIdvh/k44JCqO4wkmXpt52PfQSqCqHLe+0TWBKKeVDkqId5BRXUOE8+U15tflYwhfKXEEYA2EhTexAlDoWYnvCxtfaJjSllPIhidH2JtKsQs9r+T6V8I0IJe47SiNCm9j1TsRevN21DIqyWj84pZTyIYnu4VOacuHWpxJ+SXhqzRgeYc3paz3sUjAu2PRmK0emlFK+JTGqOuH7aQ3fFRBSk/DDm9qkA5AwCBKHam8dpVSHlxTj5wkfoKTCPbJhU5t0qg29GDJWQu6uVoxKKaV8S5fwYEICA5o0vIIPJnx3Db+5t88PvwwCQ+CLh1oxKqWU8i0iQkJ0KFn+2oYPtKxJByCmO4y9Gta+BNnbWzEypZTyLU29+coHE34Lm3QAJt8CQQ5Y2vQJt5RSyl8kNvHmKx9M+NU1/BYk/Mh4mPgz2PQGZK5rpciUUsq3JEY7OkgbfnObdKqdeiM4YuHje+qfLksppfxcYnRoTc70hO8l/PJj84+2iCMGpv0Gdn4K619phciUUsq3eDp5fDXfS/iVdiLr4NaYanD8fOgxAd6/DQoyW74/pZTyIQlR/p7wy50tb86pFhAIsx8DZwW8c7M27SilOhT/r+FXVLW8Oae2bv1gxl2w7QNY83zr7VcppbysegA1T/lewq9s5YQPcMr10HsqvPdrOLi2dfetlFJeEh4SRJTD8xYR30v4rdmkUy0gAC5ZBBHx8MpPoSS3dfevlFJekhjtebOO7yX8iirCWruGDxDRDS57HooOwWvzwOV5VyallPJVSf6c8Esrq4hoi4QP0H0MnPNn2LUUvnqsbY6hlFLt6F9XjvV4W59L+MVt0aRT25i5MHAWfHovHNrYdsdRSql24GjCQJM+l/BLW7uXTl0icP7f7V24/5sPTs9HmlNKKX/mcwm/uK0TPtj2/NmPQtYm+OA32j9fKdUp+FzCL62oavoE5s2RfpYdb2fVIlh8o17EVUp1eO2QWT1ngIoqV9tdtK3rjHvtMMqfPwjlBTDnnxDUtBsZlFLKX/hWwne3rLRJt8z6iMDpv7Pt+UvugOIcuPwFCOvSPsdXSql21KImHRG5V0TWi8haEVkiIinu5SIij4jIDvf60Z7sr7olPSK0nb+HTv25rd1nrISnz4S8Pe17fKWUagctbcN/0Bgz3BgzEngHuMu9/Bygv/sxH3jck5253Bm/zS/a1mf4D+Cnb0BRFvxjIrx6pZ1ARXvxKKU6iBYlfGNMQa2XERyrpM8GnjfW10CsiCQ3vj/7M6y5E5i3VK/T4NpPYcQPYe9y+O9ceOkyqKr0TjxKKdWKxLSwS6KI3A9cAeQD040xR0TkHeABY8yX7m0+AW43xqyq5/3zsWcBxCUkj4m66p/cNs7B4K5eSvrVTBUpBz8kffuTHEw+i23pN9g2/yYqKioiMjKyDQL0ro5YLi2T/+iI5WpJmaZPn77aGNPoLbeNNpaLyMdAUj2r7jDGvGWMuQO4Q0QWAj8H7m5KoMaYp4CnANL6phuACeNGMzrNFy6czoCPw0n58mFShk2xbf1NtGzZMqZNm9b6oXlZRyyXlsl/dMRytUeZGk34xpiZHu7rReA9bMI/APSotS7VveykXMa2MUW0Rz98T51+F+TstL14VvyfHY9n6BwYdom3I1NKqSZpaS+d/rVezga2up8vBq5w99aZAOQbYxqdY7C6cckrF20bEhAAc56Cc/8CvafA4Y3w+tV2cnSXy9vRKaWUx1palX5ARAYALmAvcL17+XvAucAOoAS4ypOdudq7H76ngsNg/LX2UeWE926FLx+G/Aw7haLerKWU8gMtSvjGmIsbWG6ABU3fn/3pU01fR9MTAAAfLElEQVQ6dQUGwayHIbYHfPIHOyTDxU/bMwGllPJhPpVZDbYjjCPYx5OnCEy+BSQQPr4buvSEmfd4OyqllDopn0r4LgPhwYFIM7o/esWkX9i7cr98GGLTYOw8b0eklFIN8qmqtDG0z0iZrUXEXsztdwa880t443ooOuLtqJRSql4+lfBd+FgPHU8EBsFl/4bTfgUbXoNHx8K6l70dlVJKncCnEr4xfpjwwfbimXk33LAcEgbDG9fZMfYrS70dmVJK1fCphO8yxj8TfrX4AXDl2/aC7prn4ekziCrY7u2olFIK8LGEb/DC0MitLTAIZtwFP3oVCjIZs+ZWeG0e5O72dmRKqU7OpxK+y3hxpMzWln4W3PQde3r+AL5/H/5vjJ00/fBmb0emlOqkfCrh+20bfkMc0ezp/WO4cQ1MuAG2vAOPT4T3f2Pv2FVKqXbkWwkfCPf3Jp36RCfDWffDLzfC+PnwzePw0g+gLN/bkSmlOhGfyq7VN151WOFxcO6DkDgU3v0VPDkFBl0AaRPt5CuOaG9HqJTqwLSG7w1jroSfvgmRifD14/DyD+GRUbD2P8cGFFJKqVbmUwkfOlgb/sn0ngxXL4GF++GKxRDXG968Hp6dZUfhVEqpVqYJ39uCw6DPVJi3BM7/O2Sug0VnQ/YOb0emlOpgfDDhd4ImnfoEBMCYuXDVu/YO3WfOhsz13o5KKdWB+GDC72Q1/LqSR8C8DyAwFJ45FzYv9nZESqkOQhO+L+rW37bvx6fDqz+FJb/TfvtKqRbzwYTfSZt06orpDle9D+OutZOn/3Ma7Pva21EppfyYz2VXreHXEhQK57knT/9gISw6CwbOAkcslOdD4jA7UFugz32MSikf5HOZQhN+PQZfAP1mwBcPwepnITAEQsJhy9uw90u49Dl7U5dSSp2EDyZ8nwvJN4REwIw77aPady/ambaemgo/fAUSB3svPqWUz/O9NvxQreF7bNSPbTu/s8I29+z4xNsRKaV8mO8l/I48lk5bSB0D134CMT3gxUth9XPejkgp5aN8KuELEBToUyH5h5hU23e/73R4+yZ4+2aoLPN2VEopH+NT2VVbc1rAEW3b8SfdDKufsU08eXu9HZVSyof4VMJPivCpcPxPYBCc8Xu4/D92SsV/zYCD33k7KqWUj9AM2xENPBeu+RiCwuCZ82D7R96OSCnlAzThd1Tx6XDNR9C1r51d69lZsPzvcGgDuFx2G1cV7P4cVjwK+76BqkrvxqyUalPa6b0ji0qCq96DL/8G2z6Aj+6yD0csdB8NhzZCcdax7UOi4JTr4PTfgYj34lZKtQlN+B1daNSxG7byD8CeL2DvcshYBWkTYOgc6HEK7F8JG1+HL/5ix+ifcqu3I1dKtTJN+J1JTHcYcbl91DXkQju/7hvXwaf3QlSyvbFLKdVhaMJXxwQEwOzHbDPP4huhosiO1hmgl3qU6gj0P1kdLygEfvBvexPX+7fBs+fZmbec5d6OTCnVQlrDVydyRMOPX4O1L8GHC+HJyXZ5WBfoOwNO++WxbQsOQmg0hEZ6J1allMc04av6idg2/H4zYfsSKDwER/fApjdh42uMih4Aqwuh8CCEd4Mr3oKkod6OWil1Eq3SpCMit4iIEZFu7tciIo+IyA4RWS8io1vjOMoLohJh9E9h6q9t+/7NG2DabwlwOaHXJDjzPjs+/3PnQ+Y6b0erlDqJFid8EekBnAnsq7X4HKC/+zEfeLylx1E+IjwOpt3O6rEPwcX/glNvhKveteP1P3c+HFzr7QiVUg1ojRr+w8BtgKm1bDbwvLG+BmJFJLkVjqV8UVwfe4NXaAy8cDHk7PR2REqpeogxpvGtGnqzyGzgdGPML0RkDzDWGJMtIu8ADxhjvnRv9wlwuzFmVT37mI89CyA+Pn7Mq6++2ux4fFFRURGRkR3vgmZ95QorOcCo735DVaCD70Y9QEVoVy9F1zwd8bPqiGWCjlmulpRp+vTpq40xYxvbrtGLtiLyMZBUz6o7gN9im3OazRjzFPAUwIABA8y0adNasjufs2zZMjpameAk5Ro5GJ6dxanrbrO9epxlkDAIxsyF/mf59ITrHfGz6ohlgo5ZrvYoU6P/fcaYmfUtF5FhQG9gndhxV1KBNSIyHjgA9Ki1eap7meroUkbZLp1fPwYBQfaC7u7P4eUfQVQKTLgexlxlu34qpdpVs6tbxpgNQEL16zpNOouBn4vIy8ApQL4xJrOlwSo/0XOifVSrctrB21Y+aQdv+/yvMOZKGPUTiB9wbDtjdNA2pdpQW51fvwecC+wASoCr2ug4yh8EBsGgWfZxYA18+TB89RiseARSRkOQA3J3QkkuJA6B7mOg/5n2ocM6KNVqWi3hG2N61XpugAWttW/VgXQfDZf9GwoPw4ZXYdMbgIH+Z9hhmw9tgA3/hVVPQ7cBMOkmGH65T7f9K+Uv9L9IeUdUou3Df+qNJ66rctovguV/h7cWwOa34JJFdqhnpVSz6fmy8j2BQTD8Urj+Czjvr7DjE1h0jh3PXynVbJrwle8SgXHXwI9fhbw98OQUWPO8nZqxvBCW/T94apq9JlBW4O1olfJ52qSjfF+/mXZ+3rd/Ycfp/+ZJO5hbSTbED4KP74EvHrZj+5QehYpC2+d/4gI7FIRSCtAavvIXCYNg3odw8dNQVWFH5rzmE1jwNcxfBv1Oh9zdIAH24u8Xf4G/DYdP77O9f5RSWsNXfkQEhl1iH7WljIJLnz1+2eFN8Nmf4fMH7RnBKdfDxJ/Zu3+V6qS0hq86psQh8IPn4IYVdvauz//srvHfD6V53o5OKa/QGr7q2BKHwA+ed9f4/59N/Cv+D8Ji7dAP4V0hbSL0mkSgU+s/qmPThK86h+rEf2gjfPeCnaDduODoPlj9DHzzOJMkCLJmwtCLYcgcvdlLdTj6F606l6ShcM4Dxy9zlkPGKg58/AQ9Dq2x4/588RCcdT/0m3Fsu/wMWPWM7SI64nI7v68O/aD8iCZ8pYJCodckdvarpMfU52HL2/DRnfDCHIjrC1HJtra/+wvAgCMGNr4GXXrDtIUw/Ac66JvyC1o9Uao2ERh8ASxYCWf90Z4RGBcUHYFTfw43rYVbttnuoWFd4I358OpPoTjb25Er1Sit4StVn6BQe+PWxAbGABx2CQy5yF4AXno/PDYeBp4HfabZbqKOWAiN1usAyqfoX6NSzRUQCKfdbO8E/uwB2PSWHfqhtqhk2/QT1/vYz9Sx0KWXV0JWnZsmfKVaKmkoXPaCHeXz4HeQvQ3KC2x//6P7IW+3HQCu6FCt9wyDgbMgcaidBL5LTwiJ8F4ZVKegCV+p1hIYBD3G2Ud9KoohdxfsWgZb3oFlDwDm2PqQKDts9MDzYNLNOg6QanWa8JVqLyERtmafNMzOA1B61H4B5O6C/P1QlGWfL38Evl1kp4GMTYOQSIiIt81AsWkQ7PB2SZSf0oSvlLeExdoZwLqPPn551lZY9ic7DWTtM4BqoTG29h+TCqnjoMd4O3Dc1ndtk1LqWDs95MBz26UYyn9owlfK1yQMtOMAOcuhvMheDyjKsjd8Hd0HxUegJMfOA7ziEXA53e8bDMMuhv3fwpI7YMkdjIoeCJHX23sF6rtGYAyU5dsvH9XhacJXylcFhdpHRFfbuyftlBO3qSiBzHW27T+uz7HlR/fBxv8RtOKf8M7NdpjoSb+AMXOhNNc2He38FDYvhqN77WiiZ/3R9jxSHZYmfKX8WUg49Jx44vLYNDjtZr6tHMG0Pg47cNxHd9pHtYBge99A2kT45gl7BnHx0xAaadcXZ9svk6LDdhiJqMR2KJBqS5rwlerIROwXwhVvwt6vYNdSiO5uzxiShh9rykkdC+/fBn8dYM8qXFVQdrTWfgKg12Q7htDgC+0XjfI7mvCV6ix6Tqz/bABg/LXQLd2OI1R9obhLL0geYe8a3vI2bPgvvHkDfPAbO6JocLi9nhAcBuPn2xFJPVFRAoWZ0LVva5RKNYEmfKWU1WeqfdQneThM/y3sXQ6rn4XvXrTzCUR0s00/q5+F9HPsZDMhEXaegT7Tj+9CWuWEtS/C0j/am9AGzoIz77NnG6pdaMJXSnlGBHqdZh8XuY4NDV2SCyv/aa8DbHv/2PYRCXDKdfa+g92fw/fv2YvFPU6BUT+Brx+Hx06xXzKOGHu/QVUlVJbYi8dRSXZoiuQRtvtpUKjdb3EO0flbYM0+yD9g35828diIpcXZdtvQKPu6qhL2r4Q9X8K+ryBrsz0jmXxLpxvlVBO+Uqrpas8DEB4H0263CbS8wE4uk73NJvRP77XbBIZC2gSY+XsYdL5NtOOutgPPHdpgt68ottsFO2ySLjoMzjL7/iCHbXLK3w+leYwG+M59/M8esENUpE2AvStsQkegaz+I7WG7qVYU2mWJQ+zyT++1sc78fesm/aIj9kut+xifHDjP9yJSSvmnwCCb/MPjbC+hfjPtTWTFWbaGHhx2/PbRKTD7sYb3Z4y932D/StjzBWRtsTepde3P+swyhk+fY+9A3vg/+PafdiaztAkw7FL7hZG5FvL22nsT+s20F53DYsHlgvduheV/t3c7n/47iEywx3S57JdKdPemJeySXHtPxDdP2jOUiAQYOgdGX+HxtQ1xVcGqRVBZCgPObZOmLk34Sqm2kzAQGNi894rYawQDzz3hruHcZcuO3Xcw5kr7MMaz2npAAJz3V9vks/xv9rrCgHPsWcTOpVCSDcERdkyktFPthe7uY+01i8JMyNlhzyT2rrD3O1SW2LMFV5UdNrv/WbBlsZ0d7ZsnoO/p9v4Hl9M2QQWH2WVxfY7Fm7me0Wt+DUU77esPf2tvpOsxHpJH2i+yhEHN+z3WoglfKdUxNKVpRgTO+D2M/DGseQ7W/cd2Pe03w56NHPke9n1th7jAgASCqar1/kBIGWmvHwSH23sXhl16rDY//FJb61/9jK31v3rFiTHEptkB8yqL4eh+QoOi7LzLySNg63uw/UPY9Ia9IA62G+2IH0JVhT3rKc2Dn7zWpFFWNeErpTqv+HQ7d/GZ99nXdb80So/a5Jqx0t6oFp0MMT3sfQvVF4UbEh5nr2tM/DlkrLKvo1NsM9WOT2wzlavKfmEMuYiVZgynDT7fvnfiz+zDGHtD3PYl9kzkw4XufXezZyLZ2+yEOx7ShK+UUg2dHYTFQvqZ9tFc7jmTazhiYHwfe+9DLc5ly+qPK6637e10ynWQs9POpJa/H/45HQoyO1bCr6ysJCMjg7KyMm+H0iwxMTFs2bLF22G0usbK5XA4SE1NJTg4uB2jUqqDq75ZrXrAvMLMJr3d5xN+RkYGUVFR9OrVC/HDPrOFhYVERTVy6ueHTlYuYww5OTlkZGTQu7feVKNUq4tMsNccmpjwAxrfxLvKysro2rWrXyb7zkpE6Nq1q9+elSnl8wICITKxfRO+iNwjIgdEZK37cW6tdQtFZIeIfC8iZ7XwOC15u/IC/cyUamNRybYNvwlao0nnYWPMX2ovEJHBwOXAECAF+FhE0o2p3a9JKaVUs0UlQ97uJr2lrZp0ZgMvG2PKjTG7gR3A+DY6llJKdT7RyV65aPtzEbkCWAXcYozJA7oDX9faJsO97AQiMh+YDxAfH8+yOl2TYmJiKCwsbIUwvaOqqqpF8d93333ceuut3H777YgIBQUFzJ8/nwkTJrRilE3nSbnKyspO+Dx9WVFRkV/F64mOWCbomOVqapnSskvpU5rH55986PlBjDEnfQAfAxvrecwGEoFA7JnC/cAi93seBX5Sax9PA5c0dqz09HRT1+bNm09Y5g1Tp041u3fvPm7ZvHnzzGOPPXbS9xUUFJywrKSkxEyZMsU4nc6Tvvfw4cPmtttuMw899JBZsmSJMcaYyspKM2vWLGOMMeXl5Wby5MmmsrKywRirTZw48aTHOpm7777bPPjggzWvAwICzLBhw8zw4cPNqFGjzPLly+t9n698dp5aunSpt0NodR2xTMZ0zHI1uUxrXjDm7mhjcnYZYJVpJL8aYxpv0jHGzDTGDK3n8ZYx5rAxpsoY4wL+ybFmmwNAj1q7SXUv6xDeeust+vTpw+LFi7nvvvsYN24cmzZt8vj9ixYtYs6cOQQGnnz+0G+//ZaxY8eyZ88eRowYAUBQUBAxMTG4XC5CQkKYMWMGr7zySqPHXLFihcfxNSYsLIzly5ezbt06/vSnP7Fw4cJW27dSykNRSfZnE5p1WtpLJ7nWy4uwNX+AxcDlIhIqIr2B/sDKlhzLV+zcuZMbbriBDz74gMsuu4w//vGPLFy4kEsuuYSqKs+uSb/44ovMnj275vXu3buZPXs2Y8eOZfz48Xz//ffAsYQ/Y8YM/vvf/wKQnZ1tv6ndw9NeeOGFvPjii40eMzLSzlO6Z88eBg0axLXXXsuQIUM488wzKS0tPWH7+++/n/T0dE477bSaeOpTUFBAly5dPCq3UqoVRafYn+2V8IE/i8gGEVkPTAd+CWCM2QS8CmwGPgAWmA7SQ2fJkiVceOGFpKen1yybM2cOAQEBbN++vdH3V1RUsGvXLnr16gXYO4mvueYaHnroIVatWsU999zDAw88ANgvgt69e7N//37+8Ic/sHTpUmbMmEFeXh75+fkADB06lG+//bZJZdi+fTsLFixg06ZNxMbG8vrrrx+3fvXq1bz88susXbuW995774T9l5aWMmnSJAYOHMg111zDnXfeiVKqnUW569tN6JrZoou2xpifnmTd/dh2/Vbz+7c3sflgQWvuksEp0dx9vodzcTbCGNPoNtnZ2cTGxta8fvPNN9m0aRMXX3wxAE6nk8mTJ/PYY4/Rr18/ABYsWMCCBQsAWLdu3XH7CwwMJCQkpEkXhnv37s3IkSMBGDNmDHv27Dlu/RdffMFFF11EeLidqPqCCy44bn11k05UVBRfffUVV1xxBRs3btS+90q1J0cMBIU1qYbv80Mr+JqZM2dy7733csstt9Qse+utt6isrDyu1t+QsLCw4+5AXbduHffffz9XX311s2MqLy/H4XA0vqFbaGhozfPAwMB6m3Q8NXHiRLKzszly5AgJCQnN3o9SqolEmtw1068SfmvVxFuif//+PProo8yYMYPS0lLeeust4uLieP311wkMDOTIkSPcdttt3Hvvvdx1113ceeed3HTTTTWvn3zySaqqqigrK8PhcJCcnMyHH37IVVddRUBAABs2bGDo0KEe15ZzcnLo1q1bqw5SNmXKFObOncvChQtxOp28/fbbXHfddfVuu3XrVqqqqujatWurHV8p5aGoZCg85PHmfpXwfcWcOXO46KKLmDdvHqNHj+bGG2+sWRcfH09aWhq33HILTz/9NMaY414HBwdz5pln8uWXXzJz5kzmzZvH0qVLGTRoEGFhYQwdOpQXXnjB41iWLl3Keeed16rlGz16NJdddhkjRowgISGBcePGHbe+ug0/ICAAYwzPPfdcoz2OlFJtICoZDqz2eHNN+M0kIjzzzDMnLC8qKmLXrl0EBQURGRlJZmbmca/Btsk//PDDzJw5k7CwMF577bVmx/HSSy/VXOQ9maKiIgB69erFxo0ba5bfeuut9W5/xx13cMcdd9S7rvqmq444CqhSfiUqqV176ahanE4nN910E/fddx8jR47k448/5rbbbqt5XX0X3ejRo5k+fbrH3TgbUlFRcUKPIaVUJxKdAk7PR6XVGr6H5s6de1zvmvoEBQWxaNEiAH79618DcMoppxAVFVXzutq8efNaHFNISAhXXHFsrkxPYlRKdSDVN195SBO+h+bOnevtEBrlDzEqpVpRVEqTNtcmHaWU8ldNrOFrwldKKX8Vldz4NrVowldKKX8V7IAwz8ey0oSvlFL+rAnt+JrwlVLKn8V73i1bE75SSvmzS5/1eFNN+Eop1UlowldKqU5CE76Pu+uuuygrK+O6667j+uuv50c/+hHLly/3dlhKKT+kCd9D06ZNO2GikKuvvpp//OMfTd5XaWkpU6dObXQsnaysLMrLy3n88ce55JJLeOKJJ3j++edrBkurqKhgypQpOJ3OBmOsduqppzY5zmr33HMPf/nLX2peBwYGMmnSJEaMGMHo0aNbdb5cpVTb0YTfDDqJuU5irpQ/0oTfRDqJ+fF0EnOl/Icm/CbSScx1EnOl/JV/jZb5/m/g0IbW3WfSMDin8QlEPKGTmOsk5kr5Mv9K+D5AJzE/nk5irpT/8K+E30o18ZZobBLzZ555hqSkJM4++2yuvvpqhg0bxsCBA2teP/bYYzqJuVLKK/wr4fuIk01iPmXKFBYtWsSBAwe47LLLSEpK4tVXX615HRYWppOYK6W8QhN+MzU0iXnfvn1Zs2YN+fn5XHPNNRQWFh73GnQSc6WUd2gvnTYQEhLCXXfd1eBrncRcKeUNmvA95MkE4fn5+dx0001ceeWVJCQkkJ+fz69//eua17XNmzevxc0gOom5UqoptEnHQ55MEB4TE8Mjjzxy3OsHH3yw3Zo+dBJzpdTJaA1fKaU6CU34SinVSWjCV0qpTsIvEr4nQxYo36KfmVK+x+cTvsPhICcnRxOIHzHGkJOT06ThHpRSbc/ne+mkpqaSkZHBkSNHvB1Ks1QPodDRNFYuh8NBampqO0aklGqMzyf84OBgevfu7e0wmm3ZsmWMGjXK22G0uo5aLqU6shY36YjIjSKyVUQ2icifay1fKCI7ROR7ETmrpcdRSinVMi2q4YvIdGA2MMIYUy4iCe7lg4HLgSFACvCxiKQbY1o2loBSSqlma2kN/wbgAWNMOYAxJsu9fDbwsjGm3BizG9gBjG/hsZRSSrVAS9vw04HJInI/UAbcaoz5FugOfF1ruwz3shOIyHxgvvtluYhsrG87P9YNyPZ2EG2gI5ZLy+Q/OmK5WlKmnp5s1GjCF5GPgaR6Vt3hfn8cMAEYB7wqIn2aECTGmKeAp9zHWmWMGduU9/u6jlgm6Jjl0jL5j45YrvYoU6MJ3xgzs6F1InID8D9jO8mvFBEX9lvqANCj1qap7mVKKaW8pKVt+G8C0wFEJB0IwZ6SLAYuF5FQEekN9AdWtvBYSimlWqClbfiLgEXudvcK4Ep3bX+TiLwKbAacwAIPe+g81cJ4fFFHLBN0zHJpmfxHRyxXm5dJdMgCpZTqHHx+LB2llFKtQxO+Ukp1El5J+CJytnvIhR0i8pt61oeKyCvu9d+ISK/2j7JpPCjTXBE5IiJr3Y9rvBFnU4jIIhHJaujeCLEecZd5vYiMbu8Ym8qDMk0Tkfxan9Nd9W3nS0Skh4gsFZHN7iFOflHPNn71WXlYJn/8rBwislJE1rnL9ft6tmm7/GeMadcHEAjsBPpge/WsAwbX2eZnwBPu55cDr7R3nG1QprnAo96OtYnlmgKMBjY2sP5c4H1AsPdifOPtmFuhTNOAd7wdZxPLlAyMdj+PArbV8/fnV5+Vh2Xyx89KgEj382DgG2BCnW3aLP95o4Y/HthhjNlljKkAXsYOxVDbbOA59/PXgBkiIu0YY1N5Uia/Y4z5HMg9ySazgeeN9TUQKyLJ7RNd83hQJr9jjMk0xqxxPy8EtnDine1+9Vl5WCa/4/79F7lfBrsfdXvOtFn+80bC7w7sr/W6vmEXarYxxjiBfKBru0TXPJ6UCeBi9+n0ayLSo571/sbTcvubie5T7vdFZIi3g2kK9+n/KGzNsTa//axOUibww89KRAJFZC2QBXxkjGnws2rt/KcXbdvP20AvY8xw4COOfYMr37IG6GmMGQH8H/bmQr8gIpHA68DNxpgCb8fTGhopk19+VsaYKmPMSOwIBONFZGh7HdsbCd+TYRdqthGRICAGyGmX6Jqn0TIZY3KMe1RR4F/AmHaKrS11uCE0jDEF1afcxpj3gGAR6eblsBolIsHYxPiiMeZ/9Wzid59VY2Xy18+qmjHmKLAUOLvOqjbLf95I+N8C/UWkt4iEYC9KLK6zzWLgSvfzS4BPjfsKho9qtEx12ksvwLZJ+rvFwBXuHiATgHxjTKa3g2oJEUmqbi8VkfHY/xFfrmzgjvdpYIsx5qEGNvOrz8qTMvnpZxUvIrHu52HAGcDWOpu1Wf5r9ykOjTFOEfk58CG2d8siY8wmEfkDsMoYsxj7Qf9bRHZgL7Bd3t5xNoWHZbpJRC7ADjWRi+2149NE5D/YnhDdRCQDuBt7kQljzBPAe9jeHzuAEuAq70TqOQ/KdAlwg4g4gVLgch+vbABMAn4KbHC3DQP8FkgDv/2sPCmTP35WycBzIhKI/YJ61RjzTnvlPx1aQSmlOgm9aKuUUp2EJnyllOokNOErpVQnoQlfKaU6CU34SinVSWjCV0qpTkITvlJKdRKa8JU6CRHpJSLviZ3rYJuILPR2TEo1lyZ8pRogIgHYsVyeMMYMAIYBY0VkvncjU6p59E5bpRogIucA1xhjLq61LBn4zBiT7r3IlGoereEr1bBB2NnLargHHIt2D5KnlF/RhK9Uw6qAyNoL3KMzhmMHwVPKr2jCV6phy4Bz60wvdwawxhjj8k5ISjWfJnylGmCMWQd8B/wBQEQSgYeww/Qq5Xc04SvVABH5DTAW+J2InA48DvQE/uGeZ1Upv6K9dJRSqpPQGr5SSnUSmvCVUqqT0ISvlFKdhCZ8pZTqJDThK6VUJ6EJXymlOglN+Eop1Un8f0Dja4S2rBJtAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.signal as sig\n", "\n", "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) - .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(Om, 20*np.log10(np.abs(.5*Pxx)), label=r'$| \\Phi_{xx}(e^{j \\Omega}) |$ in dB')\n", "plt.plot(Om, 20*np.log10(np.abs(.5*Pyy)), label=r'$| \\Phi_{yy}(e^{j \\Omega}) |$ in dB')\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", "[./uniform_white_noise.wav](./uniform_white_noise.wav)\n", "\n", "**Pink noise**\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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEYCAYAAACgDKohAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3X+UHFWd9/H3d3o6MyMEIiSgSQgJaxYNBBJIohxkn2xEiOITQFYh4i6oiB5E/PVEwXVZlkWCC0d0BQX8hagoiBKjyxrBoCC/ZCABQiAQQpQkQkIgkUD3ZHrm+/xR1TOdSXdP/6iemun6vM7JSXdVddW93TX16XvrVpe5OyIiIpVoibsAIiIycig0RESkYgoNERGpmEJDREQqptAQEZGKKTRERKRiCg1pODM7xszWxF2OYsxsrpltKDN/h5kdNJRlCrd7vZldMtTbFRmMQkNKMrP1ZpYJD5z5f1dV8Do3szfln7v73e5+cIPK2NCDq7vv6e7rGrV+kZGmNe4CyLD3f939jrgLIdEzs5S798RdDhlZ1NKQmpjZm8zsD2a23cxeNLObwul3hYs8ErZMTh3YBRS2YBaZ2aNm9qqZfdfM9jez/zWzV8zsDjN7fcHyPzOz58Nt3WVmh4TTzwZOBz4fbutX4fTxZvZzM9tiZs+a2XkF6+oIWycvm9lqYPYg9exrNYWv+2ZYzh1mdo+ZvcHMvhau70kzmzmgnheY2epw/vfNrD2cd6aZ/bHUtgZMH2tmvzazbWb2kpndbWYtg9W1yHquN7NvmdltZvYq8I9m1mZmV5jZX8zsBTO7xsw6KthuybqF8z9qZmvD1y01s/ED6vlxM3s6XPfVZmbhvKL7VTjvzWZ2e7jONWb2/nKfnTSIu+uf/hX9B6wHji0x7yfAvxJ88WgH3l4wz4E3FTyfC2wYsN77gf2BCcBm4GFgZriu5cC/Fyz/YWA00AZ8DVhZMO964JKC5y3AQ8CFwCjgIGAdcHw4/zLgbmAf4ABgVWHZitSzry7htl4Ejiwo57PAvwAp4BLgzgH1XBVuZx/gnnxZgTOBPw6yrfyyi4FrgHT47xjABqtrkbpcD2wHji743K4EloblGw38ClhcbrsV1G1e+D4dEX5m3wDuGlDPXwNjgEnAFmB+uf0K2AN4DvgQQQ/JzHAb0+L+O0naP7U0ZDBLwm+D+X8fDad3AwcC49096+5/LLOOYr7h7i+4+0aCg/gD7r7C3bPArQQHBQDc/Xvu/oq7dwEXAYeb2d4l1jsbGOfuF7v7Tg/OR3wbOC2c/37gy+7+krs/B/x3leW+1d0fKihn1t1v8KCb56bCcoeucvfn3P0l4MvAwiq3B8F7/UbgQHfv9uAckVdQ12J+6e73uHsv0AWcDXwmfD9eAS4teH2p7Q5Wt9OB77n7w+FndgFwlJlNLnjtZe6+zd3/AtwJzCjYZrH96j3Aenf/vrvn3H0F8HPgfZW+iRINhYYM5iR3H1Pw79vh9M8TfNv9k5k9bmYfrnK9LxQ8zhR5vicE/e5mdpmZPWNmfyP4hgswtsR6DwTGFwYd8EWCVg3AeIJvrHl/bkS5Cwzc1niqdzmwFvitma0zs/PD6YPVtZjC8owDXgc8VPD634TTy213sLqNp+B9dfcdwFaCVmXe8wWPX6P/fSu1Xx0IvHVAXU8H3lCmrtIAOhEuNXH354GPApjZ24E7zOwud18b8aY+AJwIHEsQGHsDLxMcWCDo6ij0HPCsu08tsb6/EnSpPB4+nxRlYYs4oODxJGBT+PhVggM2AGZW8uAXtgA+B3zOzA4FlpvZgwxe16KrK3j8IkHQHRK2+Crarrv/bpC6bSI4yOfrtgewL7DbNopss+h+RVDXP7j7OyuqpTSMWhpSEzN7n5lNDJ++THAw6g2fv0DQvx6F0QTdKFsJDrKXDpg/cFt/Al4xsy+EJ71TZnaomeVPeN8MXGBmrw/L/8mIylnKJ8xsopntQ9BXnz+x+whwiJnNCE8gX1RqBWb2nvAEsRGck+gheK8Hq2tZYRfVt4ErzWy/cFsTzOz4QbY7WN1+AnworFsbwWf2gLuvH6xMZfarXwN/b2b/bGbp8N9sM3tLJXWV6Cg0ZDC/sl2v07g1nD4beMDMdhCcSP2U91/PcBHwg7Abod4RLjcQdHVsBFYTnEAv9F1gWritJeG5hfcQ9JE/S/Bt+jsELRSA/wjX9yzwW+CHdZZvMDeG21kHPENwshx3fwq4GLgDeBood05oarjcDuA+4JvufmcFda3EFwi6oO4Pu//uAPLX1BTdbgV1uwP4N4JzDn8F/o7y51kKFd2vwlbPceF6NhF0b32F4ES7DKH8SAgRiZiZrQfO8ia8zqWZ6yblqaUhIiIVU2iIiEjF1D0lIiIVU0tDREQq1nTXaYwdO9YnT54cdzFEREaUhx566EV3HzfYck0XGpMnT6azszPuYoiIjChmVtGvI6h7SkREKqbQEBGRiik0RESkYk13TkNEmkd3dzcbNmwgm83GXZSm0d7ezsSJE0mn0zW9XqEhIsPWhg0bGD16NJMnTya8uZ/Uwd3ZunUrGzZsYMqUKTWtQ6EhMkItWbGRy5etYdO2DOPHdLDo+IM5aeaEwV84gmSzWQVGhMyMfffdly1bttS8DoWGJFKjD7hDsf4LfvEYme4eADZuy3DBLx4DaLrgUGBEq973U6EhidPoA26t668maC5ftqZv/XmZ7h4uX7am6UJDhheFhiROLQfcKA/oxdYFVBU0m7Zlim671HSpzbZt27jxxhs555xz4i7KsBHbkFsz+56ZbTazVSXmm5n9t5mtNbNHzeyIoS6jNKdqD7j5lsPGbRmc/gP6khXF715abv2l1vUfv3q8ZNAUM35MR1XTpTbbtm3jm9/85m7Tc7lcDKUZHuK8TuN6YH6Z+e8iuHPYVOBs4FtDUCZJgGoPuOVaDtWuv9S6Xn6tu+hrSgXQouMPpiOd2mVaRzrV12qRaJx//vk888wzzJgxg9mzZ3PMMcewYMECpk2bxvr16zn00EP7lr3iiiu46KKLAHjmmWeYP38+Rx55JMcccwxPPvlkTDWIXmzdU+5+l5lNLrPIicANHvx2+/1mNsbM3ujufx2SAkrTWnT8wXz+F/fybMun6bUdABjwkqfYazFkdvbQ69Bi0DEqxatdPdC++3qey8Jei1O7LQ/wWnsPhTcdyK+/1LpKaTEY+1/Fx9PvfF1vsD76t33W7S2cdXswf3TbaO468y4O2PuAyjc4jH36N59m5fMrI13njDfM4Gvzv1Zy/mWXXcaqVatYuXIlv//97znhhBNYtWoVU6ZMYf369SVfd/bZZ3PNNdcwdepUHnjgAc455xyWL18eadnjMpzPaUwAnit4viGctltomNnZBK0RJk2aNCSFk5HrpJkTeHb73nz2D5vo6JnF3ukJzJw0BoD71r1ER29v37Kp3hb2aYGuXO9u6xnV2kJPlt2WP+qgfQC495mt9Lqzx6hWZk4aw0Hj9uTnD23g1Z27d22Mam2hpxd6CtfV0sJRU/bhoHF7lqzLssefB+D4Q96wy/RNr2zi1idv5emXnm6a0BgO5syZM+j1DTt27ODee+/lfe97X9+0rq6uRhdtyAzn0KiYu18HXAcwa9Ys3VVKBvW2vxsNf4BZYz/AXed9BoCjL1vO3l27dweN6UjTRe8u3Uod6RTt1sLLXbt3K23/awf3nD+PU6+9D4CbPnZU37xj37jryKr8uhYvmA7A5295lJ09vUwoONle7iT8qc8F27jq3f3bWLJiIxeu/CVwK5+48X4Wv7s5rt8o1yIYKnvssUff49bWVnoLQj5/1Xpvby9jxoxh5cpoW0XDxXD+7amNQOFXpInhNJG6ZXPBH3iKUX3TSp0/2J7pZvF7pzMqFfy5TBjTweL3TmdblechIGjlFFvXSTMncNLMoMXz1in7cM/58/oCo5qT8PnlX3wlOJi9+OqOsstLeaNHj+aVV14pOm///fdn8+bNbN26la6uLn79618DsNdeezFlyhR+9rOfAcFV2I888siQlbnRhnNoLAX+JRxF9TZgu85nSFT6QsPa+qaVO4Fd7IBe6wimYusqpdqT8PnlLQxDt51ll5fy9t13X44++mgOPfRQFi1atMu8dDrNhRdeyJw5c3jnO9/Jm9/85r55P/7xj/nud7/L4YcfziGHHMIvf/nLoS56w8TWPWVmPwHmAmPNbAPw70AawN2vAW4D3g2sBV4DPhRPSaUZZXJBayBl/S2NRccfXLTrqNSIpGqXr0W1w4Pz0/tCg51ll5fB3XjjjSXnnXfeeZx33nm7TZ8yZQq/+c1vGlms2MQ5emrhIPMd+MQQFUcSJt/SaCkIjfw3/mLnFYqpdvlajB/TwcYiB/xyrZyN2zKY7xoaun5DojKcu6dEGqZY9xRU13VUy/LVqvZ6jPzy/d1T3bp+QyLVFKOnRKqV6Q67pwpOhA9HtbZ+Ft2yAoDR7b0sPmF6U4yekuFBoSGJVKqlMRydNHMCP/nTX4Bdh+8Otvwzz6f456PGKzAkUuqekkQaSaFRq5S19dVTJCoKDUmkTC6D0YKRGnzhESrFqL5RYiJRUWhIImVzWVpsVFPf4KfFRiWupbFkxUaOvmw5U87/H46+bPmQXtR4/fXXs2nTpr7nZ511FqtXr657vevXry877LeUM888k1tuuaXu7Q+k0JBEyuaypGjerilIXvdUtVfPR21gaHznO99h2rRpda+31tBoFIWGJFKmO7PLhX3NKGXJ6p6q9ur5Sv3oRz9izpw5zJgxg4997GP09PRw5plncuihhzJ9+nSuvPJKbrnlFjo7Ozn99NOZMWMGmUyGuXPn0tnZCcCee+7JokWLOOSQQzj22GP505/+xNy5cznooINYunQpEITDMcccwxFHHMERRxzBvffeCwQ/z3733XczY8YMrrzySnp6eli0aBGzZ8/msMMO49prrwWCnys599xzOfjggzn22GPZvHlzXfUuRaOnJJGyPdlEhEaSWhqNuJvhE088wU033cQ999xDOp3mnHPO4ZJLLmHjxo2sWhXcP27btm2MGTOGq666iiuuuIJZs2bttp5XX32VefPmcfnll3PyySfzpS99idtvv53Vq1dzxhlnsGDBAvbbbz9uv/122tvbefrpp1m4cCGdnZ1cdtllXHHFFX2/bXXdddex99578+CDD9LV1cXRRx/Ncccdx4oVK1izZg2rV6/mhRdeYNq0aXz4wx+uue6lKDQkkfLnNJpZC8kKjWqvnq/E7373Ox566CFmz54NQCaTYf78+axbt45PfvKTnHDCCRx33HGDrmfUqFHMnx/cc2769Om0tbWRTqeZPn163305uru7Offcc1m5ciWpVIqnnnqq6Lp++9vf8uijj/adr9i+fTtPP/00d911FwsXLiSVSjF+/HjmzZtXc73LUfeUJFKmO5OIcxr5ixiToBF3M3R3zjjjDFauXMnKlStZs2YNX//613nkkUeYO3cu11xzDWedddag60mn032DLlpaWmhra+t7nL917JVXXsn+++/PI488QmdnJzt37ixZpm984xt9ZXr22WcrCq6oqKUhiZTNZZv6Gg3Id08V/1nvkej5Sy+l64nSt02dCfxoRxfPbHkVd6etNcUB+3Qwdn0bfy7xmra3vJk3fPGLJdf5jne8gxNPPJHPfOYz7Lfffrz00ku88sorvP71r+eUU07h4IMP5oMf/CBQ/mfUK7F9+3YmTpxIS0sLP/jBD+jp6Sm63uOPP55vfetbzJs3j3Q6zVNPPcWECRP4h3/4B6699lrOOOMMNm/ezJ133skHPvCBmstTikJDEikIjebungpGT22JuxhDauyebezoCr65T953j0GWHty0adO45JJLOO644+jt7SWdTvPVr36Vk08+ue8GTIsXLwaCIa4f//jH6ejo4L777qt6W+eccw6nnHIKN9xwA/Pnz++74dNhhx1GKpXi8MMP58wzz+RTn/oU69ev54gjjsDdGTduHEuWLOHkk09m+fLlTJs2jUmTJnHUUYP/ekAtFBqSSJlchpSVvo1qM0jR1lSjp8q1CAodGPF2Tz31VE499dRdpj388MO7LXfKKadwyimn9D3//e9/3/d4x44dfY8vuuiiXV6Xnzd16lQeffTRvulf+cpXgKBra+D9xS+99FIuvfTS3cpw1VVXDVKb+umchiRSIk6EJ2z0lAwNhYYkUnBxX3OHRtKG3MrQUGhIIiXjRHhzXBEe3I9NolLv+6nQkERKyhXhO3t20tPbM/jCw1R7eztbt25VcETE3dm6dSvt7e01r0MnwiWRsrksLa3NHRotYfdbV08Xr2t5Xcylqc3EiRPZsGEDW7YkaxRYI7W3tzNx4sSaX6/QkMRxd7p6uhJxcR8EAfm69MgMjXQ6zZQpU+IuhhRQ95QkTv8NmJq7pZGvX5KuCpfGU2hI4iThrn2wa0tDJCoKDUkchYZI7RQakjj5q6SbvnsqPBHeTFeFS/wUGpI4+W/eLU1+cV/+ine1NCRKCg1JHHVPidROoSGJkx9N1PTdUxo9JQ2g0JDEScyQW9Q9JdFTaEji9J/TaO7uqRZ1T0kDxBoaZjbfzNaY2VozO7/I/ElmdqeZrTCzR83s3XGUU5pL/+ip5g6NfP00ekqiFFtomFkKuBp4FzANWGhm0wYs9iXgZnefCZwGfHNoSynNKDHdUxo9JQ0QZ0tjDrDW3de5+07gp8CJA5ZxYK/w8d7ApiEsnzSpxISGzmlIA8QZGhOA5wqebwinFboI+KCZbQBuAz5ZbEVmdraZdZpZp34NUwbTP3qqubunWjR6ShpguJ8IXwhc7+4TgXcDPzSz3crs7te5+yx3nzVu3LghL6SMLMm5uC9FuiWtloZEKs7Q2AgcUPB8Yjit0EeAmwHc/T6gHRg7JKWTppWU7imA9tZ2hYZEKs7QeBCYamZTzGwUwYnupQOW+QvwDgAzewtBaKj/SeqSzWVpbWmlxZr/djId6Q6FhkQqttBw9xxwLrAMeIJglNTjZnaxmS0IF/sc8FEzewT4CXCm676PUqdMLkN7a+23uxxJ2lvbNeRWIhXrVy13v43gBHfhtAsLHq8Gjh7qcklzy+ayiQoNtTQkSsP9RLhI5LK5LB2tHXEXY0h0tKp7SqKl0JDEUfeUSO0UGpI46p4SqZ1CQxInm8vSkU5I95RGT0nEFBqSOJnuhHVP6YpwiZBCQxJH3VMitVNoSOJo9JRI7RQakjgaPSVSO4WGJI66p0Rqp9CQxFH3lEjtFBqSOEkbPZXrzZHrzcVdFGkSCg1JnKR1T4Hu3ifRUWhIovT09tDd263QEKmRQkMSJX/wTNIV4aBbvkp0FBqSKPnQUEtDpDYKDUkUhYZIfRQakij5C92SNOQW0AV+EhmFhiSKWhoi9VFoSKIoNETqo9CQREnq6CmFhkRFoSGJkh96mrSWhobcSlQUGpIo6p4SqY9CQxKlr3sqYaOnFBoSFYWGJEp+6GnSWhoacitRUWhIoqh7SqQ+Cg1JlKSNnlJoSNQUGpIoSRs9ZWa0pdo0ekoio9CQREla9xTolq8SLYWGJEo2l2VUahQtlpxdvyOtW75KdGL9yzGz+Wa2xszWmtn5JZZ5v5mtNrPHzezGoS6jNJdMLjm3es1rb23X6CmJTGtcGzazFHA18E5gA/CgmS1199UFy0wFLgCOdveXzWy/eEorzSJJt3rNU/eURCnOlsYcYK27r3P3ncBPgRMHLPNR4Gp3fxnA3TcPcRmlyWRz2cRc2JfX0aruKYlOnKExAXiu4PmGcFqhvwf+3szuMbP7zWx+sRWZ2dlm1mlmnVu2bGlQcaUZqHtKpD7D/WxgKzAVmAssBL5tZmMGLuTu17n7LHefNW7cuCEuoowk6p4SqU+cobEROKDg+cRwWqENwFJ373b3Z4GnCEJEpCYKDZH6xBkaDwJTzWyKmY0CTgOWDlhmCUErAzMbS9BdtW4oCynNJdOdSczV4Hkd6Q5d3CeRiS003D0HnAssA54Abnb3x83sYjNbEC62DNhqZquBO4FF7r41nhJLM1BLQ6Q+sQ25BXD324DbBky7sOCxA58N/4nULZGhkVJoSHSG+4lwkUhlcpnkDblNd2j0lERGoSGJksiWhrqnJEIKDUmUJIdG0NsrUh+FhiRKUq8I7/Vecr25uIsiTaCqE+Fmtk8Fi/W6+7YayyPSUJnuZF4RDsH5nHQqHXNpZKSrdvTUpvCflVkmBUyquUQiDZLrzdHjPYkNjWwuy15te8VcGhnpqg2NJ9x9ZrkFzGxFHeURaZik3eo1L19fnQyXKFR7TuOoiJYRGXJJu9VrXl/3lK4KlwhUFRrungUws0sGzgvvj9G3jMhwk8RbvcKu3VMi9ap19NQEM1uYfxLeHOmOaIok0hh93VMJHD0FCg2JRq0/I/IxYJmZPQM48H3gC5GVSqQB8ldFJ7WloavCJQrVDrm9AXgYWAF8ArgRyAEnufva6IsnEh11T6mlIfWrtnvqeoLhth8CfgRMBl4GPmhm/xRpyUQiptFTCg2pX1UtDXdfDizPPzezVuAtwOHAW4FbIi2dSIQ0ekrdU1K/un4aPbwnxmPhvx9FUiKRBlH3lFoaUr+quqfM7OEolhGJg0ZPKTSkftW2NN5iZo+WmW/A3nWUR6RhNHpK3VNSv2pD480VLNNTS0FEGk3dU2ppSP2qPRH+ZwAzmwecDmwDVgGPAqvcvSvyEopEJKmhMSo1CsMUGhKJsuc0zOwQM/txkVnfA34F3A8cBFwIPB598USikx89lLQht2ZGe2u7Rk9JJAZradxB8R8g/LO7Lwkf/yzaIok0Rv6bdluqLeaSDD3d8lWiMtjoqeOAL+efmNkNZvZp4H4z+2xDSyYSsWwuS1uqDbNyt4NpTgoNiUrZ0HD3x9z99IJJ1xOMkNof+Gcz+7OZLTWz/zSz9zWwnCJ1y+ayieuayutId5DtUWhI/aK8Inw26qqSYSyTS96tXvN0TkOiEuUV4SLDWjaXTXRoqHtKolDr/TRERpxsLpu4q8HzOlo7FBoSCYWGJEbiu6d0RbhEQKEhiaHuKbU0pH4KDUmMxI+eUmhIBGINDTObb2ZrzGytmZ1fZrlTzMzNbNZQlk+aS6Y74d1TGj0lEYgtNMwsBVwNvAuYBiw0s2lFlhsNfAp4YGhLKM0m0d1TKXVPSTTibGnMAda6+zp33wn8FDixyHL/CXwF0B4vdUn06Cl1T0lE4gyNCcBzBc83hNP6mNkRwAHu/j9DWTBpTho9pe4pqd+wPRFuZi3AV4HPVbDs2WbWaWadW7ZsaXzhZERKdPdUOHrK3eMuioxwcYbGRuCAgucTw2l5o4FDgd+b2XrgbcDSYifD3f06d5/l7rPGjRvXwCLLSJbo7qmw3jt7dsZcEhnp4gyNB4GpZjbFzEYBpwFL8zPdfbu7j3X3ye4+meDeHQvcvTOe4spI5u6JHz0FuuWr1C+20Ah/t+pcYBnwBHCzuz9uZheb2YK4yiXNqbu3G8cTHxo6GS71qusHC+vl7rcBtw2YdmGJZecORZmkOSX1Vq95Cg2JyrA9ES4SpaTe6jUvX29d4Cf1UmhIIqiloZaGREOhIYmg0FBoSDQUGpII+VFDSR9yq9FTUi+FhiSCWhpqaUg0FBqSCAoNhYZEQ6EhiZA/WCZ99JRCQ+ql0JBEyA81TXpLQ0NupV4KDUkEdU+pe0qiodCQROjrnkr46CmFhtRLoSGJkB9qmvSWhobcSr0UGpIISe+eam1ppcVa1NKQuik0JBGSPnrKzOho1S1fpX4KDUmETHcGw0i3pOMuSmzaW9s1ekrqptCQRMjf6tXM4i5KbPK3fBWph0JDEiGbyya2ayqvI91BtkehIfVRaEgiZHLJvdVrnrqnJAoKDUmEfPdUkql7SqKg0JBEyOayib2wL0+jpyQKCg1JBHVPhd1TurhP6qTQkERQ95S6pyQaCg1JBIWGQkOiodCQRMh0ZzTkNt2h0VNSN4WGJIJaGtCeUktD6qfQkERQaKh7SqKh0JBEyOQyGnKb7tDoKambQkMSQS2NoKWxs2cnvd4bd1FkBFNoSCIoNPrvJdKV64q5JDKSKTSk6bm7rghHt3yVaCg0pOl19QTfrNXS0C1fpX6xhoaZzTezNWa21szOLzL/s2a22sweNbPfmdmBcZRTRrak3+o1L19/tTSkHrGFhpmlgKuBdwHTgIVmNm3AYiuAWe5+GHAL8F9DW0ppBkm/1Wtevv4KDalHnC2NOcBad1/n7juBnwInFi7g7ne6+2vh0/uBiUNcRmkC+aug1dIIu6d0VbjUIc7QmAA8V/B8QzitlI8A/1tshpmdbWadZta5ZcuWCIsozUDdUwF1T0kURsSJcDP7IDALuLzYfHe/zt1nufuscePGDW3hZNjr657S6ClAoSH1aY1x2xuBAwqeTwyn7cLMjgX+Ffg/7q4B5lK1/GghtTQ0ekrqF2dL40FgqplNMbNRwGnA0sIFzGwmcC2wwN03x1BGaQLqngqoe0qiEFtouHsOOBdYBjwB3Ozuj5vZxWa2IFzscmBP4GdmttLMlpZYnUhJGj0V0OgpiUKc3VO4+23AbQOmXVjw+NghL5Q0HY2eCmj0lERhRJwIF6mHuqcC6p6SKCg0pOlp9FRAo6ckCgoNaXoaPRVoa20DNHpK6qPQkKan7qlAa0srrS2tamlIXRQa0vQUGv10y1epl0JDml6mO0PKUqRT6biLEruO1g6NnpK6KDSk6emuff3aW9vJ9qilIbVTaEjTU2j0U/eU1EuhIU0vk8sk/mrwvI60uqekPgoNaXpqafRTS0PqpdCQpqfQ6KfQkHopNKTpZXPZxF8NntfR2qHQkLooNKTpZXIZtTRC7a3tuiJc6qLQkKan7ql+6p6Seik0pOllc1mNngp1pNU9JfVRaEjTy3SreyqvPdWuIbdSF4WGND11T/VT95TUS6EhTU+jp/qpe0rqpdCQpqfRU/3aW9vp7u2mp7cn7qLICKXQkKan7ql+uuWr1EuhIU2t13vZ2bNT3VMh3fJV6qXQkKamGzDtKv8+6AI/qZVCQ5qaQmNX6p6Seik0pKnlD466uC+Qfx8UGlIrhYY0tfyFbGppBPoNNgATAAAK0klEQVS6p3SBn9RIoSFNTd1Tu1L3lNSrNe4CSPWWrNjI5cvWsGlbhvFjOlh0/MGcNHNC3MUalhQau1JoSL0UGsNYsXAAuOAXj5HpDi7O2rgtwwW/eAyg6uAYaeFTbXmXrNjIv/3vHwG44OdP0vuuI4d1/YZCfsitRk8FhuPfwHAsUyGFxiBqOVBF8YEvWbGxaDi0p1v6puVlunu4fNkagKLbjjJ8StWv0Tt6qfejVHnzy7/UswPaYOsOrzlcm0mplkZc+3k5jd7Xyu1TUPnfUpz7eRxiDQ0zmw98HUgB33H3ywbMbwNuAI4EtgKnuvv6oSpftTsVlD4QF1u+3E5w+bI1RcNh4LS8/LYGbrvzzy/x84c2VhU+pcpV6v0otQ2oLYCqeT9KhWV+eW/ZCYD5qEHrlwTFQiPK/bzcgbWa6aW2Ueu+VkypfeqipY/TleutettRhFy5/XwoQqsSsYWGmaWAq4F3AhuAB81sqbuvLljsI8DL7v4mMzsN+ApwatRlKfVBVLtTlToQl1oeSu/om7ZV132QMiu67Z888Bw97rtNLxU+5bZb6v0otY3BWj/VHKhKlatUWOafO2FoMGrQ+iVBsSG3Ue3n+c+7moN9tV9qyu1r1R64S+0L2zLdu02rZD+vNuTy730lZdq0LTNsWiHmA96EIduw2VHARe5+fPj8AgB3X1ywzLJwmfvMrBV4HhjnZQo9a9Ys7+zsrLgc+Q9i9jNX8uB+DwbbxYrutFFrMWPPtlay4Xba06m+eTu6cvQWqaZhADi+y7TC51GUqbunl2x3L47TYkZba23vx8Cy5d/brlxvVfUbOK3wFZSZ/rdRO9m0x6t87e7/wxsye9DWmmLmpDEArP7r3wCY9sa9dnllo6fHue1tLV3MOOhmxmf3YKLvAcDfsrsfJGvVYlb0cx3sc4pCRzrVt8/2r91Ip4zuHq9in6petfWudj9vsWBesW3k/z6z3b1M33o4Gw75Uk0tEDN7yN1nDbZcnN1TE4DnCp5vAN5aahl3z5nZdmBf4MXChczsbOBsgEmTJlVViPy3LMNIef8I5O5up5UWGpqpDr056MkFT1/r6aG9NUU6Zbwu1Vr0IN0RBkt+nhm0t6bI5nqqKmupP9WO1hS9OdjZ7bRgwZIOO3eG70d1NSzYWr/ubsfcSA2YXu41QXlLLV96+uu72jnwb3sxNttBS4txwD79F/m9blSq6KsaPT3ObXf/zZm7cSJbOjK8ZsH+FtV+bgbeS5nPtdrp1W27u9tpcdttfT05+vflAVsttk9VG2O113v3eaXK1NGaItPdU3wb4d9nSziv0S2QpjgR7u7XAddB0NKo5rX55uAfJn1yl+kGXHnqjF2agxActNvTLbz82u7fzsZ0pHdpzley/I5cL5nW/uVz6RSL3zu96v7LgU3X/LZPOXLCLs3j/PTF750OFO86Ovqy5Wwp0kwuVb9S2yjVMjFg/JgONlbRVZT/PIqduyi2ngkF8zdty/DN9wTPjyx4/w4ssa1GT49r2/37yBl903JlPr9q9/PF751e8vNIme3WtVNuerX72uL3TuczN62s7mBP8X0KqPpvqdp6V1um/N/lYNt4INy9G3n+Ls7Q2AgcUPB8Yjit2DIbwu6pvQlOiEem1MFr/JiOvje80p3qogWHVLW8GUX7bfMfdv5fJUqV9aSZE5h14D4lw6fY+kv1q27PdJfcoYtto9QfUn5+NYGc/zyKlbfYevLlSvJJ74FKnbu488ktfQe+evbz/HtdzQG31PRy2yi1P1d74C63T1W77WrrXct+Xupvppbzk/WIMzQeBKaa2RSCcDgN+MCAZZYCZwD3Af8ELC93PqMWpT6I/B9NtTtV/jWVLP+Zm1YWXa7WD7tUWas9eA4WpNVso9wBHSoP2Py8Ytstth6Fxe7KnWSNaj+v5ctLtV9qqj2oljpwl9qnym2j3PRq6g3V7efltlHuy1kjxHYiHMDM3g18jWDI7ffc/ctmdjHQ6e5Lzawd+CEwE3gJOM3d15VbZ7UnwiG+i2lKNTcnjOngnvPnNXz7pZTq6sp3m1W7rmrf2+EwrLAZDdf9LUpxXUcUZVlrWU8Uf6+VngiPNTQaoZbQiEuUB+eoDcc/MqnPcN7fpD5R/L0qNEYIHZxlKGl/k1IUGiIiUrFKQ0M/jS4iIhVTaIiISMUUGiIiUjGFhoiIVEyhISIiFWu60VNmtgX4c40vH8uAH0NMiKTWG5Jbd9U7WSqp94HuPm6wFTVdaNTDzDorGXLWbJJab0hu3VXvZImy3uqeEhGRiik0RESkYgqNXV0XdwFiktR6Q3LrrnonS2T11jkNERGpmFoaIiJSMYWGiIhUTKERMrP5ZrbGzNaa2flxl6dRzOx7ZrbZzFYVTNvHzG43s6fD/18fZxkbwcwOMLM7zWy1mT1uZp8Kpzd13c2s3cz+ZGaPhPX+j3D6FDN7INzfbzKzUXGXtRHMLGVmK8zs1+Hzpq+3ma03s8fMbKWZdYbTItvPFRoEOxZwNfAuYBqw0MymxVuqhrkemD9g2vnA79x9KvC78HmzyQGfc/dpwNuAT4SfcbPXvQuY5+6HAzOA+Wb2NuArwJXu/ibgZeAjMZaxkT4FPFHwPCn1/kd3n1FwbUZk+7lCIzAHWOvu69x9J/BT4MSYy9QQ7n4Xwa1zC50I/CB8/APgpCEt1BBw97+6+8Ph41cIDiQTaPK6e2BH+DQd/nNgHnBLOL3p6g1gZhOBE4DvhM+NBNS7hMj2c4VGYALwXMHzDeG0pNjf3f8aPn4e2D/OwjSamU0muO/8AySg7mEXzUpgM3A78Aywzd1z4SLNur9/Dfg80Bs+35dk1NuB35rZQ2Z2djgtsv28td7SSXNxdzezph2HbWZ7Aj8HPu3ufwu+fAaate7u3gPMMLMxwK3Am2MuUsOZ2XuAze7+kJnNjbs8Q+zt7r7RzPYDbjezJwtn1rufq6UR2AgcUPB8YjgtKV4wszcChP9vjrk8DWFmaYLA+LG7/yKcnIi6A7j7NuBO4ChgjJnlvzQ24/5+NLDAzNYTdDfPA75O89cbd98Y/r+Z4EvCHCLczxUagQeBqeHIilHAacDSmMs0lJYCZ4SPzwB+GWNZGiLsz/4u8IS7f7VgVlPX3czGhS0MzKwDeCfB+Zw7gX8KF2u6erv7Be4+0d0nE/w9L3f302nyepvZHmY2Ov8YOA5YRYT7ua4ID5nZuwn6QFPA99z9yzEXqSHM7CfAXIKfSn4B+HdgCXAzMIngZ+Xf7+4DT5aPaGb2duBu4DH6+7i/SHBeo2nrbmaHEZz4TBF8SbzZ3S82s4MIvoHvA6wAPujuXfGVtHHC7qn/5+7vafZ6h/W7NXzaCtzo7l82s32JaD9XaIiISMXUPSUiIhVTaIiISMUUGiIiUjGFhoiIVEyhISIiFVNoiIhIxRQaIiJSMYWGSIOZ2TvM7Idxl0MkCgoNkcY7nODqY5ERT6Eh0niHAyvMrM3MrjezS63w53VFRhD9NLpI4x1G8Kuiy4DvuPuPYi6PSM3021MiDRT+HPuLBD8S9zF3vy/mIonURd1TIo31FoKf3s8BPTGXRaRuCg2RxjocuJfgng7fN7Omu52sJItCQ6SxDgdWuftTwBeAm8MuK5ERSec0RESkYmppiIhIxRQaIiJSMYWGiIhUTKEhIiIVU2iIiEjFFBoiIlIxhYaIiFTs/wP8Qa6JHITvcgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 1000 # number of samples for input signal\n", "\n", "# generate input signal\n", "np.random.seed(1)\n", "x = np.random.normal(size=N, scale=1) # normally distributed (zero-mean, unit-variance) white noise\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=.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, 2016-2018*." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 1 }