{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"nbsphinx": "hidden"
},
"source": [
"# Spectral Estimation for Random Signals\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": [
"## Introduction\n",
"\n",
"In the preceding sections various statistical measures have been introduced to characterize random processes and signals. For instance, the probability density function (PDF) $p_x(\\theta)$, the mean value $\\mu_x$, the auto-correlation function (ACF) $\\varphi_{xx}[\\kappa]$ and its Fourier transformation, the power spectral density (PSD) $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. For many random processes whose internal structure is known, these measures can be derived in closed-form. However, for practical random signals measures of interest have to be estimated from a limited number of samples. These estimated quantities can e.g. be used to fit a parametric model of the random process or as parameters in algorithms."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem Statement\n",
"\n",
"The estimation of the spectral properties of a random signal is of special interest for spectral analysis. The discrete Fourier transform (DFT) of a random signal is also random. It is not very well suited to gain insights into the average spectral structure of a random signal. We aim at estimating the PSD $\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of a wide-sense stationary (WSS) and ergodic process from a limited number of samples. This is known as [*spectral (density) estimation*](https://en.wikipedia.org/wiki/Spectral_density_estimation). Many techniques have been developed for this purpose. They can be classified into\n",
"\n",
"1. non-parametric and\n",
"2. parametric\n",
"\n",
"techniques. Non-parametric techniques estimate the PSD of the random signal without assuming any particular structure for the generating random process. In contrary, parametric techniques assume that the generating random process can be modeled by a few parameters. Their aim is to estimate these parameters in order to characterize the spectral properties of the random signal."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Evaluation\n",
"\n",
"Various measures have been introduced in order to evaluate the performance of a particular estimation technique. The estimate $\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ can be regarded as a random signal itself. The performance of an estimator is therefore evaluated in a statistical sense. For the PSD, the following metrics are of interest\n",
"\n",
"#### Bias\n",
"\n",
"The [bias of an estimator](https://en.wikipedia.org/wiki/Estimator#Bias) \n",
"\n",
"\\begin{equation}\n",
"b_{\\hat{\\Phi}_{xx}} \n",
"= E\\{ \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) - \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\}\n",
"= E\\{ \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\} - \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\n",
"\\end{equation}\n",
"\n",
"quantifies the difference between the estimated $\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ and the true $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. An estimator is \n",
"* biased if $b_{\\hat{\\Phi}_{xx}} \\neq 0$ and \n",
"* bias-free if $b_{\\hat{\\Phi}_{xx}} = 0$.\n",
"\n",
"#### Variance\n",
"\n",
"The [variance of an estimator](https://en.wikipedia.org/wiki/Estimator#Variance)\n",
"\n",
"\\begin{equation}\n",
"\\sigma^2_{\\hat{\\Phi}_{xx}} = E \\left\\{ \\left(\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) - E\\{ \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\} \\right)^2 \\right\\}\n",
"\\end{equation}\n",
"\n",
"quantifies its quadratic deviation from its mean value $E\\{ \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\}$.\n",
"\n",
"#### Consistency\n",
"\n",
"A [consistent estimator](https://en.wikipedia.org/wiki/Estimator#Consistency) is an estimator for which the following conditions hold for a large number $N$ of samples:\n",
"\n",
"1. the estimator is unbiased\n",
" $$ \\lim_{N \\to \\infty} b_{\\hat{\\Phi}_{xx}} = 0 $$\n",
"\n",
"2. its variance converges towards zero\n",
" $$ \\lim_{N \\to \\infty} \\sigma^2_{\\hat{\\Phi}_{xx}} = 0 $$\n",
" \n",
"3. it converges in probability to the true $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$\n",
" $$ \\lim_{N \\to \\infty} \\Pr \\left\\{ | \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) - \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})| > \\alpha \\right\\} = 0$$\n",
" where $\\alpha > 0$ denotes a (small) constant.\n",
"\n",
"The last condition ensures that a consistent estimator does not generate outliers. Consistency is a desired property of an estimator. It ensures that if the number of samples $N$ increases towards infinity, the resulting estimates converges towards the true PSD."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example - Discrete Fourier transform of sample functions\n",
"\n",
"The following example computes and plots the magnitude of the discrete Fourier transform (DFT) $|X_n[\\mu]|$ of an ensemble of random signals $x_n[k]$. In the plot, each color denotes one sample function and $\\Omega[\\mu] = \\frac{2 \\pi}{N} \\mu$ where $N$ denotes the length of the sample function. The magnitude spectra are plotted as continuous signals for ease of illustration."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDYwNy4xMTg3NSAyODIuOTU2ODc1IF0gL1BhcmVudCAyIDAgUiAvUmVzb3VyY2VzIDggMCBSCi9UeXBlIC9QYWdlID4+CmVuZG9iago5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTEgMCBSID4+CnN0cmVhbQp4nL1aTY8cxw29z6+Yo3NQq8hifR2tOBEQIAhsC06AIAdFXisSJCWyrPjiH5/3WN3TVb272qw/VoCAHU43q5pFvvfIHjm/Psn55TmcX+P/j2c5P8X/l6eAT29POZRFpJaET2+GT1p1aSnzzze4dPr4r9Ppu9Pjz+HkA+56ejqZLM3vinWpxf96e0qtXbO+Gayaw5Ljar54mKy+0PvzDe5TqmdVWUo+f391/uv53fnx58r96PlPeDo86VLxrPDVgjQ8VjVs9PEXV/999eLqq6dPzi8+3Lzrmzdy+vr05fn9tkZANIc1EE9Y7nAX0xK3GJye4EB+PD15dn78RzlLOD/77hTj0nLKioVt0VIFfz379vRZWMLvzs9en//wzHfw64Ti9PQ8h0JiWIoWDXXa/WC+dzjudnlXSER1KcVCjWopXotLeoC4qLQlxZRznh5iMN87Lne7vCsuKrKkqsnEWsiHuMiD5EsMdYlmXHR8iMF877jc7fLOEgphiS201krkMx7i8hD5EhvwM2mqcX6I3Xz/uNzp8s641LZIa5ZKqyKHuOiD5IvVBHiTYDKj426+P9ze6fKuuFipSwilSVbNdi0uD5EvqRjub9baTJ27+d5xudvlXXEBCC1VcixmRcohLnHPl/e4NZwfBWwiho3VJSx6fvEWHk9lySY5iLNaWCzk7mNf6bO/09eJ9+SYyypDuMw//AsugjsD1wmXP1ZPWASeHn1x9fr5Nx+/fv7uw6O/vL16+fz8xb9PX167jZBp2lzATA50cvDPN6/ef7x69PHdqxCf/L67+vI3LYz/Xzn9bA2Cv+26t5uIpS4ptdjVSigr0T4EPGxuS1pqvBaF3XrPKIzuLlGYvH0yCgXgGHsUHgIMLm5FIHzatTAM5nvGYXJ4CcTs76ZIKGpZQMCtrST6oMkAebOUdD0Mu/m+YRgd7mGY/N0aBgOvpDUMD5sNNS9FrodhN983DKPDPQyTv1vDwD5QVop80GxQjUu+jpCD+Z5hmBxewjD7uy0MVOt1BUh90GwY9MDb25rp+4VhdLiHYfJ3axgSlMKKkHHQBZQEj8i+kpgzYJO1wLowICfHDrGzHPjJWT8uMTfpHB36Fbpd8Te/AqumkGp3Ee6nKdBnB5m/+Ok2sSG2iEroy2xb/s1lA6SKpwOWTrlmLekEyZLXfwXfFDx81Zi0nb86Zsp+skCsGHKEEn57Kihu3z2MiCx6AZyrQBXlvF5atGpTWmNZ0F2eLS+SBT0pbaktWbvXjDxF3xlpLmUJPSpWsJiwHYG52WJda1ldBBnlV2soS+q2hk6kVtokL9odAIXNam60KvQHFWdCFtZajHvQGJcS3YonRPOr4laI+Oa5kkBnGqpvWE0X3JVQzCknvxAbr/0oky3VUoiMgcKZ+FaRyzGiV/LbU1qSuuhN6D8Ff/luEYbKa7PAu7D7hjFDIJu7yJC8iExwF7nygyZn0JjEij8aRMU66CGWplKSmtsRQBoVfZ3W1terOKBq3HMtC+Q4Iu/mvAS3Igi4K2hwF43K1iMkeNgWUhDfYKP0h7UpWsYU/LGBHNIPFFYcd87cdAxxMX/saoBBRDm6FRnR6xFbw46CZjrBUgSBYBCLTG435U09ifLwLLkZf6/kk5jLSDZ3jHNeewsJOHKtHPTQjnRzH8BZ9CrJoxE5I+tpJRXl3YJV30ZkfHtEmTeWPTMiYlN68yEhLP25I/KUyhZ5guVi9HjicBaLPaDC72OJ5va0Vo2gCcInNJY048CT7xlHjNMxz+6YLxUu2H6MGpvbS1y0HwqKMJcM3HIziqHnLbMdaVi6HUfePPxYE8niJmRmz1xBCLTk7LvDcffqFYCd4Cuvndg8kH1BJLriX78cm/Krm4/V4JVzVxx4sYuXHFCC4nYU7nrimYMG8xQznLj17GW+ZjxjcXN2yABpVK3FTUA46amLcsNBBacSw7OmbgUC4lPtF4PHeuYiPbQE6ZuIqCyuBWxgI2nqxnqBoookSMmfA4fd/KHNQab01bgIEQtp1QJzyLBF29AR5SkeYcMZbYCFvDdksd+OwDtgRSZvKnSQmO2+Dk4C1FB5f8JTrdiCfeCZSnCzlRUCYBXklKdgwmJejgDqyk7EprE4rKW2YuVW60rvP59ahHzSaogNDEM+Ccm69Kif5BNlqhSN4JPatjVhBRED8gn8OMeeNYiO4NA7oQzkEZcEgeOPwWItvfgRdAHyrizBwHqKwVyQvJ0SdqIAreE8O/pvRAEj7k/ueCcKchXqzw9jZQQeSZGSZkZAJsWK3JgZAWkHqI1+QCMjAEwAelFGRiAnJREHiIEREpM1acfhiREK1sh5ZgSsgaimkQ+AGejTmmPDQAe8KUpyeB7oALsgwHtSjnQAZ0if2L1sdMCWxlppRWdCEHJNLWptYgQRli7T5EAJ1MLaOst2SgAYoQrMK3NghN47MIEmRsA2QHjAgAMjlLaEnEVGQmjRT1TazAcNPbvhtpkOqOWBIc3Pf6SDwhSqnrIbISD/ACCtdcweCIHDsMKx8kQImQGtHZ4BKZBJ3DKglziuM0swgvjeJdfEEoQ6NcfhgSUgygLOymxmCTA3jjOoo+vKEsJah5oNaaYJf0MQWulL7jyBuwDTzbN/4InGhMfZ60wTTLZsHpCJJlCIkm2ljwtNQCtryJy7TjThORaCK6WNJlBtBgXgCTPwBDVLQlnEiSfgQBCYmSWMxJC61NtZgmKvFk+CgQ+Uw0ztPNX5ACkqQGi34JQ2EELMsHrzXVW9gBDyFFrRooN/3UanRgEC9pAJ52PDcZWeF4nZsI7TmFtZHcim4V7Dbovdav3FSI/MLLqCO0AfvQR0L5Dfpm8+CfqBVYNcIuinntGwVbHmIR0gHweJAFk9QD7SKID9nQmQgLaxhqAD8BQQVG9NK00WdWx642+5Wq9BHCzqQ11JKlVlrzbSsuER/Gptl5Mx1oT1NmDvDtjiRAjMNnBBJkBmD/VIBoWl2HXFQAYQG+i6woEMwBEZFV6PZIAG1UFq4AJOkIDTLp13LoArFBlcHJoDaDdEcOWurTlApCKKSQ7NQaHCaRtJ7M0Bzi6n4NA1sIGixYgQkbmObKDQyIFEkw5s4LES7ZpzZwPooASAdU02sAG4o9bKEdFOBtB/gCKvoLE/QDGhie0YemEDJRsEdLb5wAbwleGkTnRAosHVNRz4ACBMwYUNzv1BIoxJjflACIkaN3S5cmkRxDuV1n0MjEAQqdKPfGcEdGHwHXo6DpwglPK4vFPWTgpGpkbizKSQE/e9bnDnBIJgaFaPnAAprymgqRg4gS8GLTTHsbF1SDiZjZt2SmiEVyjBfGgdiP6J84eJE1CTytOzIyfwfQGaRZs4IROQpfesIyVkNmOrmt8oATWHoq3tQAmssyyxzJTAhDIL7RopCEGtHboEl9vmGTnIfwjUBoHqFyMcG9Qoa6o6tM2qPkaJrhhSblslOYpBlc64DitayOpS8Gbrr4Dr1cU7CDqemGooYzx75phqHRqhZu/A9cbRDcW8hLJlCb4PbOq6bk8rsQq+LsXbSEGsZL2yIHNcpgziHAmJZOqQs+twcDD2F0bsBaYTNmUGWTRXLZuX3AiyHNjEUqcZDNKoJEifdEDZRvjPNituegCWHkEW0JRTirGNipuMJ4Cgg+JO6OwMgZUZZBE7uJA+8xkkd8YzQx6EaQLD6YeyUtpRcmdOmcoEsYiwaahtBlgqYHQffX7SAbYil5E0XpUjwFbKX5zbjLAEu1hCKPWAsJGgICtmXhAW1YyIZTtMYJSap0HS2IywgacTV/W6IyxwraKVqwO+crCmXd2P8xdjWvXOaJi/BMYbCW5HcCXwAErSCK5AusJ2+yC4sT5SA0VcZ3BFwClcQzyAK+eUOJje12yCO/tPKXqvMoAroA6NlPSOZ0dXZV+JZLIDuiKXJOKxZUJX/t5Jgd9yQFfwNTBQuwge0ZWU0OqMrsw7qNoSD0iqlBkltxlJUe5W8mEGQ20Lb1FHdQ1pWzIhc1LXsContHJR17Cgg+3ZhfJcC5tWnJ9e17pNuxK/2forIGVaBx25nbw6V9zkGN0x00r4pALmg7sEJlJeJuaRA5DcSWxUu0TH2jloh8rAtrCUaQhO/EULGg+aFlBbqD4Pmlb4nsXSUdMCkjhnnHE1GCq6HDRtZDNXvJmb4DZz7B8PcItiA/i3GW4r4LaZHUbeiWqDA+EJbouDg/VxzIC3gfmGIKZR1Cb+HoVDvxlwfRTbQCI2Ay40MMotdTIYABeKCmmnXYlvgFvgDgccDqI2c+QQko8Ed8Cl4gM4t45pg6aFsAjUEDvkEhWThpgOkEvJbUTRGXKhDc3aiuUD5AKRAOWhj6Y3yNVIJQk9LgdRm8F4uCDOUw6If4Q69anKjrh8mJaBT3nAXApg9APWByUD7HKeDkkVZ9TlbxJDzn3ksKMuX/tB18U2zzkKwot9HofhkSkN7DlIWvhQ/iTnKGkRdmSWevKtqMv5RIUothl0gadoIENpM+ay3YAAikdFm51i1xnRBXOBog0n6U5GzGVpNCtlnnIwO1X78Y6KtlEHpjwOOdhSYHmf/I04zJF8Ta77BxxGWFpe1egAxEqmrn0EcxlzCPGoT0WHMQcqK7KhvwAxUYdja795GHMEBrOPYLcBN8FMQ+uTaP4qJK7YnmrpwjnlsNVlpPzNUQ/DDAQURGB6q3V+K02U3t9I/5IfZt/0a6lPvRq/30+u7vviPfC51jfu4k/38vLiXP3FOXTPCsqa83YvX0P/+fnLd69++Pjt1fnDf65e/PD9x7eXHxac/geWUcGbCmVuZHN0cmVhbQplbmRvYmoKMTEgMCBvYmoKMzU5NQplbmRvYmoKMTYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA5MyA+PgpzdHJlYW0KeJw9jTsOwDAIQ3dOwQUqhU8CuU/VKb3/Wpd+FvRkY3tk58abK24049DJu5BYFJ6kLy2yoexjgHrmq32RRQ6ylHIdZVCqGjlLthb3t5bwrC3SCWMqSP6R4wLiXRv3CmVuZHN0cmVhbQplbmRvYmoKMTcgMCBvYmoKPDwgL0JCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzAzCi9TdWJ0eXBlIC9Gb3JtIC9UeXBlIC9YT2JqZWN0ID4+CnN0cmVhbQp4nD1SS67EMAjb5xS+QCUgfJJjvDNUmtXM4t1/80zbN6pUUxMMJh2/UAif+31+xqGGw6TwGTod4YX3MPuPNBY621FC94ZGEJXo0DVxMjOhRbYx80bWPRlbZA0qfV6xOmEiKINZIR2Ws+EcU7WD6RtVmMVuEuRdGE2DG7Em3IW1LRRiz6Sxv5HSn2VHwt4mCDpQ9QflqmMkeTOtQdzrSXDEYMfFQp2obFoLpQjjpBvBWQikZ9zf+wKPPkQ6FBwzYnd1lLRW04nVYznXQa2Fg4c0cNRj5eACfG0iTXIFjbTPlV/MZEVyC4kucDqiOadJ9m/IPj23cHWYSe8EZ4trtdVtjN5aynh7t7jRUke6s8W1/BLX4PXwlluply/gHZ6jZ7r/h++P8x6v8TP+ANOVb3UKZW5kc3RyZWFtCmVuZG9iagoxNCAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMtT2JsaXF1ZSAvQ2hhclByb2NzIDE1IDAgUgovRW5jb2RpbmcgPDwgL0RpZmZlcmVuY2VzIFsgODggL1ggXSAvVHlwZSAvRW5jb2RpbmcgPj4gL0ZpcnN0Q2hhciAwCi9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnREZXNjcmlwdG9yIDEzIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL1N1YnR5cGUgL1R5cGUzIC9UeXBlIC9Gb250IC9XaWR0aHMgMTIgMCBSID4+CmVuZG9iagoxMyAwIG9iago8PCAvQXNjZW50IDkyOSAvQ2FwSGVpZ2h0IDAgL0Rlc2NlbnQgLTIzNiAvRmxhZ3MgOTYKL0ZvbnRCQm94IFsgLTEwMTYgLTM1MSAxNjYwIDEwNjggXSAvRm9udE5hbWUgL0RlamFWdVNhbnMtT2JsaXF1ZQovSXRhbGljQW5nbGUgMCAvTWF4V2lkdGggMTM1MCAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTIgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTUgMCBvYmoKPDwgL1ggMTYgMCBSID4+CmVuZG9iagoyMiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDg3ID4+CnN0cmVhbQp4nD2OuxHAMAhDe6ZgBPMJhn1yqZz924A/adBDOk64CTYMz9Htws6BN0HuRS+wLm+ACiGLJBn59s7FolaZxa+mReqENEmcN1HoTotK5wPj9A54Prp1HXoKZW5kc3RyZWFtCmVuZG9iagoyMyAwIG9iago8PCAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzcKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnicPVG7ccUwDOs9BUbgR/xonneXKtm/DSg5KXiAKREE5Kcs0YWfZ4jg+1nu/8gDkq1QbYQnNBWRDdPA50kRWG6kJtxe3OeEbJUj9uJcIMIQ7TwJaaQLFjsZC94XP4+rHmasuWH8vjOafVR01VEdvHsO42ZNP06U3evNrI5bm/t0764Th2tIJp/3H5yUSqeXLIM6S7iwNpoa1uO8KMZYzDj+J6qwTbK2owrB0iVIKtCAGEoSxoDFLf4iJ1oOC9qbG2nrnclOqjSKhhejDN6g9UY4inSRfJhrK4OxqZg2vvnkJTfo+2e/n69fA2ta6wplbmRzdHJlYW0KZW5kb2JqCjI0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzA0ID4+CnN0cmVhbQp4nD2SO5LDMAxDe52CF8iM+JPk82Qnlff+7T4yyVaASYkAKC91mbKmPCBpJgn/0eHhYjvld9iezczAtUQvE8spz6ErxNxF+bKZjbqyOsWqwzCdW/SonIuGTZOa5ypLGbcLnsO1ieeWfcQPNzSoB3WNS8IN3dVoWQrNcHX/O71H2Xc1PBebVOrUF48XURXm+SFPoofpSuJ8PCghXHswRhYS5FPRQI6zXK3yXkL2DrcassJBaknnsyc82HV6Ty5uF80QD2S5VPhOUezt0DO+7EoJPRK24VjufTuasekamzjsfu9G1sqMrmghfshXJ+slYNxTJkUSZE62WG6L1Z7uoSimc4ZzGSDq2YqGUuZiV6t/DDtvLC/ZLMiUzAsyRqdNnjh4yH6NmvR5led4/QFs83M7CmVuZHN0cmVhbQplbmRvYmoKMjUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA1NCA+PgpzdHJlYW0KeJwzNjZXMFAwNDJX0DUyNlUwMjRQMDczUUgx5IIxc8EssGwOF1whhAmSz4GrzOFKAwBMkA8VCmVuZHN0cmVhbQplbmRvYmoKMjYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA3MSA+PgpzdHJlYW0KeJwztjRQMFCwMFPQNTQ2VDCyNFYwNzNQSDHkAgqBWLlcMLEcMMvMEsQyNDdDYumaGUJlkVgg43K4YAbnwMzL4UoDAPG0FiMKZW5kc3RyZWFtCmVuZG9iagoyNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY3ID4+CnN0cmVhbQp4nDO2NFAwULA0V9A1NDZUMDYwUTA3M1BIMeSCMXPBLLBsDhdMHYRlBmIYGZogscyAxoEl4QyQGTlw03K40gDOgxXTCmVuZHN0cmVhbQplbmRvYmoKMjggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzAgPj4Kc3RyZWFtCnicNVFJbsMwDLzrFfOBAOIuv8dBT+3/rx3SCWBgaEuczREbGxF4icHPQeTGW9aMmvibyV3xuzwVHgm3gidRBF6Ge9kJLm8Yl/04zHzwXlo5kxpPMiAX2fTwRMhgl0DowOwa1GGbaSf6hoTPjkg1G1lOX0vQS6sQKE/ZfqcLSrSt6s/tsy607WtPONntqSeVTyCeW7ICl41XTBZjGfRE5S7F9EGqs4WehPKifA6y+aghEl2inIEnBgejQDuw57afiVeFoHV1n7aNoRopHU//NjQ1SSLkEyWc2dK4W/j+nnv9/AOmVFOfCmVuZHN0cmVhbQplbmRvYmoKMjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMjcgPj4Kc3RyZWFtCnicNU87sgMhDOs5hS6QGYxtYM+zmVQv92+fZLINEv5I8vRERyZe5sgIrNnxthYZiBn4FlPxrz3tw4TqPbiHCOXiQphhJJw167ibp+PFv13lM9bBuw2+YpYXBLYwk/WVxZnLdsFYGidxTrIbY9dEbGNd6+kU1hFMKAMhne0wJcgcFSl9sqOMOTpO5InnYqrFLr/vYX3BpjGiwhxXBU/QZFCWPe8moB0X9N/Vjd9JNIteAjKRYGGdJObOWU741WtHx1GLIjEnpBnkMhHSnK5iCqEJxTo7CioVBZfqc8rdPv9oXVtNCmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDUgPj4Kc3RyZWFtCnicRVC7jUMxDOs9BRcIYP0se553SJXbvz1KRnCFIVo/kloSmIjASwyxlG/iR0ZBPQu/F4XiM8TPF4VBzoSkQJz1GRCZeIbaRm7odnDOvMMzjDkCF8VacKbTmfZc2OScBycQzm2U8YxCuklUFXFUn3FM8aqyz43XgaW1bLPTkewhjYRLSSUml35TKv+0KVsq6NpFE7BI5IGTTTThLD9DkmLMoJRR9zC1jvRxspFHddDJ2Zw5LZnZ7qftTHwPWCaZUeUpnecyPiep81xOfe6zHdHkoqVV+5z93pGW8iK126HV6VclUZmN1aeQuDz/jJ/x/gOOoFk+CmVuZHN0cmVhbQplbmRvYmoKMzEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDcgPj4Kc3RyZWFtCnicTVG7bUQxDOvfFFzgAOtreZ4LUl32b0PJCJDCIKEvKaclFvbGSwzhB1sPvuSRVUN/Hj8x7DMsPcnk1D/muclUFL4VqpuYUBdi4f1oBLwWdC8iK8oH349lDHPO9+CjEJdgJjRgrG9JJhfVvDNkwomhjsNBm1QYd00ULK4VzTPI7VY3sjqzIGx4JRPixgBEBNkXkM1go4yxlZDFch6oCpIFWmDX6RtRi4IrlNYJdKLWxLrM4Kvn9nY3Qy/y4Ki6eH0M60uwwuileyx8rkIfzPRMO3dJI73wphMRZg8FUpmdkZU6PWJ9t0D/n2Ur+PvJz/P9CxUoXCoKZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzOCA+PgpzdHJlYW0KeJxFUktyxTAI2+cUXCAz5mfj87xOV+n9t5VwOt089AwICTI9ZUim3DaWZITkHPKlV2SI1ZCfRo5ExBDfKaHArvK5vJbEXMhuiUrxoR0/l6U3Ms2u0Kq3R6c2i0Y1KyPnIEOEelbozO5R22TD63Yh6TpTFodwLP9DBbKUdcoplARtQd/YI+hvFjwR3Aaz5nKzuUxu9b/uWwue1zpbsW0HQAmWc95gBgDEwwnaAMTc2t4WKSgfVbqKScKt8lwnO1C20Kp0vDeAGQcYOWDDkq0O12hvAMM+D/SiRsX2FaCoLCD+ztlmwd4xyUiwJ+YGTj1xOsWRcEk4xgJAiq3iFLrxHdjiLxeuiJrwCXU6ZU28wp7a4sdCkwjvUnEC8CIbbl0dRbVsT+cJtD8qkjNipB7E0QmR1JLOERSXBvXQGvu4iRmvjcTmnr7dP8I5n+v7Fxa4g+AKZW5kc3RyZWFtCmVuZG9iagozMyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiFtCNEGUglgQpWYmZhBJOAMilwYAybQV5QplbmRzdHJlYW0KZW5kb2JqCjM0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjU1ID4+CnN0cmVhbQp4nEWRS5IDIAhE956CI4D85DyZmlVy/+00mEw2dpeo/YRKI6YSLOcUeTD9yPLNZLbptRyrnY0CiiIUzOQq9FiB1Z0p4sy1RLX1sTJy3Okdg+IN566cVLK4UcY6qjoVOKbnyvqq7vy4LMq+I4cyBWzWOQ42cOW2YYwTo81Wd4f7RJCnk6mj4naQbPiDk8a+ytUVuE42++olGAeCfqEJTPJNoHWGQOPmKXpyCfbxcbvzQLC3vAmkbAjkyBCMDkG7Tq5/cev83v86w53n2gxXjnfxO0xru+MvMcmKuYBF7hTU8z0XresMHe/JmWNy031D51ywy91Bps/8H+v3D1CKZogKZW5kc3RyZWFtCmVuZG9iagozNSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MSA+PgpzdHJlYW0KeJxFkEsSwyAMQ/ecQkfwRwZ8nnS6Su+/rSFNs4CnsUAGdycEqbUFE9EFL21Lugs+WwnOxnjoNm41EuQEdYBWpONolFJ9ucVplXTxaDZzKwutEx1mDnqUoxmgEDoV3u2i5HKm7s75R3D1X/VHse6czcTAZOUOhGb1Ke58mx1RXd1kf9JjbtZrfxX2qrC0rKXlhNvOXTOgBO6pHO39BalzOoQKZW5kc3RyZWFtCmVuZG9iagozNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDgwID4+CnN0cmVhbQp4nEWMuw3AMAhEe6ZgBH4mZp8olbN/GyBK3HBPunu4OhIyU95hhocEngwshlPxBpmjYDW4RlKNneyjsG5fdYHmelOr9fcHKk92dnE9zcsZ9AplbmRzdHJlYW0KZW5kb2JqCjM3IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjM2ID4+CnN0cmVhbQp4nE1QS25EIQzbc4pc4EkkIQHOQ9VV5/7bscNU7SqGGH9ID+myVR7rU2J1iezypU2XyjJ5FajlT9v/UQwCbv/QyEG0t4ydYuYS1sXCJDzlNCMbJ9csH487TxtmhcbEjeOdLhlgnxYBNVuVzYE5bTo3QLqQGreqs95kUAwi6kLNB5MunKfRl4g5nqhgSncmtZAbXD7VoQNxWr0KuWOLk2/EHFmhwGHQTHHWXwHWqMmyWcggSYYhzn2je5QKjajKeSsVwg+ToRH1htWgBpW5haKp5ZL8HdoCMAW2jHXpDEqBqgDB3yqnfb8BJI1dUwplbmRzdHJlYW0KZW5kb2JqCjM4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNDkgPj4Kc3RyZWFtCnicMza0UDBQMDQwB5JGhkCWkYlCiiEXSADEzOWCCeaAWQZAGqI4B64mhysNAMboDSYKZW5kc3RyZWFtCmVuZG9iagozOSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE1NyA+PgpzdHJlYW0KeJxFkLkRQzEIRHNVQQkSsAjqscfRd/+pF/lKtG8ALYevJVOqHyciptzXaPQweQ6fTSVWLNgmtpMachsWQUoxmHhOMaujt6GZh9TruKiquHVmldNpy8rFf/NoVzOTPcI16ifwTej4nzy0qehboK8LlH1AtTidSVAxfa9igaOcdn8inBjgPhlHmSkjcWJuCuz3GQBmvle4xuMF3QE3eQplbmRzdHJlYW0KZW5kb2JqCjQwIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzMyID4+CnN0cmVhbQp4nC1SOY4kMQzL/Qp+YADr8vGeHkzU+/90SVUFBapsyzzkcsNEJX4skNtRa+LXRmagwvCvq8yF70jbyDqIa8hFXMmWwmdELOQxxDzEgu/b+Bke+azMybMHxi/Z9xlW7KkJy0LGizO0wyqOwyrIsWDrIqp7eFOkw6kk2OOL/z7FcxeCFr4jaMAv+eerI3i+pEXaPWbbtFsPlmlHlRSWg+1pzsvkS+ssV8fj+SDZ3hU7QmpXgKIwd8Z5Lo4ybWVEa2Fng6TGxfbm2I+lBF3oxmWkOAL5mSrCA0qazGyiIP7I6SGnMhCmrulKJ7dRFXfqyVyzubydSTJb90WKzRTO68KZ9XeYMqvNO3mWE6VORfgZe7YEDZ3j6tlrmYVGtznBKyV8NnZ6cvK9mlkPyalISBXTugpOo8gUS9iW+JqKmtLUy/Dfl/cZf/8BM+J8AQplbmRzdHJlYW0KZW5kb2JqCjQxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTcgPj4Kc3RyZWFtCnicMza0UDCAwxRDLgAalALsCmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMzEgPj4Kc3RyZWFtCnicRY/LDQQhDEPvVOES8hk+qYfVntj+r+swmkFC+EEiO/EwCKzz8jbQxfDRosM3/jbVq2OVLB+6elJWD+mQh7zyFVBpMFHEhVlMHUNhzpjKyJYytxvhtk2DrGyVVK2DdjwGD7anZasIfqltYeos8QzCVV64xw0/kEutd71Vvn9CUzCXCmVuZHN0cmVhbQplbmRvYmoKNDMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzggPj4Kc3RyZWFtCnicNVI5rt1ADOt9Cl0ggHbNnOcFqX7u34aUXwpDtFaKmo4WlWn5ZSFVLZMuv+1JbYkb8vfJCokTklcl2qUMkVD5PIVUv2fLvL7WnBEgS5UKk5OSxyUL/gyX3i4c52NrP48jdz16YFWMhBIByxQTo2tZOrvDmo38PKYBP+IRcq5YtxxjFUgNunHaFe9D83nIGiBmmJaKCl1WiRZ+QfGgR61991hUWCDR7RxJcIyNUJGAdoHaSAw5sxa7qC/6WZSYCXTtiyLuosASScycYl06+g8+dCyovzbjy6+OSvpIK2tM2nejSWnMIpOul0VvN299PbhA8y7Kf17NIEFT1ihpfNCqnWMomhllhXccmgw0xxyHzBM8hzMSlPR9KH5fSya6KJE/Dg2hf18eo4ycBm8Bc9GftooDF/HZYa8cYIXSxZrkfUAqE3pg+v/X+Hn+/AMctoBUCmVuZHN0cmVhbQplbmRvYmoKNDQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDggPj4Kc3RyZWFtCnicLVE5kgNBCMvnFXpCc9PvscuR9//pCsoBg4ZDIDotcVDGTxCWK97yyFW04e+ZGMF3waHfynUbFjkQFUjSGFRNqF28Hr0HdhxmAvOkNSyDGesDP2MKN3pxeEzG2e11GTUEe9drT2ZQMisXccnEBVN12MiZw0+mjAvtXM8NyLkR1mUYpJuVxoyEI00hUkih6iapM0GQBKOrUaONHMV+6csjnWFVI2oM+1xL29dzE84aNDsWqzw5pUdXnMvJxQsrB/28zcBFVBqrPBAScL/bQ/2c7OQ33tK5s8X0+F5zsrwwFVjx5rUbkE21+Dcv4vg94+v5/AOopVsWCmVuZHN0cmVhbQplbmRvYmoKNDUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNzEgPj4Kc3RyZWFtCnicTZBNDkIhEIP3nKIXMKHzA4/zaFzp/bd28PnigvRLIUOnwwMdR+JGR4bO6HiwyTEOvAsyJl6N85+M6ySOCeoVbcG6tDvuzSwxJywTI2BrlNybRxT44ZgLQYLs8sMXGESka5hvNZ91k35+u9Nd1KV199MjCpzIjlAMG3AF2NM9DtwSzu+aJr9UKRmbOJQPVBeRstkJhailYpdTVWiM4lY974te7fkBwfY7+wplbmRzdHJlYW0KZW5kb2JqCjQ2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjEwID4+CnN0cmVhbQp4nDVQyw1DMQi7ZwoWqBQCgWSeVr11/2tt0DthEf9CWMiUCHmpyc4p6Us+OkwPti6/sSILrXUl7MqaIJ4r76GZsrHR2OJgcBomXoAWN2DoaY0aNXThgqYulUKBxSXwmXx1e+i+Txl4ahlydgQRQ8lgCWq6Fk1YtDyfkE4B4v9+w+4t5KGS88qeG/kbnO3wO7Nu4SdqdiLRchUy1LM0xxgIE0UePHlFpnDis9Z31TQS1GYLTpYBrk4/jA4AYCJeWYDsrkQ5S9KOpZ9vvMf3D0AAU7QKZW5kc3RyZWFtCmVuZG9iagoyMCAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMgL0NoYXJQcm9jcyAyMSAwIFIKL0VuY29kaW5nIDw8Ci9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0NiAvcGVyaW9kIDQ4IC96ZXJvIC9vbmUgL3R3byAvdGhyZWUgNTMgL2ZpdmUgNzcgL00gOTEKL2JyYWNrZXRsZWZ0IDkzIC9icmFja2V0cmlnaHQgOTcgL2EgOTkgL2MgL2QgL2UgMTAzIC9nIDEwNSAvaSAxMDkgL20gL24gMTEyCi9wIDExNCAvciAvcyAvdCAvdSAxMjQgL2JhciBdCi9UeXBlIC9FbmNvZGluZyA+PgovRmlyc3RDaGFyIDAgL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRm9udERlc2NyaXB0b3IgMTkgMCBSCi9Gb250TWF0cml4IFsgMC4wMDEgMCAwIDAuMDAxIDAgMCBdIC9MYXN0Q2hhciAyNTUgL05hbWUgL0RlamFWdVNhbnMKL1N1YnR5cGUgL1R5cGUzIC9UeXBlIC9Gb250IC9XaWR0aHMgMTggMCBSID4+CmVuZG9iagoxOSAwIG9iago8PCAvQXNjZW50IDkyOSAvQ2FwSGVpZ2h0IDAgL0Rlc2NlbnQgLTIzNiAvRmxhZ3MgMzIKL0ZvbnRCQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRm9udE5hbWUgL0RlamFWdVNhbnMgL0l0YWxpY0FuZ2xlIDAKL01heFdpZHRoIDEzNDIgL1N0ZW1WIDAgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9YSGVpZ2h0IDAgPj4KZW5kb2JqCjE4IDAgb2JqClsgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAKNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCAzMTggNDAxIDQ2MCA4MzggNjM2Cjk1MCA3ODAgMjc1IDM5MCAzOTAgNTAwIDgzOCAzMTggMzYxIDMxOCAzMzcgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNgo2MzYgNjM2IDMzNyAzMzcgODM4IDgzOCA4MzggNTMxIDEwMDAgNjg0IDY4NiA2OTggNzcwIDYzMiA1NzUgNzc1IDc1MiAyOTUKMjk1IDY1NiA1NTcgODYzIDc0OCA3ODcgNjAzIDc4NyA2OTUgNjM1IDYxMSA3MzIgNjg0IDk4OSA2ODUgNjExIDY4NSAzOTAgMzM3CjM5MCA4MzggNTAwIDUwMCA2MTMgNjM1IDU1MCA2MzUgNjE1IDM1MiA2MzUgNjM0IDI3OCAyNzggNTc5IDI3OCA5NzQgNjM0IDYxMgo2MzUgNjM1IDQxMSA1MjEgMzkyIDYzNCA1OTIgODE4IDU5MiA1OTIgNTI1IDYzNiAzMzcgNjM2IDgzOCA2MDAgNjM2IDYwMCAzMTgKMzUyIDUxOCAxMDAwIDUwMCA1MDAgNTAwIDEzNDIgNjM1IDQwMCAxMDcwIDYwMCA2ODUgNjAwIDYwMCAzMTggMzE4IDUxOCA1MTgKNTkwIDUwMCAxMDAwIDUwMCAxMDAwIDUyMSA0MDAgMTAyMyA2MDAgNTI1IDYxMSAzMTggNDAxIDYzNiA2MzYgNjM2IDYzNiAzMzcKNTAwIDUwMCAxMDAwIDQ3MSA2MTIgODM4IDM2MSAxMDAwIDUwMCA1MDAgODM4IDQwMSA0MDEgNTAwIDYzNiA2MzYgMzE4IDUwMAo0MDEgNDcxIDYxMiA5NjkgOTY5IDk2OSA1MzEgNjg0IDY4NCA2ODQgNjg0IDY4NCA2ODQgOTc0IDY5OCA2MzIgNjMyIDYzMiA2MzIKMjk1IDI5NSAyOTUgMjk1IDc3NSA3NDggNzg3IDc4NyA3ODcgNzg3IDc4NyA4MzggNzg3IDczMiA3MzIgNzMyIDczMiA2MTEgNjA1CjYzMCA2MTMgNjEzIDYxMyA2MTMgNjEzIDYxMyA5ODIgNTUwIDYxNSA2MTUgNjE1IDYxNSAyNzggMjc4IDI3OCAyNzggNjEyIDYzNAo2MTIgNjEyIDYxMiA2MTIgNjEyIDgzOCA2MTIgNjM0IDYzNCA2MzQgNjM0IDU5MiA2MzUgNTkyIF0KZW5kb2JqCjIxIDAgb2JqCjw8IC9NIDIyIDAgUiAvYSAyNCAwIFIgL2JhciAyNSAwIFIgL2JyYWNrZXRsZWZ0IDI2IDAgUgovYnJhY2tldHJpZ2h0IDI3IDAgUiAvYyAyOCAwIFIgL2QgMjkgMCBSIC9lIDMwIDAgUiAvZml2ZSAzMSAwIFIgL2cgMzIgMCBSCi9pIDMzIDAgUiAvbSAzNCAwIFIgL24gMzUgMCBSIC9vbmUgMzYgMCBSIC9wIDM3IDAgUiAvcGVyaW9kIDM4IDAgUgovciAzOSAwIFIgL3MgNDAgMCBSIC9zcGFjZSA0MSAwIFIgL3QgNDIgMCBSIC90aHJlZSA0MyAwIFIgL3R3byA0NCAwIFIKL3UgNDUgMCBSIC96ZXJvIDQ2IDAgUiA+PgplbmRvYmoKMyAwIG9iago8PCAvRjEgMjAgMCBSIC9GMiAxNCAwIFIgPj4KZW5kb2JqCjQgMCBvYmoKPDwgL0ExIDw8IC9DQSAwIC9UeXBlIC9FeHRHU3RhdGUgL2NhIDEgPj4KL0EyIDw8IC9DQSAxIC9UeXBlIC9FeHRHU3RhdGUgL2NhIDEgPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgL0YxLURlamFWdVNhbnMtT21lZ2EgMjMgMCBSIC9GMi1EZWphVnVTYW5zLU9ibGlxdWUtdW5pMDNCQyAxNyAwIFIgPj4KZW5kb2JqCjIgMCBvYmoKPDwgL0NvdW50IDEgL0tpZHMgWyAxMCAwIFIgXSAvVHlwZSAvUGFnZXMgPj4KZW5kb2JqCjQ3IDAgb2JqCjw8IC9DcmVhdGlvbkRhdGUgKEQ6MjAyMTAxMTUxNDM1MzQrMDInMDAnKQovQ3JlYXRvciAoTWF0cGxvdGxpYiB2My4zLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcpCi9Qcm9kdWNlciAoTWF0cGxvdGxpYiBwZGYgYmFja2VuZCB2My4zLjIpID4+CmVuZG9iagp4cmVmCjAgNDgKMDAwMDAwMDAwMCA2NTUzNSBmIAowMDAwMDAwMDE2IDAwMDAwIG4gCjAwMDAwMTUwMzEgMDAwMDAgbiAKMDAwMDAxNDc2MCAwMDAwMCBuIAowMDAwMDE0ODAzIDAwMDAwIG4gCjAwMDAwMTQ5MDIgMDAwMDAgbiAKMDAwMDAxNDkyMyAwMDAwMCBuIAowMDAwMDE0OTQ0IDAwMDAwIG4gCjAwMDAwMDAwNjUgMDAwMDAgbiAKMDAwMDAwMDM5OCAwMDAwMCBuIAowMDAwMDAwMjA4IDAwMDAwIG4gCjAwMDAwMDQwNjggMDAwMDAgbiAKMDAwMDAwNTIxNCAwMDAwMCBuIAowMDAwMDA1MDA2IDAwMDAwIG4gCjAwMDAwMDQ2OTEgMDAwMDAgbiAKMDAwMDAwNjI2NyAwMDAwMCBuIAowMDAwMDA0MDg5IDAwMDAwIG4gCjAwMDAwMDQyNTQgMDAwMDAgbiAKMDAwMDAxMzM5OSAwMDAwMCBuIAowMDAwMDEzMTk5IDAwMDAwIG4gCjAwMDAwMTI3MzcgMDAwMDAgbiAKMDAwMDAxNDQ1MiAwMDAwMCBuIAowMDAwMDA2Mjk5IDAwMDAwIG4gCjAwMDAwMDY0NTggMDAwMDAgbiAKMDAwMDAwNjgyOSAwMDAwMCBuIAowMDAwMDA3MjA2IDAwMDAwIG4gCjAwMDAwMDczMzIgMDAwMDAgbiAKMDAwMDAwNzQ3NSAwMDAwMCBuIAowMDAwMDA3NjE0IDAwMDAwIG4gCjAwMDAwMDc5MTcgMDAwMDAgbiAKMDAwMDAwODIxNyAwMDAwMCBuIAowMDAwMDA4NTM1IDAwMDAwIG4gCjAwMDAwMDg4NTUgMDAwMDAgbiAKMDAwMDAwOTI2NiAwMDAwMCBuIAowMDAwMDA5NDA2IDAwMDAwIG4gCjAwMDAwMDk3MzQgMDAwMDAgbiAKMDAwMDAwOTk2OCAwMDAwMCBuIAowMDAwMDEwMTIwIDAwMDAwIG4gCjAwMDAwMTA0MjkgMDAwMDAgbiAKMDAwMDAxMDU1MCAwMDAwMCBuIAowMDAwMDEwNzgwIDAwMDAwIG4gCjAwMDAwMTExODUgMDAwMDAgbiAKMDAwMDAxMTI3NCAwMDAwMCBuIAowMDAwMDExNDc4IDAwMDAwIG4gCjAwMDAwMTE4ODkgMDAwMDAgbiAKMDAwMDAxMjIxMCAwMDAwMCBuIAowMDAwMDEyNDU0IDAwMDAwIG4gCjAwMDAwMTUwOTEgMDAwMDAgbiAKdHJhaWxlcgo8PCAvSW5mbyA0NyAwIFIgL1Jvb3QgMSAwIFIgL1NpemUgNDggPj4Kc3RhcnR4cmVmCjE1MjQ4CiUlRU9GCg==\n",
"image/svg+xml": [
"\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",
"N = 256 # number of samples\n",
"M = 5 # number of sample functions\n",
"\n",
"# generate random signal\n",
"np.random.seed(1)\n",
"x = np.random.normal(size=(M, N))\n",
"h = sig.firwin2(N, [0, 0.4, 0.42, 0.65, 0.67, 1], [0, 0, 1, 1, 0, 0])\n",
"x = [np.convolve(xi, h, mode=\"same\") for xi in x]\n",
"\n",
"# DFT of signal\n",
"X = np.fft.rfft(x, axis=1)\n",
"Om = np.linspace(0, np.pi, X.shape[1])\n",
"\n",
"# plot signal and its spectrum\n",
"plt.figure(figsize=(10, 4))\n",
"plt.plot(Om, np.abs(X.T))\n",
"plt.title(\"Magnitude spectrum\")\n",
"plt.xlabel(r\"$\\Omega[\\mu]$\")\n",
"plt.ylabel(r\"$|X[\\mu]|$\")\n",
"plt.axis([0, np.pi, 0, 30])\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise**\n",
"\n",
"* Increase the number `N` of samples. What changes? What does not change with respect to the evaluation criteria introduced above?\n",
"* Is the DFT of a single sample function a consistent estimator for the spectral properties of a random process?\n",
"\n",
"Solution: Increasing the number of samples does only lead to an increase in the number of discrete frequencies $\\mu$. The amplitude of the fluctuations (variance) of the spectra within $1.3 < \\Omega < 2$ is not decreased when increasing the number of samples. The DFT of a single sample function is hence not a consistent estimator since at least the second condition is violated."
]
},
{
"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
}