{
"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 Welch Method\n",
"\n",
"In the previous section it has been shown that the [periodogram](periodogram.ipynb) as a non-parametric estimator of the power spectral density (PSD) $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of a random signal $x[k]$ is not consistent. This is due to the fact that its variance does not converge towards zero even when the length of the random signal is increased towards infinity. In order to overcome this problem, the [Bartlett method](https://en.wikipedia.org/wiki/Bartlett's_method) and [Welch method](https://en.wikipedia.org/wiki/Welch's_method) \n",
"\n",
"1. split the random signal into segments, \n",
"2. estimate the PSD for each segment, and \n",
"3. average over these local estimates. \n",
"\n",
"The averaging reduces the variance of the estimated PSD. While Barlett's method uses non-overlapping segments, Welch's is a generalization using windowed overlapping segments. For the discussion of Welch's method we assume a wide-sense ergodic real-valued random process."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Derivation\n",
"\n",
"The random signal $x[k]$ is split into into $L$ overlapping segments of length $N$, starting at multiples of the step size $M \\in {1,2, \\dots, N}$. These segments are windowed by the window $w[k]$ of length $N$, resulting in a windowed $l$-th segment $x_l[k]$ with $0\\leq l\\leq L-1$. The discrete-time Fourier transformation (DTFT) $X_l(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the windowed $l$-th segment is then given as \n",
"\n",
"\\begin{equation}\n",
"X_l(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sum_{k = 0}^{N-1} x[k + l \\cdot M] \\, w[k] \\; \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega\\,k}\n",
"\\end{equation}\n",
"\n",
"where the window $w[k]$ defined within $0\\leq k\\leq N-1$ should be normalized as $\\frac{1}{N} \\sum\\limits_{k=0}^{N-1} | w[k] |^2 = 1$. The latter condition ensures that the power of the signal is maintained in the estimate. The stepsize $M$ determines the overlap between the segments. In general, $N-M$ number of samples overlap between adjacent segments. For $M = N$ no overlap occurs. The overlap is sometimes given as ratio $\\frac{N-M}{N}\\cdot 100\\%$. \n",
"\n",
"Introducing $X_l(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ into the definition of the periodogram yields the periodogram of the $l$-th segment\n",
"\n",
"\\begin{equation}\n",
"\\hat{\\Phi}_{xx,l}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{1}{N} \\,| X_l(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2\n",
"\\end{equation}\n",
"\n",
"The estimated PSD is then given by averaging over the segment's periodograms $\\hat{\\Phi}_{xx,l}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ \n",
"\n",
"\\begin{equation}\n",
"\\hat{\\Phi}_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{1}{L} \\sum_{l = 0}^{L-1} \\hat{\\Phi}_{xx,l}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\n",
"\\end{equation}\n",
"\n",
"Note, that the total number $L$ of segments has to be chosen such that the last required sample $(L-1)\\cdot M + N - 1$ does not exceed the total length of the random signal. Otherwise the last segment $x_{L-1}[k]$ may be zeropadded towards length $N$.\n",
"\n",
"The Bartlett method uses a rectangular window $w[k] = \\text{rect}_N[k]$ and non-overlapping segments $M=N$. The Welch method uses overlapping segments and a window that must be chosen according to the intended spectral analysis task."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example\n",
"\n",
"The following example is equivalent to the previous [periodogram example](periodogram.ipynb#Example---Periodogram). We aim at estimating the PSD of a random process which draws samples from normally distributed white noise with zero-mean and unit variance. The true PSD is consequently given as $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = 1$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Bias of the Welch estimate: \t\t -0.0114\n",
"Variance of the Welch estimate: \t 0.0255\n"
]
},
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDYwMS42NjU2MjUgMjgxLjYyODc1IF0gL1BhcmVudCAyIDAgUiAvUmVzb3VyY2VzIDggMCBSCi9UeXBlIC9QYWdlID4+CmVuZG9iago5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTEgMCBSID4+CnN0cmVhbQp4nOVaSW8exxG9f79ijvKBo67q/WhBtoAAQSxbSC65CBRNSxAVWEucn59Xs75ZviaNHCNBAPlU3dXVtde0dB8u0t13rvuAf3900r3Cv/uLw28Pl+SkTykmjfj1I/+qBT9qyRGw2/z22+Xy6+X599jkCxa9ulx86sO4yOc+Bg+qh0us4YB+ZFRj7b0b4XUHRidGOjK6x7EhQl8gBNgacs6ZwNC7ifHlBW7hj8uLN93zH6UT17359aKlj2I8ta+xe/Pu8sz17rvuzYfuhzeXgdEF/91rDq545sBom4VI7r26GouG6Dec4paTCnYqGoMwJ0YfEcalXrQGSaG6RJxkL5N32pe6YbNAj/CA9qrPJYfsY9nw2EnjK+yl1gAqZkNom5Mvvs8hea+1iBAn3UsTCq7S5SqJOTHa5hSy9jEGkKqmsOG0kykmmKYkn8PGzgltc4pJep+8hhqyZOLkV5l+xzLX3TjYuheZ7bjC2opZ6u0DtrwYHRzBGalbfujTSACGNy/vPrz9+7df3n76cvO3h7v7t93Lf11e4++knFOvwXEO7nomRu6184N2LBxMXrM3sYVDgqekPYcFvM4hmfDEQfcGtnCo2LbsOSzgdQ6l9lKIQ7wqg6jrs9uzWNHrPERSD29cmeSrYkgofdADkwVtMAkBul+YSEMbkhEkDgpf0QYTuIkPxOS6QtT5vhx0vqLXmaiDKSVicl0n6vG/B7WvaIOJ5j4XYnJdJ2sy2jChFHWVSYw4/MJESSe/dyf+F4dwirOk7vNd94/uUwdBLNEhpktMJSmInLnE+CdfzPhiKeqj1u7nV93zl3f/fn979/OrF93tl8cSY8RPXpHtTYRfLiHCFQYiDqELKFF7i8WiI3n0MGovu3y/gBKkr1mrCyN5UpinlO3uKygIJsllnGskz4iSwe2OvoIyxB51UkdyBHoX6+7sK2hnz756lyfy0ucxljL5AkrMPSSt883U3IeadmdfQcG5JKjGkVqQg01d290JlYKbjDFBzeMCsYS6v3lCJYa+SAGPaYHisuL+7gmVCP9HZZKmyxfYrc/72ydUguu997ClaQG05+r+/gkVTX1ROJDMC2pf3F4DhCK04HpzqjOHiLpL9zogVMyWUqluUrEkuFTYa4FQ05pHaq2z0DkektkDowLjLsjHZT5SQQg9eAChAp9POcY866H6Phx8gFBTXCihxknTigpLDl5AKOo2eI24Mgmtgmh08ANCUXMYszILrVJRyOz1QChMsMctV52EVtS/evAFQiXA0B1qqFkGBJB68AZCYabgoNlNt4Sat88HfyDUnFlQ4c0M7NIO7kCopboM55TJ4VBFgOFeDYTi7HCl6HQyVoXDloM7EGqW4UIuYb4kKyQP7kCouVLAgUKeF1h2OKhhRc1D1fD5SLX07uAOhFoIA7XO8dQjnuSDOxA6CF2lTgy8xLNksKJmeqFEFDPTAg0n+YDQgQFspEx69pYcD95AqCAdwvRKmfTmg5zkBELtkioqgflE0Z0kBUKtIyu+hJDmBfUkLRCKVqn3NekiM6LJMTEQaraNiKRz6PYomo6pgVC71YRmbJGhxJPkQKjpWVWrzreKaHJMDoRaYEWEiTIZRkA0OSYHQsWjk/AOVzAtQItxTA6Emhoko8iYvAF59yQ5ECohod/V6HVeUE6SA6GWJyCCmy0jGNnBGwi1KIbD6ew+UPlJciDUUpeiBIh+LnbCSXIgVDzECeai0wKEk2NyINSiGIp0Wcop1MrH5EDokNU1+DRfa5GT5ECoJIT9Wtyc1QPCyTE5EGrJAYrLi+IQTo7JgVC71orUMEUx9M4nuYFQU3TMKHEnd4gIJ8fcQKh4/BhinWuliHByzA2EClpAdPA+ylx1+pPkQKgpOpk0k3XHIRsd6tQVNYeD+LIcCfHkmBwINaHhn4uxYvFJciDUoowLcNB5QSonyYFQS7nqtZRJ0THnk+RA6JARU8hz+RYRT47JgVDrKFGa62x7Z3O6zfSuWvgIcS6VXndPbGCe2Lp8vr9O2e2bHKa8XKf8fD8NUoY5yqHrcQiZ3e1D9/yvzgYjM2HpRy3X1Dv0C6jCMyrQ63Q3qMVCtWLbhxBahFAvInDA2VODDMJho5isb2yQCfSH9qrkUlq7IcKYuyHuxdzazaNPhQClgrBFh5yA42kKRVuXd6PWQGVIiyzTOt/N0BQosmsNTTox1eaC7VCENwlNFBhjtpq1JQoiHfQmyAm5xRmMUUOh4kZEadF5aDcPtX7zpm9snlecOrR/2hIEchR0EgXHq02zQk2FMGMlPaqYphyhJNSpCRmnecAeR5MagpOWAd5Y3xWG6ky0bfY54IqRrF2TLyqgmEIKOGLzAmH4CERBag5NeSEHql/E7RhbdIpEEGxYLjDspkmjO0AGyN7KiOYFIlSiDUWAyU1PhwsX1EAIVrjCFl0vglv2ITVd2MZ/bgh/rnl/aiMpUPoqLa4oU6oGtFmueX1mz7BTRXp3oXkt1uPB7EPU1LQXtS8YVsU7aekDbo4EiHZFm2wtHKCUQ6QsTfewpIUNcYGPmCm0G1DqJjfUJG27t69DKJTaASH0apk1hna8RyeEUjBX+EhTH4j4Jgc03DQDlK1eY4TxpZa4qLKTF8SNkptKszYyBC9WXjxmLqiREZtL05oRTRHroeCSW0YArcFMFPEqNoMQcnn2CXcnj9g8WvoEzwhtaWGhJUKAIm1TtpiL+tuhxcyttGqESStKV9+WA/0t2gzkXz9+B7q+YQVnKxAUUrezAnJ+Rr2R2lkGJW/KJkipbZGHcZwVHKU2w1XBHWZ0Jgj4j4QrmBSkhSm0DhhwMTU6JJrUvJhoHaEFK7FZ50r4+uklrXZ/mebyYnVqNf2iHrWKFuENJUVNrhyH8Ve+8tCAnb7y8Bx9RWlcTihNxQml4Tej64ybUBplE8oTa4ZpLs0wTZ8Z5hkzwTRJ3sA0LyaYpsIM0+yXYZ7wEkxzXIZpWkswz2QZpsnrBl7nqwzzFJVgnpUSTBNRhnnwSTCNNxnmISbBNKrcwOtAkmEaOxLMw0WGaYbIME0KGaZ5IMM09WOYhnsbuJ5Z92ZQRzCN4ximoRvDPFpbYR6gMUxjMoZ5GMZwObNvHmwxzOMrgmlIxTCNohimgRPDNFbawPXMvnlExFGKJ0EE07yHYZrqbEKdP7NvntAwzHMYhuuZfW9mKgTz5ITg7XxkgdcpyOZD/DgFsWzglidJf+r76+YT8hNmMDv6P/HgakPN37Ov7+4GucbnVjI8trpfPnDr8Axo+KA3WhYq2Gmtfef+4cvX9w9vv969695+etd9/fztrvvpl5fzl+/n33t7yrW8RsO9jW/RYDfrGfHjcsaSj7BpC7mslmr5eYFRdTBtzc7Kiw1aJtLbC8P2dWXZmGCrk7M98WJ2cMCZ3Xq0BbzdCLLAo18gH2pV3kLCGTdGlTYOR0E+buFVamK3XtDZHd/aI7sX8yO7ob4xrZ/UOF27xok41Lw7zKOUwVFtPnuG7+xsfdBH75/QFq0rkQvGleMLqGwfh/IYaG5cn8fgBwt8/qN22Wz02X/M6i7BPs/HMVC5wURnPEhyU0jq0xj3hvXzK45n/3w2UFo5Jkgv44hu5jHR3A0kqU+SAnp5/EH/Usb4sjnNh4Ew24DFT/9/s6Vc+X430F555FXi+TOvn357PxazwyI0gNMq/ITMU9EDIxWkkLqlFHrkqdgj708Oo9mdOeiqdtVZeQ9bPPdS1zD1RFafJ4N56nT5suUZZp52lorWBtVrYDza1x/7NjNSR3EBjaih6PtqMDRaT+ni4JbeHoDUqGnAE5pqb3vIyhGuVguaky26yD7uQTgiTS2h2+yMAITfJW7PYQ8P49DlXTanRiLPxT5UsoTeLRLSfawo9rDRbXZO0pZ6fMJWdjtDLLT9bn8OXAK8LWvYnFrjmYRblM4Rz26Pd15vms/BeuFTr1okCUnjp/YxRUY3JKtmcFrt+yQ42Szy/yE2Zftmonn6PncSXL59em8UFKR2kc2N9df/ENpw6Mfj2uvLfwHTkHDvCmVuZHN0cmVhbQplbmRvYmoKMTEgMCBvYmoKMzAyMgplbmRvYmoKMTcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMjcgPj4Kc3RyZWFtCnicPZJLjsQwCET3OQUXaMn8/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+jZr0eZXneP0BbPNzOwplbmRzdHJlYW0KZW5kb2JqCjMxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjI3ID4+CnN0cmVhbQp4nDVPO7IDIQzrOYUukBmMbWDPs5lUL/dvn2SyDRL+SPL0REcmXubICKzZ8bYWGYgZ+BZT8a897cOE6j24hwjl4kKYYSScNeu4m6fjxb9d5TPWwbsNvmKWFwS2MJP1lcWZy3bBWBoncU6yG2PXRGxjXevpFNYRTCgDIZ3tMCXIHBUpfbKjjDk6TuSJ52KqxS6/72F9waYxosIcVwVP0GRQlj3vJqAdF/Tf1Y3fSTSLXgIykWBhnSTmzllO+NVrR8dRiyIxJ6QZ5DIR0pyuYgqhCcU6OwoqFQWX6nPK3T7/aF1bTQplbmRzdHJlYW0KZW5kb2JqCjMyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ1ID4+CnN0cmVhbQp4nEVQu41DMQzrPQUXCGD9LHued0iV2789SkZwhSFaP5JaEpiIwEsMsZRv4kdGQT0LvxeF4jPEzxeFQc6EpECc9RkQmXiG2kZu6HZwzrzDM4w5AhfFWnCm05n2XNjknAcnEM5tlPGMQrpJVBVxVJ9xTPGqss+N14GltWyz05HsIY2ES0klJpd+Uyr/tClbKujaRROwSOSBk0004Sw/Q5JizKCUUfcwtY70cbKRR3XQydmcOS2Z2e6n7Ux8D1gmmVHlKZ3nMj4nqfNcTn3usx3R5KKlVfuc/d6RlvIitduh1elXJVGZjdWnkLg8/4yf8f4DjqBZPgplbmRzdHJlYW0KZW5kb2JqCjMzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ3ID4+CnN0cmVhbQp4nE1Ru21EMQzr3xRc4ADra3meC1Jd9m9DyQiQwiChLymnJRb2xksM4QdbD77kkVVDfx4/MewzLD3J5NQ/5rnJVBS+FaqbmFAXYuH9aAS8FnQvIivKB9+PZQxzzvfgoxCXYCY0YKxvSSYX1bwzZMKJoY7DQZtUGHdNFCyuFc0zyO1WN7I6syBseCUT4sYARATZF5DNYKOMsZWQxXIeqAqSBVpg1+kbUYuCK5TWCXSi1sS6zOCr5/Z2N0Mv8uCounh9DOtLsMLopXssfK5CH8z0TDt3SSO98KYTEWYPBVKZnZGVOj1ifbdA/59lK/j7yc/z/QsVKFwqCmVuZHN0cmVhbQplbmRvYmoKMzQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA2OCA+PgpzdHJlYW0KeJwzMrdQMFCwNAEShhYmCuZmBgophlxAvqmJuUIuF0gMxMoBswyAtCWcgohbQjRBlIJYEKVmJmYQSTgDIpcGAMm0FeUKZW5kc3RyZWFtCmVuZG9iagozNSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI1NSA+PgpzdHJlYW0KeJxFkUuSAyAIRPeegiOA/OQ8mZpVcv/tNJhMNnaXqP2ESiOmEiznFHkw/cjyzWS26bUcq52NAooiFMzkKvRYgdWdKeLMtUS19bEyctzpHYPiDeeunFSyuFHGOqo6FTim58r6qu78uCzKviOHMgVs1jkONnDltmGME6PNVneH+0SQp5Opo+J2kGz4g5PGvsrVFbhONvvqJRgHgn6hCUzyTaB1hkDj5il6cgn28XG780Cwt7wJpGwI5MgQjA5Bu06uf3Hr/N7/OsOd59oMV4538TtMa7vjLzHJirmARe4U1PM9F63rDB3vyZljctN9Q+dcsMvdQabP/B/r9w9QimaICmVuZHN0cmVhbQplbmRvYmoKMzYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjEgPj4Kc3RyZWFtCnicRZBLEsMgDEP3nEJH8EcGfJ50ukrvv60hTbOAp7FABncnBKm1BRPRBS9tS7oLPlsJzsZ46DZuNRLkBHWAVqTjaJRSfbnFaZV08Wg2cysLrRMdZg56lKMZoBA6Fd7touRypu7O+Udw9V/1R7HunM3EwGTlDoRm9SnufJsdUV3dZH/SY27Wa38V9qqwtKyl5YTbzl0zoATuqRzt/QWpczqECmVuZHN0cmVhbQplbmRvYmoKMzcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4MCA+PgpzdHJlYW0KeJxFjLsNwDAIRHumYAR+JmafKJWzfxsgStxwT7p7uDoSMlPeYYaHBJ4MLIZT8QaZo2A1uEZSjZ3so7BuX3WB5npTq/X3BypPdnZxPc3LGfQKZW5kc3RyZWFtCmVuZG9iagozOCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE0NyA+PgpzdHJlYW0KeJw9T7kNAzEM6z0FFzjAeixb81yQ6rJ/G8pGUggiQPGRZUfHClxiApOOORIvaT/4aRqBWAY1R/SEimFY4G6SAg+DLEpXni1eDJHaQl1I+NYQ3q1MZKI8rxE7cCcXowc+VBtZHnpAO0QVWa5Jw1jVVl1qnbACHLLOwnU9zKoE5dEnaykfUFRCvXT/n3va+wsAby/rCmVuZHN0cmVhbQplbmRvYmoKMzkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNDkgPj4Kc3RyZWFtCnicNY9LDgMhDEP3OYUvMFJ+hHAeqq6m9982YVoJCQvbL8EWg5GMS0xg7Jhj4SVUT60+JCOPukk5EKlQNwRPaEwMM2zSJfDKdN8ynlu8nFbqgk5I5OmsNhqijGZew9FTzgqb/svcJGplRpkDMutUtxOysmAF5gW1PPcz7qhc6ISHncqw6E4xotxmRhp+/9v0/gJ7MjBjCmVuZHN0cmVhbQplbmRvYmoKNDAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA0OSA+PgpzdHJlYW0KeJwzNrRQMFAwNDAHkkaGQJaRiUKKIRdIAMTM5YIJ5oBZBkAaojgHriaHKw0AxugNJgplbmRzdHJlYW0KZW5kb2JqCjQxIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTU3ID4+CnN0cmVhbQp4nEWQuRFDMQhEc1VBCRKwCOqxx9F3/6kX+Uq0bwAth68lU6ofJyKm3Ndo9DB5Dp9NJVYs2Ca2kxpyGxZBSjGYeE4xq6O3oZmH1Ou4qKq4dWaV02nLysV/82hXM5M9wjXqJ/BN6PifPLSp6FugrwuUfUC1OJ1JUDF9r2KBo5x2fyKcGOA+GUeZKSNxYm4K7PcZAGa+V7jG4wXdATd5CmVuZHN0cmVhbQplbmRvYmoKNDIgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzIgPj4Kc3RyZWFtCnicLVI5jiQxDMv9Cn5gAOvy8Z4eTNT7/3RJVQUFqmzLPORyw0QlfiyQ21Fr4tdGZqDC8K+rzIXvSNvIOohryEVcyZbCZ0Qs5DHEPMSC79v4GR75rMzJswfGL9n3GVbsqQnLQsaLM7TDKo7DKsixYOsiqnt4U6TDqSTY44v/PsVzF4IWviNowC/556sjeL6kRdo9Ztu0Ww+WaUeVFJaD7WnOy+RL6yxXx+P5INneFTtCaleAojB3xnkujjJtZURrYWeDpMbF9ubYj6UEXejGZaQ4AvmZKsIDSprMbKIg/sjpIacyEKau6Uont1EVd+rJXLO5vJ1JMlv3RYrNFM7rwpn1d5gyq807eZYTpU5F+Bl7tgQNnePq2WuZhUa3OcErJXw2dnpy8r2aWQ/JqUhIFdO6Ck6jyBRL2Jb4moqa0tTL8N+X9xl//wEz4nwBCmVuZHN0cmVhbQplbmRvYmoKNDMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA2OCA+PgpzdHJlYW0KeJwzMzZTMFCwMAISpqaGCuZGlgophlxAPoiVywUTywGzzCzMgSwjC5CWHC5DC2MwbWJspGBmYgZkWSAxILrSAHL4EpEKZW5kc3RyZWFtCmVuZG9iago0NCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE3ID4+CnN0cmVhbQp4nDM2tFAwgMMUQy4AGpQC7AplbmRzdHJlYW0KZW5kb2JqCjQ1IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTMxID4+CnN0cmVhbQp4nEWPyw0EIQxD71ThEvIZPqmH1Z7Y/q/rMJpBQvhBIjvxMAis8/I20MXw0aLDN/421atjlSwfunpSVg/pkIe88hVQaTBRxIVZTB1DYc6YysiWMrcb4bZNg6xslVStg3Y8Bg+2p2WrCH6pbWHqLPEMwlVeuMcNP5BLrXe9Vb5/QlMwlwplbmRzdHJlYW0KZW5kb2JqCjQ2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMzM4ID4+CnN0cmVhbQp4nDVSOa7dQAzrfQpdIIB2zZznBal+7t+GlF8KQ7RWipqOFpVp+WUhVS2TLr/tSW2JG/L3yQqJE5JXJdqlDJFQ+TyFVL9ny7y+1pwRIEuVCpOTksclC/4Ml94uHOdjaz+PI3c9emBVjIQSAcsUE6NrWTq7w5qN/DymAT/iEXKuWLccYxVIDbpx2hXvQ/N5yBogZpiWigpdVokWfkHxoEetffdYVFgg0e0cSXCMjVCRgHaB2kgMObMWu6gv+lmUmAl07Ysi7qLAEknMnGJdOvoPPnQsqL8248uvjkr6SCtrTNp3o0lpzCKTrpdFbzdvfT24QPMuyn9ezSBBU9YoaXzQqp1jKJoZZYV3HJoMNMcch8wTPIczEpT0fSh+X0smuiiRPw4NoX9fHqOMnAZvAXPRn7aKAxfx2WGvHGCF0sWa5H1AKhN6YPr/1/h5/vwDHLaAVAplbmRzdHJlYW0KZW5kb2JqCjQ3IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ4ID4+CnN0cmVhbQp4nC1ROZIDQQjL5xV6QnPT77HLkff/6QrKAYOGQyA6LXFQxk8Qlive8shVtOHvmRjBd8Gh38p1GxY5EBVI0hhUTahdvB69B3YcZgLzpDUsgxnrAz9jCjd6cXhMxtntdRk1BHvXa09mUDIrF3HJxAVTddjImcNPpowL7VzPDci5EdZlGKSblcaMhCNNIVJIoeomqTNBkASjq1GjjRzFfunLI51hVSNqDPtcS9vXcxPOGjQ7Fqs8OaVHV5zLycULKwf9vM3ARVQaqzwQEnC/20P9nOzkN97SubPF9Phec7K8MBVY8ea1G5BNtfg3L+L4PePr+fwDqKVbFgplbmRzdHJlYW0KZW5kb2JqCjQ4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMTcxID4+CnN0cmVhbQp4nE2QTQ5CIRCD95yiFzCh8wOP82hc6f23dvD54oL0SyFDp8MDHUfiRkeGzuh4sMkxDrwLMiZejfOfjOskjgnqFW3BurQ77s0sMScsEyNga5Tcm0cU+OGYC0GC7PLDFxhEpGuYbzWfdZN+frvTXdSldffTIwqcyI5QDBtwBdjTPQ7cEs7vmia/VCkZmziUD1QXkbLZCYWopWKXU1VojOJWPe+LXu35AcH2O/sKZW5kc3RyZWFtCmVuZG9iago0OSAwIG9iago8PCAvQkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA3MgovU3VidHlwZSAvRm9ybSAvVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJzj0jWyMFWwMDBQyOXSNTI0BjNzuHQtjRXMDM1ALEMzQxjTyMRSwdwYzDQ2NoeJmhiYwhVAzYKqNTWDGAtl5nClAQCTgxVOCmVuZHN0cmVhbQplbmRvYmoKNTAgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyMTAgPj4Kc3RyZWFtCnicNVDLDUMxCLtnChaoFAKBZJ5WvXX/a23QO2ER/0JYyJQIeanJzinpSz46TA+2Lr+xIgutdSXsypognivvoZmysdHY4mBwGiZegBY3YOhpjRo1dOGCpi6VQoHFJfCZfHV76L5PGXhqGXJ2BBFDyWAJaroWTVi0PJ+QTgHi/37D7i3koZLzyp4b+Ruc7fA7s27hJ2p2ItFyFTLUszTHGAgTRR48eUWmcOKz1nfVNBLUZgtOlgGuTj+MDgBgIl5ZgOyuRDlL0o6ln2+8x/cPQABTtAplbmRzdHJlYW0KZW5kb2JqCjIyIDAgb2JqCjw8IC9CYXNlRm9udCAvRGVqYVZ1U2FucyAvQ2hhclByb2NzIDIzIDAgUgovRW5jb2RpbmcgPDwKL0RpZmZlcmVuY2VzIFsgMzIgL3NwYWNlIDQwIC9wYXJlbmxlZnQgL3BhcmVucmlnaHQgNDYgL3BlcmlvZCA0OCAvemVybyAvb25lIC90d28gL3RocmVlCjUzIC9maXZlIDU1IC9zZXZlbiA2OCAvRCAvRSA4MCAvUCA4MyAvUyA5NyAvYSAxMDAgL2QgL2UgMTA1IC9pIDEwOSAvbSAvbgoxMTQgL3IgL3MgL3QgL3UgXQovVHlwZSAvRW5jb2RpbmcgPj4KL0ZpcnN0Q2hhciAwIC9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnREZXNjcmlwdG9yIDIxIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zCi9TdWJ0eXBlIC9UeXBlMyAvVHlwZSAvRm9udCAvV2lkdGhzIDIwIDAgUiA+PgplbmRvYmoKMjEgMCBvYmoKPDwgL0FzY2VudCA5MjkgL0NhcEhlaWdodCAwIC9EZXNjZW50IC0yMzYgL0ZsYWdzIDMyCi9Gb250QkJveCBbIC0xMDIxIC00NjMgMTc5NCAxMjMzIF0gL0ZvbnROYW1lIC9EZWphVnVTYW5zIC9JdGFsaWNBbmdsZSAwCi9NYXhXaWR0aCAxMzQyIC9TdGVtViAwIC9UeXBlIC9Gb250RGVzY3JpcHRvciAvWEhlaWdodCAwID4+CmVuZG9iagoyMCAwIG9iagpbIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwCjYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgMzE4IDQwMSA0NjAgODM4IDYzNgo5NTAgNzgwIDI3NSAzOTAgMzkwIDUwMCA4MzggMzE4IDM2MSAzMTggMzM3IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYKNjM2IDYzNiAzMzcgMzM3IDgzOCA4MzggODM4IDUzMSAxMDAwIDY4NCA2ODYgNjk4IDc3MCA2MzIgNTc1IDc3NSA3NTIgMjk1CjI5NSA2NTYgNTU3IDg2MyA3NDggNzg3IDYwMyA3ODcgNjk1IDYzNSA2MTEgNzMyIDY4NCA5ODkgNjg1IDYxMSA2ODUgMzkwIDMzNwozOTAgODM4IDUwMCA1MDAgNjEzIDYzNSA1NTAgNjM1IDYxNSAzNTIgNjM1IDYzNCAyNzggMjc4IDU3OSAyNzggOTc0IDYzNCA2MTIKNjM1IDYzNSA0MTEgNTIxIDM5MiA2MzQgNTkyIDgxOCA1OTIgNTkyIDUyNSA2MzYgMzM3IDYzNiA4MzggNjAwIDYzNiA2MDAgMzE4CjM1MiA1MTggMTAwMCA1MDAgNTAwIDUwMCAxMzQyIDYzNSA0MDAgMTA3MCA2MDAgNjg1IDYwMCA2MDAgMzE4IDMxOCA1MTggNTE4CjU5MCA1MDAgMTAwMCA1MDAgMTAwMCA1MjEgNDAwIDEwMjMgNjAwIDUyNSA2MTEgMzE4IDQwMSA2MzYgNjM2IDYzNiA2MzYgMzM3CjUwMCA1MDAgMTAwMCA0NzEgNjEyIDgzOCAzNjEgMTAwMCA1MDAgNTAwIDgzOCA0MDEgNDAxIDUwMCA2MzYgNjM2IDMxOCA1MDAKNDAxIDQ3MSA2MTIgOTY5IDk2OSA5NjkgNTMxIDY4NCA2ODQgNjg0IDY4NCA2ODQgNjg0IDk3NCA2OTggNjMyIDYzMiA2MzIgNjMyCjI5NSAyOTUgMjk1IDI5NSA3NzUgNzQ4IDc4NyA3ODcgNzg3IDc4NyA3ODcgODM4IDc4NyA3MzIgNzMyIDczMiA3MzIgNjExIDYwNQo2MzAgNjEzIDYxMyA2MTMgNjEzIDYxMyA2MTMgOTgyIDU1MCA2MTUgNjE1IDYxNSA2MTUgMjc4IDI3OCAyNzggMjc4IDYxMiA2MzQKNjEyIDYxMiA2MTIgNjEyIDYxMiA4MzggNjEyIDYzNCA2MzQgNjM0IDYzNCA1OTIgNjM1IDU5MiBdCmVuZG9iagoyMyAwIG9iago8PCAvRCAyNCAwIFIgL0UgMjUgMCBSIC9QIDI3IDAgUiAvUyAyOSAwIFIgL2EgMzAgMCBSIC9kIDMxIDAgUiAvZSAzMiAwIFIKL2ZpdmUgMzMgMCBSIC9pIDM0IDAgUiAvbSAzNSAwIFIgL24gMzYgMCBSIC9vbmUgMzcgMCBSIC9wYXJlbmxlZnQgMzggMCBSCi9wYXJlbnJpZ2h0IDM5IDAgUiAvcGVyaW9kIDQwIDAgUiAvciA0MSAwIFIgL3MgNDIgMCBSIC9zZXZlbiA0MyAwIFIKL3NwYWNlIDQ0IDAgUiAvdCA0NSAwIFIgL3RocmVlIDQ2IDAgUiAvdHdvIDQ3IDAgUiAvdSA0OCAwIFIgL3plcm8gNTAgMCBSID4+CmVuZG9iagozIDAgb2JqCjw8IC9GMSAyMiAwIFIgL0YyIDE1IDAgUiA+PgplbmRvYmoKNCAwIG9iago8PCAvQTEgPDwgL0NBIDAgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTIgPDwgL0NBIDEgL1R5cGUgL0V4dEdTdGF0ZSAvY2EgMSA+PgovQTMgPDwgL0NBIDAuOCAvVHlwZSAvRXh0R1N0YXRlIC9jYSAwLjggPj4gPj4KZW5kb2JqCjUgMCBvYmoKPDwgPj4KZW5kb2JqCjYgMCBvYmoKPDwgPj4KZW5kb2JqCjcgMCBvYmoKPDwgL0YxLURlamFWdVNhbnMtT21lZ2EgMjYgMCBSIC9GMS1EZWphVnVTYW5zLVBoaSAyOCAwIFIKL0YxLURlamFWdVNhbnMtdW5pMDMwMiA0OSAwIFIgL00wIDEyIDAgUiA+PgplbmRvYmoKMTIgMCBvYmoKPDwgL0JCb3ggWyAtOCAtOCA4IDggXSAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMSAvU3VidHlwZSAvRm9ybQovVHlwZSAvWE9iamVjdCA+PgpzdHJlYW0KeJxtkEEOhCAMRfc9RS/wSUtFZevSa7iZTOL9twNxQEzdNNC+PH5R/pLwTqXA+CQJS06z5HrTkNK6TIwY5tWyKMegUS3WznU4qM/QcGN0i7EUptTW6Hijm+k23pM/+rBZIUY/HA6vhHsWQyZcKTEGh98LL9vD/xGeXtTAH6KNfmNaQ/0KZW5kc3RyZWFtCmVuZG9iagoyIDAgb2JqCjw8IC9Db3VudCAxIC9LaWRzIFsgMTAgMCBSIF0gL1R5cGUgL1BhZ2VzID4+CmVuZG9iago1MSAwIG9iago8PCAvQ3JlYXRpb25EYXRlIChEOjIwMjEwMTE1MTQ1MTUyKzAyJzAwJykKL0NyZWF0b3IgKE1hdHBsb3RsaWIgdjMuMy4yLCBodHRwczovL21hdHBsb3RsaWIub3JnKQovUHJvZHVjZXIgKE1hdHBsb3RsaWIgcGRmIGJhY2tlbmQgdjMuMy4yKSA+PgplbmRvYmoKeHJlZgowIDUyCjAwMDAwMDAwMDAgNjU1MzUgZiAKMDAwMDAwMDAxNiAwMDAwMCBuIAowMDAwMDE1NTk5IDAwMDAwIG4gCjAwMDAwMTUwMDIgMDAwMDAgbiAKMDAwMDAxNTA0NSAwMDAwMCBuIAowMDAwMDE1MTg3IDAwMDAwIG4gCjAwMDAwMTUyMDggMDAwMDAgbiAKMDAwMDAxNTIyOSAwMDAwMCBuIAowMDAwMDAwMDY1IDAwMDAwIG4gCjAwMDAwMDAzOTggMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDAzNDk1IDAwMDAwIG4gCjAwMDAwMTUzNDUgMDAwMDAgbiAKMDAwMDAwNDg0MyAwMDAwMCBuIAowMDAwMDA0NjM1IDAwMDAwIG4gCjAwMDAwMDQzMDUgMDAwMDAgbiAKMDAwMDAwNTg5NiAwMDAwMCBuIAowMDAwMDAzNTE2IDAwMDAwIG4gCjAwMDAwMDM5MTYgMDAwMDAgbiAKMDAwMDAwNDE0MSAwMDAwMCBuIAowMDAwMDEzNjQzIDAwMDAwIG4gCjAwMDAwMTM0NDMgMDAwMDAgbiAKMDAwMDAxMjk4OCAwMDAwMCBuIAowMDAwMDE0Njk2IDAwMDAwIG4gCjAwMDAwMDU5NDggMDAwMDAgbiAKMDAwMDAwNjE4MSAwMDAwMCBuIAowMDAwMDA2MzMyIDAwMDAwIG4gCjAwMDAwMDY3MDMgMDAwMDAgbiAKMDAwMDAwNjk0MSAwMDAwMCBuIAowMDAwMDA3MzEzIDAwMDAwIG4gCjAwMDAwMDc3MjQgMDAwMDAgbiAKMDAwMDAwODEwMSAwMDAwMCBuIAowMDAwMDA4NDAxIDAwMDAwIG4gCjAwMDAwMDg3MTkgMDAwMDAgbiAKMDAwMDAwOTAzOSAwMDAwMCBuIAowMDAwMDA5MTc5IDAwMDAwIG4gCjAwMDAwMDk1MDcgMDAwMDAgbiAKMDAwMDAwOTc0MSAwMDAwMCBuIAowMDAwMDA5ODkzIDAwMDAwIG4gCjAwMDAwMTAxMTMgMDAwMDAgbiAKMDAwMDAxMDMzNSAwMDAwMCBuIAowMDAwMDEwNDU2IDAwMDAwIG4gCjAwMDAwMTA2ODYgMDAwMDAgbiAKMDAwMDAxMTA5MSAwMDAwMCBuIAowMDAwMDExMjMxIDAwMDAwIG4gCjAwMDAwMTEzMjAgMDAwMDAgbiAKMDAwMDAxMTUyNCAwMDAwMCBuIAowMDAwMDExOTM1IDAwMDAwIG4gCjAwMDAwMTIyNTYgMDAwMDAgbiAKMDAwMDAxMjUwMCAwMDAwMCBuIAowMDAwMDEyNzA1IDAwMDAwIG4gCjAwMDAwMTU2NTkgMDAwMDAgbiAKdHJhaWxlcgo8PCAvSW5mbyA1MSAwIFIgL1Jvb3QgMSAwIFIgL1NpemUgNTIgPj4Kc3RhcnR4cmVmCjE1ODE2CiUlRU9GCg==\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 = 128 # length of segment\n",
"M = 64 # stepsize\n",
"L = 100 # total number of segments\n",
"\n",
"# generate random signal\n",
"np.random.seed(5)\n",
"x = np.random.normal(size=L * M)\n",
"\n",
"# estimate PSD by Welch's method\n",
"nf, Pxx = sig.welch(x, window=\"hamming\", nperseg=N, noverlap=(N - M))\n",
"Pxx = 0.5 * Pxx # due to normalization in scipy.signal\n",
"Om = 2 * np.pi * nf\n",
"\n",
"# plot results\n",
"plt.figure(figsize=(10, 4))\n",
"plt.stem(Om, Pxx, \"C0\", label=r\"$\\hat{\\Phi}_{xx}(e^{j \\Omega})$\", basefmt=\" \")\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, 2])\n",
"plt.legend()\n",
"\n",
"# compute bias/variance of the estimator\n",
"print(\"Bias of the Welch estimate: \\t\\t {0:1.4f}\".format(np.mean(Pxx - 1)))\n",
"print(\"Variance of the Welch estimate: \\t {0:1.4f}\".format(np.var(Pxx)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise**\n",
"\n",
"* Compare the results to the periodogram example. Is the variance of the estimator lower?\n",
"* Change the number of segments `L`. What changes?\n",
"* Change the segment length `N` and stepsize `M`. What changes?\n",
"\n",
"Solution: When comparing both the estimates of the PSD in the previous periodogram and above example, it is obvious that the variance of the Welch estimator is lower. Increasing the number of segments `L` lowers the variance further. Increasing the segment length `N` increases the total number of discrete frequencies in the estimated PSD. Since in above example the total number of segments is kept constant, the variance increases. Lowering the stepsize `M` has the same result, since the total number of samples is reduced for a fixed number of segments."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Evaluation\n",
"\n",
"It is shown in [[Stoica et al.](../index.ipynb#Literature)] that Welch's method is asymptotically unbiased. Under the assumption of a wide-sense stationary (WSS) random process, the periodograms $\\hat{\\Phi}_{xx,l}(e^{j \\Omega})$ of the segments can be assumed to be approximately uncorrelated. Hence, averaging over these reduces the overall variance of the estimator. It can be shown formally that in the limiting case of an infinite number of segments (infintely long signal) the variance tends towards zero. As a result Welch's method is an asymptotically consistent estimator of the PSD. \n",
"\n",
"Note, that for a finite segment length $N$ the properties of the estimated PSD $\\hat{\\Phi}_{xx}(e^{j \\Omega})$ depend on the length $N$ of the segments and the window function $w[k]$ due to the [leakage effect](../spectral_analysis_deterministic_signals/leakage_effect.ipynb)."
]
},
{
"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
}