{
"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": [
"## The Periodogram\n",
"\n",
"The [periodogram](https://en.wikipedia.org/wiki/Spectral_density_estimation#Periodogram) is an estimator for the power spectral density (PSD) $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of a random signal $x[k]$. For the following it is assumed that $x[k]$ is drawn from a wide-sense ergodic real-valued random process."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Definition\n",
"\n",
"The PSD is defined as the [discrete-time Fourier transformation (DTFT) of the auto-correlation function (ACF)](../random_signals/power_spectral_densities.ipynb#Definition)\n",
"\n",
"\\begin{equation}\n",
"\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\mathcal{F}_* \\{ \\varphi_{xx}[\\kappa] \\}\n",
"\\end{equation}\n",
"\n",
"Hence, the PSD can be computed from an estimate of the ACF. Let's assume that we want to estimate the PSD from $N$ samples of the random signal $x[k]$ by way of the ACF. The signal $x[k]$ is truncated to $N$ samples by multiplication (windowing) with the rectangular signal $\\text{rect}_N[k]$ of length $N$\n",
"\n",
"\\begin{equation}\n",
"x_N[k] = x[k] \\cdot \\text{rect}_N[k]\n",
"\\end{equation}\n",
"\n",
"where $x_N[k]$ denotes the truncated signal.\n",
"The ACF is estimated by applying its definition in a straightforward manner. For an ergodic random signal $x_N[k]$ of finite length, the estimated ACF $\\hat{\\varphi}_{xx}[\\kappa]$ can be expressed [in terms of a convolution](../random_signals/correlation_functions.ipynb#Definition)\n",
"\n",
"\\begin{equation}\n",
"\\hat{\\varphi}_{xx}[\\kappa] = \\frac{1}{N} \\cdot x_N[k] * x_N[-k]\n",
"\\end{equation}\n",
"\n",
"Applying the DTFT to both sides and rearranging the terms yields\n",
"\n",
"\\begin{equation}\n",
"\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{1}{N} \\, X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\, X_N(\\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega}) = \n",
"\\frac{1}{N} \\, | X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2\n",
"\\end{equation}\n",
"\n",
"where the intermediate and last equalities result from the symmetry relations of the DTFT. This estimate of the PSD is known as the periodogram. It can be computed directly from the DTFT\n",
"\n",
"\\begin{equation}\n",
"X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sum_{k=0}^{N-1} x_N[k] \\, \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,k}\n",
"\\end{equation}\n",
"\n",
"of the truncated random signal."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example - Periodogram\n",
"\n",
"The following example estimates the PSD of an ergodic random process which draws samples from normal distributed white noise with zero mean and unit variance. The true PSD is hence given as $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = 1$. In order to compute the periodogram by the discrete Fourier transformation (DFT), the signal $x[k]$ has to be zero-padded to ensure that squaring (multiplying) the spectra does not result in a circular convolution."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Bias of the periodogram: \t 0.0236\n",
"Variance of the periodogram: \t 0.7911\n"
]
},
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDU5NS4zMDYyNSAyODEuNjI4NzUgXSAvUGFyZW50IDIgMCBSIC9SZXNvdXJjZXMgOCAwIFIKL1R5cGUgL1BhZ2UgPj4KZW5kb2JqCjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMSAwIFIgPj4Kc3RyZWFtCnic3VtNrxy3Ebzvr5ijfFiKzebn0YJsAQGCWLGQXHIR5GdZgqUgkh3nkB+fqtkhhzvcx5WusSDjvVIPP7urqzkcWd6fZHm72OU9/v6xyPICf9+eLH77cAolGLXRBfz2a/eby2KiyykAtVe//XI6/Xx6+i2a+IxnXpxOao2sz2gywSuM0GzOR/DXDnShoJ8L2h7vwa0Td+nkLQaMwZuM4aNLIjd73TFv7Nbp6Rlm/8fp2avl6feyiF1e/XxyzkhQYZ/OlLC8+un0xBr7zfLq/fLdq9Pa0UkkmyI+2dx30aPzTkSsKVFCVueDXvUUrntyNpnkVGLse+rROz2VYlKy1osvNnY9yXFOrkSDiR+WrUfvLBx2MKTsS0ka8lVPhzlphid5a7NebVCHznvSlIzmWEIqWaTryR3n5JM34nPw0vfUo/OefIxGCgyjc9Ff9XSYU4j8t1hKuXL1Dp33FEIw1mpO3idJXU+6z+lfeMwuZwuPVxs2dy7GpUx3ffMBLZ5ohmiwtLTtBxMvBujv/Pzh/eu//f7j64+fz3/58PD29fL8n6eX+LPtzY3QcdHEY7zemkMyblE1yUeywc3Aac3HYFw8NN+wx5uPYrAne/NHz6rNF6xfPjTfsMebz8VIbs0PIdKaF2dNsof2d/DxDkSiQRzuPTw2Afin8e7YQwMnPXiP3W49DCGx95Dg28cd3sFJD8khULseHpuDs2rycZN38PEenIXj7Nusj84BFGjCcZ93cNKDA4HmrodH59AST99Dl40e7QEBHfad9l0YL2OQhZUxxaS4fHpY/r58XDAHZjQjTkLM0cHI0vcv/6UTfQ1J02lwZfnri+Xp84d/v3vz8NcXz5Y3n+cZ0GMe6tU6Dv/HE5jJx3A0biBoHdkjWt2ss8kH0wsCeeBjcj5d7LwilFaLnnkb6JMRpMocN+tqc0XT249JjUP+ApWttkEQQmuG7pm2gaLRpFLQ9mYewSoH4w3y5GTxW7uRNCn5ehA7mOHoOWksm7UHmxwG0bDgTCwubpOLhYntMIYdFESEOKazizmmm3rL9fecjNcIpt9ssCmhHBZsB7HOpUQm/NUaQ/fHiTUsZpOQfWLYbLF66TixHRQPt+ZgN+cpYuQwjAqhB5ddituQ6R8lHkaxg8y0CeSw7RpkSz76ZMMCs4AiMC62Yr1hlFxbdygIC65pxVd7KpWD8QXCKLy6pNVSkDP06Jgdmgp2HEzQ7BH6x23ZQQkJvCYxbLGE2DYxHF25Q8XB9cHH1Twinx+NN8x5iyjRUqNEEO6ajrveoRBRLmhyzd6bNAylgVSs0TOvVPNibDk6dYcyG6qLIW+OglkYf23ckIIlsaXaJZPtsIY7CkEaQ5HKXoJYk6Nb7aCoMzmDh6uvIPaDO/p3h2Y1KnCzuqMggHxsvWK+IMwR57VtChg/jGVHwQ8x2sZMAg6IwxI2UAScFqGm6kxBBGUg6g5F/JeEkqK6C37Xo/EF8gxRxaObJUI1DmTdobC3KDZsnSjYyI7mFUTmgy7XUscNQvADYXeokIIU424PBCTCo/mGOdC7Ws8fV2NqBRlou0PpZgjWYHN9AKnt6AA7CEdDPEslYwdiyAN5dyiWJkGyu9q6KFVwZ7oCHsnDZ43NCgYDf3coWg0JK755IaVBGebYQHh7kOLrioAS3EDhHQoWCsmVskW+AyfE41AqtoaxhXWsxp6p5TiUHUVRiydzDK11JLqj0+5gAtdJkDZRbFYaqLxDA7wdvdlmj5x3GMsGgSEMFhzEv5lilfzA5R0q8NuSoQy1PkAhMdhXEIIEy6+uDR21iQxk3qEYuvXgjeoD4IUwWFcMeYWFYGsbqW8g8w5FXMes2da5osJ1w1AaKIUKvnNaPB4GMu9QaipFZVHq0uS1Pu2NLwgThhSUWbVlkIIObN6hpDtbQpUTDqwQB+9qoGcSghauvoVcYAcy71AyC1KFrUlOQQt6aL1hkqGpRYNUpQtOSAOddyjoP0Boamu83BLdFVxzkoeK3cJUQQw3dPeOQlv4EmwV1KDrQXlvEDw+WvxQBT3PsAY679BkKZZD3U348Ki+dxDMAoUIMt2cUdXeEOAdyujO1lvbHgiDBG8YCAAVTEtFimdHFd6h1Ozopnq6UuoePWAHkXGxgi0lKlZqlOIdygOAmFLY2A7F1bUWvwDUqjaXOgSquoHMO5QpBeWg32JCeYgxzLCBULme9Fw3E5QwCvIOjfABH33VT4rNPSryhoGCvFdtLgVGGDV5hzKTq8ux2edRle9ghF7Q0lSCZr0hyzsUm4MJSG6VZTrK8gplymxI1WoJQhhleYdismja2dY0KGHQ5TuIbFtQe1RJCc+9Ics7dD2oKgj9ah8GXd4wuGKAL0izLTdkeYfiGTSetEabByEMunwHWVwEX1LJ1Tzf0OUdygzjkd5rmQ5WuJblGxI8GD+Vtiagg1GWd6isu5mLrxW9yqjLdzDBsyLUrVbrcEOWdyiLZIeJVvMyqPKG0RSsnKtC9KCMUZZ3qDikI8gXV+e68ulxCRtIJQWdXTnFgxJGWd6hTLZM03XsIAY9Gl8gdFGCtNLGgxNGWd6h2EuL4qZ6CmhpUOU7iASdvM1Vl3lwwqjKOxTEbJO06OSh+VGUN0xWVi5tikhdoybv0IL5IpG3oAAnDJJ8B9cjB/BtHUlONyR5h+ZgJIV9g7ABV5L8AqznHin7Oga+6BhYvEMF64485HRrNlg7avIdZLoHG0HbVPNwQ5V3aMC2JpZN1b4MqrxhngUonKvait4Q5R0KCWQyCEtcfSCPonwHWQ2VhLluaxNACaMq71BBfWNFs9QzN0TUQZZXCA6G+tq20zkwwqjKOxQtZymxLQsoYRDlO8j0DPlUaSVgoUZN3qHc3hS9pGrvB03esKJkWqkchHx+Q5N3KPwFXlgPTAJy7yDJdxARgVIhthVnjTbQeIciZSGrS4r1kFMOgnxDEL+r6qtnSCHGG4K8QyF/EMrIJ/UBbNigyHcwUVSqb54IXh8VeYcmHlf5WnoGUMJRj/cYCbTW7SHrDTneoWLX8xLIt/rA8J65f/nMhKXWWdmW8OXyZQfzX3gk/+nt45bL8fC+tzw9bvnp7fYO8PIK8HCcL8lavuVbnv7Z8qVeNcQ026GCtwmG+FVmhsl41LmeDjZt8EzhHJBFbcHoZ5ZImeIiCL/wQHBiGBFkAfQQw53JnEHiCC8bU4B2m1kqylgfCgLHljJt0kFBoOzGbOBGM8tgSU+uJKjAqeFZgomFcSpQGfO1RLpTqFzOyU/XkuVYjiFjLefbyFN7sHJy8wbPkS9YvBTk2DIfpGAfUXtBAlGMTy0xSkvH9Xk+G2E9l0pxKc6XEsyBNc/cxWmLQp5OyBjQolPDBKYDKaUCLTjtWRXSImlQT8E661pY6qEEQoiVO+vDqwHwjGzJXTP3VdBhLjzzjlPDsh4BleD4OmLW9frqBxzvYpgvUIBEAhMUN48bRR0F+spQ/24+aZ6YwHPwx/t5KGLWUdBxzhqnLgk/KxEtqs63xvPIKZbsxdm5nynvCRQIMRSc0/g6KyVnAKnldGc64GrsXoQAn684o4avLkmVc0O+g4KYZr09p7/IEpqn8cWGeXDzggR4ShA6Ls8tKRciz9R1bgh1C9p1BUJnagh5imoZSw5qubOQIN1kUd356S7ygpVPKYcc5tEAZQL5CJ0ewzS6EK6JWRYTmXsF0o2qaIA8vkMrfGFCfoykyTlZQAoV1DQFuWHaOS+I+ABZkO9MJzC+NFv4ZL5D+VghpDBwebqT61iBQTXkeC91g6GFcgU5R6dLpHwxFKOiOJ4nJtRNZrWzcIz5WgbmEQfFBC10Rw7wikKEoi9zcaPe0CQjy99JYeK4P6ir8ipTZxOCJS9D+YL5T0NHBLGjaBTSbh47wtMX8ch5YZ68+eoEZQ8v0SA3T9uE/o9YoMgwm/eurFkjk22euwfKShAGS+dUps7Ot+PZZyRGOPI8LBjlgjoy6DzZnvnGwCWfVayXObehoA0uWjhIsNO9dNYES4GPZu8Iu0glrzJn9DPmjaKA5GbnUX6WAokD1eQi4mLKB3wNJgmb7u7tI9/YgQ2QcPOdzj3PknjhocyVCw/ISpCSeAQ742BDva8pJT93DAOvCAnJaV5sOB5COkuVMQ9c3lYCUSd7JzuVYASFC2l1XmHxXgfcG5JO79AVa09sILg3zJNEpqFE3iOy87yTTOL5WOQ1hzlVOh4F+GKhU+dKCCTEdGex7vOSUWzh5QDLF/JzwXZGpU7tmcCWMo8vXgcOfHE719wWAhkS1Qtc6M7meLgP3XE+7VJ44zliv8uco898L8pLK9B/qDju14FoM/GFwHQpUYtpLJg6T9NmhjxIcZmnBXfkmvBFa4mRSmc6ysTDTcQD6t87zsbTCZBKipbvJ6fehkTmeSHJ8RLhLN+i0glIjNideYJguYrdLmvSmc688HKfsM7CXs69EpwWeNQ5159nnkahSobyBRXcWXVrQKUBWdfPRUHkfW6UwIpJzQnQeoPAxQKhOpl3zjIU9arF/tw5JeBVDZALWvXxDmclg5qIr2DmnReUMHBgONL8uAUhkR0b9HFeDJ6RGEEYMYFc50NEFuPLLfz188TIV1CwQ40Q50nnzJelBWUbarI5X4mHp+N/kIBunvDW85YUQQcx31nzAhHqQFtKYTCdkIG48TyBv7OUPNvLPALUMq+M8BOFcs4yFwRwNaW8KQq37AxffvHZq1v+1C5GY5kvZ6Z64nG+K16wVGT57SQVlfmXXYy++UHSV5wI76MSHvPyQiiHxgNhC56yggDLk6F0HyvsQ+m+VrgJbuNj37Z9APUVl8Cv7rDfPSw/WH/NKt68Sv9oy3adz+WjLlk/6Xrb7ta79TMjtO7FbqfKYXuWV+y/+/zbuw+vf3v4aXn98aflt0+/Pyw//Pi83rp/+q3yg7H2tRvW6/Kt2ynw/vVlLBK6EcqAciuokwsqow7mO4jdFLKB75N6MG2Gb04dCnG1t7qjyjd8vLjQd8WX8VtXbVQNe9PPoKFoFRkIDEMy6hqwcqOrDnRdqzLM4NcrdJ9s19W+LDfW9Q2/3ntWv95bA4cbfSN4lnnwhPYBFC8Hknouu+ZuwAfH2r8T7D6oCtoG63kfTC6nsvykKhkqosuL1DNTtd887un3bkn0ySf/oZed+FWZC5e3uXZ1yYrz+4ntTSnU0NrN+nz9YOTJP56slmowYY2X94e29rHZPKwmkFMSWZrjP3L55bXa1Wjer4YIQl4jufz7+dpy7/eb1faRr8Yu13rG78Z++OXdhbfXh1D8bk/hJ5fBULz8xff7oP7tpsKdb8/ufOwyvC87+IK0TXduD+AO5ccm+41hekKXPpZZ+ugdrW/cjejuZl82j0+bM37Zu7/DW8Jpu6eryftt8pcl4SuHUDo4QGPE6MsF4105XwAmFHQppBUshTc/wQhAoZiYvR0WADUSfuK1gdpVNOosdVIP1sV/c7pGi6wjuWo0m5g4lL7/YvjJW1qf34fKK4brQ/2k1LZJtfk3jM9TqqwD7E3X90Vr912jzLSc9XX/oGS3LlA/VHdrUu7m/N2tpeoa3Rd1779f/32o3Vbtk+o2ddz/jXvtmgFnBNg8+sJ/a2i37/068vjvyh3CazFVrl6x2/8LO479blO/QZuo5fjiUbb7STeI7/eP7yy/vdsJtHtcDSr67eG26l9FwIrNCxgqb2NFLsF99n15+h+uZ8UrCmVuZHN0cmVhbQplbmRvYmoKMTEgMCBvYmoKNDI4NgplbmRvYmoKMTcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMjcgPj4Kc3RyZWFtCnicPZJLjsQwCET3OQUXaMn8/DlPRrPquf+2H3Z6FlGhAEVRuGtKk+jyUpccKtmb/OgVU8XN5O+JhsQ0cTfQwSn3taMI/gS4DmbuDNFKiUG9dYnK8pGx89fX05cH78vbYbBRepaYV5+SsQYL8nR08QHm3Nruf5XvK5OOLL1KT0XvS71YlqgPMfti9SncxuYb23ownkzxazZRq5lT1toiugzURUo3sdULUO1RgeNcR28VbQMQgTtBB5UJtWpswUhItEXWL8xpQvfE/+0Bul/axHXsg9i0jWd8RRpD0N24R1nDzDHqfGpzw9rT2SbrlOJMLteyOFYJNErLWGpL8Kx6XRRjgxlopRDN0WpsZWh26OtleD/IDG87irGopUrPnjen4Fx97NcUevA8ix3Bm5zYLejnkKxCBi02YeosXa/wOed9/X4AKj19RgplbmRzdHJlYW0KZW5kb2JqCjE4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTUyID4+CnN0cmVhbQp4nD1PyxFDIQi8W8U2wIwggtbzMjmZ/q8BTTyxsrgf8YEKYhaQVIe4w63ixYW1o6vjU6QdtAqLg+YGlr8SsYK8gevW6Rg9Zpt4iufGGDpjhrBwzJEMWdrFM+62L0WODYK7YVah6SmWPuR6YRsHUnqztF2hpnAupiJjhnHbaZ9bJdKO0y9K/ZquIr3D1JK1i8affX8BvPc2ZwplbmRzdHJlYW0KZW5kb2JqCjE5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggOTIgPj4Kc3RyZWFtCnicPYyxDcAwCAR7pvgFImGMbdgnSuXs3+YtJ2ng9A/X0qA4rHF2VTQfOIt8eEv1hI3ElKaVR1Oc3doWDiuDFLvYFhZeYRGk8mqY8XlT1cCSUpTlzfp/dz3Hqxu6CmVuZHN0cmVhbQplbmRvYmoKMTUgMCBvYmoKPDwgL0Jhc2VGb250IC9EZWphVnVTYW5zLU9ibGlxdWUgL0NoYXJQcm9jcyAxNiAwIFIKL0VuY29kaW5nIDw8IC9EaWZmZXJlbmNlcyBbIDEwMSAvZSAxMDYgL2ogMTIwIC94IF0gL1R5cGUgL0VuY29kaW5nID4+Ci9GaXJzdENoYXIgMCAvRm9udEJCb3ggWyAtMTAxNiAtMzUxIDE2NjAgMTA2OCBdIC9Gb250RGVzY3JpcHRvciAxNCAwIFIKL0ZvbnRNYXRyaXggWyAwLjAwMSAwIDAgMC4wMDEgMCAwIF0gL0xhc3RDaGFyIDI1NSAvTmFtZSAvRGVqYVZ1U2Fucy1PYmxpcXVlCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDEzIDAgUiA+PgplbmRvYmoKMTQgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDk2Ci9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL0l0YWxpY0FuZ2xlIDAgL01heFdpZHRoIDEzNTAgL1N0ZW1WIDAgL1R5cGUgL0ZvbnREZXNjcmlwdG9yIC9YSGVpZ2h0IDAgPj4KZW5kb2JqCjEzIDAgb2JqClsgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAKNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCAzMTggNDAxIDQ2MCA4MzggNjM2Cjk1MCA3ODAgMjc1IDM5MCAzOTAgNTAwIDgzOCAzMTggMzYxIDMxOCAzMzcgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNgo2MzYgNjM2IDMzNyAzMzcgODM4IDgzOCA4MzggNTMxIDEwMDAgNjg0IDY4NiA2OTggNzcwIDYzMiA1NzUgNzc1IDc1MiAyOTUKMjk1IDY1NiA1NTcgODYzIDc0OCA3ODcgNjAzIDc4NyA2OTUgNjM1IDYxMSA3MzIgNjg0IDk4OSA2ODUgNjExIDY4NSAzOTAgMzM3CjM5MCA4MzggNTAwIDUwMCA2MTMgNjM1IDU1MCA2MzUgNjE1IDM1MiA2MzUgNjM0IDI3OCAyNzggNTc5IDI3OCA5NzQgNjM0IDYxMgo2MzUgNjM1IDQxMSA1MjEgMzkyIDYzNCA1OTIgODE4IDU5MiA1OTIgNTI1IDYzNiAzMzcgNjM2IDgzOCA2MDAgNjM2IDYwMCAzMTgKMzUyIDUxOCAxMDAwIDUwMCA1MDAgNTAwIDEzNTAgNjM1IDQwMCAxMDcwIDYwMCA2ODUgNjAwIDYwMCAzMTggMzE4IDUxOCA1MTgKNTkwIDUwMCAxMDAwIDUwMCAxMDAwIDUyMSA0MDAgMTAyOCA2MDAgNTI1IDYxMSAzMTggNDAxIDYzNiA2MzYgNjM2IDYzNiAzMzcKNTAwIDUwMCAxMDAwIDQ3MSA2MTcgODM4IDM2MSAxMDAwIDUwMCA1MDAgODM4IDQwMSA0MDEgNTAwIDYzNiA2MzYgMzE4IDUwMAo0MDEgNDcxIDYxNyA5NjkgOTY5IDk2OSA1MzEgNjg0IDY4NCA2ODQgNjg0IDY4NCA2ODQgOTc0IDY5OCA2MzIgNjMyIDYzMiA2MzIKMjk1IDI5NSAyOTUgMjk1IDc3NSA3NDggNzg3IDc4NyA3ODcgNzg3IDc4NyA4MzggNzg3IDczMiA3MzIgNzMyIDczMiA2MTEgNjA4CjYzMCA2MTMgNjEzIDYxMyA2MTMgNjEzIDYxMyA5OTUgNTUwIDYxNSA2MTUgNjE1IDYxNSAyNzggMjc4IDI3OCAyNzggNjEyIDYzNAo2MTIgNjEyIDYxMiA2MTIgNjEyIDgzOCA2MTIgNjM0IDYzNCA2MzQgNjM0IDU5MiA2MzUgNTkyIF0KZW5kb2JqCjE2IDAgb2JqCjw8IC9lIDE3IDAgUiAvaiAxOCAwIFIgL3ggMTkgMCBSID4+CmVuZG9iagoyNCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MCA+PgpzdHJlYW0KeJw9kEsSwyAMQ/ecQkfA+H+edLpK7r+tDZ1ssBiE9MB9YiKjFieCr8SHBqXDJPBsFYR7MNkRcoTkBE2GsoMkcQ0NBqXCpmOZ78mmddJKrLzRftl3NGaddIotRYd2If/n9SLco+Aa6xk8D2AxyNpKpeyZMFplpq7yqOi1H9PhPQ9Eq8Xl9Qau8NpHN6koKkvq/kR3NNj+kbf7Ht8fmWU4JAplbmRzdHJlYW0KZW5kb2JqCjI1IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNzkgPj4Kc3RyZWFtCnicTc27DcAgDATQnik8AuD/PlGqsH8bGyJCYz/pTjrBDhXc4rAYaHe4WvGlUZh96pkSklBzPURYMyU6hKRf+ssww5jYyLbvt1buF94bHBkKZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzcKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnicPVG7ccUwDOs9BUbgR/xonneXKtm/DSg5KXiAKREE5Kcs0YWfZ4jg+1nu/8gDkq1QbYQnNBWRDdPA50kRWG6kJtxe3OeEbJUj9uJcIMIQ7TwJaaQLFjsZC94XP4+rHmasuWH8vjOafVR01VEdvHsO42ZNP06U3evNrI5bm/t0764Th2tIJp/3H5yUSqeXLIM6S7iwNpoa1uO8KMZYzDj+J6qwTbK2owrB0iVIKtCAGEoSxoDFLf4iJ1oOC9qbG2nrnclOqjSKhhejDN6g9UY4inSRfJhrK4OxqZg2vvnkJTfo+2e/n69fA2ta6wplbmRzdHJlYW0KZW5kb2JqCjI3IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTY1ID4+CnN0cmVhbQp4nEWPOxIDIQxDe06hI4B/wHk2k4q9fxvLO0kaLIwlP6IrOvbKw2NjysZrtLEnwhbuUjoNp6mMr4qnZ12gy2EyU29czVxgqrDIbk6x+hh8ofLs5oSvVZ4YwpdMCQ0wlTu5h/X6UZyWfCS7C4LqlI3KwjBH0vdATE2bp4WB/I8veWpBUJnmjWuWlUdrFVM0Z5gqWwuC9YGgOqX6A9P/TKe9P9z0PYAKZW5kc3RyZWFtCmVuZG9iagoyOCAwIG9iago8PCAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMzgKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnicRZBNssMwCIP3PgVHMJgf+zydeav0/ttKJH1dZGAMSPoylodEHnkPO+xMLFDxauYSWvIaulW8jmhOcdtdVzon7OYUbpifvjFVTKiiewlVNUPoo3Ew+Tpew33123v4sd6KZX0V+dXJma2cVu3Utb3ZMQ03mI43d16qkICqJKIPHV//jlenqCI2bDfcSDHB6OCLhd0N48moYtXltkXjYYJpmPV6YIxfcRBS8Qugl+4NnKEPcNkX+Ndxeo0wdAy4gaL41nkUC5HoUUvbkvXOwM4qhRvMyJs7NFWIQVVQtWE94PPBnj/8a/x9AKh+YIAKZW5kc3RyZWFtCmVuZG9iagoyOSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzOCA+PgpzdHJlYW0KeJw1UkuSW0EI279T6AKuav7NeZya1eT+2wjsrKCBFhKQFjjIxEsMUY1yxR95gvE6gb/r5Wn8Pt6F1IKnIv3AtWkb78eaNVGwNGIpzD72/Sghx1Pj3xDouUgTZmQyciAZiPu1Pn/Wm0w5/AakaXP6KEl6EC3Y3Rp2fFmQQdKTGpbs5Id1LbC6CE2YG2siGTm1MjXPx57hMp4YI0HVLCBJn7hPFYxIMx47Zy15kOF4qhcvfr2N1zKPqZdVBTK2CeZgO5kJpygiEL+gJLmJu2jqKI5mxprbhYaSIvfdPZyc9Lq/nEQFXgnhLNYSjhl6yjInOw1KoGrlBJhhvfaFcZo2SrhT0+1dsa/fZyZh3Oaws1IyDc5xcC+bzBEke90xYRMeh5j37hGMxLz5XWwRXLnMuSbTj/0o2kgfFNfnXE2ZrSjhH6rkiRXX+P/83s/PP5A3fbEKZW5kc3RyZWFtCmVuZG9iagozMCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMwNCA+PgpzdHJlYW0KeJw9kjuSwzAMQ3udghfIjPiT5PNkJ5X3/u0+MslWgEmJACgvdZmypjwgaSYJ/9Hh4WI75XfYns3MwLVELxPLKc+hK8TcRfmymY26sjrFqsMwnVv0qJyLhk2TmucqSxm3C57DtYnnln3EDzc0qAd1jUvCDd3VaFkKzXB1/zu9R9l3NTwXm1Tq1BePF1EV5vkhT6KH6UrifDwoIVx7MEYWEuRT0UCOs1yt8l5C9g63GrLCQWpJ57MnPNh1ek8ubhfNEA9kuVT4TlHs7dAzvuxKCT0StuFY7n07mrHpGps47H7vRtbKjK5oIX7IVyfrJWDcUyZFEmROtlhui9We7qEopnOGcxkg6tmKhlLmYlerfww7bywv2SzIlMwLMkanTZ44eMh+jZr0eZXneP0BbPNzOwplbmRzdHJlYW0KZW5kb2JqCjMxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNTQgPj4Kc3RyZWFtCnicMzY2VzBQMDQyV9A1MjZVMDI0UDA3M1FIMeSCMXPBLLBsDhdcIYQJks+Bq8zhSgMATJAPFQplbmRzdHJlYW0KZW5kb2JqCjMyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjI3ID4+CnN0cmVhbQp4nDVPO7IDIQzrOYUukBmMbWDPs5lUL/dvn2SyDRL+SPL0REcmXubICKzZ8bYWGYgZ+BZT8a897cOE6j24hwjl4kKYYSScNeu4m6fjxb9d5TPWwbsNvmKWFwS2MJP1lcWZy3bBWBoncU6yG2PXRGxjXevpFNYRTCgDIZ3tMCXIHBUpfbKjjDk6TuSJ52KqxS6/72F9waYxosIcVwVP0GRQlj3vJqAdF/Tf1Y3fSTSLXgIykWBhnSTmzllO+NVrR8dRiyIxJ6QZ5DIR0pyuYgqhCcU6OwoqFQWX6nPK3T7/aF1bTQplbmRzdHJlYW0KZW5kb2JqCjMzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ1ID4+CnN0cmVhbQp4nEVQu41DMQzrPQUXCGD9LHued0iV2789SkZwhSFaP5JaEpiIwEsMsZRv4kdGQT0LvxeF4jPEzxeFQc6EpECc9RkQmXiG2kZu6HZwzrzDM4w5AhfFWnCm05n2XNjknAcnEM5tlPGMQrpJVBVxVJ9xTPGqss+N14GltWyz05HsIY2ES0klJpd+Uyr/tClbKujaRROwSOSBk0004Sw/Q5JizKCUUfcwtY70cbKRR3XQydmcOS2Z2e6n7Ux8D1gmmVHlKZ3nMj4nqfNcTn3usx3R5KKlVfuc/d6RlvIitduh1elXJVGZjdWnkLg8/4yf8f4DjqBZPgplbmRzdHJlYW0KZW5kb2JqCjM0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ3ID4+CnN0cmVhbQp4nE1Ru21EMQzr3xRc4ADra3meC1Jd9m9DyQiQwiChLymnJRb2xksM4QdbD77kkVVDfx4/MewzLD3J5NQ/5rnJVBS+FaqbmFAXYuH9aAS8FnQvIivKB9+PZQxzzvfgoxCXYCY0YKxvSSYX1bwzZMKJoY7DQZtUGHdNFCyuFc0zyO1WN7I6syBseCUT4sYARATZF5DNYKOMsZWQxXIeqAqSBVpg1+kbUYuCK5TWCXSi1sS6zOCr5/Z2N0Mv8uCounh9DOtLsMLopXssfK5CH8z0TDt3SSO98KYTEWYPBVKZnZGVOj1ifbdA/59lK/j7yc/z/QsVKFwqCmVuZHN0cmVhbQplbmRvYmoKMzUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA5MCA+PgpzdHJlYW0KeJxNjUESwCAIA++8Ik9QRND/dHrS/1+r1A69wE4CiRZFgvQ1aksw7rgyFWtQKZiUl8BVMFwL2u6iyv4ySUydhtN7twODsvFxg9JJ+/ZxegCr/XoG3Q/SHCJYCmVuZHN0cmVhbQplbmRvYmoKMzYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA2OCA+PgpzdHJlYW0KeJwzMrdQMFCwNAEShhYmCuZmBgophlxAvqmJuUIuF0gMxMoBswyAtCWcgohbQjRBlIJYEKVmJmYQSTgDIpcGAMm0FeUKZW5kc3RyZWFtCmVuZG9iagozNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI1NSA+PgpzdHJlYW0KeJxFkUuSAyAIRPeegiOA/OQ8mZpVcv/tNJhMNnaXqP2ESiOmEiznFHkw/cjyzWS26bUcq52NAooiFMzkKvRYgdWdKeLMtUS19bEyctzpHYPiDeeunFSyuFHGOqo6FTim58r6qu78uCzKviOHMgVs1jkONnDltmGME6PNVneH+0SQp5Opo+J2kGz4g5PGvsrVFbhONvvqJRgHgn6hCUzyTaB1hkDj5il6cgn28XG780Cwt7wJpGwI5MgQjA5Bu06uf3Hr/N7/OsOd59oMV4538TtMa7vjLzHJirmARe4U1PM9F63rDB3vyZljctN9Q+dcsMvdQabP/B/r9w9QimaICmVuZHN0cmVhbQplbmRvYmoKMzggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjEgPj4Kc3RyZWFtCnicRZBLEsMgDEP3nEJH8EcGfJ50ukrvv60hTbOAp7FABncnBKm1BRPRBS9tS7oLPlsJzsZ46DZuNRLkBHWAVqTjaJRSfbnFaZV08Wg2cysLrRMdZg56lKMZoBA6Fd7touRypu7O+Udw9V/1R7HunM3EwGTlDoRm9SnufJsdUV3dZH/SY27Wa38V9qqwtKyl5YTbzl0zoATuqRzt/QWpczqECmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4MCA+PgpzdHJlYW0KeJxFjLsNwDAIRHumYAR+JmafKJWzfxsgStxwT7p7uDoSMlPeYYaHBJ4MLIZT8QaZo2A1uEZSjZ3so7BuX3WB5npTq/X3BypPdnZxPc3LGfQKZW5kc3RyZWFtCmVuZG9iago0MCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE0NyA+PgpzdHJlYW0KeJw9T7kNAzEM6z0FFzjAeixb81yQ6rJ/G8pGUggiQPGRZUfHClxiApOOORIvaT/4aRqBWAY1R/SEimFY4G6SAg+DLEpXni1eDJHaQl1I+NYQ3q1MZKI8rxE7cCcXowc+VBtZHnpAO0QVWa5Jw1jVVl1qnbACHLLOwnU9zKoE5dEnaykfUFRCvXT/n3va+wsAby/rCmVuZHN0cmVhbQplbmRvYmoKNDEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNDkgPj4Kc3RyZWFtCnicNY9LDgMhDEP3OYUvMFJ+hHAeqq6m9982YVoJCQvbL8EWg5GMS0xg7Jhj4SVUT60+JCOPukk5EKlQNwRPaEwMM2zSJfDKdN8ynlu8nFbqgk5I5OmsNhqijGZew9FTzgqb/svcJGplRpkDMutUtxOysmAF5gW1PPcz7qhc6ISHncqw6E4xotxmRhp+/9v0/gJ7MjBjCmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA0OSA+PgpzdHJlYW0KeJwzNrRQMFAwNDAHkkaGQJaRiUKKIRdIAMTM5YIJ5oBZBkAaojgHriaHKw0AxugNJgplbmRzdHJlYW0KZW5kb2JqCjQzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTU3ID4+CnN0cmVhbQp4nEWQuRFDMQhEc1VBCRKwCOqxx9F3/6kX+Uq0bwAth68lU6ofJyKm3Ndo9DB5Dp9NJVYs2Ca2kxpyGxZBSjGYeE4xq6O3oZmH1Ou4qKq4dWaV02nLysV/82hXM5M9wjXqJ/BN6PifPLSp6FugrwuUfUC1OJ1JUDF9r2KBo5x2fyKcGOA+GUeZKSNxYm4K7PcZAGa+V7jG4wXdATd5CmVuZHN0cmVhbQplbmRvYmoKNDQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzIgPj4Kc3RyZWFtCnicLVI5jiQxDMv9Cn5gAOvy8Z4eTNT7/3RJVQUFqmzLPORyw0QlfiyQ21Fr4tdGZqDC8K+rzIXvSNvIOohryEVcyZbCZ0Qs5DHEPMSC79v4GR75rMzJswfGL9n3GVbsqQnLQsaLM7TDKo7DKsixYOsiqnt4U6TDqSTY44v/PsVzF4IWviNowC/556sjeL6kRdo9Ztu0Ww+WaUeVFJaD7WnOy+RL6yxXx+P5INneFTtCaleAojB3xnkujjJtZURrYWeDpMbF9ubYj6UEXejGZaQ4AvmZKsIDSprMbKIg/sjpIacyEKau6Uont1EVd+rJXLO5vJ1JMlv3RYrNFM7rwpn1d5gyq807eZYTpU5F+Bl7tgQNnePq2WuZhUa3OcErJXw2dnpy8r2aWQ/JqUhIFdO6Ck6jyBRL2Jb4moqa0tTL8N+X9xl//wEz4nwBCmVuZHN0cmVhbQplbmRvYmoKNDUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNyA+PgpzdHJlYW0KeJwzNrRQMIDDFEMuABqUAuwKZW5kc3RyZWFtCmVuZG9iago0NiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMSA+PgpzdHJlYW0KeJxFj8sNBCEMQ+9U4RLyGT6ph9We2P6v6zCaQUL4QSI78TAIrPPyNtDF8NGiwzf+NtWrY5UsH7p6UlYP6ZCHvPIVUGkwUcSFWUwdQ2HOmMrIljK3G+G2TYOsbJVUrYN2PAYPtqdlqwh+qW1h6izxDMJVXrjHDT+QS613vVW+f0JTMJcKZW5kc3RyZWFtCmVuZG9iago0NyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzOCA+PgpzdHJlYW0KeJw1Ujmu3UAM630KXSCAds2c5wWpfu7fhpRfCkO0VoqajhaVafllIVUtky6/7UltiRvy98kKiROSVyXapQyRUPk8hVS/Z8u8vtacESBLlQqTk5LHJQv+DJfeLhznY2s/jyN3PXpgVYyEEgHLFBOja1k6u8Oajfw8pgE/4hFyrli3HGMVSA26cdoV70PzecgaIGaYlooKXVaJFn5B8aBHrX33WFRYINHtHElwjI1QkYB2gdpIDDmzFruoL/pZlJgJdO2LIu6iwBJJzJxiXTr6Dz50LKi/NuPLr45K+kgra0zad6NJacwik66XRW83b309uEDzLsp/Xs0gQVPWKGl80KqdYyiaGWWFdxyaDDTHHIfMEzyHMxKU9H0ofl9LJrookT8ODaF/Xx6jjJwGbwFz0Z+2igMX8dlhrxxghdLFmuR9QCoTemD6/9f4ef78Axy2gFQKZW5kc3RyZWFtCmVuZG9iago0OCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI0OCA+PgpzdHJlYW0KeJwtUTmSA0EIy+cVekJz0++xy5H3/+kKygGDhkMgOi1xUMZPEJYr3vLIVbTh75kYwXfBod/KdRsWORAVSNIYVE2oXbwevQd2HGYC86Q1LIMZ6wM/Ywo3enF4TMbZ7XUZNQR712tPZlAyKxdxycQFU3XYyJnDT6aMC+1czw3IuRHWZRikm5XGjIQjTSFSSKHqJqkzQZAEo6tRo40cxX7pyyOdYVUjagz7XEvb13MTzho0OxarPDmlR1ecy8nFCysH/bzNwEVUGqs8EBJwv9tD/Zzs5Dfe0rmzxfT4XnOyvDAVWPHmtRuQTbX4Ny/i+D3j6/n8A6ilWxYKZW5kc3RyZWFtCmVuZG9iago0OSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE3MSA+PgpzdHJlYW0KeJxNkE0OQiEQg/ecohcwofMDj/NoXOn9t3bw+eKC9EshQ6fDAx1H4kZHhs7oeLDJMQ68CzImXo3zn4zrJI4J6hVtwbq0O+7NLDEnLBMjYGuU3JtHFPjhmAtBguzywxcYRKRrmG81n3WTfn67013UpXX30yMKnMiOUAwbcAXY0z0O3BLO75omv1QpGZs4lA9UF5Gy2QmFqKVil1NVaIziVj3vi17t+QHB9jv7CmVuZHN0cmVhbQplbmRvYmoKNTAgMCBvYmoKPDwgL0JCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNzIKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnic49I1sjBVsDAwUMjl0jUyNAYzc7h0LY0VzAzNQCxDM0MY08jEUsHcGMw0NjaHiZoYmMIVQM2CqjU1gxgLZeZwpQEAk4MVTgplbmRzdHJlYW0KZW5kb2JqCjUxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjEwID4+CnN0cmVhbQp4nDVQyw1DMQi7ZwoWqBQCgWSeVr11/2tt0DthEf9CWMiUCHmpyc4p6Us+OkwPti6/sSILrXUl7MqaIJ4r76GZsrHR2OJgcBomXoAWN2DoaY0aNXThgqYulUKBxSXwmXx1e+i+Txl4ahlydgQRQ8lgCWq6Fk1YtDyfkE4B4v9+w+4t5KGS88qeG/kbnO3wO7Nu4SdqdiLRchUy1LM0xxgIE0UePHlFpnDis9Z31TQS1GYLTpYBrk4/jA4AYCJeWYDsrkQ5S9KOpZ9vvMf3D0AAU7QKZW5kc3RyZWFtCmVuZG9iagoyMiAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMgL0NoYXJQcm9jcyAyMyAwIFIKL0VuY29kaW5nIDw8Ci9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0MCAvcGFyZW5sZWZ0IC9wYXJlbnJpZ2h0IDQ2IC9wZXJpb2QgNDggL3plcm8gL29uZSAvdHdvIC90aHJlZQovZm91ciAvZml2ZSA2OCAvRCAvRSA4MCAvUCA4MyAvUyA5NyAvYSAxMDAgL2QgL2UgMTA1IC9pIDEwOSAvbSAvbiAxMTQgL3IgL3MKL3QgL3UgMTI0IC9iYXIgXQovVHlwZSAvRW5jb2RpbmcgPj4KL0ZpcnN0Q2hhciAwIC9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnREZXNjcmlwdG9yIDIxIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDIwIDAgUiA+PgplbmRvYmoKMjEgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zIC9JdGFsaWNBbmdsZSAwCi9NYXhXaWR0aCAxMzQyIC9TdGVtViAwIC9UeXBlIC9Gb250RGVzY3JpcHRvciAvWEhlaWdodCAwID4+CmVuZG9iagoyMCAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzQyIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjMgNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjEyIDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTIgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwNQo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTgyIDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoyMyAwIG9iago8PCAvRCAyNCAwIFIgL0UgMjUgMCBSIC9QIDI3IDAgUiAvUyAyOSAwIFIgL2EgMzAgMCBSIC9iYXIgMzEgMCBSIC9kIDMyIDAgUgovZSAzMyAwIFIgL2ZpdmUgMzQgMCBSIC9mb3VyIDM1IDAgUiAvaSAzNiAwIFIgL20gMzcgMCBSIC9uIDM4IDAgUgovb25lIDM5IDAgUiAvcGFyZW5sZWZ0IDQwIDAgUiAvcGFyZW5yaWdodCA0MSAwIFIgL3BlcmlvZCA0MiAwIFIgL3IgNDMgMCBSCi9zIDQ0IDAgUiAvc3BhY2UgNDUgMCBSIC90IDQ2IDAgUiAvdGhyZWUgNDcgMCBSIC90d28gNDggMCBSIC91IDQ5IDAgUgovemVybyA1MSAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YxIDIyIDAgUiAvRjIgMTUgMCBSID4+CmVuZG9iago0IDAgb2JqCjw8IC9BMSA8PCAvQ0EgMCAvVHlwZSAvRXh0R1N0YXRlIC9jYSAxID4+Ci9BMiA8PCAvQ0EgMSAvVHlwZSAvRXh0R1N0YXRlIC9jYSAxID4+Ci9BMyA8PCAvQ0EgMC44IC9UeXBlIC9FeHRHU3RhdGUgL2NhIDAuOCA+PiA+PgplbmRvYmoKNSAwIG9iago8PCA+PgplbmRvYmoKNiAwIG9iago8PCA+PgplbmRvYmoKNyAwIG9iago8PCAvRjEtRGVqYVZ1U2Fucy1PbWVnYSAyNiAwIFIgL0YxLURlamFWdVNhbnMtUGhpIDI4IDAgUgovRjEtRGVqYVZ1U2Fucy11bmkwMzAyIDUwIDAgUiAvTTAgMTIgMCBSID4+CmVuZG9iagoxMiAwIG9iago8PCAvQkJveCBbIC04IC04IDggOCBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTMxIC9TdWJ0eXBlIC9Gb3JtCi9UeXBlIC9YT2JqZWN0ID4+CnN0cmVhbQp4nG2QQQ6EIAxF9z1FL/BJS0Vl69JruJlM4v23A3FATN000L48flH+kvBOpcD4JAlLTrPketOQ0rpMjBjm1bIox6BRLdbOdTioz9BwY3SLsRSm1NboeKOb6Tbekz/6sFkhRj8cDq+EexZDJlwpMQaH3wsv28P/EZ5e1MAfoo1+Y1pD/QplbmRzdHJlYW0KZW5kb2JqCjIgMCBvYmoKPDwgL0NvdW50IDEgL0tpZHMgWyAxMCAwIFIgXSAvVHlwZSAvUGFnZXMgPj4KZW5kb2JqCjUyIDAgb2JqCjw8IC9DcmVhdGlvbkRhdGUgKEQ6MjAyMTAxMTUxNDUyNDArMDInMDAnKQovQ3JlYXRvciAoTWF0cGxvdGxpYiB2My4zLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcpCi9Qcm9kdWNlciAoTWF0cGxvdGxpYiBwZGYgYmFja2VuZCB2My4zLjIpID4+CmVuZG9iagp4cmVmCjAgNTMKMDAwMDAwMDAwMCA2NTUzNSBmIAowMDAwMDAwMDE2IDAwMDAwIG4gCjAwMDAwMTcwMjMgMDAwMDAgbiAKMDAwMDAxNjQyNiAwMDAwMCBuIAowMDAwMDE2NDY5IDAwMDAwIG4gCjAwMDAwMTY2MTEgMDAwMDAgbiAKMDAwMDAxNjYzMiAwMDAwMCBuIAowMDAwMDE2NjUzIDAwMDAwIG4gCjAwMDAwMDAwNjUgMDAwMDAgbiAKMDAwMDAwMDM5NyAwMDAwMCBuIAowMDAwMDAwMjA4IDAwMDAwIG4gCjAwMDAwMDQ3NTggMDAwMDAgbiAKMDAwMDAxNjc2OSAwMDAwMCBuIAowMDAwMDA2MTA2IDAwMDAwIG4gCjAwMDAwMDU4OTggMDAwMDAgbiAKMDAwMDAwNTU2OCAwMDAwMCBuIAowMDAwMDA3MTU5IDAwMDAwIG4gCjAwMDAwMDQ3NzkgMDAwMDAgbiAKMDAwMDAwNTE3OSAwMDAwMCBuIAowMDAwMDA1NDA0IDAwMDAwIG4gCjAwMDAwMTUwNTYgMDAwMDAgbiAKMDAwMDAxNDg1NiAwMDAwMCBuIAowMDAwMDE0Mzk5IDAwMDAwIG4gCjAwMDAwMTYxMDkgMDAwMDAgbiAKMDAwMDAwNzIxMSAwMDAwMCBuIAowMDAwMDA3NDQ0IDAwMDAwIG4gCjAwMDAwMDc1OTUgMDAwMDAgbiAKMDAwMDAwNzk2NiAwMDAwMCBuIAowMDAwMDA4MjA0IDAwMDAwIG4gCjAwMDAwMDg1NzYgMDAwMDAgbiAKMDAwMDAwODk4NyAwMDAwMCBuIAowMDAwMDA5MzY0IDAwMDAwIG4gCjAwMDAwMDk0OTAgMDAwMDAgbiAKMDAwMDAwOTc5MCAwMDAwMCBuIAowMDAwMDEwMTA4IDAwMDAwIG4gCjAwMDAwMTA0MjggMDAwMDAgbiAKMDAwMDAxMDU5MCAwMDAwMCBuIAowMDAwMDEwNzMwIDAwMDAwIG4gCjAwMDAwMTEwNTggMDAwMDAgbiAKMDAwMDAxMTI5MiAwMDAwMCBuIAowMDAwMDExNDQ0IDAwMDAwIG4gCjAwMDAwMTE2NjQgMDAwMDAgbiAKMDAwMDAxMTg4NiAwMDAwMCBuIAowMDAwMDEyMDA3IDAwMDAwIG4gCjAwMDAwMTIyMzcgMDAwMDAgbiAKMDAwMDAxMjY0MiAwMDAwMCBuIAowMDAwMDEyNzMxIDAwMDAwIG4gCjAwMDAwMTI5MzUgMDAwMDAgbiAKMDAwMDAxMzM0NiAwMDAwMCBuIAowMDAwMDEzNjY3IDAwMDAwIG4gCjAwMDAwMTM5MTEgMDAwMDAgbiAKMDAwMDAxNDExNiAwMDAwMCBuIAowMDAwMDE3MDgzIDAwMDAwIG4gCnRyYWlsZXIKPDwgL0luZm8gNTIgMCBSIC9Sb290IDEgMCBSIC9TaXplIDUzID4+CnN0YXJ0eHJlZgoxNzI0MAolJUVPRgo=\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",
"\n",
"N = 128 # number of samples\n",
"\n",
"# generate random signal\n",
"np.random.seed(5)\n",
"x = np.random.normal(size=N)\n",
"\n",
"# compute magnitude of the periodogram\n",
"x = np.concatenate((x, np.zeros_like(x)))\n",
"X = np.fft.rfft(x)\n",
"Om = np.linspace(0, np.pi, len(X))\n",
"Pxx = 1 / N * abs(X) ** 2\n",
"\n",
"# plot results\n",
"plt.figure(figsize=(10, 4))\n",
"plt.stem(Om, Pxx, \"C0\", label=r\"$|\\hat{\\Phi}_{xx}(e^{j \\Omega})|$\")\n",
"plt.plot(Om, np.ones_like(Pxx), \"C1\", label=r\"$\\Phi_{xx}(e^{j \\Omega})$\")\n",
"plt.title(\"Estimated and true PSD\")\n",
"plt.xlabel(r\"$\\Omega$\")\n",
"plt.axis([0, np.pi, 0, 4])\n",
"plt.legend()\n",
"\n",
"# compute bias/variance of the periodogram\n",
"print(\"Bias of the periodogram: \\t {0:1.4f}\".format(np.mean(Pxx - 1)))\n",
"print(\"Variance of the periodogram: \\t {0:1.4f}\".format(np.var(Pxx)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise**\n",
"\n",
"* What do you have to change to evaluate experimentally if the periodogram is a consistent estimator?\n",
"* Based on the results, is the periodogram a consistent estimator?\n",
"\n",
"Solution: The conditions for consistency have to be checked for the limiting case of a infinitely long signal ($N \\to \\infty$). Increasing the length `N` of the random signal in above example reveals that the periodogram can be assumed to be bias free, $b_{\\hat{\\Phi}_{xx}} = 0$. However, its variance $\\sigma^2_{\\hat{\\Phi}_{xx}}$ does not tend to zero. The reason for this is that by increasing the length $N$ of the random signal also the number of spectral coefficients is increased by the same amount."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Evaluation\n",
"\n",
"From above numerical example it should have become clear that the periodogram is no consistent estimator for the PSD $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. It can be shown that the estimator is asymptotically bias free for $N \\to \\infty$, hence\n",
"\n",
"\\begin{equation}\n",
"\\lim_{N \\to \\infty} E\\{ \\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\} = \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\n",
"\\end{equation}\n",
"\n",
"This is due to the [leakage effect](../spectral_analysis_deterministic_signals/leakage_effect.ipynb) which limits the spectral resolution for signals of finite length.\n",
"\n",
"The variance of the estimator does not converge towards zero\n",
"\n",
"\\begin{equation}\n",
"\\lim_{N \\to \\infty} \\sigma^2_{\\hat{\\Phi}_{xx}} \\neq 0\n",
"\\end{equation}\n",
"\n",
"This is due to the fact that by increasing $N$ also the number of independent frequencies $\\Omega = \\frac{2 \\pi}{N} \\mu$ for $\\mu = 0,1,\\dots,N-1$ increases.\n",
"\n",
"The periodogram is the basis for a variety of advanced estimation techniques for the PSD. These techniques rely on averaging or smoothing of (overlapping) periodograms."
]
},
{
"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
}