{ "cells": [ { "cell_type": "markdown", "metadata": { "nbsphinx": "hidden" }, "source": [ "# Spectral Estimation for Random Signals\n", "\n", "*This jupyter/Python notebook is part of a [collection of notebooks](../index.ipynb) in the masters module [Digital Signal Processing](http://www.int.uni-rostock.de/Digitale-Signalverarbeitung.48.0.html), Comunications Engineering, Universität Rostock. Please direct questions and suggestions to .*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parametric Methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Motivation\n", "\n", "Non-parametric methods for the estimation of the power spectral density (PSD), like the [periodogram](periodogram.ipynb) or [Welch's method](welch_method.ipynb), don't rely on a-priori information about the process generating the random signal. Often some a-priori information is available that can be used to formulate a parametric model of the random process. The goal is then to estimate these parameters in order to characterize the random signal. Such techniques are known as *[parametric methods](https://en.wikipedia.org/wiki/Spectral_density_estimation#Parametric_estimation)* or *model-based methods*. The incorporation of a-priori knowledge can improve the estimation of the PSD significantly, as long as the underlying model is a valid description of the random process. The parametric model of the random process can also be used to generate random signals with a desired PSD." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Process Models\n", "\n", "For the remainder we assume wide-sense stationary real-valued random processes. For many applications the process can be modeled by a linear-time invariant (LTI) system\n", "\n", "![Synthesis process model](synthesis_model.png)\n", "\n", "where $n[k]$ is [white noise](../random_signals/white_noise.ipynb) and $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ denotes the transfer function of the system. In general, the random signal $x[k]$ will be correlated as a result of the processing of the uncorrelated input signal $n[k]$ by the system $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. Due to the white noise assumption $\\Phi_{nn}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = N_0$, the PSD of the random process is given as\n", "\n", "\\begin{equation}\n", "\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = N_0 \\cdot | H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2\n", "\\end{equation}\n", "\n", "Parametric methods model the system $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ by a limited number of parameters. These parameters are then estimated from $x[k]$, providing an estimate $\\hat{H}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the transfer function. This estimate is then used to calculate the desired estimate $\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the PSD." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Autoregressive model\n", "\n", "The [autoregressive](https://en.wikipedia.org/wiki/Autoregressive_model) (AR) model assumes a recursive system with a direct path. Its output relation is given as\n", "\n", "\\begin{equation}\n", "x[k] = \\sum_{n=1}^{N} a_n \\cdot x[k-n] + n[k]\n", "\\end{equation}\n", "\n", "where $a_n$ denote the coefficients of the recursive path and $N$ the order of the model. Its system function $H(z)$ is derived by $z$-transformation of the output relation\n", "\n", "\\begin{equation}\n", "H(z) = \\frac{1}{1 - \\sum_{n=1}^{N} a_n z^{-n}}\n", "\\end{equation}\n", "\n", "Hence, the AR model is a pole-only model of the system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Moving average model\n", "\n", "The [moving average](https://en.wikipedia.org/wiki/Moving-average_model) (MA) model assumes a non-recursive system. The output relation is given as\n", "\n", "\\begin{equation}\n", "x[k] = \\sum_{m=0}^{M-1} b_m \\cdot n[k-m] = h[k] * n[k]\n", "\\end{equation}\n", "\n", "with the impulse response of the system $h[k] = [ b_0, b_1, \\dots, b_{M-1} ]$. The MA model is a finite impulse response (FIR) model of the random process. Its system function is given as\n", "\n", "\\begin{equation}\n", "H(z) = \\mathcal{Z} \\{ h[k] \\} = \\sum_{m=0}^{M-1} b_m \\; z^{-m}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Autoregressive moving average model\n", "\n", "The [autoregressive moving average](https://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model) (ARMA) model is a combination of the AR and MA model. It constitutes a general linear process model. Its output relation is given as\n", "\n", "\\begin{equation}\n", "x[k] = \\sum_{n=1}^{N} a_n \\cdot x[k-n] + \\sum_{m=0}^{M-1} b_m \\cdot n[k-m]\n", "\\end{equation}\n", "\n", "Its system function reads\n", "\n", "\\begin{equation}\n", "H(z) = \\frac{\\sum_{m=0}^{M-1} b_m \\; z^{-m}}{1 - \\sum_{n=1}^{N} a_n z^{-n}}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parametric Spectral Estimation\n", "\n", "The models above describe the synthesis of the samples $x[k]$ from the white noise $n[k]$. For spectral estimation only the random signal $x[k]$ is known and we are aiming at estimating the parameters of the model. This can be achieved by determining an analyzing system $G(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ such to decorrelate the signal $x[k]$\n", "\n", "![Analysis process model](analysis_model.png)\n", "\n", "where $e[k]$ should be white noise. Due to its desired operation, the filter $G(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ is also denoted as *whitening filter*. The optimal filter $G(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ is given by the inverse system $\\frac{1}{H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})}$. However, $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ is in general not known. But this nevertheless implies that our linear process model of $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ also applies to $G(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. Various techniques have been developed to estimate the parameters of the filter $G(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ such that $e[k]$ becomes decorrelated. For instance, by expressing the auto-correlation function (ACF) $\\varphi_{xx}[\\kappa]$ in terms of the model parameters and solving with respect to these. The underlying set of equations are known as [Yule-Walker equations](https://en.wikipedia.org/wiki/Autoregressive_model#Yule-Walker_equations).\n", "\n", "Once the model parameters have been estimated, these can be used to calculate an estimate $\\hat{G}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the analysis system. The desired estimate of the PSD is then given as\n", "\n", "\\begin{equation}\n", "\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})}{|\\hat{G}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})|^2}\n", "\\end{equation}\n", "\n", "where if $e[k]$ is white noise, $\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = N_0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example\n", "\n", "In the following example $n[k]$ is drawn from normal distributed white noise with $N_0 = 1$. The Yule-Walker equations are used to estimate the parameters of an AR model of $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. The implementation provided by `statsmodels.api.regression.yule_walker` returns the estimated AR coefficients of the system $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$. These parameters are then used to numerically evaluate the estimated transfer function, resulting in $\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = 1 \\cdot \\vert \\hat{H}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\vert^2$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bias of the Welch estimate: \t\t 2.8932\n", "Variance of the Welch estimate: \t 62.0748\n", "Bias of the parametric estimate: \t -0.0291\n", "Variance of the parametric estimate: \t 15.3476\n" ] }, { "data": { "application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDYwMS42NjU2MjUgMzI1Ljc0MDYyNSBdIC9QYXJlbnQgMiAwIFIgL1Jlc291cmNlcyA4IDAgUgovVHlwZSAvUGFnZSA+PgplbmRvYmoKOSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDExIDAgUiA+PgpzdHJlYW0KeJztnEuTZMdthff1K+5ytGBNvh9LMmgxwhubEsPaeMMYj0hOzNCmaJn++f4OcG9V1nT3NMu0eqEwFQx2o7LzgQRwDoAsxe3dKW7fbWF7x7+/bHH7in+/OwV++3BqIZ5bqy1Vfn2//ppTPfcS9ON7Bt/8+v3p9OfT68+Z5mf+7KvTKbdz2f+sn2vJvTJ3neWB9P2NNBbmdPEywyK1hX7aHpm+1rGl3s9p+8vb7U/bj9vrz5O2k7Z/5Hgc9Tw4bDi3GeJsoY/CPl9/+fa/fnjz9g9ffbG9+fnxTT++j9MfT19vPx1rBNS5rIFCkTw7XTmHY7YvuJFfTl98s73+fdxi2L758ymNc40MS+k86/bNv51ehXP43fbNu+0fvrG1/090cPpqu9UBS54T9zryuutVeqcWnp/w03qIsZ9zCrOOVGq+UUf926sjcdowUi1x3f0qvVMdz0/4jFmEdo5pltjKDG1RR3wJ68ghnce82fpFdK93fGKqZ1RAwJi5j1561pEWFbyAReRJRJyzMNO69UV6ryKenfDT6sgjs17LOc0R46KO9BIWUQZ+GPqMbd39Kr1THc9P+Gl1lJ7OtRaGptTKjTpewDpqm+ccW+7lBvgW6Z3qeH7CT6ujtnjOLacyS499UUe+WsdP/FnYPgubvDgeapkE6CEEevOBKU8ax96ChobLD+fmA1jwsy/fvvv2X/76x29//Pmzf/rw9rtvty///fS1/ve388VfzzL+94CNxh7O9hhQZexu6vNsDjn6C8L2dVbMP84HmrhI79XEMt1VE+tsn9REG+c8L5p4EQ+8zhqByvHQKK7ie3WxTnhVxs18n9RGDDjd1TDqyxpGLLj2Q8u4iu9WxzLhoo51vk+rg3O2q3X0F7aODpF6xDou4rvVsUy4qGOd7xF16DQRf8pXu4gvHDESC/aHhnEV36uJdcKrJm7me0ITYrfzahLxhSMGGca5PLSJq/huTSwTLppY53tKE7hRW2zihYNFgm+FR2ziIr5bE8uEiybW+Z7SBB6UFpt44ThxJWwfnipV3MktlgkXcrHO94QmMh40rzaRfnucIGMyTcQUaxst9XoKTLL/0/mkhzpGyjXN7Q8fK2k51OQgY6fqHNwlo4ZGRvL+hL1Hh0ikpbYakqRz3yWpzohtkN2/P9V8niGLdk4Oy/7TkHSck5YZk09LahrZynkEtp23MQSrc2qpDrwkJkBaz6mEZkvNfI5aqY9z76k0jYyBZZ3OYGBp9DFMKspSNUEv55rncHFUzmm7ZYVQZ7TdorczWWjJY+uR6Aaw2dSJHekMHQWSi0TdaoSp12obhr1bAtslxtHDtDqaSH1BbDMXwNHMXiRqDkbYJLWco1kA0tHYUpG06aLE1jVHCLH7Nho37/aCGF3EbmfhBH2w1aH9xRiJyiYeIBUnT7btNpOraRCEjbXoiHhA7dPEsJ6m+p/YIE7Uiu2Pu0xJipKeRs5gm4nrWVepues5hzh6MzGKKClq7o5WU0jaSQpFY6rmnlxc9AtLXE1rQfseBMcQ0tDhE4wi7FeuGyu5JRNz/xu7TKEGs8KkSlJHpcMsJs05u4nnwUgwuZAq/0jM1fasi5HRAl6h2eiEgUaplFjQyLyqSxuXKBPfJjrPHGSaeBxJFH+PqhliO8l4WeWDbJxQiXuw8cSD5kyxybijz4LhjNilpUgEwC9btIPnC1iaeRZs13aeUVqXXyKf59axB9sklpadlCSGB+7ANlNwmWzmENkXxtqtpiuU6tHHE1CwPFMrdjmmWUks5ZwGDmKTV3SnnVRyg9HnaCZEc7VI4zJRlgxtmJw46t4UsbWC35hLp9rOnC3wSwQZJsvv84ydOkTughy1uhobbl0tWOhu8bpgW8HsvXRaiEGjuHknXACblha1/EgkxKZdrHp0u9HUiFOV9W1y3D9UM8LUFd1qbqZ1RQtfs/s6rhe8hwuQFvWjXbUNV+HC3VszEgx6telxq94V6ViWe2wo1eVFGTpWiIy4YFscOpwZro6UsZZmasHZiiJD5vzcLRdiU+Bso5gxSotp5jZsi1NW79bIkkGWJs1kTDN59NDt4RbZgq7KcdmsC8Phnrr5bMY5a3NrxFwI5uzH5O3IehUVUUBOPp7wlXZr5MyFT03OjmPwBIqQoLuzsj/enBSc8CQF61pTMnF9ADWSEiSLBWtUNBouOk08DYJGZt6cTYc5JQKjh7GJTogDNgOesENTl7kDIj66nfMlulV+HS4eO2TJmIaFtiyf8YiH3WE802JYFrQ7kAmoag8WH3OuKvTK8xvGGbJjQ8bBDODkJTWrrCkhQdujY2XLg/hnihMzTGawtStOx2Z7w12nXS6mE8vYd4xbZu/S6Gp7dtDJeACzRnm1r8UKwYMgsMz0CKwJk7FqLVWGqvDJB3PQ7mGNaMC25pSpZ+w7e5xiLwOdFdsvNug4S7zohB4DITbiOJs24g+KmBbT8iRdiqZ4OS/R1Vy6YHQjm9oQC2zsUgtW0U3HCm0FJJPqC7c0nFXCRbD9ZJG7lHagIWJoDQYqcRt7AS2Lj+Dyzk3agRWKpoTOHo2d5CNwZRZJYxeDR+YpCJtkJhRsOMxK3Kur9AnxziF/HYfLj3C4B92shZ1hJznkgHagnCxKwOUSjAB1kLOYhoXW3IdhBgZo2MjVAB7EEsFazhDB4XQuhWhgJ+uqrnUGECOnsJVoIgd2iieSkuXMk5+61sLuODYBjS0ksYhplyzql2Fn03AYtpJNX+p3RExP5A2ccL8FELKiM07O2GwQhO8QDQFBEcoeW7aoh5u1TnCKIiTwkG6hc/dJbE6arUIRSSEhWBbmwl+Vmee0GXA+dmW0CIocDTVBjRKAviLSUQqIV3f2OfHUbgwDbRg2KLL01o129NBm1Q5YV7lHNqLKrmKWlNMQxaHiRnEaXqGTgRUYP78ZgeUmzWTEdrB93xj4W0xfMB7piGgzFIEgbLodVh4jlDksnoXaii4SPYVWCsBjoIa3dSfLiSsPdjSCY7fVAJVCNIXuIGW3ziUxi4am48bHBMFoVgOrG7Czln2C7v4ncwMUFf1RSMjs0jhgMEarIMbWa4LXuRgHrFy8ccDCiFB2eg4QB+fyORtxVWEhiQ6KFJPRGeOMFn+JiKKtEM7cfKziL9ZZ5BQY5DSuoyYjNysqo8ie7AYkZmn8ryiGK5cwwil44y5EK0SwZ69OfWGGOAgYKDFaGeYsUdSQiDXt8rjEaEalFAGy3Kox4lKwFBsNmuD0LVRNQubgg7Mlh5W4Is4xULMdRgQQVoFHoDqYZNoTjSlWT9zWtokA2TVKjIQagDiCk6AswqRENI6S7IyjzGnsh0RA5kfcEYqCt+Qjm2clGAQI14U2ffhgcT/wfKAQPA8f8TxD7A8+hWE3xUV0bZkD2INngRrySGTN+ABMi5PBLshhpBDuxsVYfEN7VT5JdCjpyHhmDRMzZkko/rB4E8UIc9EZuPUsAuliGAa3W4smmZFIZnOLKAZjjYoSDce1DRq4w7J1BT2Trdpg44mhKzHsooPFr9doolJmKYolqivKHGq0aLeOVSSDFeVYJSUYtp0dUmLsRjQN8lKT3Tq+0CwmK/Xq5LvdDRvd+JWJIvZS8Lom6t6KW7xYInqQtRJwZk/7aGWVOLflb4Do7gj4Ps7GgZHiqtN9FN/vBGUFbOPz2MPm6dsgv1DEFhEj+00+nCgHEVDIxrlBbU/fgshM9JCtWnB3cdK1Ens2wLXU0C0GKdnDQgOAPsSnkRu3VQaEJYuBsVOuNPtocATXi2byFQ2ZqYn+xkFOOjQ3lODIAWGAjcQz6jjsOVuUVmrIjF3BW+iys0nlEMTjniyJI66m4mKSVO5zCMPIUNLw1I48ByfKRj0x+GZMgL8RuittkQXAPaplAUlYPHd47AFrsEMSBODPwhaDLKKYJTBZYMwNixvrGQj0yIZnZbooZdg1VIVQO2YWIEMkAO80BSAl+nhh8sRuLFcLwXcuntMzRmijAb/hiU0RLHfpKprHGxNWdsjNc9BtKlJAbz2XlOIISFiEShLY4z6JgBlLmKRqQvHiB63K8kgoSAKwWVwtZM8ahc5tKmloyvFq8VRFGRrpDzERaMbtExdvciF0GUT1FKxASn5gcmH0GLiN7nyqTGC7JziQcHQlffJ8YCl4MiigbtaPZcZW3ebkN6SoBPY4LUM9UkdhNQdU6pgFD8MzZ5VP+NnmrujPMzhV81Fv1VsZyYdFZcv+ogUTGcnw3r4njiSpxRLBISsp1ZM+gTOwpeQWI5DKbJOEGdRGOGOeqdQ1ejVgCKDFTzb9CKwmt+qhQhKaRzVFjtmL5S2YMShdVBpQLMgZTLcrsdsX07ZHKU1E2cTcMvZJjJdJwCuG34iIMniVZQjRg++RbGLukz8WMCofNrmSTfBDxEiwgWtGT1SC8LqTBW/CL0zLqxBwnXOYhDxVJ4qIgJelYLMgdhUAqZpRSEmrywXZwmRzeW6zHdlmg7EG820RN0vSoiBbjFDVOyKROWuOIDaWrRxNzRvCu2UsUYhdq0YrEPCzzZGE2CTvU9VJXKw3F3P1hHFFeqyHY5prK9/kAkryWuTUfJvnm5PdEMwAVg6arE6S5c04AJFexBBsMVVlAflU9QskyJzFV1RuovhuaEHiPPcsUkgOkYgGzlz7PDLOgf+WojRQNV5LWrJSI5JGeFxRtQDHtdFFUA6bFGvCharVNnIRK8pyqSqiD/S5WEhuNBf3hk6XPW9VPm8ljSqTx4Q9FRWSMyRtMICJ/ZuZVAH5kOtoewC6wX6uKhUkBWMiBNEvZRcLyJtiF3Mwdc0+iUo+QW7uyW/28kATkgOaZMB4M796xt6UH/KnWRkp8OkKkU8mOwKnJZRki39ZuUA0yi0lwFqjKbsLsZEmqc9Cgq2ICWC4+ITSV35oRtUzLl4xfOIzk+g9U3KxQLiIumgjIWajnWRXRA098dFltOgxKE/5NMs35cAgW/IVVWgHvInDiKEle9FFRcvMh5Z1q7qyTyKw7ZYaK5vBA00j04rJYkUq9RH3gqeKAlsySMvRUVvvnqMDSWASKbOWnNWBsgShKhl7VNo9evdSJ2mWvIL0xPJuWPH0jF6waiRKwUprVxNbqdHy/AAwJBc2K77jhuo572VsaC+YqpepmoEPZ/aiAIsXK2tYgwJstE0nQaqXIVhPsGgNCz1PDCpbaWpIeDOeRkiFVEKz7YhcooFbSQLUrlKAtpexQ1NTFheOMlzEeEf0RDsLUGE53dSEMZt3ELpV48OVitqk0CfPv4FTxMPEhK3sabVONo0hSwz4GvkoMi7CMnrPorA7sy4yLiyYkKv94aZ+MUVYaomryiRYp+ECVqMCtsiZJkEFVtwhSVTwJWWWGJzZGztZhXimlLRDeKetWAWkcL1uc0wFGombghgazhJP0TUrJ5B41CTir4IXzmuNnKK0DeDvdpiCIu3OezYOlpMvWJ3kFsuhMFWv4rCffQ5hH7dutZYuV3KxoG+EaWfpIVRXH+5IWIVf2IqwVr9HJbHkwIRE1XyYwjYyhHuKpTbzXqaEmahqo2vMKqhV78+AiUTePqrfF2zGokVRKRP8C3ZEJkvWeShTkAclSbYgKzafRIinlEajVdRyL5jSDl7oZ1RJ0WosRo+Yw85IOlddqqYMkcWOqDBuoY/YCeJVkXoNDp6/1SDAU4rFxBU8MepUybaJUJixjSVBt1hbRXWaLECjM5nS9NECtpqm3ReJV7aAA1qKEohg6yzcmxkZUUqFKriCm8IIho6wHX7poleuEM9zsSTII8BabIPEtuiTCNnALbsvWL6lb8RZXHOqE2TS4I0bgAJVQmJsDhCiWhgijAl0Rvejg9OmEA3BMdCZfuL6zHVrEYINlcd16WSIuxgII2IQPW3q4MUzQhs8MLS8VwtL9/tS2VWU0LRNZtFd21VYNexhsPJIAoDtr4qa6VHf5vW82F1sJU4825aUUboYDMuuJpTaXNiEYCVXOwxZtZfxqrJs4k1orqfibLTio1PZjB+megeqquINmY92uwHrdO2pkYkBNq+Soj6/dCVyqtLaEbHI6pOoHYkTdL9d2ITFIfCFgAeT6/ulD0PBqo5Gm9OvAPPsXg5UXxm+sItzcwtWqQpSlZtbSK7Dy7VNxLMPP0wYu1Roh2uY25VmvqN6rdCOYGuHUXZoNa5KNpR7bh4mydjyPvpR8V3V3Uc69FFt+akqeRItUPm6+OOF8em2PL5KXB9LW16oCle2qHlpyysXkIYNF462vAoHBJPopn90M4hw4G8z8q+KqxoUQzkwXutEVO0saw6osKnuq4SqgapCEtRFuGnJw45h4FalXFryDIc5Wjy8bclrme41u2tLPuEkOQcvrN205NVmtRLG0pIfonlWYLm25K0TXrw0ctORBxi9abt25GXfIl3bTUceKSmh9WBvO/JEo1aPitPeg+i6xDC8KrQ25NXG8jmWhrx13nv4uCFvP0Wvgi4N+SxuOiyurA159YjJIo+KztGQb6JfTtbWhryCwp76LA15Vaqavay+bcir1bZzKj30Gt5zGmrbFIfXtU+PgQ5R7e3Spx8qmxQvOVz79LJWubWLr316vD30vfN+bdMXVQp78crKpU3fDeeirXVt01s/fmJ3Rx3m0qe3pmEI3mFf+vTN6pBenzna9FEN0LZXldY2vbaSvcF006bXM4y+b2Zp0yfZYQiemF/b9Nnq+vVBl54bIz2ZLr926TVEcc2OtPTpq5pjo7XbPr3eoADYD/r02iJJrw9e+vQwpW4dk9s+vfUQYvNm96VPzz1bstw/7tMnPVyHNLj8aNQrrSMSZtPj0qgnlHTQyGtIS6O+N3Nmr/EsjXpcA19M0wsrl0a9OiVq7Jjal0a9PGmoOLJ91Khvuqbp13TbqC+qUdj0R6NeAd/Beu3TC62nc9Y0VfGy+1ecTTvArW36bqY7XYlrm16t0x497V/a9MV6bMET6EubHrsQudvrINc2PVbfbL7bLr2a5+LmllwuXfog/VsZ9LZLj9GTSngp5Nqkb0oLHA7WJn0Qs3XqtXbpmxpspXpLf+/SqyoBKngd5Nqlx36HIG+77dJjkdxAPKomR5eeoK8iu/f0L116cQHyEK+DXPv0Vplt/gDgpk8P7fae/tqlR38xf9ylz2XWfey1Sy9SaMmaqh3Zn6elzZ5yFCenKoJEhUyvPbTU9zZ98mdrSeUO+16c1zW4UY+EVY/0srfO1fHwkKeSClzVIiH/wYCHgE8Uk5zLyxp6rG+98KbC0iFU69muu6iOQv5gU+jNXDx69UHdBK8a5ANpVfrjfjzvlUW4W1t5iJNbnmIvvO3kSsvht56nyArdOZThtrwndFae8J68rmp6RhgPiFP+2LzmXFTJdN1ZAjSd2Kv+5pwKE8ZWi9PYfN5zMx3amONHjffevOv8uPj2TafY2PU952/51uhj30751MPS+77icu+z1defF3+yKrKpZ6u/2DH9C8Hww3B5vlvTdTru4YFY/FnBf6jqeBWLVF/G6i2JKs630mNDb06LWOXd68SLWAlWVxV6WU7SY7m+Dj0mXg5yFXupIxJylJxexES4R1a7kfbrxFfx9SDvb8WXUy/LXRX0qI7f6HvOX5yefay7/crHuid9UY/Zj4WmmKoh6Qf7vtxDuVmHzCIcZrG/Y16+K6bv5F3+Ug9+7C/922JqmLTurOezoGJQ9lfLr39PsNaD5lf/rbfLp6IcuTr3Cvas+ZADK2EnNufmdNL+/ngR/epfX9lIL5E2ZwbhWGMf89aGiK2AQCTPqlrYu5300W7e2cDuNQz//LPbkdd1f2djn/hCnOcvD78S98/f/+BfiLM/Isbsf8VPSQ98B8CoVu5WVBTI/dmv1a3PgI4nQB9/OXO5cz1ReuzOr/I9IqQjIvzKqxdQP3n1ylv+7m7+4bqvtv0DFdGjvz3zD57c6qs/2QdKvHubDqP+wbHzWJODrovfm1hkdOQUlvFv7INqBQKndi7/fp8H2jKc2O5betJ+O+l5SH1/4PaIAf71xx80YjHkj6x/38BvMH9U/7ztX40+bndURlZ/aHrllB+4w0V8bwSsF7j5fy+4wwv+4zDS0o//1xGJvz18IPW5f/nHP/jLobUYc1q09uQffNjdzB4er3t6ws3+89ip3rWO5YOnVv7hXr/8O/C/r0//A1zhT00KZW5kc3RyZWFtCmVuZG9iagoxMSAwIG9iago2MTEwCmVuZG9iagoxNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMyNyA+PgpzdHJlYW0KeJw9kkuOxDAIRPc5BRdoyfz8OU9Gs+q5/7YfdnoWUaEARVG4a0qT6PJSlxwq2Zv86BVTxc3k74mGxDRxN9DBKfe1owj+BLgOZu4M0UqJQb11icrykbHz19fTlwfvy9thsFF6lphXn5KxBgvydHTxAebc2u5/le8rk44svUpPRe9LvViWqA8x+2L1KdzG5hvbejCeTPFrNlGrmVPW2iK6DNRFSjex1QtQ7VGB41xHbxVtAxCBO0EHlQm1amzBSEi0RdYvzGlC98T/7QG6X9rEdeyD2LSNZ3xFGkPQ3bhHWcPMMep8anPD2tPZJuuU4kwu17I4Vgk0SstYakvwrHpdFGODGWilEM3RamxlaHbo62V4P8gMbzuKsailSs+eN6fgXH3s1xR68DyLHcGbnNgt6OeQrEIGLTZh6ixdr/A55339fgAqPX1GCmVuZHN0cmVhbQplbmRvYmoKMTcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNTIgPj4Kc3RyZWFtCnicPU/LEUMhCLxbxTbAjCCC1vMyOZn+rwFNPLGyuB/xgQpiFpBUh7jDreLFhbWjq+NTpB20CouD5gaWvxKxgryB69bpGD1mm3iK58YYOmOGsHDMkQxZ2sUz7rYvRY4NgrthVqHpKZY+5HphGwdSerO0XaGmcC6mImOGcdtpn1sl0o7TL0r9mq4ivcPUkrWLxp99fwG89zZnCmVuZHN0cmVhbQplbmRvYmoKMTggMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA5MiA+PgpzdHJlYW0KeJw9jLENwDAIBHum+AUiYYxt2CdK5ezf5i0naeD0D9fSoDiscXZVNB84i3x4S/WEjcSUppVHU5zd2hYOK4MUu9gWFl5hEaTyapjxeVPVwJJSlOXN+n93PcerG7oKZW5kc3RyZWFtCmVuZG9iagoxNCAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMtT2JsaXF1ZSAvQ2hhclByb2NzIDE1IDAgUgovRW5jb2RpbmcgPDwgL0RpZmZlcmVuY2VzIFsgMTAxIC9lIDEwNiAvaiAxMjAgL3ggXSAvVHlwZSAvRW5jb2RpbmcgPj4KL0ZpcnN0Q2hhciAwIC9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnREZXNjcmlwdG9yIDEzIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL1N1YnR5cGUgL1R5cGUzIC9UeXBlIC9Gb250IC9XaWR0aHMgMTIgMCBSID4+CmVuZG9iagoxMyAwIG9iago8PCAvQXNjZW50IDkyOSAvQ2FwSGVpZ2h0IDAgL0Rlc2NlbnQgLTIzNiAvRmxhZ3MgOTYKL0ZvbnRCQm94IFsgLTEwMTYgLTM1MSAxNjYwIDEwNjggXSAvRm9udE5hbWUgL0RlamFWdVNhbnMtT2JsaXF1ZQovSXRhbGljQW5nbGUgMCAvTWF4V2lkdGggMTM1MCAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTIgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTUgMCBvYmoKPDwgL2UgMTYgMCBSIC9qIDE3IDAgUiAveCAxOCAwIFIgPj4KZW5kb2JqCjIzIDAgb2JqCjw8IC9CQm94IFsgLTEwMjEgLTQ2MyAxNzk0IDEyMzMgXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIzNwovU3VidHlwZSAvRm9ybSAvVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJw9UbtxxTAM6z0FRuBH/Gied5cq2b8NKDkpeIApEQTkpyzRhZ9niOD7We7/yAOSrVBthCc0FZEN08DnSRFYbqQm3F7c54RslSP24lwgwhDtPAlppAsWOxkL3hc/j6seZqy5Yfy+M5p9VHTVUR28ew7jZk0/TpTd682sjlub+3TvrhOHa0gmn/cfnJRKp5csgzpLuLA2mhrW47woxljMOP4nqrBNsrajCsHSJUgq0IAYShLGgMUt/iInWg4L2psbaeudyU6qNIqGF6MM3qD1RjiKdJF8mGsrg7GpmDa++eQlN+j7Z7+fr18Da1rrCmVuZHN0cmVhbQplbmRvYmoKMjQgMCBvYmoKPDwgL0JCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjM4Ci9TdWJ0eXBlIC9Gb3JtIC9UeXBlIC9YT2JqZWN0ID4+CnN0cmVhbQp4nEWQTbLDMAiD9z4FRzCYH/s8nXmr9P7bSiR9XWRgDEj6MpaHRB55DzvsTCxQ8WrmElryGrpVvI5oTnHbXVc6J+zmFG6Yn74xVUyoonsJVTVD6KNxMPk6XsN99dt7+LHeimV9FfnVyZmtnFbt1LW92TENN5iON3deqpCAqiSiDx1f/45Xp6giNmw33Egxwejgi4XdDePJqGLV5bZF42GCaZj1emCMX3EQUvELoJfuDZyhD3DZF/jXcXqNMHQMuIGi+NZ5FAuR6FFL25L1zsDOKoUbzMibOzRViEFVULVhPeDzwZ4//Gv8fQCofmCACmVuZHN0cmVhbQplbmRvYmoKMjUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4OSA+PgpzdHJlYW0KeJw9jbsRwDAIQ3umYARjftY+uVTO/m1sH+cGPQQnYYAbq64BD84OfoTWvukjKZrUR7CILjJFee6tKJvWdVTKJFzK3PmTwnDU4By6v9XacbrE0Vv3/tOJHZkKZW5kc3RyZWFtCmVuZG9iagoyNiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMwNCA+PgpzdHJlYW0KeJw9kjuSwzAMQ3udghfIjPiT5PNkJ5X3/u0+MslWgEmJACgvdZmypjwgaSYJ/9Hh4WI75XfYns3MwLVELxPLKc+hK8TcRfmymY26sjrFqsMwnVv0qJyLhk2TmucqSxm3C57DtYnnln3EDzc0qAd1jUvCDd3VaFkKzXB1/zu9R9l3NTwXm1Tq1BePF1EV5vkhT6KH6UrifDwoIVx7MEYWEuRT0UCOs1yt8l5C9g63GrLCQWpJ57MnPNh1ek8ubhfNEA9kuVT4TlHs7dAzvuxKCT0StuFY7n07mrHpGps47H7vRtbKjK5oIX7IVyfrJWDcUyZFEmROtlhui9We7qEopnOGcxkg6tmKhlLmYlerfww7bywv2SzIlMwLMkanTZ44eMh+jZr0eZXneP0BbPNzOwplbmRzdHJlYW0KZW5kb2JqCjI3IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjMwID4+CnN0cmVhbQp4nDVRSW7DMAy86xXzgQDiLr/HQU/t/68d0glgYGhLnM0RGxsReInBz0HkxlvWjJr4m8ld8bs8FR4Jt4InUQRehnvZCS5vGJf9OMx88F5aOZMaTzIgF9n08ETIYJdA6MDsGtRhm2kn+oaEz45INRtZTl9L0EurEChP2X6nC0q0rerP7bMutO1rTzjZ7aknlU8gnluyApeNV0wWYxn0ROUuxfRBqrOFnoTyonwOsvmoIRJdopyBJwYHo0A7sOe2n4lXhaB1dZ+2jaEaKR1P/zY0NUki5BMlnNnSuFv4/p57/fwDplRTnwplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ1ID4+CnN0cmVhbQp4nEVQu41DMQzrPQUXCGD9LHued0iV2789SkZwhSFaP5JaEpiIwEsMsZRv4kdGQT0LvxeF4jPEzxeFQc6EpECc9RkQmXiG2kZu6HZwzrzDM4w5AhfFWnCm05n2XNjknAcnEM5tlPGMQrpJVBVxVJ9xTPGqss+N14GltWyz05HsIY2ES0klJpd+Uyr/tClbKujaRROwSOSBk0004Sw/Q5JizKCUUfcwtY70cbKRR3XQydmcOS2Z2e6n7Ux8D1gmmVHlKZ3nMj4nqfNcTn3usx3R5KKlVfuc/d6RlvIitduh1elXJVGZjdWnkLg8/4yf8f4DjqBZPgplbmRzdHJlYW0KZW5kb2JqCjI5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ3ID4+CnN0cmVhbQp4nE1Ru21EMQzr3xRc4ADra3meC1Jd9m9DyQiQwiChLymnJRb2xksM4QdbD77kkVVDfx4/MewzLD3J5NQ/5rnJVBS+FaqbmFAXYuH9aAS8FnQvIivKB9+PZQxzzvfgoxCXYCY0YKxvSSYX1bwzZMKJoY7DQZtUGHdNFCyuFc0zyO1WN7I6syBseCUT4sYARATZF5DNYKOMsZWQxXIeqAqSBVpg1+kbUYuCK5TWCXSi1sS6zOCr5/Z2N0Mv8uCounh9DOtLsMLopXssfK5CH8z0TDt3SSO98KYTEWYPBVKZnZGVOj1ifbdA/59lK/j7yc/z/QsVKFwqCmVuZHN0cmVhbQplbmRvYmoKMzAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjMgPj4Kc3RyZWFtCnicRZC5dQQxDENzVYESeIA66hk/R7P9pwtpvN5A+niEeIg9CcNyXcWF0Q0/3rbMNLyOMtyN9WXG+KixQE7QBxgiE1ejSfXtijNU6eHVYq6jolwvOiISzJLjq0AjfDqyx0Nb25l+Oq9/7CHvE/8qKuduYQEuqu5A+VIf8dSP2VHqmqGPKitrHmravwi7IpS2fVxOZZy6ewe0wmcrV/t9A6jnOoAKZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiFtCNEGUglgQpWYmZhBJOAMilwYAybQV5QplbmRzdHJlYW0KZW5kb2JqCjMyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNDUgPj4Kc3RyZWFtCnicMzK3UDBQsDQBEoYWJgrmZgYKKYZclhBWLhdMLAfMAtGWcAoingYAn30MtQplbmRzdHJlYW0KZW5kb2JqCjMzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjU1ID4+CnN0cmVhbQp4nEWRS5IDIAhE956CI4D85DyZmlVy/+00mEw2dpeo/YRKI6YSLOcUeTD9yPLNZLbptRyrnY0CiiIUzOQq9FiB1Z0p4sy1RLX1sTJy3Okdg+IN566cVLK4UcY6qjoVOKbnyvqq7vy4LMq+I4cyBWzWOQ42cOW2YYwTo81Wd4f7RJCnk6mj4naQbPiDk8a+ytUVuE42++olGAeCfqEJTPJNoHWGQOPmKXpyCfbxcbvzQLC3vAmkbAjkyBCMDkG7Tq5/cev83v86w53n2gxXjnfxO0xru+MvMcmKuYBF7hTU8z0XresMHe/JmWNy031D51ywy91Bps/8H+v3D1CKZogKZW5kc3RyZWFtCmVuZG9iagozNCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDgwID4+CnN0cmVhbQp4nEWMuw3AMAhEe6ZgBH4mZp8olbN/GyBK3HBPunu4OhIyU95hhocEngwshlPxBpmjYDW4RlKNneyjsG5fdYHmelOr9fcHKk92dnE9zcsZ9AplbmRzdHJlYW0KZW5kb2JqCjM1IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjM2ID4+CnN0cmVhbQp4nE1QS25EIQzbc4pc4EkkIQHOQ9VV5/7bscNU7SqGGH9ID+myVR7rU2J1iezypU2XyjJ5FajlT9v/UQwCbv/QyEG0t4ydYuYS1sXCJDzlNCMbJ9csH487TxtmhcbEjeOdLhlgnxYBNVuVzYE5bTo3QLqQGreqs95kUAwi6kLNB5MunKfRl4g5nqhgSncmtZAbXD7VoQNxWr0KuWOLk2/EHFmhwGHQTHHWXwHWqMmyWcggSYYhzn2je5QKjajKeSsVwg+ToRH1htWgBpW5haKp5ZL8HdoCMAW2jHXpDEqBqgDB3yqnfb8BJI1dUwplbmRzdHJlYW0KZW5kb2JqCjM2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTQ3ID4+CnN0cmVhbQp4nD1PuQ0DMQzrPQUXOMB6LFvzXJDqsn8bykZSCCJA8ZFlR8cKXGICk445Ei9pP/hpGoFYBjVH9ISKYVjgbpICD4MsSleeLV4MkdpCXUj41hDerUxkojyvETtwJxejBz5UG1keekA7RBVZrknDWNVWXWqdsAIcss7CdT3MqgTl0SdrKR9QVEK9dP+fe9r7CwBvL+sKZW5kc3RyZWFtCmVuZG9iagozNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE0OSA+PgpzdHJlYW0KeJw1j0sOAyEMQ/c5hS8wUn6EcB6qrqb33zZhWgkJC9svwRaDkYxLTGDsmGPhJVRPrT4kI4+6STkQqVA3BE9oTAwzbNIl8Mp03zKeW7ycVuqCTkjk6aw2GqKMZl7D0VPOCpv+y9wkamVGmQMy61S3E7KyYAXmBbU89zPuqFzohIedyrDoTjGi3GZGGn7/2/T+AnsyMGMKZW5kc3RyZWFtCmVuZG9iagozOCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDQ5ID4+CnN0cmVhbQp4nDM2tFAwUDA0MAeSRoZAlpGJQoohF0gAxMzlggnmgFkGQBqiOAeuJocrDQDG6A0mCmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNTcgPj4Kc3RyZWFtCnicRZC5EUMxCERzVUEJErAI6rHH0Xf/qRf5SrRvAC2HryVTqh8nIqbc12j0MHkOn00lVizYJraTGnIbFkFKMZh4TjGro7ehmYfU67ioqrh1ZpXTacvKxX/zaFczkz3CNeon8E3o+J88tKnoW6CvC5R9QLU4nUlQMX2vYoGjnHZ/IpwY4D4ZR5kpI3Fibgrs9xkAZr5XuMbjBd0BN3kKZW5kc3RyZWFtCmVuZG9iago0MCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMzNlMwULAwAhKmpoYK5kaWCimGXEA+iJXLBRPLAbPMLMyBLCMLkJYcLkMLYzBtYmykYGZiBmRZIDEgutIAcvgSkQplbmRzdHJlYW0KZW5kb2JqCjQxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTcgPj4Kc3RyZWFtCnicMza0UDCAwxRDLgAalALsCmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMzEgPj4Kc3RyZWFtCnicRY/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/AOopVsWCmVuZHN0cmVhbQplbmRvYmoKNDUgMCBvYmoKPDwgL0JCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNzIKL1N1YnR5cGUgL0Zvcm0gL1R5cGUgL1hPYmplY3QgPj4Kc3RyZWFtCnic49I1sjBVsDAwUMjl0jUyNAYzc7h0LY0VzAzNQCxDM0MY08jEUsHcGMw0NjaHiZoYmMIVQM2CqjU1gxgLZeZwpQEAk4MVTgplbmRzdHJlYW0KZW5kb2JqCjQ2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjEwID4+CnN0cmVhbQp4nDVQyw1DMQi7ZwoWqBQCgWSeVr11/2tt0DthEf9CWMiUCHmpyc4p6Us+OkwPti6/sSILrXUl7MqaIJ4r76GZsrHR2OJgcBomXoAWN2DoaY0aNXThgqYulUKBxSXwmXx1e+i+Txl4ahlydgQRQ8lgCWq6Fk1YtDyfkE4B4v9+w+4t5KGS88qeG/kbnO3wO7Nu4SdqdiLRchUy1LM0xxgIE0UePHlFpnDis9Z31TQS1GYLTpYBrk4/jA4AYCJeWYDsrkQ5S9KOpZ9vvMf3D0AAU7QKZW5kc3RyZWFtCmVuZG9iagoyMSAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMgL0NoYXJQcm9jcyAyMiAwIFIKL0VuY29kaW5nIDw8Ci9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0MCAvcGFyZW5sZWZ0IC9wYXJlbnJpZ2h0IDQ2IC9wZXJpb2QgNDggL3plcm8gL29uZSAvdHdvIC90aHJlZQo1MyAvZml2ZSA1NSAvc2V2ZW4gODcgL1cgOTcgL2EgOTkgL2MgMTAxIC9lIDEwNCAvaCAvaSAxMDggL2wgL20gMTEyIC9wIDExNAovciAxMTYgL3QgXQovVHlwZSAvRW5jb2RpbmcgPj4KL0ZpcnN0Q2hhciAwIC9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnREZXNjcmlwdG9yIDIwIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDE5IDAgUiA+PgplbmRvYmoKMjAgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zIC9JdGFsaWNBbmdsZSAwCi9NYXhXaWR0aCAxMzQyIC9TdGVtViAwIC9UeXBlIC9Gb250RGVzY3JpcHRvciAvWEhlaWdodCAwID4+CmVuZG9iagoxOSAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzQyIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjMgNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjEyIDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTIgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwNQo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTgyIDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoyMiAwIG9iago8PCAvVyAyNSAwIFIgL2EgMjYgMCBSIC9jIDI3IDAgUiAvZSAyOCAwIFIgL2ZpdmUgMjkgMCBSIC9oIDMwIDAgUgovaSAzMSAwIFIgL2wgMzIgMCBSIC9tIDMzIDAgUiAvb25lIDM0IDAgUiAvcCAzNSAwIFIgL3BhcmVubGVmdCAzNiAwIFIKL3BhcmVucmlnaHQgMzcgMCBSIC9wZXJpb2QgMzggMCBSIC9yIDM5IDAgUiAvc2V2ZW4gNDAgMCBSIC9zcGFjZSA0MSAwIFIKL3QgNDIgMCBSIC90aHJlZSA0MyAwIFIgL3R3byA0NCAwIFIgL3plcm8gNDYgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAyMSAwIFIgL0YyIDE0IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL0NBIDAgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTIgPDwgL0NBIDEgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTMgPDwgL0NBIDAuMyAvVHlwZSAvRXh0R1N0YXRlIC9jYSAxID4+Ci9BNCA8PCAvQ0EgMC44IC9UeXBlIC9FeHRHU3RhdGUgL2NhIDAuOCA+PiA+PgplbmRvYmoKNSAwIG9iago8PCA+PgplbmRvYmoKNiAwIG9iago8PCA+PgplbmRvYmoKNyAwIG9iago8PCAvRjEtRGVqYVZ1U2Fucy1PbWVnYSAyMyAwIFIgL0YxLURlamFWdVNhbnMtUGhpIDI0IDAgUgovRjEtRGVqYVZ1U2Fucy11bmkwMzAyIDQ1IDAgUiA+PgplbmRvYmoKMiAwIG9iago8PCAvQ291bnQgMSAvS2lkcyBbIDEwIDAgUiBdIC9UeXBlIC9QYWdlcyA+PgplbmRvYmoKNDcgMCBvYmoKPDwgL0NyZWF0aW9uRGF0ZSAoRDoyMDIxMDExNTE1MDA0NyswMicwMCcpCi9DcmVhdG9yIChNYXRwbG90bGliIHYzLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZykKL1Byb2R1Y2VyIChNYXRwbG90bGliIHBkZiBiYWNrZW5kIHYzLjMuMikgPj4KZW5kb2JqCnhyZWYKMCA0OAowMDAwMDAwMDAwIDY1NTM1IGYgCjAwMDAwMDAwMTYgMDAwMDAgbiAKMDAwMDAxNzM0MCAwMDAwMCBuIAowMDAwMDE2OTY3IDAwMDAwIG4gCjAwMDAwMTcwMTAgMDAwMDAgbiAKMDAwMDAxNzE5MyAwMDAwMCBuIAowMDAwMDE3MjE0IDAwMDAwIG4gCjAwMDAwMTcyMzUgMDAwMDAgbiAKMDAwMDAwMDA2NSAwMDAwMCBuIAowMDAwMDAwMzk5IDAwMDAwIG4gCjAwMDAwMDAyMDggMDAwMDAgbiAKMDAwMDAwNjU4NCAwMDAwMCBuIAowMDAwMDA3OTMyIDAwMDAwIG4gCjAwMDAwMDc3MjQgMDAwMDAgbiAKMDAwMDAwNzM5NCAwMDAwMCBuIAowMDAwMDA4OTg1IDAwMDAwIG4gCjAwMDAwMDY2MDUgMDAwMDAgbiAKMDAwMDAwNzAwNSAwMDAwMCBuIAowMDAwMDA3MjMwIDAwMDAwIG4gCjAwMDAwMTU2MzggMDAwMDAgbiAKMDAwMDAxNTQzOCAwMDAwMCBuIAowMDAwMDE0OTg3IDAwMDAwIG4gCjAwMDAwMTY2OTEgMDAwMDAgbiAKMDAwMDAwOTAzNyAwMDAwMCBuIAowMDAwMDA5NDA4IDAwMDAwIG4gCjAwMDAwMDk3ODAgMDAwMDAgbiAKMDAwMDAwOTk0MSAwMDAwMCBuIAowMDAwMDEwMzE4IDAwMDAwIG4gCjAwMDAwMTA2MjEgMDAwMDAgbiAKMDAwMDAxMDkzOSAwMDAwMCBuIAowMDAwMDExMjU5IDAwMDAwIG4gCjAwMDAwMTE0OTUgMDAwMDAgbiAKMDAwMDAxMTYzNSAwMDAwMCBuIAowMDAwMDExNzUyIDAwMDAwIG4gCjAwMDAwMTIwODAgMDAwMDAgbiAKMDAwMDAxMjIzMiAwMDAwMCBuIAowMDAwMDEyNTQxIDAwMDAwIG4gCjAwMDAwMTI3NjEgMDAwMDAgbiAKMDAwMDAxMjk4MyAwMDAwMCBuIAowMDAwMDEzMTA0IDAwMDAwIG4gCjAwMDAwMTMzMzQgMDAwMDAgbiAKMDAwMDAxMzQ3NCAwMDAwMCBuIAowMDAwMDEzNTYzIDAwMDAwIG4gCjAwMDAwMTM3NjcgMDAwMDAgbiAKMDAwMDAxNDE3OCAwMDAwMCBuIAowMDAwMDE0NDk5IDAwMDAwIG4gCjAwMDAwMTQ3MDQgMDAwMDAgbiAKMDAwMDAxNzQwMCAwMDAwMCBuIAp0cmFpbGVyCjw8IC9JbmZvIDQ3IDAgUiAvUm9vdCAxIDAgUiAvU2l6ZSA0OCA+PgpzdGFydHhyZWYKMTc1NTcKJSVFT0YK\n", "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2021-01-15T15:00:47.197924\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.3.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import scipy.signal as sig\n", "\n", "\n", "K = 4096 # length of random signal\n", "N = 3 # order of AR model\n", "a = np.array((1, -1, 0.5)) # coefficients of AR model\n", "\n", "# generate random signal n[k]\n", "np.random.seed(2)\n", "n = np.random.normal(size=K)\n", "\n", "# AR model for random signal x[k]\n", "x = np.zeros(K)\n", "for k in np.arange(3, K):\n", " x[k] = a[0] * x[k - 1] + a[1] * x[k - 2] + a[2] * x[k - 3] + n[k]\n", "\n", "# estimate AR parameters by Yule-Walker method\n", "rho, sigma = sm.regression.yule_walker(x, order=N, method=\"mle\")\n", "\n", "# compute true and estimated transfer function\n", "Om, H = sig.freqz(1, np.insert(-a, 0, 1), worN=256)\n", "Om, He = sig.freqz(1, np.insert(-rho, 0, 1), worN=256)\n", "# compute PSD by Welch method\n", "Om2, Pxx = sig.welch(x, nperseg=511, return_onesided=True)\n", "\n", "# plot PSDs\n", "plt.figure(figsize=(10, 5))\n", "plt.plot(Om, np.abs(H) ** 2, label=r\"$\\Phi_{xx}(e^{j\\Omega})$\")\n", "plt.plot(\n", " Om2 * 2 * np.pi,\n", " 0.5 * Pxx,\n", " \"k-\",\n", " alpha=0.3,\n", " label=r\"$\\hat{\\Phi}_{xx}(e^{j\\Omega})$ (Welch)\",\n", ")\n", "plt.plot(Om, np.abs(He) ** 2, label=r\"$\\hat{\\Phi}_{xx}(e^{j\\Omega})$ (parametric)\")\n", "\n", "plt.xlabel(r\"$\\Omega$\")\n", "plt.axis([0, np.pi, 0, 20])\n", "plt.legend()\n", "plt.grid()\n", "\n", "# compute bias/variance of the estimators\n", "print(\"Bias of the Welch estimate: \\t\\t {0:1.4f}\".format(np.mean(Pxx - np.abs(H) ** 2)))\n", "print(\"Variance of the Welch estimate: \\t {0:1.4f}\".format(np.var(Pxx)))\n", "print(\n", " \"Bias of the parametric estimate: \\t {0:1.4f}\".format(\n", " np.mean(np.abs(H) ** 2 - np.abs(He) ** 2)\n", " )\n", ")\n", "print(\n", " \"Variance of the parametric estimate: \\t {0:1.4f}\".format(np.var(np.abs(He) ** 2))\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**\n", "\n", "* Change the order `N` of the AR model used for estimation by the Yule-Walker equations. What happens if the order is smaller or higher than the order of the true system? Why?\n", "\n", "* Change the number of samples `K`. Is the estimator consistent?\n", "\n", "Solution: Choosing the order of the estimated AR model differently from the true process results in a mismatch of the model with the consequence of potentially large deviations between the estimated PSD and true PSD. These deviations are typically larger when choosing the order smaller since the model does definitely not fit to the true process. However, choosing the order larger is typically not that problematic since some of the AR coefficients are estimated as approximately zero due to the lower order of the process generating the random signal. Increasing the number of samples seems to lower the bias and variance of the estimated PSD, the estimator can therefore be assumed to be consistent." ] }, { "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": { "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 }