matlab/sparsity_9b_music.ipynb

{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Sparse Spikes Deconvolution with MUSIC Algorithm\n", "================================================\n", "\n*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_matlab/) for details about how to install the toolboxes.\n", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n", "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n", "$\\newcommand{\\CC}{\\mathbb{C}}$\n", "$\\newcommand{\\RR}{\\mathbb{R}}$\n", "$\\newcommand{\\EE}{\\mathbb{E}}$\n", "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n", "$\\newcommand{\\Ww}{\\mathcal{W}}$\n", "$\\newcommand{\\Vv}{\\mathcal{V}}$\n", "$\\newcommand{\\Nn}{\\mathcal{N}}$\n", "$\\newcommand{\\NN}{\\mathcal{N}}$\n", "$\\newcommand{\\Hh}{\\mathcal{H}}$\n", "$\\newcommand{\\Bb}{\\mathcal{B}}$\n", "$\\newcommand{\\Ee}{\\mathcal{E}}$\n", "$\\newcommand{\\Cc}{\\mathcal{C}}$\n", "$\\newcommand{\\Gg}{\\mathcal{G}}$\n", "$\\newcommand{\\Ss}{\\mathcal{S}}$\n", "$\\newcommand{\\Pp}{\\mathcal{P}}$\n", "$\\newcommand{\\Ff}{\\mathcal{F}}$\n", "$\\newcommand{\\Xx}{\\mathcal{X}}$\n", "$\\newcommand{\\Mm}{\\mathcal{M}}$\n", "$\\newcommand{\\Ii}{\\mathcal{I}}$\n", "$\\newcommand{\\Dd}{\\mathcal{D}}$\n", "$\\newcommand{\\Ll}{\\mathcal{L}}$\n", "$\\newcommand{\\Tt}{\\mathcal{T}}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\al}{\\alpha}$\n", "$\\newcommand{\\la}{\\lambda}$\n", "$\\newcommand{\\ga}{\\gamma}$\n", "$\\newcommand{\\Ga}{\\Gamma}$\n", "$\\newcommand{\\La}{\\Lambda}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\Si}{\\Sigma}$\n", "$\\newcommand{\\be}{\\beta}$\n", "$\\newcommand{\\de}{\\delta}$\n", "$\\newcommand{\\De}{\\Delta}$\n", "$\\newcommand{\\phi}{\\varphi}$\n", "$\\newcommand{\\th}{\\theta}$\n", "$\\newcommand{\\om}{\\omega}$\n", "$\\newcommand{\\Om}{\\Omega}$\n"], "metadata": {}, "cell_type": "markdown"}, {"source": ["This numerical tour explores the use of the\n", "MUSIC algorithm to perform sparse deconvolution.\n", "\n", "\n", "The MUSIC algorithm was introduced by [Schmidt](#biblio).\n", "It is closely related to Prony's method [Prony](#biblio), and is very popular\n", "in signal processing [KrimViberg](#biblio), where it is mostly used for\n", "line spectral estimation (i.e. find locations\n", "of Diracs in the Fourier spectrum), and we re-formulate here as a problem\n", "of finding Diracs' over the temporal domain.\n", "\n", "\n", "We follow here the exposition of [LiaoFannjiang](#biblio) who propose a\n", "theoretical analysis of the performances of this method.\n", "Several related algorithms exists, see for instance\n", "[RoyKailathHuaSarkar,DemNeedNg](#biblio)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 2, "cell_type": "code", "language": "python", "metadata": {}, "input": ["addpath('toolbox_signal')\n", "addpath('toolbox_general')\n", "addpath('solutions/sparsity_9b_music')"]}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ms = 20; % markersize"]}, {"source": ["MUSIC Algorithm\n", "---------------\n", "We consider here the problem of recovering a Radon measure $m_0 \\in \\Mm(\\mathbb{T})$ defined on the torus\n", "$\\mathbb{T}=\\RR/\\ZZ$ from\n", "of possibly noisy low-pass measurements\n", "$$ y = y_0 + w \\qwhereq y_0 = \\Phi m_0 $$\n", "where $w$ is the noise term,\n", "and so that $ \\Phi m $ computes the first $M+1$ low frequencies of\n", "the Fourier transform of $m$, i.e.\n", "$$ \\forall \\ell \\in \\{0,\\ldots,M\\}, \\quad (\\Phi m)_\\ell = \\int_{\\mathbb{T}} e^{-2\\imath\\pi x \\ell} d m(x). $$\n", "We only consider positive frequency because we assume $m$ is a real measure.\n", "\n", "\n", "In the following, we consider a discrete measure of the form\n", "$m_0=m_{a_0,x_0}$ where we used the following notation\n", "$$ m_{a,x} = \\sum_{i=1}^N a_i \\de_{x_i} $$\n", "where $a \\in \\RR^N$ and $x \\in \\mathbb{T}^N$.\n", "\n", "\n", "We thus have $\\Phi m_{a,x} = \\Phi_x^M a$ where\n", "$$ \\Phi_x^{M} = ( e^{-2\\imath\\pi x_j \\ell} )_{0 \\leq\\ell \\leq M, 1 \\leq j \\leq N}$\n", "  \\in \\RR^{(M+1) \\times N}.\n", "$$\n", "\n", "\n", "Note that $\\Phi_x^{M}$ is a (rectangular) Vandermonde matrix, since\n", "$$ \\Phi_x^M = (V_{\\ga(x)}^M)^* \\qwhereq \\ga(x) = ( e^{2\\imath\\pi x_j} )_{j=1}^N \\in \\CC^N. $$\n", "\n", "\n", "For $y \\in \\RR^M$ and for $0<L<M$ we denote the Hankel matrix\n", "$$\n", "\tH_y =\n", "\t\\begin{pmatrix}$\n", "\t\ty_0 & y_1 & \\ldots & y_{M-L} \\\\\n", "\t\ty_1 & y_2 & \\ldots & y_{M-L+1} \\\\\n", "\t\t\\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "\t\ty_L & y_{L+1} & \\ldots & y_M\n", "\t\\end{pmatrix}$\n", "\t\\in \\RR^{(L+1) \\times (M-L+1) }.\n", "$$\n", "\n", "\n", "Since $y=\\Phi_x^M a+w$, one can check that one has the following factorization\n", "$$\n", "      H_y = \\Phi_x^{L} \\diag(a) ( \\Phi_x^{M-L} )^* + H_{w}.\n", "$$\n", "\n", "\n", "We suppose that $\\min(L,M-L+1) \\geq N$. The fundamental property, which\n", "is based on the above factorization, is that\n", "$$\n", "\ts \\in \\{x_{0,1},\\ldots,x_{0,N}\\}$\n", "\t\\quad\\Longleftrightarrow\\quad\n", "\t\\phi_L(s) \\in \\text{Im}(H_{y_0})\n", "\t\\qwhereq\n", "\t\\phi_L(s) = (e^{-2\\imath\\pi \\ell s })_{0 \\leq s \\leq L} \\in \\CC^{L+1}.\n", "$$\n", "\n", "\n", "From observations $y$, we define the following detection function\n", "$$\n", "  \\forall s \\in \\mathbb{T}, \\quad\n", "      d_y(s) = \\norm{ U^{\\bot,*} \\phi^L(s) }^2,\n", "$$\n", "where we considered the following SVD factorization\n", "$$\n", "      H_y = [U, U^\\bot] \\diag_j(\\si_j) [V, V^\\bot]^*\n", "$$\n", "$$\n", "\t\\qwhereq\n", "\t\\choice{\n", "\t\tU \\in \\CC^{(L+1) \\times N}, U^\\bot \\in \\CC^{(L+1) \\times (L+1-N)} \\\\\n", "\t\tV \\in \\CC^{(M-L+1) \\times N}, V^\\bot \\in \\CC^{(M-L+1) \\times (M-L+1-N)}.\n", "\t}$\n", "$$\n", "\n", "\n", "The detection property can hence be conveniently re-written using this\n", "detection function as\n", "$$\n", "\t\\{x_{0,1},\\ldots,x_{0,N}\\} = \\enscond{ s }{ d_{y_0}(s)=0 }.\n", "$$\n", "\n", "Sparse Spikes Deconvolution\n", "---------------------------\n", "We compute here some low frequency measurements.\n", "\n", "\n", "Number $M$ of recorded frequency."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["M = 18;"]}, {"source": ["Imaging grid, for display purpose."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Q = 1024*4;\n", "z = (0:Q)'/Q;"]}, {"source": ["Input positions $x_0 \\in \\mathbb{T}^N$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x0 = [.1 .4 .5 .7 .9]';"]}, {"source": ["Number $N$ of spikes."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["N = 5;"]}, {"source": ["Input amplitudes $a_0 \\in \\RR^N$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["a0 = [.7 -.8 .9 1 -.9]';"]}, {"source": ["Observation operator at position $x$, i.e. $\\Phi_x$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Phi = @(x)exp(-2i*pi*(0:M)'*x(:)');"]}, {"source": ["Observations $y=y_0+w$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = .15;\n", "w = (randn(M+1,1)+1i*randn(M+1,1)); w = w/norm(w);\n", "y0 = Phi(x0)*a0; \n", "y = y0 + sigma*norm(y0)*w;"]}, {"source": ["Compute the observed low frequency signal\n", "$$\n", "     f = \\frac{1}{2M+1} \\sum_{\\ell=-P}^P y_\\ell e^{2\\imath\\pi \\ell \\cdot}$\n", "$$\n", "where we completed the observation using $y_{-\\ell} = \\bar y_{\\ell}$\n", "since the input measure is real."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["PhiT = @(x)exp(2i*pi*x(:)*(-M:M)) / (2*M+1);\n", "f = real( PhiT(z) * [conj(y(end:-1:2)); y]);"]}, {"source": ["Display the observed signal $f$, together with the spikes locations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARgUlEQVR4nO3d25ab\nOAIFUJiV//9l5oEUcdsum6ukI/Z+mJVJVadkkHR0gxqnaRoAIM3/ahcAAPYQYABEEmAARBJgAEQS\nYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJg\nAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAET6U7sAO43jOE3T618uf379KgA9\nyQuwx5R6JbcAbiJvCXGapg8pNY7j54QDoA95M7DP5mxbFhhfw8wUDaAPXQXY23CSWABdyltC/I2V\nQ4Bb6WEGNi8YTtO0ZJhZF0D33hxG78nb0/YAdKCfJUQAbkWAARBJgAEQSYABEEmAARBJgAEQSYAB\nEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQ\nSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJ\ngAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmA\nARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAENGMfa\nJSBPhwE2agkQZByHcRym6e8fYLU/tQtwJtEFYebo+jsBm6bp39/AV10F2PS3Jfwnxp7+76RtQGPG\ncZjm/4EtugqwtyQWpBBibNLhHhgQY976+vnjNFg/ZAMBBlQ2DeMwOIjIZgIMqGkcpr+zrmmygsgm\nHQaYTS+AO+gwwAC4AwEGtOLhSAd8J8AAiCTAAIgkwACIJMCgnttv+HjxIUf0/yopaNEcXcuhBb04\nbGcGBsXN846fp3drlwZSCTCo52keBmwhwACIJMCgNicZHpiOsp5DHFDcYyett4a9BBjUME+5zL3g\nAAEGFfw9vVG7GBDNHhjUYeoFB5mBXW7Z49BhMbPtNbOAykFmYCVopTxRJeA4AQY1mY3BbgKsHF0V\nwIkE2LWWVX5LRgDnEmBQlIk4nEWAQWmm4595mxQrCTAAIgkwACIJsHIsjPBKrYDdBBgAkQTYhbwp\nB35j3slxAgyow/COgwQYlGPaAScSYFCUaQecRYABEEmAARBJgEF99sZgBwFWlKdWeWVXDPb5U7sA\nvXvKK30VrDAP9TQXPhNgVxnHYRrGYXgIrXEcxnEYNEqAE1hCvNjjGNJ4EuA8AqwIG18AZxNgAESy\nB3alx0OH489+mMnYXTmVsLAkwSkE2MXmHkvXBf+lQXCcJUQAIgkwaIJVNdhKgEF91tNgBwFWmrdJ\nAZxCgAEQSYBdwhwL4GoC7Cp2NQAuJcAAiCTAgBY57sRXAgyASAIMgEgCDErwOkw4nQADIJIAgybc\n58yCyShnEWAARBJgAEQSYBXcZ7EI4DoCDIBIAgyASALsfA5ZARQgwACI9Kd2ATYbf84/TC/TnPHh\naMTrVwHoSViAjeO4JNPjnxdyC+AmeltCHMdxdEQduuCBEz4Lm4F9Nc/AniZqr98AQLquAuxtOEks\ngC71s4Ro5ZBmraybVsxgk7AZ2DRNr6cQ5wXDt1+CRqiScLqwABvehdPyN3IL4D76WUIE2uc9NZxI\ngNVhtwPgIAEGQCQBBkAkAXYyC4MAZQiw89mjBihAgAEQSYABEEmAARBJgAEQSYBBWxxkfeSRfz4Q\nYNAQR1hhPQEGQCQBBteyAgYXEWBwOQuDcAUBVo3dae7G71LhXAIMgEgCDIBIAgyASAIMgEgC7Ez2\nqAGKEWAARBJgAEQSYNAWDwjCSgIMgEgCDIBIAgy+s6YHDRJg8MWcXjKsFpuCn43jfa+PAIPvPN5H\n4+6ZYQKsJkPL7rm/XGquYLcdYAkw+OT421Vu27k8keUXuXMFE2Cwli74oDt3tQXcsH4KMFhF50vL\n7lk/BRhApBtOuZ4IMIBU95x4LQQY/KrWCNfx1J7Mz2m5oVcQYKdRQWu5tHe4+QiXU6hFFxFgZ1JN\ny1uiywCCxqmipxNgxJumQkMHK3tstVSYAlX0hvVTgAEdamfnKWhhpoXLtYkAI9hTe4trflxqTg61\nYqXEl1YLMLItw9ugcW6ulN6t5MJdLcdfcvb0rw2BkS/AgFUqjtB37O40GF2N71HNV6zB6/aBAIvU\nyOI+97GM0BvvhV9VLHDWhUokwCrb0bocHJ+du4RS/t+Pc+Rq3PZiRnzq17uT0rcIsEjFDo7D8NKd\nNT4Ja7lsVzvl1gT1LQIszJ0bZyPW34Ij047WQiKoUxvSSstuAizP47m7pvq46gpcED1jIi2lVwIM\n+OTtPFIkdOO3dYKI+yvAkkRUqTJcChr3WxWNqLopKw3dB9gUUV3We61YFT/g8raeWs8GwZOm2vvb\nmSsn6j7AClXnWqeEK7aHxweDhsY6Ds6SeFuFxCkiLmP3AcYlHl88M4TUdfb5cHOz4m1fadt8aUCD\nRapCgLFTU4uZlJc1atlX2pZfGvDbJ9p0vqbBz7WJAIvxYZVydy3ct4n19ptLdmfprY4g3b80IPrT\n3SLA+u7vdte/I5tYp0fpVqcMPz/ou86s1NT7n8of3G/q4/NW/wHWfhWs+EjN0ybWmmJ8+J72L/Um\nV2fkmgL0kaPFPsV1kbP1Xgi/MvoPsK3a3LM93etnXN/etExm67v1uDpzh06gAwLsP241btrx+P3B\nKRpB+ruPW2dR9+kKcgmwfx4rd2ut93N5zlprWtNiP3+PNt8Td/NJa90CdwmwTQsdbbbbc399+I6H\nexppvbeaJdOOs2pdO1ubHTSlWwTYmpvUwb08xSkTrEba5w3VfblXN9KvXjsZebVbBNgO96kBb71+\n9vVXwzhgVn5IdOLLvY4/F1jeucVQjSN0FWDjj9oFadrXjvW37k+Tbtnbl3sdaQorb3dTtaKRwuiB\niuknwMZxnH5szbCDg+VG6uuJxXi6Gjv+5Uauya083bVGevOmHFxZ6fixgVDj1MuVngPs6c9mYwD7\ntJ8Of2oX4HIPqfbpBUi/PRRV5gDI12OBZxXjaaHpOmt+0OP3nHupv3/bt+84fkfO2gM7+LP2FWPT\nf7XmUOtZl+LDP3VW3Shw6798w7HKub4MQ/5MsZ8lxN1aOH9Y8sxImQ+7NUpvfmrmUlsv7Okt4sQD\n6Jyog+spwKjgt9FfSoY1Fbef86ZAJ9XU1eBW+gmw+ezG7O3S7dtm1lnDa+3jbHrz78pJWwGtXUa+\n2rrmuVtQWqeU84h+AmwYhuUU4sb/6tNXgypBsV2Zlfb9ypWvL80qI2V1ZU393NTnBlX4fQ7tTuXo\n41N81VWAne4mleBSTx3i561jF3yH0y+au3BENxEY4V4B9jQUjRhsRreHp8dp1xx8ClqiKeOsC+Kq\n0p97Bdir3GxIsWTYwWO70UFe3dWXzrBjk+qXq5vWdMcAm6tOl+3tQ8OoWGXnt/Mt7+hbo8u7c7qt\n9/TrVe2mX/vq7aVY+fGrxw+L2wXY46LWTdpqHPflClWu6hWJeDw8VLBu3C7AhmHbbODraOussZhh\n3ZO43cpHBaYy+67Jh/8q4iL/9sYc7umOAXa667qqszqpuBb+OFHeMV3ubDTw28fZWvG+fn/u1KSR\nkjdV6zprBW8JsNad1Uk10sLXs9h7kes6tazu8u0DHsepqyUJMNq19ehHSXWHt7t/9L6ny4/8y216\nW9r1HyF6ctPTUR0B1qfoDSTeerqPR/qgUxYkqaKn+DlOgHVox2sGadzjHTw4Inl6uvz4P9gBVyCU\nAPuu2HLBFS8K0TJ7stzNgyOSx6fLr9tlvG6u8NRM9v2Up+ZmkJdIgDVtd6OaG6cTED15nTkd/9eG\n4cxdxjJDvXPr89JM9v23T/+3wbbW9xC2/9/IfFsNtiUOavP3TOZaEnfHpYg4x/FayPbLvIkAA+7r\nhine00e2hHjIFcOZjl/VGKHNhaAgJSuwm3VzAuyo61Z1tEziFKvALYzwHhfoRGkVAmyVwuvdLTRO\n2GduLJf25ueeZ6miZJfy9Av5emIPrDnGcSl+66n76ya2KlaHW2gsjd/uvp8WMAOD83XWTTSokXeM\nPT6NQHlmYPSpwEIWDCHR1fg0cTczMICedTxNNAMD6Fx/0TUzA1ur+2faAbIIsEN6HdcAtE+AAWRw\nLumJAAMgkgDbxr4XjyJeSQ69EmAbmLzzlTyDYgTYThaj+Y2KAWUIsM0MsQFaIMC26eA12AB9EGA7\nWSYCqEuAbdbIa7C5grk1xTjCepwAg//YOjrRDUEtAoxulY8WSQYlCTA4k+VlKEaAARBJgMEJLB5C\neQIMjvJ0IFQhwOA0NsCgpD+1CwA9EF1QnhkYAJEEGEAAvwHjlQADIJIAAyCSAAMgkgADIJIAAyCS\nAAOoxgtcjhBgAHU4Fn+QAIO/PGcDWQQYPfPrkqFjAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyA\nSAIMoHWesn9LgAEQ6U/tAmw2/rxZYXoZkIwPL114/SoAPQkLsHEcl2R6/PNCbgHcRG9LiOM4jl5+\nB3ADYTOwr+YZ2NNE7fUbAEjXdIBtzZ633yCxgGbNvzBBL7VP0wG2KXvebokB0KumA+zVNE2vpxDn\n6Hr7JRjWPUNj5xTihAXY8C6clr+RW7xa/zstVR/I0tspRABuQoABEEmAARBJgAEQSYABEEmAARBJ\ngAEQSYABNM1T9r8RYACt85T9WwIMgEgCDIBIAgygMrtc+wgwgJrsb+0mwACIJMDo39ffqGIBBxIJ\nMBgGyzgQSIABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAAbRrHD2k+CsBBkAkAQZAJAEGQCQB\nBlDZ1/dN85YA4xZ0ENAfAQZAJAHG3TmmDKEEGACRBBgAkQQYAJEEGPczn0d0KhHCCTBuZjmz4WQ9\nhBNg3MnTiUMZRtsckf1MgAEQSYBxS4a1kE+AcSdPa4YWaGiJ9eyt/tQuABQyh9f0mGHSi2bYkN1B\ngHE/cgu6YAkRgEgCjFuzCwa5BBgAkQQYAJEEGACRBBhAi2zQfiXAAJrhVyVsIcAAmjAND4/Yy7AV\nBBhAK8bhZ9Fwmv7lGb8QYNyIcS3tsuW1nQADaMAyvBJjqwkw7suQF6J5mS9AG6ZpGsdl58v46isB\nBtCK+Rf+DJYH1rGECNAK54w2EWAARBJg3IsRLu1TRVcSYNyUPoKW2QBbwyEO7ksfQYNUy/XMwEoY\n2xjtt1CMFsowKEZjZRgUo7EyDM0U4zMBxu1M0zAMRrkQr8MAixg4UNtooQbSdbUHJroA7mOcuhuI\njuO/DyXSAPZpPx26moG9av8GALBPaoA9Ta0EFcDdpAaYxAK4uQ5PIQJwBx0e4gDgDlKXENdY9snK\nh/TnH/14TrJKGUpemd9+VuG7k3JH3n61TDFauCOF97a/1ooCZfhQjIo92FKAxmc43QbY02H6krfh\nw48udqz/68ef/+bqK/O5GGXK8LUYZW7KyjtSqxiPnXXFO1KysXytFXXbSMUebMh5AMkeWFHTNLUw\nolGGRTtjzHEcq/cacxlauCCNFOO2Gumpvup2BsZX1fuI6v11U4rNfj4XoG4ZmlJsCXFNMXhLgN1R\nIy1z6bJrFWD+0cv/Vrwg1e9FU6onaN3lu9k0TUs7lWG/sYR4U9U7iIo/fTH9GKpekEauBk2Zg9PI\n5rOe1wraOYX4epSjytGm+efWPej1uDdepgCfi7F8teIdefulYsVo7Y6UnPG4I18L1nhAtF4+AHjL\nEiIAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEE\nGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQYAJEEGACRBBgAkQQY\nAJEEGACRBBgAkQQYAJEEGACRBBgAkf4P7Bgx3JYerwQAAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "plot(z, f);\n", "stem(x0,a0, 'r');"]}, {"source": ["Noiseless Recovery \n", "-------------------\n", "We first check the exact recovery when using $y_0$.\n", "\n", "\n", "We set here $L=M/2$, which is a standard choice."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["L = M/2;"]}, {"source": ["Build Hankel matrix,\n", "one must have min(size(H))>N."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["MusicHankel = @(y)hankel(y(1:L),y(L:M));"]}, {"source": ["Compute SVD of $H_{y_0}$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["[U,S,V] = svd(MusicHankel(y0),0); \n", "S = diag(S);"]}, {"source": ["Display the decay of the eigenvalues of $H_{y_0}$. Since there is no\n", "noise, only the $N$ largest one are non-zero. This is a convenient way\n", "to detect $N$ if it is unknown."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAM9ElEQVR4nO3d0VLr\nOBaGUWuK939lzYUPbpOEEIIT699aq6a6OGf6QhXAX7Ql6NZ7XwAgzf/OXgAAPEPAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRg\nAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIn2cvYAXau2/j3s/bx0AvEDZHdi+\nXtd/BCBdzYDJFUB5lUeIF5qsATysD3/0UjNgvd/YhI38yWitjby8jXUeyzoPFLHIJWqdZy/hZzVH\niDe1ZrQIUEfNHdhytQnb3vGsf5nwBgiAe8oGbPnXsMvd+vonGQNIlzGNfdr9cbOMAdwUcVZXeQf2\nI7sxgFxTB2wlYwCJBOwfGQPIImBfyBhACgG7QcYAxidg35IxgJEJ2A9kDGBMAvYQGQMYjYD9gowB\njEPAfk3GAEYgYE+SMYBzCdifyBjAWSb674G9Tu+3/xOaALyOgB1mbZiMAbyHEeKRTBQB3kbAjidj\nAG8gYK8iYwAvJWCvJWMALyJg7yBjAIcrFbD2eQWwD1mJ7zK2XVwcctUAg6oTsNba1q39x6O5yNj+\n2r0tGsDj6gQsyz5jADyhYMDa1yxc/HHYnRnAuVrae+o6Aeu9r6/+9sH29+ctCiDG/mkZEbM6ARv5\n3OtX3OkAeESdgO03Xiklu/4VwG4nAjyoTsCWnG7t3Vny9n8pGcC1UgGrSskArglYEiUD2AhYpIuS\nyRgwIQHLdvED0UoGzEPAKjBaBCYkYKUoGTAPAatJyYDyBKw4JQOqErBZKBlQjIBNxxV8oAYBm5cr\n+EA0AZud0SIQSsD4537J7vzWfIBTCBiXrkt2rTUNA04mYHzrkZIBnOV/Zy8AAJ4hYABEEjB+dvO4\ny1wROJeA8ZB9w3r/9z8NA07kEgePut6HbQ1zIxF4PwHjT/a/zkPGgHcyQuQAJorA+wkYh1kbJmPA\newgYR7IVA95GwDierRjwBi5x8BIudwCvZgfGC5koAq8jYLyciSLwCtUC1lprnpTjsRUDDlcqYK21\n3nvvXcPGZCsGHMglDt7K5Q7gKNUCtu69+u7ReLEb656aA9gy5rMB44ibXdUJ2Do/vP5YsYbldwHD\nUO689R9TqTMw4rjcATytzg5sf3fDriuLrRjwhDoBW3QrmcsdwG8ZITIQE0XgcQLGcPy4GPAIAWNE\ntmLAjwSMcdmKAXeUusRBPS53AN+xAyOAiSJwTcCIYaII7AkYSWzFgI2AkcdWDFhc4iCUyx2AHRjB\nTBRhZnZgxNv/LuCLmNmcQWF2YFTw3VbM5gwKEzAAIgkYAJEEDIBIAkYd11c2XOKAwgSMUvbFcsMe\nanONnmquG2YfBiXZgVGcfRhUJWAARBIw6rMJg5IEjCloGNQjYMxCw6AYAWMiGgaVCBhz0TAoQ8CY\njoZBDQLGjDQMCqjzmzja1wdS99sXuMsv6YB0dQK2L1bz7hqguoIjxNaa7RePMEiEaHV2YN8xWuQO\ng0TYxM2uqgXsevulWNynYbCKO4gpOEKE3zJLhESlAub0i6dpGMQpFTD14i80DLKUChj8kYZBEAGD\nLzQMUggYAJEEDC7ZhEEEAYMbNAzGJ2Bwm4bB4AQMvqVhMDIBg3s0DIYlYPADDYMxCRj8TMNgQAIG\nD9EwGI2AARBJwOBRNmEwFAGDX9AwGIeAwe9oGAxCwODXNAxGIGDwDA2D0wkYPEnD4FwCBs/TMDiR\ngMGfaBicRcAAiCRg8Fc2YXAKAYMDaBi8n4DBMTQM3kzA4DAaBu8kYHAkDYO3ETA4mIbBewgYHE/D\n4A0+zl7AkdrnM6P3fu5KYG2Yr0R4nTo7sLVevffee/PulwHYh8FLldqBLbuMnb0QAF6rVMC2brXW\n9h/f/HfgDQwSCRI3uyoVsJsUi3NpGCn2T8uImNU5A4NhOQyDV6izA9vf3bDrYjT2YXC4OgFbdIux\naRgcywgR3scsEQ4kYPBWW8NaEzP4EwGDE2zpkjF4moDB+TQMniBgAEQSMAAiCRicz916eIKAwbvJ\nFRyi1A8yQ4qtYX60GZ5mBwZAJAEDIJKAARBJwOBMfjsiPE3AAIgkYABEEjA4mSkiPEfAAIgkYABE\nEjAAIgkYnM8xGDxBwACIJGAARBIwACIJGAzBMRj8loABEEnAAIgkYABEEjAYhWMw+JWPsxdwpLb7\n7u+9n7gSAF6tVMAW3QKYRrURYmutmcIQyxQRHldzB9Za27ZiFz2zRQO4Ke7df6mA3YyTYgE8Yv+0\njIhZnRFixMsNwFHq7MB671vD7LrItR6D+RKGH9UJ2KJbADOpM0IEYCoCBkAkAYPh+GkweISAARBJ\nwACIJGAwIlNE+JGAARBJwACIJGAARBIwGJRjMLhPwACIJGAARBIwACIJGIzLMRjcIWAARBIwACIJ\nGAzNFBG+I2AARBIwACIJGACRBAxG5xgMbhIwACIJGACRBAyASAIGARyDwTUBAyCSgAEQScAAiFQw\nYM1ZARU5BoML1QKmXgCTKBWw1lrv/exVAPAOH2cv4OUu9mQKR651iuhLmBeJm2DVCdj60m//3EKl\nWACP2D8tI2JWJ2DbS2+QCDCDUmdgAMyjYMBsvyjMZXrYFAwYADMQMAAiCRgAkQQMwjgGg5WAARBJ\nwACIJGCQxxQRFgEDIJSAARBJwACIJGAQyTEYCBgAkQQMgEgCBkAkAYNUjsGYnIABEEnAAIgkYABE\nEjAI5hiMmQkYAJEEDIBIAgbZTBGZloABEEnAAIgkYABEEjCI5xiMOQkYAJEEDIBIH2cv4Ejtc4zS\nez93JQC8WrUdWO+9994cCDAZx2BMqFTAbLwA5lFqhLjspojf/Y3IAdwUN7uqFrC1T/tPg2IxiXWK\n6Oudp+2flhExqzNCjHi5AThKnR3Y/u6GXRdAeXUCtugWwEzqjBABl+mZioABEEnAAIgkYABEEjAo\nxTEY8xAwACIJGACRBAyqMUVkEgIGQCQBAyCSgAEQScCgIMdgzEDAAIgkYABEEjAAIgkY1OQYjPIE\nDIBIAgZAJAEDIJKAQVmOwahNwACIJGAARBIwqMwUkcIEDIBIAgZAJAEDIJKAQXGOwahKwACI9HH2\nAo7UPt9n9t7PXQkAr1ZtB9Z77703ExOA6koFzMYLbnIMRkmlRoir1tq+ZBe7MZEDuCludlUqYOur\nf5EoxQJ4xJ23/mMqNUJc5Aq+YYpIPXV2YOv7BRcRASZRJ2CKBTCVaiNEACYhYDALx2AUI2AARBIw\nACIJGACRBAwm4hiMSgQMgEgCBkAkAYO5mCJShoABEEnAAIgkYABEEjCYjmMwahAwACIJGACRBAyA\nSAIGM3IMRgECBkAkAQMgkoABEEnAYFKOwUgnYABEEjAAIgkYzMsUkWgCBkAkAQMgkoABEEnAYGqO\nwchVMGDNtyPABD7OXsCRpAtgHqV2YL333vvZqwDgHUrtwG662JYpHFxYj8F8ZxA3xKofMMUCeMT+\naRkRs1IjRADmIWCAy/REKhgwM0OAGRQMGAAzEDAAIgkYsCyOwQgkYABEEjAAIgkYAJEEDPjHMRhZ\nBAyASAIGQCQBAyCSgAH/cQxGEAEDIJKAARBJwIAvTBFJIWAARBIwACIJGACRBAy45BiMCAIGQCQB\nAyCSgAEQScCAGxyDMT4BAyCSgAEQScCA20wRGZyAARBJwACIJGBDaCGTGus8lnUeKGKRS846IwgY\ncE9ry7L0wZ+6EYtcrPNoAgbcdvH8GvZxtl9Ya9b5VxfrHJmAAcEGf8Juctc58spb7/3sNbyQcTM8\nq/KTgV8ZthIfZy/gtWrnGV7q+u3fgN9PN9+jWufTst7zGyECwQZsAG8jYMBtF21IScWY67xelXX+\nXfEzMACqsgMDIFLxSxzLsrQ2+i5zuyppnQca/PO+vx8bsc6IRa5GXuryudqIRS7Dr7NywILu0K9f\nJYM/c5ecdUZ86gd/DZevj9qRP+n7hQ3+qd9expFfz/3aRl7nUnuE2Hsf+aXfRCxyyVnn4N9ym9ba\n4E/b1brOlJc0Yp0cpfIOLEvE917EAzfF+O/El927lsHXGSRihLga//tdwM4X9AW9PXPPXsi31rVt\n/xz2VR12YblG/nSvUkZzvfftoTTyN/tSe4QYZNgv5c3gX8eb/mkZ+FVNeTGZ0xrXYb999sZ9F3CU\nkd/pLFEXqFIuJq1SPu8jL3KJWufgK1wFvZ7rB6Ovc/D1AcBNRogARBIwACIJGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAi/R/M\nPkLBFnaMPgAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "plot(S, '.-', 'MarkerSize', ms);\n", "axis tight;"]}, {"source": ["Imaging function $d_{y_0}$, evaluated on a thin grid $z \\in \\mathbb{T}^Q$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Ubot = U(:,N+1:end);\n", "d = Ubot'*exp(-2i*pi*(0:L-1)'*z(:)');\n", "d = sum(abs(d).^2) / L;"]}, {"source": ["Display $d$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAASb0lEQVR4nO3d63Kr\nuLYGUHQq7//KnB902CwgGNuANKfGqF1dWfbaHQWEvqmL02UcxwEAovm/2g0AgE8IMABCEmAAhCTA\nAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJ\nMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAgpYYCVUmo3AYDb/dRu\nwJVEF0A/Us3AxnEcx7F2K/ZI1ta0cEdKaaIZjXApZjrGaaXREf8Lpfzvh2phTjZf3/pNYRiGZu7I\n2EAbGtHIHWlEOx2j/XRItYS4q4l7UEoL7VhGe89tmNrRwh0ZGrga7si/rWjgajRzKWo34bVUS4gA\n9EOAARBSwgBrYPINwO0SBhh/aSHaW2gDS+7IkqsRiwADICQBBkBIAgyAkAQYACEJMABCEmAAhCTA\nAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACH9\n1G4AqZSyfmUca7QD6IAA4zJTei0Tq5ShFBkG3MISItfYptf8x+20DOB7ZmB8aze6ZuMowIBbmIHx\nleP0msgw4A4CjM+dSa/VXwa4SrwlxPI7EI6bgfPgLS73VnqZhAGXCxZgpZQ5nJZfH7/F5d5KL4A7\nBAswrrI7H3prMfDd9JomYTIPuErCACv/js2rP5qZ/RU/02e2Jn9dpJd/AYirRFvozxNg4zhOV3/+\nYn69XqOaczB5ml/cJtmyV395OU3CoFmrTZmKLTkpT4DZ9zrp5UU6SLLvv3WEhwKIIViALWdXqyMb\nu2+x9G54uIpAy4IF2LAXTvMre7k1Oi+35DqwYl9z5lLMpt2Y2q14Lf0HmYtfxzdpZPPJKmJT5vLO\nfVlWup1fil8BrkL6ABsGnRL+tvqvB/RsuhSGi0aK3TO6CLCh+zWBnp9G/rLqFT0/Iy5FUL0EGE09\nkwK1EU31irp2PxnZoUDTr6GrAOt2lb+1nzrQ49GbPp+Rb34rDXV1FGA98zSyEqvQvptLMQlXvggw\ngH3hBvTvxcryvgKswxWS3n5evtTbM3Lww8YayvvUV4D1qcHnsLdRsjUu/lKDDwgnCTDokVGblYjb\noj0GmPITIIHuAixcifENUc0H+lngfflj9nMpguouwHrTcmAbGqiu5QeElwQYdRg4aom41cHdgvYK\nAQZASD0GWCfr2kFLKmhKJ8NFUD0GGHCsh1FbhZeAAAMgJAFGNT2U+a0x7WArbq8QYDkJBiC9fgMs\n/RAftKSCBqUfLoLqNMAM7vCSUXtiuGhWpwEGHMs9asfd9WFJgFGZMh8qCp3lAoya4j45EYUeqmBL\ngCVknAJ60HWAWbwCzvCZxTb1G2DmKPBSylE75Q/Vp5/aDXhb+e19414ETe/uvgW8JfG0wwiRQ7AA\nK6XM4bT8evXK9i1aZtMOqoj+6AULsA+Uf2vI9MEWrkcmLvOb4iLzUonWS7IF2HYJMX1iwUkeBY4t\nR8sQYZbnEMe0bDg5f+kj3COgCYaL1uQJsA8oSIGTDBcNCraEuJxdrU5z7L4FfCncruqxZD9O54IF\n2LAXTvMrcguu5YgNLet6CTGfuNWlURIeFne4mAkw6ov+FLVPfUBKAgy6oEogHwGmOAUIqfcAU5YC\nJznS0preAwzoR4JjCywJsDzUhtzBtINmCbBUVJfAGTmKEgFGE5T58LAE9a4Ag+Rs/JCVAAMgJAFm\n8QogJAEGcJZ6tykCLAn7HHDMM5KPAANeMO2gTQIMoC9pJqMCjFYo84G3CDAAQhJgkJlJLYkJsGGw\neEVqOXY7YEuAAbxHvdsIAZaBxwmOXXjuzoy2HQIsCQ8V0BsBBrxmnziNNB8CGwQYTTFKAucJMABC\nEmAAhPRTuwFvK79rTOO/67jl37WnMc0qL3wq024HbAULsFLKnEzLr4d/E6vYSAHILuES4irYTop7\nfECVDQ+LO1wkE2wG9gFLi9A5YXNSuLWrbAG2nX5JLMAwcEa4jZiES4hEF+HB6ZF1swSS7TgEm4GN\n47g9hTjPuj7b/aIpRkngpGABNuwtCc6vSC+AflhCBCAkAQY5WYklPQEWW7ItWa6lb5CbAPsfxweA\nkwwXLRBgQGZWKRITYACEJMBojsUZuEO+yagAA85SW9AUAQZASAIMgJAEGAAhCbDA8m3JchU7VfRA\ngP3DHjVpKG5IT4ABad1aj6p3qxNgQGZmookJMID8Uk4WBRgtsjgDl8s3GRVgwBvUFrRDgAEQkgCL\nShUMdE6ABZZvRRvgPAEG2Zid0wkBtmaPmgTMzumBAANyeqASVe/WJcCAtMxEcxNgAMllnSYKMBpl\ncQYulHIy+lO7AW8rv6PauLkhB28BV5lqCw8Z1QULsFLKHE7Lr4ff9JpeWb0FQD7BAuylZYwlpv4F\nSBVgu5Oz8u9GSvpsA/hMibbtnCrAdkksgDO2mzKNcwpxh/NvxGV5mX4Em4GN47g9ajgtGO6+BUBW\nwQJs2Aun+RW5BUwem4n6UEFFlhABMku8ISLAaJfNSLhE1gmiAAM+obagOgEGvC1rRU8sAgyAkARY\nPI48AQwCDICgBBgAIQmwfQ5wE5FOS1cEGKRif3SosU+sdKhCgAF8peWiIXeyCjCaZi0XvtRyvn5J\ngAEQkgADICQBFoz1NBphdZfqBFg8iVe0Ac4TYACEJMAACEmAARCSAPuTPWpi0V3pjQCDPBzwGSr9\n94bUu1UIMFpnaIDPpP9vBwowAEISYACEJMAACEmARZJ+RZtwbE9SkQADPqScoq6f2g14W/kt+cbN\n01MW1eD2XQAyCRZgpZQ5mZZfz+QWQCeyLSGWUopVeYAOBJuBvTTNwFYTte1fOP1vc26CGJRtE9fh\nG+Gq/1QBthtOFhVzUEm85PpMal2HBPXucrQMEWZ5lhBDXG4+E3pQgCp6GBGDzcDGcdyeQpwWDHff\nAuhW+oEwWIANe+E0vyK3APqRZwkxvejL6wDXEmDAV3rYa6FNAgz4nFUBKhJgAIQkwAAISYABEJIA\ne2H6dD20TBedOKnbGwEGGRi4W6DefZgAIwZDA5zXyWRUgAEQkgCLweQDYEWAhdHDggDAeQIM+JYV\nAqoQYMBXrA1QiwADICQBBkBIAgyAkATYaz5CC+3r5KO7LAkwwlBJ7DJws9LPYyLAAC7TSJnVSU0j\nwAAISYAFYI0IYEuAARCSAAMu0MLGD70RYMC3LHFThQADICQBBkBIAgyAkOIFWPl18Bcu/6aNfDgR\n2OWjJn36qd2A95RSxt9+uvx6+RcebxTPmSoJQ9XM1WClqyEwWIAdmyJtlWGrP24zD+BC1cusj791\nuAlAqgDbJbEAzliOliHCLN4e2F+my738JwCJ5ZmBHe+NxWWTA2BXsABbbnFlTSwISrHFw4IF2LC3\np7V6RZjB83zUhOfl2QMDoCsCDICQBBgAIQkwgrHXwor+0C0BdpZxk9bokDMntya9HQQVYBBYV6NV\nIOrdZwgwAEISYE1TxAH8RYC1zhoRwC4BBkBIAgy4hpMLPEyAARCSAAMgJAFGPJaqmOkJPRNgQGxO\n6nZLgL1B4U87dMX2PXyPevs9UoMAg7h6G61icXceIMAACEmAtavDBQGA8wQYACEJMABCEmDAZZzU\n5UkCDICQBBhRqfShcwKMkJzPZHBSt3sCDICQBNh77FHTAjOPEJ4cLvrsEj+1G/C28tsjxs3tOngL\ngGSCBVgpZQ6n5dez6ZXdt2IxzwM4lmoJMXporeT6aQAuFmwG9lLZzFxWryQLOYCrbMfPxmULsHkJ\ncfUK8Jg+DxQksBwtQ4RZniXEEJcb0hNdPCbYDGwcx+1Rw+nIxu5bAGQVLMCGvXCaX5FbXZk+ZOOe\nd8vdJ88SIgBdEWBv88s4qMvMI5YHhotuu4QAA7hLn7nyGAHWIjM8gJcEWKMUbgDHBBgAIQkwAEIS\nYMDFnNTlGQIMgJAEGIGp9LvV7SefWBJgAIH1XMMJsE8o/KnFzIOtbruEAAO4kXr3PgIMgJAEWHOs\nEQGcIcAACEmAAbew8cPdBBhwPcvgPECAEZsjXtAtAQYE46ATEwEGEFXnyw8C7ENWrniemQdbPXcJ\nAQZwO/XuHQQYwL16niTdSoC1xRoRwEkCDICQBBhwCweduNtP7Qa8rfw+E+Nmre3gLRKbBkr3HHoT\nLMBKKXM4Lb+eTa/svgUkYFbHLNUSotCCHnjQJ7I82AzsjNX0q/x7k4UcQRmt2Lp2PCvROlmqAJuu\n/iqi7kssWy88TGcLrf3h4qD0b1OqJcTBBAtokpHpDsFmYOM4bo8aTmuG0+uhDyK2X6ABtCNYgA17\nyTS9EjGxAPhYtiVE+uQzs9AhAQbcRWHBrQQYEIY4ZEmAAZHY7J448zUIMAjBaAVbAuxb1jSAkwwX\n1xJgX1EUAycZLi4nwFphjQjgLQKMJJzYht4IMOBeCgtuIsCAG124MC4IWRFgQBj2iVkSYNA6B3xY\nMRmdCDCAeNQ0gwD7nsNvwEmGi2sJsCZYIwJ4lwAjFeUt9EOAkYdZLHRFgAH3umTjx9yaLQEGTTNw\nz8ywWRFg0DoDN0tqmpkAA3jU9wmkppkIsAv4bAdwkuy5kACrz4fArqWYgE4IMFJRCkA/BBhwO8vs\n3EGAAa2zzM6ueAFWfh38hSfbA/cxcLNieFv6qd2A95RSxt8Hevn1/EqNRgE8R00zizcDOzCO41jv\n3kpPgCcFm4F9YDUtuynhPt6jFnuXm+7F3ZXMdOPUwnzgmS76gXCLWPkDrOKc7KTmG8g/ls+4GCOT\n5WgZIsxSLSHC3ebEmv83mEaf4yQ9lwsWYOM4zqcQl6c56raKTuzOt2TY3dpcbaMF8ZYQt0uCq1fa\nXzMkooPVwmlucfk4a+BmRZ20EmwGBhUdxImk4Rl62pIAu4wl/nZcfi9O/tt0AM7TW74nwGrSg0M4\nedRQacx5esslBFhl+nEI52+TogQeI8DgMheWIymD8IOl3ZTXgasIMDjywVHAq8Zcs/OJ6zCR5VsC\njLSqPPBGW+6jd60IMHK65FH/+JNYimXO0E++JMCu9NYSv76blTKZM/ST7wmwmvTgln1ZYShQ4G4C\nDP70cYXxZWmS+JdIWaXgQgKMzD4eAQ2djcga5O9KXNN8Q4CR1pcP/Pf/dykItxJg1zNsAWeocr4k\nwC52smy3INCyC+/OB8OTEQ1OEmBwl1pnQBp3ctohyHlJgJFc9V++ZyD+WO4gP08X+osAI7PPRsAL\nx01DMC+dyScdaZcAu56NWT6m5/RGMn1DgFXgBEezbro1b8VSJ33j+JoIcs4QYCRXfULcSSC95cw1\ncd0msvyAAIP/3DpSGIY4cNw9ZPlfBNhd/uqRBrKW3TRS+HQgB9z0jwmwWxz3SP31YWdWER8oLNQu\nSwc3xYXiJAEG/7m1sFC1vMXlmsjyYwLsRjpfUw5ux2N36rgNRu2e/dU39IoDAuwuu92u7iBVGkjU\nWm04uOxTix64L22ORBV7xe4qYt1O2mD/5IAAoyMtFLm7bWigtGhIt6P5qhuYlL+UKsDKr9oN+c+q\nwGymXZ36a078fBt2v2m3Q5VnZNJtB/hGngArpYy/2smw4d9nUh+tbjtcPnxTthnWc6G9/cG7vRST\nuWP03CvOK2OWizQF2OrrppIMIJD20+GndgPu1f4NAOAzeZYQAeiKAAMgpDx7YMNixyvTDwXArlQB\nBkA/Mh/iqDghO/7WywOTVdrw5JX563s9fHei3JHdd59pRgt3ZHVsuNalWL7bzzOy9cxz8Y20AbZ7\nqr76t37sWP/LH3/+mMGtV+a4Gc+04WUznrkpJ+9IrWYsB+uKd+TJh+Vlr6j7jFQcwYY4H0ByiONR\n0+esa7eiiT3CFtowtFRjtvBLZKY2tHBBGmlGtxoZqV5KOwPjpepjRPXxuimPzX6OG1C3DU15bAnx\nTDPYJcB61MiTWf23pUzfev5nxQtS/V40pXqC1l2+m8y/Eq+1343XFEuInao+QFT87rP5l2cOVS9I\nI1eDpsy/3LV2Q5qWea2gnVOI26McVY42Td+37kGv7a+prHLmrak7svvWY81o7Y48OeNxR142rPGA\naL19ALDLEiIAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmA\nARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEIS\nYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEL6fy7BALdDr6L3AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "plot(z, d, 'b');\n", "stem(x0, 1+x0*0, 'r.');\n", "axis([0 1 0 1]); box on;"]}, {"source": ["Computing the Diracs' Locations by Root Finding\n", "-----------------------------------------------\n", "Instead of evaluating the function $d_{y_0}$ on a thin grid, it is\n", "possible to compute directly its zeros using root finding.\n", "\n", "\n", "Indeed, one has\n", "$$\n", "  d_y(s) =  \\norm{ U^{\\bot,*} \\phi^L(s) }^2\n", "       = \\sum_j \\sum_{\\ell,\\ell'} \\bar U^{\\bot}_{\\ell,j}$\n", "       U^{\\bot}_{\\ell',j} z^{\\ell-\\ell'}$\n", "$$\n", "where we denoted $z=e^{2\\imath\\pi s}$.\n", "\n", "\n", "So the zeros of $d_y$ are equal to the zeros of the polynomial $P_y$\n", "that are on the unit circle, where we introduced\n", "$$\n", "      P_y(z) = \\sum_\\ell C_\\ell z^\\ell\n", "      \\qwhereq \\choice{\n", "           C_\\ell = \\sum_{j} B_{j,\\ell}, \\\\\n", "           B_{j,\\cdot} = U_{\\cdot,j} \\star \\bar U_{-\\cdot,j},  \\\\\n", "          }$\n", "$$\n", "where $\\star$ is the product of convolution.\n", "\n", "\n", "This means that\n", "$$\n", "    \\enscond{s}{d_y(s)=0} = \\frac{1}{2\\pi} \\text{Angle}\\pa{\n", "          \\enscond{z}{P_y(z)=0 \\qandq \\abs{z}=1}$\n", "      }$\n", "$$\n", "where Angle returns the angle of a complex number.\n", "\n", "\n", "Compute the coefficients $C$ of the polynomial $P_{y_0}$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["B = [];\n", "for j=1:size(Ubot,2)\n", "    u = Ubot(:,j);\n", "    v = flipud(conj(u)); \n", "    B(:,j) = conv(u,v);\n", "end\n", "C = sum(B,2);"]}, {"source": ["Find its roots."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["R = roots(C(end:-1:1));"]}, {"source": ["Locate those that are along the circle."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["R1 = R(abs( abs(R)-1 )<1e-4);"]}, {"source": ["Display roots."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAOPElEQVR4nO3d0XLi\nuAJF0ehW/v+XdR/cQ7uBELAtS0deq+ahZzLVcWzQloQhpdb6BQBp/tf7AABgCwEDIJKAARBJwACI\nJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACI9N37AKZVSln+UGvteyQAU7ICa6jW\nWmu9lQyAAwlYKxZeAE3ZQmyrlLIumdUY8IKJ70cErJWlVY8Px+4P0LumOgbH4BjGOYa+BxDHFmJD\n3Z8PABOzAmtimUm5ERGgnf6r5ksZYZsCGJPx4VO2EAGIJGAARBIwACIJGACRBAyASAIGQCQBAyCS\ngAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZA\nJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQKTv3gcAV1fKP/9aa6fjgDRWYDCWu54BPxEw\n6EmuYDMBAyCSgAEQScCgJ7dswGbuQjxAKaU+jENl9eLG41fhpta/r4R5pMD7BGyX8vIleN3iTR4p\nsIEtxF1qrS8qVUp5XTgANrMCa2hp290G4y1p1meAOe4eAtbKT33SLeDm6eyWN9lCbMIDEaA1K7CD\nLRuGtVZbhQBNCdgB1om6/Vm3AJqyhQhAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgbZVSeh8CwJy+\nex/AtKQLoCkrsFZqrbXW3kcBMC0rsLPdVmbydhGb1+IeIVdgq2YPATubUWliTwejbVe8lHLg38aw\n1hdUzD4lYLDd44hzVGCe/j13VRMzLk7A4DN30TqzInffS8y4OAFry7AyjQFrcTuM9cpskGODEwgY\nvDJgtx4pGdckYPBcYgmUjEsRMPjHHOP+Y8mifxx4SsDgjznSdWf5WWSMKQkYzJmuNRljSgLGpU2f\nrrX1vuLXNX5k5iZgXNSl0nWn1mo1xgQEjCsydttUZAICxrUYr9dkjGh+nQpXYZj+ye1X//gwWbII\nGJdwS5d6/eRuNQbjEzAmZ+H1PksxsggYM7Pw2kDDSCFgTMvCazPbiUQQMOakXjvZTmR8AsZslnWD\nbcNDaBgjEzCmYuF1OA1jWALGPNSrEQ1jTD6Jg0moV1PrhjnJDMIKjBkYWM9hKcZQBIx46nUmDWMc\nAkY29TqfhjEIASOYevWiYYxAwEilXn1pGN25C5FIc9RrPfgn/ijLb3bufRRclxUYeear1+O/ptAw\nOhIwwsxRr8loGF0IGEmmqdfTAT+0Al4MoxcBY3il/Pnn6+trinp9/fCKV+5P9rdh/14saErAGNt6\nKCwldoSf3/2H/2sY7QkY9HG33spdfv3xWCwNozG30TOw2UfA+GhBV1ZgAEQSMAZmhQL8TMCIImnD\nerw0LhaNeQ2MoZVZ7pu/hNWVWm6pd+1oSsB2KT+/OWn9vk5P4228NzaXj5jiBAK2XVlNMMuzyaZu\n7THNh25c1tIwV5B2vAbWUCnFJHQPY98EPAVoxwqsodvn69S71wZWX+Upo94cbCT+yvnZQ8Ba+alP\nuvUmJ2oaNhJfeDq75U22EJvwQNzD2ZuJdNGOFdh26+2Ru7s5nn6J9zlpk7EIowUB2+XFnYeerttY\nfs3HK2E0YguR4Wj/lDSMwwkYA7HRNCuXlRYEDDiJRRjHEjBGYXSbm0UYhxMwBmKMm55pCgcSMIZg\nXLsCExSOJWCMwuh2ESYrHEXAgPOYpnAgAaM/d88DGwgYcDa7iBxCwOjMWHY1VtscRcDoz4gGbCBg\nQAdW3uwnYPRkFLsma24OIWB0ZiwDthEwoA/rb3YSMKADK2/2EzC6MQEH9hAwejINBzYTMKAbq3D2\nEDCgD+tvdhIw+jD1BnYSMLoxAQf2EDAAIgkY0JPNZDYTMKAb28jsIWAARBIwOiilmHoDOwkYAJEE\nDIBIAgZAJAEDIJKAAZ15KxjbCBjQk/tR2UzAAIgkYABE+u59ANO6bevbIflHKV9fX3X5gzMD7CBg\nTaw/acKnTvx191q9hgE72EIEIJIV2Nmuu7X49FZpi7DjrE+wk5rCWwj2ELCzXa5bN7U+bxhHsDsb\naj0giNmnbCHSlVG2GYMh07MCa6LWet2twhfuFmHOzEG0imsSsFZ067m63ELvzswjPd2ddYKZni1E\nACIJGMzgbr1l+cUV2EKESYgWV2MFBvTkBVE2EzAAIgkYAJEEDIBIAgZAJAGjDx/7BuwkYHTgrjNg\nPwEDurEQZw8BA3qyHGczAQMgkoDRje0jYA8Bow8bR8BOAgb0YQnOTgIGdGMhzh4CBkAkAaObWqtN\nJGAzAQM6MHdhP7+RGdq6G6i96HPjBTB2sgKjs7ln4o8/3NQ/LpxKwOjJHPya5p61cBoBAzowd2E/\nAaM/83FgAwGjs7ln4o8/3NQ/7lvMVziKgEFbt2LVql5/zD1r4TQCxhDmnpVLF7QgYPRnPn4dc89U\nOJmAMQpD20WYr3AUAWMIBrUrMEfhWAIGnMdMhQMJGAMxQwfeJ2CMwtx8bqUUl5hjCRhjsQgD3iRg\nDMQMfVaWX7Tg94HtclsuPD451ysJT92PGOyAdwjYdutx9umYaxTeoNZqF3EyZiQ0YguxoVKKsXgD\nDZuJS0k7VmANLbPOu+nni11H1kzbp+E6viDwewjYu+4eZ78+J3/6HzyZ32ERNgcX8VdPZ7e8ScDe\n9VF4rB4O4TROwBWkHQHbbr1KuLub4+mX+IhFWDqXj9YEbJcXdx7q1iEswkIt9XLtaErAGNeyCNOw\nGOv3Pn59+SWetOY2eoYmXTEeNwxtIdKYgBHAqynAIwFjdLe30/U+EGAsAkYADQMeCRgZvBg2uscL\n5JLRmICRxCJsWPef+6letCdgxLCROKy/7/q6/QPtCRhJNGxY9ng5n4ARRsNG41rQi4CRR8PG4SOj\n6EjAiKRhI1Av+hIwUhk3+1IvuhMwgvmVK72oFyMQMLJp2PnUi0EIGPFuv3Wl94FcgnoxDgFjBu7p\nOId6MRQBYxIa1pp6MRoBYx4a1o56MSABYyoadrjl9cVaq3oxGgFjNhp2IAsvRvbd+wDgeOuGGXw3\ncwIZnBUY07IU2+z2tgT1YmQCxsw0bINbutSLwdlCZHK2Ez/iRBHECoxLsBT7lbsNiWMFxlVYiv3k\n1nWnhSxWYFzLLWNWYwuveJFLwLic22B98Ya51ZB0thC5qCvvKNozZA4CxqWtfxXLFUZz6WImAsbV\n3YbyuQf3uX86rknA4I+7+zv+GejXr5alBUC6mJWAwT+eZOzuXo9SIhqmW0xPwOCJv4N+2p2K61sr\npYu5CRjMQLe4IAFra/lsnt5HwQ6PW4gj7c7pFlcmYK1c/E2yc3v6cR7n9OPucSVaXJmAteKzHuZx\ntwj7rxnreDz9bKr9dXn6+BEtWAjY2cbZfeIDv12sp1fzkOmLx8nczHH3ELADfLSrYzy6DteaX92t\n4zseSSIBO4BxCuB8Po0egEgC1pbFGUAjAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiPTd+wBmUEqp\ntT7+x9ufH78KwE4Ctsu6Uo90C6AdW4i71FpfVKqU8rpwAGxmBdbQ0ra7DcZb0qzPAHPcPQTsXXeP\ns1/z89P/oFvAzdPZLW8SsHd9FJ6nt3UAcCABO9iSrlqrrUKApgTsAOtE3f6sWwBNuQsRgEgCBkAk\nAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwC6nlNL7EByDY3AMHOC79wFM6/Z8qLX2PRKAKVmBNVRrrbWa2QG0\nUKwPWivl70kWM+AFA/JHBKytdb0AOJDXwA5wt65airX8R/UCaETADvBTpdQLoB0bXE08XZMBcCAB\nAyCSLcST9H1b2OO9JOs14vmHdP69LS/Of5dT0evxMNp5WH/3Qb5pr/PgnaMbCNh5bjd3dHlW/HQ8\nJ+vyRoK7dzI8/uDnV+T18fT6vpd6SLz4aq+EdBkionkj80k6PiV++tallPPHjhfH01GXUzGg6zwk\nXn/TXufh5O84AQE71VBzKx8UcuNULJyHRcfzMNQQMT5biMfr9bawj259POdJ0uVuzE+/qfFi4Tws\nep0H7xzdQMCO1+ttYe///afN8ro8Gz/6pia8C+dh0fc8uASfErAzLHOrQe4yWp6igxzMOZ7+vMt5\n6HIqep3/0c7DUPqeh6GGiCCmXQBEchMHAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACR\nBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASP8Hm7jZeE/81LQAAAAASUVO\nRK5CYII=\n", "output_type": "display_data"}], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "plot(exp(2i*pi*z), 'k');\n", "plot(R, 'b.', 'MarkerSize', ms);\n", "plot(R1, 'r.', 'MarkerSize', ms);\n", "axis equal; box on;\n", "axis([-1 1 -1 1]*2);"]}, {"source": ["Sort according to angle, and discard one out of two (because of double roots)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["[x1,I] = sort(mod(angle(R1),2*pi)/(2*pi));\n", "x1 = x1(1:2:end);"]}, {"source": ["Compare the found locations with the ground trust."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Error (should be 0): 3.1113e-09\n"], "output_type": "display_data"}], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["disp( ['Error (should be 0): ' num2str(norm(x1-x0))] );"]}, {"source": ["Noisy Recovery \n", "---------------\n", "In practice, one has access to noisy observation $y=y_0+w$, and one\n", "approximates the support $x$ using the best $N$ local minima\n", "$ \\tilde x =  \\{\\tilde x_1,\\ldots,\\tilde x_N\\} \\in \\mathbb{T}^N $ of\n", "$d_{y}$."], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 1__\n", "\n", "Compute and display $d_y$\n", "VD\n", "_y\n", "isp"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAASuklEQVR4nO3d25Kj\nOhIFUDRR///LzANdHAqwjW1AytRaMXHCbfd0yVy0pZTsKuM4DgAQzf9qNwAAPiHAAAhJgAEQkgAD\nICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTA\nAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADIKSEAVZKqd0EAC73U7sB\nZxJdAP1INQMbx3Ecx9qt2CNZW9PCGSmliWY0wqGYuTAOK432+F8o5b831cKcbD6+9ZvCMAzNnJGx\ngTY0opEz0oh2Loz20yFVCXFXE+eglBbasYz2ntswtaOFMzI0cDSckb+taOBoNHMoajfhtVQlRAD6\nIcAACClhgDUw+QbgcgkDjEdaiPYW2sCSM7LkaMQiwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAA\nQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABENJP7QaQTSn/Hvjl\n7MClBBhnmtJriq7lY4DTKSFyjlKGUoZx/C+xljEGcDoBxgkeTbZkGHAdAca3npcKlRCBiwgwvnJw\nocskDDidAONbL9PLJAy4ggDjc9OujeN/GeBE8bbRl9+OcNz0nU9e4nRvBdI4CjDgZMECrJQyh9Py\n8fOXOJ3PeJHPcozl2g5BCbE70we2pv99/C8M79/hJmG0b/4go2s1hGAzsCPK30tv9cfOZ2arb8r4\n4GufzL1IaTX96jPASrS3nSfAxnGcjv78YH6+XqPastpzMT+ek+zgbnhHlJRWF/Zbe5RyWC3KVGzJ\nQXkCzLrXc0/uxtWEbPevnfIVvdPA1llqii9fHvYKht1OwmIJFmDL2dVqy8buS0yO3IqP6oo6uMTm\nIcv8VZbd6vm9xxUswIa9cJqfkVtPHDw227ri8f8vsSwn3DJsS8GgffECjHd9Vgm57r7VKbRjtSba\nZ9Gsz3edg230XWgnMNppCcxclkEJsOSMLtn16MJwwaw4IC0TYPkZXbLr0e9v68rzfOrwgMQiwDJr\ndvDYbMPokJSKS4Al1+DN2WCTevN8AGF4QRQCDHrkN2iTgABLyzgavtftpwtCEGCZGU3DEz6SGJ0A\now6j2pY5O4QgwKjAsLeil+Hk7BCFAMtJbYQnXBvvMiVtkwADeEbeNyt9gI2GTsCWKsUT0+8Grt2K\n19IHWBlM/wEO++0wA/Sb6QNs/VsZexDizfp4TRUHpx3OTree/Fr2BuUPsCHOyThRh28ZOEWg3qOL\nAJsYUQKf6WRKGu499hJggcYUwNXC9dS3idVV9hJg/XBnwhGxemp29RVgnXTuge7MTs5II9492s5O\nVyKe7o4CLFC33gln5H7Hj7mz06FwJ72jAAMgEwEG8FonGxFj6S7Acl+Cud8dcJGgX6zVV4BFPEPv\n6uE9wjeCdtZs9RVgNMisEfiMAKMmA+HbfDDtsOpD435qN+Bt5feWGvdux+nV3ZcWf0e/CfBP3GFK\nsAArpczhtHy8emb70sygEmAl6Jg+WIB9oPzNq+eTs9DMLIFvlGij+2wBti0hJk4s4F3fdNFT/SZx\nj7LsLUOEWZ4Ae15dBJjoG9LodBdihLFFR5wOqCX03RdsBjaO43YX4jTf2n3pwT8S+5wl43Tc4OPC\nV/qiGUPkKWmwABv2wml+RtkQoB+dlhDzMUwGeiPAAAhJgAEcZcm2KQIM6IVK+0r0A9JpgBlGtcbp\nAN7VaYDRlNBjwPSM9miWAMtA/8IT0ctE8IgAS0IPBfRGgAEQkgAD6FGCpYd+A8zSNPCBTF1H9KWH\nfgMM6IrNLPkIMABCEmDh5ahmZCrLpOTs0CABloHCCJdygdEmAQaZWfghMQEG0J0cNeGuA8y6C9Ct\nBFPzrgMMgLgEGJDfuWuBijeNEGAAhCTAYsu0x8yoFniLAAMgJAEGaZ04ozU/ziTNqRRgec4lbKWp\nMHOuHBdG7wGW4ywCdKj3AAMgqJ/aDXhb+S35jX9nT+VvKXA0twJILViAlVLmZFo+Hv4mVuljXauP\ndwktmna1GCfXlbCEuAq23PK9UanM6SRNVsFmYB9QWgzEXm2oKFztKluAbadfEgtg9iSkwi3EJCwh\nfiDCmYLKzI/TSDOqDzYDG8dxuwtxnnV9tvrltiQlCz+kFyzAhr2S4PyMaiFAP5QQAQhJgAEQkgCL\nygoH1GX5vDoBRnN0CpzIUC8xAUZb9DXAQQIMoBfJ5qMCbBjUsuEwdwrtEGCQ0EUxk2nwTgICDHIS\nNqQnwAAISYABEJIACynZViKADwgwIC17JnMTYDTHpxo40aW1CtdqXQIMoAv5lh4EGAAhCbB/lALg\nIHcKjRBgkM2lAZOsBkVoAgwSEjP0QIABEJIAi8cKBMAgwIJSIAIQYEBOahXpCTBa5FMNnOKGWkWU\nazVEI98lwAC6kG/pQYD9kXKQApCSAPtPvuEJQGI/tRvwtvI7Sxo3gfPkJeAs06qPm4zqggVYKWUO\np+Xj4Te9pmdWL0E/RAv9CBZgLy1jDIDEUgXY7uSs/N2YET3bjK+Bi5Ro29hSBdiu6IkFcI/tokzj\n7EIEIKRgM7BxHLdbDaeC4e5LxGWrG9+4+eJp/FptvHkfCxZgw144zc/ILeB+Ub5NKh8lxD9ciABR\nCDAAQhJgwCfUKqhOgAFvs9xMCwRYJMa8ADMBFoyRL09k3S0NuwQYACEJMIDMEi89CDCalvje4zpK\nqStZj4YAo11Z7zpSMti6nwBb82UcwLsMtqoQYACEJMAACEmAARCSAAvDyhxNsVpMdQIsEgvFADMB\nBkn48BO9EWAAaeUu8wowgMwSz8sFGE2zU4B3uWD6IcD2uQcgrlpzDv3GzQTYjsQzbuAi+o37CTAA\nQhJgAIQkwGJQWwdYEWBhqLDTGntEqUuAARDST+0GvK38DvnGzZSkLEaD21chMd8jRYeCBVgpZU6m\n5eOZ3AKYpB/WZCshllKKqnw6TimwFWwG9tI0A1tN1LZ/gUDsFOA4l8o3wo3+UwXYbjhJLOhK3Ts+\ndNVu2VuGCLM8JcRzD7dRP/CuuNEVVLAZ2DiO212IU8Fw9yUAsgoWYMNeOM3PyC2AfuQpISYWuqoO\ncBEBBnzOajEVCTAAQhJgEJ4iM30SYAAJ9VDaFWAEYKEFPpB+Xi7AgCSUUnsjwB4y6gdomQADOI2B\n750EGAAhCTAAQhJgrVOOANglwAKwsYqWWfWhFgEGQEgCDCCbTubEAgxi8+lddvVwVQgwYrDQAqwI\nMCADM9EOCbAXjPoB2iTAnjGgA2iWAAM4k/Xa2wgwAEISYE0zjgN4RIC1zjoc7VM0owoBBkBIAgwg\nlX5mwwIMAvPpXXZ1clUIMMKw0AIsxQuw8uvJX7izPUB1ZqJ9+qndgPeUUsbf63T5ePkXzv2J06jf\nvQHQmmAB9twUaasMW/1xm3kADAHLV6kCbJfEAm4WtHKz7C1DhFm8NbBHpsO9/C8AieWZgT1fG4so\n4ggO4DbBAmy5xJUvsSCuoEUzQgsWYMPemtbqGWEGdKur9ZM8a2AADN18DccgwCCursbasCXAiMS3\nSa30M9aGLQEGxGZM0y0B9ppRPzTOTLRPAgyAkAQYwPlUbm4gwBrl0gd4ToC1S1kf4AkBBkBIAgw4\nh1Wf6nr7OkoBBkBIAoxgDPMnDgIIMIiqq2LRI4K8ZwLsKPcJtEmQd0uAHeIOAWiNAAO4isrNpQQY\nwCVUbq4mwFrU24c5AD4gwAAISYABEJIAA07jY+YVdbj0IMAISS8JCDDi6W2YudXhWBu2BBgQlSDv\nnAADICQBdpTVaYCm/NRuwNvKb4yMm9rBk5cA7jcNfHVIFwkWYKWUOZyWj2fTM7svReFyBzgiVQkx\nbmgB8K5gM7CXymadavWMkAPYte0/G5ctwOYS4uoZ4B5Wfao4JXqWvWWIMMtTQgxxuDmREw5LHQ4a\ngs3AxnHcbjWctmzsvkRWPX+qods3DivBAmzYC6f5GblFJ1zpg/26ZCohAtAVAfaGnstWAK0RYACE\nJMAALqRycx0B1hYXOsBBAqw5NlYBHCHAgPOpJdyp26MtwICTqSLcr89jLsCIyto4dE6AQSS+fmLi\nODAIMACCEmBvU7YCaIEAe4+qBUAjBBgAIQmwhihOQkp2zF5EgLVFiRLgIAEGnM+c4zY9f6JAgBGb\nXhK6JcAIrLeBZ89jbdgSYEAwgpyJAAMgJAH2NqvTAC0QYACEJMBaYVYHubnHTyfAGmJdGrJyd19B\ngAFXMee4WudHWIARXuf3cLPMOe7R83H+qd2At5Xf7mrcnLcnL5GVTaHQrWABVkqZw2n5eDY9s/sS\nhCanJ44Ds1QlRKFFbi7wiePAJNgM7IjV9Kv8HbCdFXK+zAZIpkSb3qYKsOnoryLqimmZdRcgnydD\n/zalKiEOqogA3Qg2AxvHcbvVcKoZTs8H3YioIAnpTZUbd/qJggXYsJdM0zOxEgvgSxGKfNfKVkKk\nQ5Ykm+XUXK3zcbsAAyAkAQYBWDuZmM+xJMCASAQ5MwH2OYNBgIoE2IcMAwHqEmAA91G5OZEAq8/6\nPHTCnX4uAUYGPm9Ebwx8BwEGXMrYgusIMGidAJg4DqwIMAhAsWjiOLAkwAAISYABEJIA+5zVaYCK\nBFhl9sJCVwx8TyTASEK/AL0RYMC1jC1Op3IzEWDQNF3VRASyJcCAGAQ5KwIMgJAEGAAhCbCvWJ0G\nqEWA1ST8oE/u/VMIsMqsS5/IhLhlTs3sy7vekZwJMGhXmq5Kl30uA9/JT+0GQDDLzvSGfkRXNXEc\n2IoXYOW3/xgfXNGllEcvwZemq2++vlZ/BO4ULMCW4bQNqqLQwJW2cTWtuvmyDKgi1RrYOI73z71s\nHOjEo8mW6IJags3APrCalqkupnfdfOjRPzsNYlxZRBeuiJU/wJpNLF3eFS6aEB85WU7oSw7R7ONB\nz6URs+wtQ4RZqhIiXOHIjXxFvxyhA3mD6DqLIzkLNgMbx3G7C9G2Q6528Po6fYbhuh7M23gsWIAN\neyXB1TPCjBMdnwbZzgM3U0KEF94aEckwuI0AO4Ghd2tqnQ6Tf7iTAKtD4F3nxBRpcJ8YaRj4fk+A\nVWO0ntJZpzVr15b1fd3DfpYVAQb7qne1+boq01nOJcDgoY8jRLd7onxBzlkEGOz4JoF0uHAPAXYO\n67FNOeV0yCFonACrQNSlZ0DDEa6TLwmwOozuW1a9T0m82UyX/bHEV8XHBBjsOKWn0FN/SZfNcwIM\nLqHnhasJsDMZcbfj41KVUT9EIcBOo9djxXoPR7hIPibA7maA37IWupL0V4hcX6r1q+ZyEGDwx+nd\nhM76M7psXhJgpNXCSF8XDNcRYPBP9bQD3iLAzvRyyK+LbNxFE6bj572fupl7YdZCqSAoAXa3Trqn\ncK7rQZzxrZfHpJ8gP8LReESAkdzxZLq0jzDEhtMJMDJrZNx6fKs0cJwAO5lydkTtlGgaacY9Ht0p\nfd5BjsYHBNh92ukle9NIF9BIMxrx/F7o7U5xND4jwEjuyH6BfpoBmQiw86kihnPbCPf5hdHhQHt7\nQHq+dxyNdwmwm7RQPywN3A1V2vBkSHFnc6pfALsqXhWPDkjFA+VoxCLArtJAWPDCdI7u7CAeRWkL\n45ta3ClLy6PhyLyUKsDKr9oN6bczataj5KhyplYtaeCCrWZ1/HsO8mHvauz5aByRJ8BKKeOvFjJs\n+O2Y7h/m80gLw9vpSlj99M4vj+WdwsTROKKMWW6dKcBWjxtJMoBw2k+Hn9oNuFb7JwCAz+QpIQLQ\nFQEGQEh51sCGxYpXpjcFwK5UAQZAPzJv4qg4IXv+o5cbJqu04c4j8+hn3Xx2opyR3VfvaUYLZ2S1\nbbjWoVi+2s89snXPffGNtAG2u6u++o++bVv/y7c/f8zg0iPzvBn3tOFlM+45KQfPSK1mLDvrimfk\nzpvl5VVR9x6p2IMNcT6AZBPHrabPWdduRRNrhC20YWhpjNnCl8hMbWjhgDTSjG410lO9lHYGxkvV\n+4jq/XVTbpv9PG9A3TY05bYS4pFmsEuA9aiRO7P6t6VMP3r+b8UDUv1cNKV6gtYt303mr8Rr57vx\nGqSE2KnqHUTFnz6bvzxzqHpAGjkaNGX+ctfaDWla5lpBO7sQt1s5qmxtmn5u3Y1e26+prLLnrakz\nsvvSbc1o7YzcOeNxRl42rPGAaL19ALBLCRGAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBI\nAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQ\nkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACH9H6LxHL1b\n8kh7AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["__Exercise 2__\n", "\n", "Display the roots of $P_y$ that are inside the unit disk\n", "oefficients\n", "oots\n", "eep those inside\n", "isplay"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAANLElEQVR4nO3d25Ia\nNxhG0W6X3/+VyQXOGAPDcFBL+n6tVblIpSoJhpa2pG5m9tPptAFAml+jXwAAvEPAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRg\nAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJ\nwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkCk36Nf\nAFSz73vz/+bpdGr+34R0u4EBb3hcqbbDquf/C4IIGPzsbkJmGDvTvjDoQMDgjtswpIyU3FcOrxIw\n+ONq6q8xNEr+oeBMwFjaavP75Z+3/B+W8gSM5awWre+IGekEjFWYr7/jnSGUgFGc2fklX2+X94r5\nCRg16daHlIz5CRjVmHkbsg5gZgJGEbp1KG8vExIw4plbe/JuMw8BI5WZdCBvPjMQMPKYPSfhg2As\nASOJGXNO58/Fh0JnAkYG6Zqfz4jOBIzZmRaz+LzoRsCYl6kwl0NFOhAwZiRdNcgYhxIw5iJd9cgY\nBxEwZiFdhflwOYKAMQWL9BXIGG0JGIOZ1FZjsUIrAsZI5rJl+ej5nIAxhvkL1wAfEjAGMHPxxcXA\n2wSMrsxW3HIflPcIGJ2YpHjM4oZXCRg9mJt4huuElwgYhzMr8RIXDE8SMA5kJuI9rhyeIWAcxRzE\nJ9w05UcCxiHUiyZcSDwgYDRm4UxbGsZ3BIyWzDUcwXXFXQJGM2YZjmNnzy0Bow31ogOXGZd+jX4B\nVGBaoY/zNfa1G2NxdmB8RLroz1XHmR0Y7zOPMIR9GGcCxpvUi4E0jE3AeI96MZyG4R4YL1MvprLv\n5rFF+eB5jXoxIZflmhwh8gLTBHNynLgmAeNZ6sXMXJkLEjCeol7M73Q62YQtxT0wfqZeBHG5rsMO\njB+YDsjiftg6BIxH1ItErthFCBjfUi9yuR+2AgHjPvUinYaVJ2DcoV6UoWGFCRjX1GsUM21zHuio\nzWP0/EO9hriaYL39bbmqq7ID45px3tnt9sCGoS37sKoEjL/8VG+q0rCSBIw/jG1qszirR8DYNjcJ\nWIMH64sRMNRrsNs33kdxKA0rwz2Ppf0/kF0GUzh/HD6Ko1mxlWHmWpdHt1mWhtXgCHFRHt1mZdJV\ng4ABi3IzLJ2AASvyzbACBAxYlIalE7AV7fu+bdeD1k0BFuRmWDRPIS7H81dw5kHcdHZgK1Iv8CBu\nAQK2Fsf9QBkCthCHh0AlArYW9YIHHFFkEbBVGJlwyc9QLkDAluDwEG5dDojz3/t9K1l+j34BdKJe\ncOvusPCryVPYgdVnRQnPk64gArYEYxJeYtkXQcCKcxgCrzJkUghYZVaR8B5Pc0QQsOKsJeFtGjY5\nASvL2INPWPzNT8AqMwLhQxaCMxOwmjy7AZ8ziCYnYACP2IRNS8AKsv2CVgylmQkYwCMeqZ+WgFVj\n+wUsQsAAfmATNicBK8X2C1iHgAH8zCZsQgJWh+0XsBQBA3iKTdhsBKwI2y9gNQIG8AKbsHkIWAVG\nFPThnGMqAlaEcQXdWDJOQsDiGUvQk8XiPASsAiMKOrNwnIGAAbzGknESApbN0/PAsgQM4B1OEYcT\nsGDGD4zi5GMGApbNKIKBLCLHEjCAd1g+DidgqTy+ASxOwADe5xRxIAGLZMzADJyCjCVgqYwcYHEC\nBvARJyKjCFgeowXm4SxkIAGLZMwACBjAp5yLDCFgYYwTmI0TkVEELI/RArAJGAChBAzgU6fTyfF+\nfwKWxAgB+CJgYdwAAzgTMIA2nJF0JmAADTgd6U/AYvgFYACXBAyASAIG0IzbYD0JGEAbDvk7E7AM\nlnUAVwQshsUdwCUBA2jJeUk3AgbQjJOSngQMgEgCFsBXmAFuCRgAkQQMgEgCBtCSX27ZjYABEEnA\nZmcpB3CXgAXwCCLALQEDIJKAARBJwADac/e6AwEDaMx96z4EDIBIAjY1pxAA3xGw2TmLALhLwACI\nJGAARBIwACIJGACRBAyASAIGcAhfgzmagAG05wswHQgYAJEEDIBIAgZAJAEDIJKAARBJwACIJGDz\n2vdt2077vvkyCcAtAZvUVbQ0DOCKgMXQMIBLAgZAJAEDINLv0S+AZ/nJahzh6mjaZUYQO7BJXc0j\nphWOcHtj1a1WggjYvE6nbdv200m9AO4QMAAiCRgAkQQMgEgCBuu6vb3qhitBPEYPS1MsctmBARBJ\nwADa232l7ngCBnCIk/PZgwkYAJEEDIBIAgZAJAEDIJKAARBJwGbnYVyAuwRsah7DBfiOgAEQScAA\nGnPy34eAAbTn/L8DAQMgkoABEEnAAjhPB7glYLNzkg5wl4ABEEnAAFpy5t+NgAE05uS/DwEDIJKA\nZXAoAXBFwAI4jgC4JWAAzTgs6UnAAFpyZNKNgAEQScAynE4nRxMAlwQMgEgCBtDGvu9ugPUkYABE\nErAkboMBfBGwGI4mYGbWl/0JGEAbVpmdCRgAkQQsiW+DAXwRMIBPWVkOIWAADbgB1p+A5bHWA9i2\n7ffoF8Br3Aab3NWHY1G+AkNyFDswaOZ2HjOzLcL54RACFsmKD0DA8ljrAWwCBvAJP4F+IAFL5RRx\nQrfzmJkNjiNgkaz4pnX5yfiUyrOOHMtj9NCYbi3FanIgO7BgVn/AygQslXUfjGUFOZyAZTOEYCDr\nyLEELJjBA6xMwABe5utfMxCweE4RgTUJWDZrQOjPqnESAlaB4QSdWTvOQMDiGUjQk/XiPAQM4DVW\njZMQsAr8mmZgQQIG8CxPz09FwOqwCQOWImBFWBXC0awRZyNgpRhgcCgrxakIWB2GFhzH6nBCAlaN\nYQYHsUacjYCVYoDBEawL5yRgBRls0JzV4YQErBrDDNqyIpyWgNVkyEFD1oVzErCCDDZoxVpwZgJW\nloEHHzoPIivCaQlYTYYcNGEozUzAyvIj6uEThs/8BKw4gxDeZvs1OQGrzPCD91j5RRCw+gxFeIP1\n3/wErDiDEF7lt1amELD6PM0BzzNYggjYKgxLeJLtVwoBW4IBCc+wzssiYKtwkAiX9v3PXxf/xM/d\nCCNga9Ew2Larbv39e/XKImALMThh+7dYX//E2i6RgC3HQIW7rPDiCNhazkNUw+CKeiUSsOUYqEAN\nArYiTySystslnEVdKAFbl4axrNPpz1/btqtXLgFblJth4GcephOwdRm6rMzqrQABW5qbYazJD92o\nQcCwFGUt6lWGgK3OzTAWpF41CBgaxkI8uFGJgLFtFqSswSqtGAHjDw90UJtbX/UIGP/QMEpSr5IE\njL/cDKMk9apKwPiHhlGSepUkYFzTMCrx2GFhAsYdGkYNruHaBIz7NIx0bn2VJ2B8S8PIpV4rEDAe\n0TASqdciBIwfmAXIol7rEDB+5od0kEK9liJgPEXDmJ96rcY3JHiBCYJpuTgXZAfGCzzTwZzUa00C\nxms0jNmo17IEjJdpGPNQr5W5B8abTBwM5+ccLs7Hz/s0jFFce2yOEPmEs0SGUC/OBIyPaBidqRdf\nHCHSwFfAXE4cSr24JGA0Y3LhUC4wrggYLZliOIItPncJGI1pGG25oviOgNGe9TKtqBcPCBhHMfXw\nIZcQjwkYBzIB8R5XDs/wPTAO5FtivEG9eJIdGD2YkniSS4XnCRidmJh4zBXCqwSMrkxS3OXC4A0C\nRm+mKi65HnibgDGGaQvfF+RDAsYwGrYynz6fEzAGM5GtxsaLVgSM8cxo67BeoSEBYxYyVpt00ZyA\nMRfTXD2WJhxEwJiRjNUgXRxKwJiUuS+aj48OBIypmQfj+MjoRsAIYE6M4GOiMwEjhvlxWj4ahhAw\n8rw3XZ7/Jdd7W9LFQAJGquenzqtfqOmSb0K6GE7AyPbjNHr310G76t+mW8zj9+gXAB/5mkZNrEfz\nDjMbAaOI86y677svQbelW0xLwCjldkO2babdd1y8gdLFpNwDo7L/Z+F/LnKX/AO6RRA7MCq7tyHb\n7Mlu6RaJ7MBYjsn6i7eCaALGuv7dlq0yg4sWZQgYbFv1mNX+07EsAYNr+82Xn+OGiWKxAgGDH8zf\ns9tXuM33IqE5AYOX3Q3G2aEDatT/F+YkYNDMg8C0YsDCFwEDINKv0S8AAN4hYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAi\nCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnA\nAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQKT/ANar\nlWZchgT8AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["__Exercise 3__\n", "\n", "Keep only the best $N$ ones, i.e. those that are the closest from the\n", "unit circle.\n", "We denote those as $\\tilde x \\in \\mathbb{T}^N$.\n", "Recover an approximation $\\tilde a$ of the amplitudes $a_0$\n", "By solving in least squares (using backslash operator) the system\n", "$ \\Phi_{\\tilde x} \\tilde a = y. $\n", "Display the recovered measure $ m_{\\tilde a, \\tilde x} $.\n", "osition\n", "eep only the best N ones.\n", "ompute amplitude by solving a least square.\n", "isplay the recovered measure."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAS7ElEQVR4nO3d4ZKj\nKgIGUN2a939l94fTTm6SNmpU+OCc2tqa29PTTRD4BAkZp2kaACDN/0oXAACOEGAARBJgAEQSYABE\nEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQS\nYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABE+lO6AAeN4zhN0+sXlz+//i0A\nLckLsMeUeiW3ADqRt4Q4TdNKSo3juJ5wALQhbwa2bs62ZYHxNcxM0QDa0FSAvQ0niQXQpLwlxN9Y\nOQToSgszsHnBcJqmJcPMugCa92Yzekve7rYHoAHtLCEC0BUBBkAkAQZAJAEGQCQBBkAkAQZAJAEG\nQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZA\nJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAk\nAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQB\nBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZApAYD\nbBzH0kUA4HJ/ShfgTKILQo3jOAzTNJUuB1HG9prMOP57Ua+R1t7rhWhzJ52GYe6reijbNTUDe0t/\ngGr9vd38l2HT4w0orGvwGRiQZRymYRgGucVOAgwoz+NrDhBgQDHTNK8ZLl+wfsgODQaYDgBZ5n0c\ndhGzV4MBBkSZ5pvOaZqGwd0nOwgwoBY/GxJhEwEGQCQBBkAkAQaU59kXBwgwoJhx9PZljhNgAEQS\nYABEEmBQyDjaMw7fEGBAeZKcAwQYFDCO4zgMo/OTXngvM9sJMChj+jm1U4bBMQIM7vb4mY3OnobD\nBBgUMG/gmGde80eKFC4QBBJgUMYy9ZJecMyf0gVo3zI6WSviwfh4fFKfC4mO4eBLZmB30Ev5r2l4\n+BTHPtPriSrgADMwKOBvaI3jaOiGo8zA7uNJB8CJBNi1llV+q0QA5xJgcCsTcTiLAIO7/dtAX7QY\nVXmsCqdJsZEAAyCSAAMgkgC7j4URXmkVcJgAAyCSALuQk3LgN+adfE+AAWW4veNLjpK6ys9Jd3//\n03l3DC/TDm1ioSo4QIBd4ucDMqbHj8yQYQymHXAeS4hXEVcAlxJgAEQSYFfxMbtsp7HAAQLsEj/r\nh+Ni/op3rfLo50lp4WLUQLfgAAF2lWma5q1V0zR5Hga7uNVjCwF2Fd0P4FIC7ELLOuLf/wFwHgEG\nQCQBdr15IuYxGMCpBBjcxOHOCwvqnEKAXcv+Q1b03DieekbPVcFhAgyASAIMqmBVDfYSYFCelWY4\nQIDdzREDAKcQYNdypC8rNI6FquAAAXYJsQVwNQF2FU81AC4lwACIJMCAGtnuxEcCDIBIAgyASALs\nWs5CZPb2JF+NY6EqOECAARBJgEEV+tmz4GNlOIsAAyCSAAMgkgC71tuzEPtZLGKdVrBQFRwgwACI\nJMAAiCTAzmeTFcANBBgAkf6ULsBuy7aI10MuHndMOAIDoG1hATaO45JMj39eyC2ATrS2hDiO49ud\n66UIVFZoHIu3VeENJ6wLm4F9NAfG00Tt9RsASNdUgL0NJ4kF0KR2lhCrWjmERxvbphUz2CVsBjZN\n0+suxHnB8O1fQSU0SThdWIAN78Jp+UqFufV2qyTMRvs4fqgKDmhnCRGon3NqOJEAK8PTDoAvCTAA\nIgkwACIJsJNZGAS4hwA7n2fUADcQYNeyh54VGsdCVXCAAAMgkgADIJIAAyCSAIO62Mj6yFv+WSHA\nruWMfFa8No5uN/3oJxwgwACIJMDgWibhcBEBBpfrdmEQLiXAivF0mt74LBXOJcAAiCTAAIgkwK7l\nLERWaBwLVcEBAgyASALsTJ5RA9xGgAEQSYABEEmAXctZiKx42zj6fINgf6+YEwgwACIJMAAiCTD4\nrMM1PaifAIMP5vSSYaX0+VBwu3Hst34EGHzm7X1Urs8ME2AlubVsnuvLpeYG1u0NlgC7lrMQ031/\nusrKP++qcaxneVdVca6exxgBBluZTn2p56H2Bh22TwEGmxh8qVmf7VOAAUTqcMr1RIABpOpz4rUQ\nYNdyFmK0q6/ebz++w+2pDb/c+X1avV3Qewiw02igpVw6OnR+h8sptKKLCLAzaab3W6LLDQSV00RP\nJ8CIN0033Tp0uLLHl5YGc0MT7bB9CjCgQfU8eQpamKmhunYRYAR76m9x3Y9LzcmhVWyUeGi1ACPb\ncnsbdJ+bK2V0u3PhrpTvDzl7+mlDYOQLsGs5C5EVWY3j0jv09ao48HSnwp5X+TOqucYqrLcVAixS\nJYv79GO5Q698FH5VsMBZFZVIgBV2oHfZOD47dwnl/p8f55va6LYyI17169VJGVsEWKTbNo7D8DKc\nVT4Jq7lsVzvl0gSNLQIsTM+dsxLbL8E3047aQiJoUBvSSsthAuxaV5yF+LjvrqoxrrgbKuTckdHV\nW1xaFXpKqwQYsObtPFIkNOO3dYKI6yvAkkQ0qXuoCir3WxONaLopa7DNB9gU0Vy2e21YBV/gclpP\nkTKk9DHuVFV/fztz5UTNB9hNzbnULuGC/eHxjUFDZQMHZ0m8rELiFBHV2HyAcYnHg2eGkLbOMSsX\nNyvejpW2zkMDKixSEQKMg6pazOR+WXctx0pb86EBv72iXftrKnxduwiwa514FuLKKuXhVnjsIdbb\nb75zOEvvdYuoCLhWtVXR/KEB0a+uiwBrZrx763D7++Yh1ulRutcpt58r2m4zG1V1/tP9G/erevm8\n1X6A1d8EC76l5ukh1pZirHxP/VW9y9UZuaUAbeToba/iusjZey2E3z3aD7C96nxme7rX17i9v+mZ\nzLYP63FtpodBoAEC7D+6um868Pb7L6doBGnvOu6dRfUzFOQSYP88Nu6zeu9ZZyGu/5iz1pq29Nj1\n79Hnd6k8I27dmHPfrzquvVBP10uA7VroqHMUPvfjww+8uaeS3tvVLJl6nNXq6nm02UBX6iLAtlyk\nBq7lKU6ZYFXSPztU9nCvZqTXXj0ZebUuAuyAflrAW6+vfXttuA+Y3X9LdOLhXt+/L/B+5xZDM47Q\nVICNP0oXpGofB9bfhj9dumZvD/f6pitsvNxVtYpKCmMEuk07ATaO4/Rjb4Z9ebNcSXs9sRhPtXHg\nJ1dSJ115umqVjOZV+XJlpeG3DYQaTzzrqKw5wJ7+bDYGcEz96fCndAEu95Bqawcg/famqHs2gHzc\nFnhWMZ4Wmq6z5Rc9fs+5Vf2pPsfpU+G+vyJnPQP78ncdK8auf7VlU+tZVbHyo85qQjdc+g/fMI7j\nanicVYYhf6bYzhLiYTXsP7xzz8g9L3ZvlHa+a+ZSeyv29B5x4gZ0TtRAfQowCvjt7i8lw6qK2/W8\nuWGQqqo26Eo7ATbv3Zi9nX2/7WaNdbzaXs6uk383TtpuUFs18tHeNc/DgtI6pZzfaCfAhmFYdiHu\n/FdrfxvUCG57KrPRsY9c+Xho1j1SVle2tM9dY25Qgz/mq6dTOdp4FR81FWCn66QRXOppQFx/dKzC\nDzi90lyFbzQTgRH6CrCnW9GIm83o/vD0dtotG59uWKKpf3Pwo7MqJKK1M8tqogX1FWCvtJOrLRn2\n5bbd6CAv7uqqC3oyVIPi1dVMb+oxwOam02R/W+kYBZvsfDrfckbfFk1endPtvaYfa7WZce2jt1Wx\n8eUXjx8W3QXY46JWJ301jutyhSK1ekUifh8eGlgzuguwYdg3G/h4t3XWvZjbuidxTysf3TCVOVYn\nK/8qopJ/OzGHPvUYYKe7bqg6a5CK6+GPE+UD0+X1u4G4EzJ/ezl7G97H78+dmlRS8rNa1ilNtId7\nYgFWu7MGqUp6+HYWey9y3aCWNVy+fYPH97TVOwkw6rV368edyt7eHv7Vx95d/s1PrtPb0m5/CdGT\nm5a26giwNkU/QOKtp+v4zRh0yoIkRbQUP98TYA06cMwglXu8gl/ekTy9u/z7H9gANRBKgH1223LB\nFQeF6JktWa7ml3ckj+8uv+4p43Vzhaducuy3PHU3N3mJBFjVDnequXPaAdGS15nT9z9tGM58ynjP\nrd657XnpJsf+7dN/VtjX2r6Fbf8TmbtVYV+qxBT7CL7Oz5nMtTSEA1VxaSM66yzE10JmNvxfCTCg\nXx2meEsv2RLiV664nWn4qMYIdS4EBbmzAbtYnRNg37puVUfPJM5tDbiGO7zHBTpRWoQA2+TmhyY1\ndE44Zu4sl47m5+5nKeLOIeXpA/laIsCqU/PxE2046yzE38ag9oaJvW5rvTV0k2/2Mf7+M0/7iU9V\nVEONnUiAwfkaGyYqVMlN3uO7EbifXYi06YaFLBhCoqvVVQEzMICWNTxNNAMDaFx70TUzA9uq+fe0\nA2QRYF9p9b4GoH4CjO6cddAcXOS3Jmpf0hMBBkAkAbaP5148ij3XHlogwHYweecjeQa3EWAHWYzm\nNxoG3EOA7eYWO92JB83BFTTRjQTYPg0cgw3QBgF2kGUigLIE2G6VHIPNFcytuY0trN8TYPAfe+9O\nDENQigCjWfdHiySDOwkwOJPlZbiNAKM7zkKkcproRgIMTmDxEO4nwOBb3h0IRQgwOI2FH7jTn9IF\ngBaILrifGRjdcdAcldNENxJgAAF8AsYrAQZAJAEGQCQBBkAkAQZAJAEGQCQBRnccNEc93m6Y10Q3\nEmAAZcipLwkw+Mv7bCCLAKNlPi4ZGibAAIgkwOiOg+aonCa6kQADIJIAAyCSAAMgkgADIJIAAyCS\nAAOonXfZvyXA6I6D5qicJrrRn9IF2G15h8TrNX5884QWANC2sAAbx3FJpsc/L+QWQCdaW0Icx9Gb\n2AF6EDYD+2iegT1N1F6/AYB0VQfY3ux5+w0SiyfjOGoTVGL+wISnUert8xFeVR1guy6hSw7QlaoD\n7NU0Ta+7EOfoevtXMGx7D40npxAnLMCGd+G0fEVu8Wr7Z1pqPpCltV2IAHRCgAEQSYABEEmA0R3P\nSqmcJrqRAAMgkgADIJIAAyCSAAOomnfZ/0aA0R2fV0DlXpuoXR1vCTAAIgkwACIJMIDCrGofI8AA\nSvJ86zABBkAkAUb7Pn6iigUcSCTA6M7bg+Ys41APZyFuJMAAiCTAAIgkwACIJMAAiCTA6I6zEKmc\nJrqRAAMgkgADIJIAA6jXOHqT4q8EGACRBBgAkQQYAJEEGN1x0By1eTpvWhPdSIDRhY8H0gNxBBgA\nkQQYvbNNGUIJMAAiCTC646A5KqeJbiTAAIgkwACIJMAAiCTAACpli+w6AUavPCeHcAKMXrmzhXAC\njO5M0yS9qNCyKOAsxI0EGL1wHCI1k1kHCDAAIgkwgLo4iWOjP6ULACXZpkxV5uiapuknw7TONQKM\nroyek1OxcZqmYRzHcZzb5ziOMmyFJUR6MY8F87jwcIcLpBJgdEqGUS3LAxsJMIBK/OemallI5Dee\ngdEdgwI1GwfnnG1lBkYvpmmaB4f/fgVq8ff57Pyodpps3/jIDIyu/FuikV6QzgyMjvzsQJykF1Wz\nhLiNAAMgkiVE+uUYDmqkUW5mBgZQHYuIWwgwgBqZiX0kwABq4VPrdhFgAEQSYPTFHS7100Q3EmB0\nyhhBzTwA28I2evpljKBCmuV2ZmB3qORjO2ooRg1lGBSjsjIMilFZGYZqirFOgNGd+UCp0qUAvtVg\ngEXcOFDaaKEG0jX1DEx0AfSjwU/8fPwYU5EGcEz96dDUDOxV/RcAgGNSA+xpaiWoAHqTGmASC6Bz\nDe5CBKAHDW7iAKAHqUuIWyzPye4P6fVf/bhPskgZ7qyZ337XzVcn5Yq8/dt7ilHDFbn52fbHVnFD\nGVaKUXAEWwpQ+Qyn2QB72kx/52VY+dW3bev/+PLnr1xdM+vFuKcMH4txz0XZeEVKFeNxsC54Re7s\nLB9bRdk+UnAEG3LegOQZ2K2maarhjkYZFvXcY47jWHzUmMtQQ4VUUoxuVTJSfdTsDIyPio8Rxcfr\nqtw2+1kvQNkyVOW2JcQtxeAtAdajSnrmMmSXKsD8q5f/L1ghxa9FVYonaNnlu9k0TUs/lWG/sYTY\nqeIDRMHfvph+DEUrpJLaoCpzcLqzWdfyWkE9uxBft3IU2do0/96yG70en43fU4D1Yix/W/CKvP2r\n24pR2xW5c8bjinwsWOUBUXv5AOAtS4gARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABE\nEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQS\nYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARPo/amzZrMvqN6gAAAAA\nSUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo3()"]}, {"collapsed": false, "outputs": [], "prompt_number": 30, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["__Exercise 4__\n", "\n", "Display the evolution of roots as the noise level $\\sigma$ increases."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3debBk\nZ32f8eft9W4zd0YzSCPNaEFCLDIgY4xdyEBCCCYkxk5wxXalkipTjrN4wYQ4gcSETMWhyjg2pjBx\nuRyXE4pyUrbjNTY23oihgKQsGwHBxggJSWiZkWY0y917O/nj3NPTd5nlzu3uc95znk/dGmmWe/vt\n7vO+3/N73/ecDkmSIElSbGp5N0CSpOthgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmS\nomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJk\ngEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJ\nkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSomSASZKi\nZIBJkqJkgEmSomSASZKiZIBJkqJkgEmSotTIuwFS2YQQxv4zkyQZ+8+UYhfsGNJ1uHJKjbdbTfOx\npIgYYNLV7RohReg7hW2YNAUGmLSLncEQS0+Jt+XSXhlg0qZtQ385ukYpn5SUMsBUaVUb30efb+mf\nrErPAFPlVC20LscwU+wMMFWF4/Xl+MooUgaYSs7ReU+GL5evlYrPAFM5mVv7ZJKp+AwwlY0j7xh5\nHqAiM8BUEubWRPnyqoAMMEXPsXWafLVVHAaYYuVImiNffBWBAab4OHoWhG+E8mWAKSaOmMWUvi++\nKZoyA0xxMLqKz/dIU2aAqegcFuPi+6WpMcBUXEUYCk/u+HSSk3aZa+CkoqagkXcDpF0UIbq0H+kb\nZ4xpomp5N0DaIoQwHPWKOfD95G4fgqxdpW/i8D2VxssKTEVRwKpr5/whsDT9dkRuNMOK8+aqBKzA\nVAjFrLp2Xe6amX474jd8Z63GNEYGmHIW3bl5G/6nQ/B1GY2xvNuiMjDAlKdiFl6jDo50khoswCz0\n82xR9FwY07i4jV75iKjq+ukQepBADdowCwtwGF4XQ+OLLKJjQMXkJg7lIK6R6xBsQAIBmtCEBsTR\n9GJzq732ySlETVU6cVTkOcOdjsICzMM8zMAsNCHAp50BGwc3d+i6OYWoKSngLvlr95FsFjFAHRrQ\nhjbMwb0RPp1ishTTXhlgmoYSjE2/H0IaYAEa0IAWtOEA3BXz8yqUEhwnmibXwDRx5RiV5qEHQIAa\n1KEGAbo5t6tUXBXTnliBaYJKNhL93xAGQDaRWB+pw55bludYECU7cjQhBpgmpZRj0J+PTCTWsrnE\nJrTgRLmeae6iXjTVdLgLURNRyvQCZqGV5VZjZCIROO0murHyth26KiswjVnpT5y/NDKRWMvWw4al\n2JGSPusclfVkSPtngGmcKjLWPJxNJLI1w+rQhMWyP/3pq8hxpb0ywDQ2lRplHssKzWEdVssyrAEL\n1XgRpqn0lb2ugwGm8ahUeqWe3C3DhtOJM1V6KaamgoeZrsAA0xhUdlh5+jIZlsZYs3ovyBRU9mDT\nTgaY9sXR5NmRjfWj1zinGVarQb+6L86EeNQpZYDp+jmOpC6EwNYMGy3FaEC36i/R2HnsCa8D03Vz\nBBlaTJLhNWG1kSTblEDbK5nGzEvEhBWYro/ptdN6VoeRxViaZ5uFWDqleNFXbJw8DivOANOeOWpc\nTm/Hetj2ZbH0SrEzvnTjlH7CXN6tUA5847U3ptdVhJC+NJeWxWo7MqwOLXjC13BsPCyryQDTHjhM\nXKt62JxMDFvnE+tZkg3vPfUVX8zx8OCsIANM18oBYm/SjRuj84k7d9nXs1vZ/6Wv6hg4l1g1vt+6\nJqbX9VgIDOcTR/fX17bu7Eg/VWwWDsAf+wrvixlWKW6j19WZXtdpOWEloTkyc7h9i31mAD1Yhde7\nL3xfkiQJIbi9viIMMF2F6bVf5xMuJJsThvWRImxnkiXQhW918N0XLxGrDgNMV2J6jc2Z5NJHYe5M\nsrQjpkPuAL7HwXdfPGIrwgDTZZleY/ZkwlMJTyY0tybZ8Gs4wdiHD5lh+5LOJebdCk1WI+8GqKBM\nrwl6JAG4J8DIrvp0R2Iz215/Id8mlkGaYR7DJWaAaRem1zT8RQLwmkAv25fYhhmYh1lo5928sjDD\nSswA03am11R9PAH4h4GNbD/9IhyFBSCA78K+DDclejyXkgGmLUyv6Qvhr5JfgK/COjRgDo7CzXk3\nqyzMsBIzwLSd/XyaQvg0L3xlaCTJSwMr0Ic2HIUTnOHo0bybVw5mWFkZYLrEHj5lIfwC/C1ugpvg\nlfBYVoQdhuM8xF0G2LiYYaXkNnptcs9xHubTxS5O800HPnn2a44MXl4b3Ftbuu3A/fWv/wSvzrt5\npWJ0lY8VmGDn0tdM+OIGfwoPbf1nJx0CxmwFznL6IF/kU1+67x89/8Pfwm8f5cwqc5/nJT/J2/Nu\nXtm4sb5kfC9FCOFdyY9+Hz9z88ee4n/Af+XBHp+F/7fbPzbDxiiE/wivhtfwhsB98M3wcuabKyvd\n+aSZd+NKym1KZeIUYqWF93GQpfck//btvO/mjz3F78IfcqHHs3Au77ZVQZK8C/4P/CEfvcgH4G3w\napYxvSbIOyWWiVOIFRWOwVvgB3kdf/RNfPLw6XN8FZ6Bi3RgA9bzbmFFJMk78m5C5bihozQMsMoJ\n4UPwHbx+ljvhJo7zxCIXWIYVWINZWtnWgrN5N1WakCzD3g9N6MN6kvzrvBulPXMKsYJexMIsR+Ag\nwLPccIFFFmAeDsI8i3AYjkFrx3e6AKYS+Tm4Fe6Au+CuEP5z3u3Rnhlg1RLCB2EOYAmegdM8yS0P\nc+fqTXPcCifgTriV2+EeeAUcGfle00ulEcIH4Ha4HW6F2+FOuDvvRmnPnEKsmiassgLn4GH4HJ+6\n5b7DnGuz8R2v/eUGPZpwkPoT3L7C7fO8/hOGlkrpOXAL3Agz0IT1zRkJRcVlzGoJ4f3wYngBN93K\nC+FeuA/uhhcBJLN5t0+aihB+D76RxmEOQAs6cAEGn0uSl+bdNO2BAVY5IXwI7oA7CbdyCNokT+Xd\nJmm6Qrif9tdzCyzCDHTgIny5lyROSsXEAKucEAL8NCxAG55Jkrfm3SIpB+EOeC7pORw9WIfHSD6X\nc6u0J55uVEt2G4IfyLshUt5OwAk4AunM+ZqrYPExwCrEm+hIlxyD2+EozEIdOgZYfAywajG9pE1H\n4BY4BovZPo7DfqJQZAywqvDmb9IWB+E43AaHYBY68Bz4sOkVEwOsEpw8lLZbhONwJ43FXosOsHpi\njjmLsJj4VlWCfVLaKZzi8LFzN/L0DOt1+h1aj3HbxbCIZ3uRsAIrPycPpV3Vj/Wfy1eO88Q8K202\nVpg/ypmPJYldJhYGWCV4OintdJCLt/PoHTxyhLNtNjZon+UI/HWctIiEAVZy9kPpcm7k6Tt45G4e\nvInT86xs0D7PobfzU4lFWCQMsDKzE0pX8EVe8EE+eA9/cZwnZlgHznNogWX4F2mGefJXcAZYydkD\npSu4mwe/hi8c3TjDMsCJhcePts8M/9YMKzgDrLQsv6Srei5fOfrEGR4nDTAWOHbiFMch+9TmXFun\nqzDAysyTR+nK7uZBHoCvwDo0YAbOkwZYyiKsyAywcrLXSdcinE74U3gcVqANC7AMXxu4KcEirPAM\nMEkV9rvwZ3AKOtCCQ7AOn4VvvvRPPB0sLAOshOxv0rX6Dfg8XIQuzMA89Lf8vUVYkTnSlZABJl2r\nO0LnUS5CD2bgUAOOwyPbu499qph8V8rGniZdu6UQHsgKsDm4FV7Ugg0DLA5OIUqqrj+CB0Z+ezcc\n6HBixz/zuuZiMsBKxT4m7ckDW3/7IBxilwBTMdXyboAkFciDl/lzd3MUkAFWHpZf0p6c3C2Qzk+/\nHbpeBpikijq5xxM+i7CiMcBKwvJLGou9pppyZIBJqq5tcXUt6WURVhyetpdB2qN8K6UpcLajOKzA\nSsIeJU2NRVhBGGDRsy9J0+TJYnEYYGVgj5KmzBPHIjDAJGlvPGUsCAMsbq4nS6osA0ySroeziLkz\nwCJm/5Hy4sxHERhgcbMXSTnyJDJfBpgkXQ9PH3NngMXK7RuSKs4Ak6Tr5yxijgywKNlnpCJwFiRf\nBlis7DmSKs4Ak6R9cUYkLwZYfOwtUnE4F5IjAyxK9hlJMsAkab+cF8mFARYZ+4lUNM6I5MUAi4+9\nRZIwwCRJkTLAJGm/kiRxen/6DLCY2EMkacgAi4wLYJKUMsAkaTycI5kyA0ySxsDZkekzwKLhB4BJ\n0igDTJIUJQNMksbGZbBpMsAkaTyc5J8yAywOntZJ0jYGWDQ8uZOkUQaYJI2T8yVTY4BJ0tg4UzJN\nBpgkKUoGWAS8hFmSdjLAJElRMsAkTcSqexk0Yc5NRcApREWkH8IqdGADejCAGtxWsQPYPjsdjbwb\nIKlEamEDOrAOG9CFDRjAsyF8rQO6xs3ThKJLrynxbVIE5sNglRVYhTVYh/WsFOtDHxrw16pxJFuB\nTYdrYBGwJygO/c3/JjCADnRgDdZgCc7DU/CLLoxpfDxNKDpP5RSHGwLrsMYKLMMSrLBZja3CCqzB\nBqxBF7rwrlIf1Xbb6fBVLjp7giJwR2B5c8aw12UZzsFFuAirWZilqba+9ftOlvTYtttOh1OIkvat\nDwEa0KTRZBFOQB1qMIA+dLMlsW1OlndG0TsiToEBJmnfatCAOjRhhjBPc4EX38Ac1AHowyDnJk6V\n5dd0uI1e0v7cE6hl6VWDdOiuwdPJywD49RAcaDQJVmCF5iyEiu5VgTo0NucPaUEbZmDh0j/5e0ky\nB8382qiycqWx0FwKVqG9KXAGVqCT3XVjWH49tMtx+2MhVGQTB3beqfAlLjT7gIro3YFn4Bycg+Vs\na3wvW+YK0ITPX/a4HW7cKHF6YeedCl/iQrMPqCg+HDb3xS9v/XU1u+lhf6T8esCD1s47Da6tSrqM\nT4XNlFrOLkhO7w01yIqtdPPhAEK23TDAbJ5NVqUYYJK2ejRwIcut4U15hzVWutAVsvRK92bUsz+s\nw8csOzQlBpikzFLgQjYr2M2u3hpkuZWqsblpvpXd/HC4db4OHzG9ND0GmCTWmJtdX2Mdetl2jGFi\nhay0GkZXWnilM4dpcQb8qtG1nctgk2aASZX2l9xzjFOLbFyqt8hmCMPW0GpAe+Qf1LMS7ecco3eR\nJInXcU6aASZV16e57wTLdfq1NJfCyE2h0jKrB+1shpAszNL1sH9nbilnBphUUb/PG25iNY2uAbVa\nbbB5O6getLZWWk2YAeDbDS0ViAEmVdGH+O7b6dQYDKj1qXdptlsbW67lakMHgBcbWiooA0yqolnW\ngDS6OrQa9EI9ac11LhVbC+aWis4AkyrnffzLE9Cn3qG1QbvGICF0aDVqPdosts/n3UDpmhhgUhV1\naa4x26QLdGitMF+nv0H7Xh7Iu2nStTLACiqE3+QV38Z7kvAbAMnfzbtBKot/ys89l9Y6M8P0atAD\nljjwJn4r79ZJe2CAFVEIf8Drvo3/xA0ve/YWnjzIxedx00OfuWvz8wGlfejRWGN2mQVG0muF+e/h\n5/NumrQ3BlgxvZ7v4vjLnngNHz/GqRnW15l59GW3w5vzbpiit8L8Mgt1+h1aTbqBpEPrX/HjebdL\n2jMDrHBC+AA3v5Xn8Ro+/jI+c5hzs6x1aN3BI99K87d4U94NVNwucnCelT71Nhs1Bj0a5zmUd6Ok\n62GAFU6SvDUEeD4v4fN38dAhzqeTPGvMttnAANP+nOVIk+5w8nCD9qPcHt5P8ra8WybtkQFWTJ/i\n1+677Qceu5mnFlgenikf5GLeDVP0nuJmYIkDDXoDaivMP8RdfCPh++B/weNPwV/BV+DJJPmRvBsr\nXYkBVkRJcl8In73xB54+wtkDLM2wXqffp36Apbybpuh9tXtrv1k/wFJ6G44LLD7TfQ4J3Awvh/7N\nPDUL8zAXwk8kyQ/n3V7psgywgkqSez/J6kEuLnJhhvUwSIDF2oVHueN2Hsm7dYpY0iR0b5ltrtUY\n9Kmvd2d4EpahATNwCE4dIjkKF/GEScVmgBXXAZYOsDTbyz6laUCtPrh9/lHftL06OfKpFif9fCbg\nz1k7Prv5iZRrcBpOwxok0IAWbLSgBc0Q/kOSvDvv5kq78/PWCq1Po77UZxW62ecHNuBm37I9OLnj\nM5nMMCD8NhwEYAPOwGl4GB6GB+Fh6D0GD8NXkuQtOTc0Zn6g5aR5Ml9o9X6fNViDjexzb2vQDdxm\nr7gmO9NLqeRbCO+FWejBKjwDT8NZWILeBqzBOqzn3UzpSgywYuvBcDDpZUXYOtyWd8MUv+QdAOGN\n0IIlOAen4akBPAlPw+kk+ed5t1G6EivcQgshJA/BBVjPZhETqMFBeJlv3NXtWoG14d942GvynEKc\ntFreDdDV9KAD61kptgqrcDbvVkVi1+Wu+vTbIWkCnEIsvBZ0oZstg6Uf9O6JxzVrZKuHw9960GsK\nguuvk+dAWHh3JAygk9Vh61kR9st2j2syk+4Hhya0oQ3NvJukinD+cNI8GY1BDfrQgS70IO0UnZwb\nFYtZqGeFa3oZggEmlYMVWAxeldDPZhE7sAHrsAb/xSLs6lpZ4ZV+pRfoSioBK7B49KA3shcR6Obc\noii0swoWCFDzoJfKwgosEm9OGGQB1hmpw37KIuwq0tWv4TJY0ylEqSw8GY3HICvCeiNFWO8q36Q2\nkK2BATUDTCoLK7Ciu7QZ95/sVoStwY9ahF1JY2v51TLApLKwAiu0JEm2XE0ygH5WhCXZnaUswq4o\nPcTrWcnqGphUGlZgUamPBFg3u0nHBrzTIuyyGiNfTS9klkrEAIvK25PNIqyfZVj6tZp3wwqsDvUs\nt+oGmKbC23BMhwEWm/rWicRhEfaDdpjd1bdGV8N7IWoqvA3HFBhgsTmZkIwUYcPpRIuwy6hvzbC6\nFZhUFgZYhGrZDo7R9bAuvMUibBf1HV8e9FI52JcjsH0+/SeylbBtGbaRU/uKrba1CKs7hSiVhQFW\ndLvPpKdFWD/76mYXh32nRdgWD4QwWngZYFKZGGBx+pmEJJtI7I3cpGMD3myGXVIb+aobYFK5GGDR\n+m/Zbo5tG+v9mJUR26Ir/X/vxKGJcg/91BhgMQtcqsNGY+yN9p9NYUcRVoMXu79ZE+Ye+ukwwGL2\niwmMBNhgpA77m2YYZIm17UtSOdid43DZSYlfTbYUYaMx9loz7FIFNpxC9CIwqTQMsAhcZToizald\n5xJfXfUMq3kRmFReduf4/WayS4YN67BXBl5V3RgLTiFqutzBMU1251L4nWRzQwe7xdgGfENFO9XO\n9PKI16S5g2Nq7M5l8dHk0qbEnatiPfi6KmaYASaVmN05Dts/2XJXf5TNJbJbjPXhpYF7qxVjIZtF\nDAaYVDp253L53yMZxm6lWBfuCdxTlRgbTa9ggEnlYncunU9s3dMxujA2ev/f5weeX/4YcwpR0xRC\ncAFsmuzOZfTJZLPiSA3rsG0rZH24K3BXmWMs7CjCyvxspYrxfCEa6RrY3t6vl4fNqiuNq2GMbQuz\n1KMlPBKeCmH0KaYxdsxjXpNhBTZlVmDRuJ6O8Wc7SjG2TiqOVmMnAsfLVp9YgWlqvAJs+jxfiMn1\nn9+9JGx+6srgMnXYYGu8AafKcGA8HcJg5LdpgN3oMa8JsPyaPl/xmOy3h7wwbNmROJphOycV01/P\nxH14PDMyhZieHtfgqMe8JsAAmz5vbRqT9Gqw6+8kX0wA7g7ApTt3sPV/tsXY4UAC52PtlmHkV5xC\nlMrFNbDqeTChvmNoZ7cYG9ZqBwMHAgfjG/xrI0tfwQDTxLgAlgsrsEp6MAur23b0ustVY+mfzwUC\nrERTkIUd/+8wowlx/nD6rMDiM85zvccSHk82j4IrV2Ojn9XSDszEEQSjN+AY7kictBAuhjlCOD35\nh5IqzVXH+ExqrfjmsH1/x85titt2eQToFvr4ubjjOrAAByd5zIdXwO2wAOfgK/D59yTJj0zu4VQE\n13ONpsbBAIvPxDc7HQ3XFGOMbO/rF/QoWg5hsCPAajA/mRcwvApeC4ehAevwBHyE5MFJPJQKxP2H\neXENLEqT7TDp1vl0/yEj+xV3XAG9+QcJIYRQyA68kCRLW6dcJ7uP42/Di2AO6tCBm6A7uQeTqs4A\ni881fbTK/p1LABYDgx0ZxmZ6DSuxARBCo5AZVhtpNZMOsG+A4zADwDoswMrkHkyqOgNMV3QhG/zT\n/YfhUmrtXBHrhBCgWbAYCzt+O8EAeyEcgTZ0YQPmYG1yD6ZCcP4wRwZYrKbdbVYTgPb2i6C3xRjQ\nD2GmSP25NhKxk95J3z6xscByg16v3Vhtz601Zzk7sQeTKs8Ai9KUZhF32kgA6psPvWNRbNNqCHOF\nybBaOsM5kl4BehOY8KwzuIvH2mwEkh6NZRaenb1hZXF+vI+iQvH65XwZYNq7fkItDCcS2S3DVkKY\n0E6/vWolSSeEbRsRJzHqvITPH+J8k26feo/GDOsDaisYYCXn/GGODLCI5Tn5PkhqkJaB224ITDZZ\ntxTCgWL07eHi3USnEG/hyQMs1el3aa4z06Tbpfn4Xccn82iSDLBo5TaLOKKeJEm4NJ248wqxgqiN\npBeTqcDeyO/dyuMNekCf+hqzNQYXWBz346hAcu+AMsDilvsOqEaSDLLP3Np268QA50I4XIAirLZ1\nenPzhlLNMMbbiDyPL6flVzp/2KSbEJ7mxnH9fBWT84f5MsAiVoQiDGglSQuW0rvpjMwoUpgMCwGS\nrVVXGHOReAtPtugAfeobtAPJBu1V5sb5GJK2MsA0HjXob73zVKoQJ6iDJGS7TmD8Gzl+ircfYaVB\nb7T8WmH+s7x0bI+hgsl98kMYYCVQkI40nyQXsjvnDkZmEfNvWWqSN+S4idOj5VeNQc+eJU2e3Sxu\nBZlFTC0myZkQBjtmEU+FcCz3iK1tbRNjC7Bf49uPsFJjMKDWpTksvz7I94/nAVQ8xel0FWeAlUFB\nirDU6EeyUJwiLL2eeXQn4pg+GewYp2oMEkKfepdmjUGXZp3+eH66iqo4Pa7KDLDoFaoIO5okj4fQ\nL2aApZJxroF9hq87SDeQDKj1aGzQTghrzP4Q7x/DT1chFae7yQDTmCXZbo5iBVgd2Hox874rsC/w\n4iMspZOHPRodWkCHVijGM9bkWH4VhAFWBmkRVpBONYA+9EcCrA6PhnB7vs3beS3Yvk+jj3GqTn9A\nLV39Svdu1Om/mV/d74+WdA0MMI1ZHzpZBTa8h1M971ZttmCQ/XbfAbbG3CJdYLj0lRAa9AbjWltT\nIRXnTFEYYGVSkK7Vh+5IgKWVT/7H2Xi30SdhtrP5YxuNXj30E0IaY6/m4/tqp6Rrlv/AorEozlaO\n1awCI5s/ZKTyyU1t64Vp+wmwlcB69jMb0KTR7jVDt8agX4BSU5NTkC6mIQOsVIpQhK1CL9uyPvws\nrhfk3apLa2DDfRzXN9V3ent6pf9faw0Swov4i/G0VkWVe//SKAOsPApShK1BD4Aa1CEZ2wVX+5M2\nYnQN7Dqa9UhgBYA6NNi81qsGTRLCjZweQztVVEXoXNrGACubfIuwj4awkYVWLTu8mnm1ZujOsDm3\nF0b2lux1tu+BwEaWfA1owMxmegHtzbpMZWb5VTQGWKnkXoSdhUFWe2Wza8VYANtWgdXgq3sZjP4k\n0B+J5fTndLMnaTcqO8uvYrLnlVBeRdjPhtCHOtShldVhTbgv9/PW4YRhbWRz5LX7nbD5XenMYfps\n6tlWywY08n6CMXsmhCYMv0Ib2jAHi3AIDsMhWIQDcBAW4AD84xxecMuvAjLAyiavIuw9IQyybQ2t\nLCaaxbkNx3DvRrKX+cNf2lp4NUa+PS3CGnCoEE8xIh8JIZ1/bcIMzEI7+5qB9gaNbZfBp6//8O1L\n4GcDC1mYLcLXT/YtsPwqLAOsnKZZhJ0c6d592ID5bAoxgTcU4bw1PcxHb2/Vuobv+vkA2XxoPfve\nWnav4vTHnijAs4vEyRDmoAnzMA8LMLv1ji3DV7rRg/RrAzrQzX7twga0oAkb2SlFHe4Pm/XZcyZV\nDVt+FZMBVkK5r4R1RqYQC2F4Mdpw9u/KAfa+sHnKX89GSbLfDjfit+AFDmrXZHiKswpAd2QeNy3Z\n08XE0c8xCKMfK7fzazDy1YMOrEMDzsLB0J1tNumMq/GWX0VmgJVWjtsR+1l6vbkg563NbP/hcOBs\nwRsCh+CXRlr4rrBltrAxUh0ML2pLf0Ib7i3GUyu8kzsCYBVa0Ibe1vs+bxeyiqy+YyIxbJ1RTEZu\nwdmh2equ1ufPcfgUx17O/ftpfJpell+FZYCV09SKsJ3DE9lGxO8uTrdPA2x4d5BaFk4d+Adhc/ll\nuIugPjKghpHCa/gnry/M84rWts8rGP2qQQhZaTaXLotlb1Aj+3UYaaPfyaUwq9f7M6zfwLP384pn\neM4b+ch1N9X0KrKizPFo7KaTYSd3695NeGuhuv3sSD41s0WU4WpWP1tx6WU3EUn/ZHQWC0hH1TcV\n6XnFrDYyQdvK3pY0qpiFBbgRjmS7EA9kS2fzMAOt7F83s0irbYm0QBJIagyadOdZ+XXe/AF+aK8t\ndPKw+KzASi6XicQfKVR6AX+c8NqweeOMsHWw23V9JW3+aNVVhzl4c8GeV7Ra2Z7DOTgAC3D8Go+Z\n/x42v7+VFWfbSrQ61BhQSwhAjUEaY4tceC/vfAc/tqd2Wn4VXP63ztNETS3A0rnEXQuyonh9djnX\ncHdGMxsKmyMzVKOFWvo1C99b4OdVeNvmmRfgBHzXWA6VL2V5lm7Mn4E2vWZjlbn06yIHlzhwkYPn\nOHyWI6e56SHu+iW+46o/2NWvKBhgJWc/3OJN4dJthuvZ7FXj8uk1A2/zpRuDNMMa8K4JHoqhT71H\nY52Z9GuVuRXmlziwzMI5Dp/n0BmOnuamR7jjo3zzVX5WAbCCSJYAAAnASURBVO6LravyTSo/u+J2\nfz/A8Jqj7Jq1YYylBdksvNMXLT4P8vw6/R6NDq00w1aZW2bhIgcvsHiWI2c4eopjj3DHA9x7uR9i\nl4mFa2Dll+7msENe8isJwPeGS/dtHJZfP+6rFLe7+RLwEf5Oi076YdmrzK0xu8Zsh1afevqZbQ16\ni1y8wMGdP8G9GxFxXKsEJxJVQR/kB9MYS6uxZRYusHiBxTMcfYqbH+O2s6eOJMe2f5dnexHxraoK\nu6Wq6Yf5yXRSMc2wJQ48yw1nOfIEx8+dO8z9JK+/9I891YuLg1qFmGGqpvv49DFOAR1ayywss/AM\nz3maG9fPzPCX8DGSd4PpFSHXwKrFDFMFfYpXAi/iiw16G7TXmL3IwfXODKvQg7lL/9LeEReHs2ox\nwFRlNZI6/T71ZDVwDp6Gr8IX4GPwB5Zf8bECqxwzTJU1IITPNZiFHizDs3AWVjb/1n4RHQOsWtIt\n9WaYqutX4MXQhA4swTPwNCwBp/JumfbMgayKDDBVWXgLHIMGrMN5OAVfhC+fSnbuqVexOZBVlBmm\nKgt/DQ4AsAyn4auw8oUk+Zqcm6U9chSrKHcMq8pCOMORozSgB+dgcBo+myRXuUGiisY1sIpyMUxV\nliRHQ/gs3ADAGjwBD3K1O/yqaAyw6prapzZLhfS7cDu0oAPPJsn3590e7Zkn4FVnEaYKcgq9HGp5\nN0D5sw5TpZhepWGAVV3ajc0wVYrpVQ4GmMwwVYhz5mVigAk8IVU1eJZWMgaYNrkpUeXm0lf5GGDa\nwgxTKZlepWSA6RIXw1RKpldZGWDawgxTKZlepWSAaTszTGXitsMSM8C0CzNsi2/0dYiVx3C5eW6i\ny6ruysF3B87DEjwNp2EDBtCABfhq9V6NaFX3AK4MA0xXUrkh4PsCS7AKy3Aevgwr0GGQkEA9wK3w\naGVejZhV7tCtJO9Gryup1qeufH9gDXrQgz5swCr9DTagA33oJzQf43DezdRVmV4VYYDpKqpygfPb\nAn1IYAADSIe+Dh1Yh3U2Y6wDj4bwtY6MBWZ6VYcBpqtLM6zMI8I7wqXQqmVfTZIBPejABizDKqzC\nCnw5hFVYhG8r8WsSJ9OrUgwwXZMyZ9i/DwwgjERXA5owTwj0k83ZxHVYgQvwLJyDZQA+k9WmJ0v5\nysTG9Koat9HrWg3Xw/JuyFi9NzCAWhZgdWjCDMzBAfgb1KAPPejCOizBhSy9VCimVwUZYNqDEl4f\nNhipuhrQgjbMwkE4Bn+YzGbTFP0syTq7/ZiTZXpNImR6VZNTiNqbUu1L/EDYjK50ASxAHVoQ4H2b\nz66ZJLfAcggNCJBAL982awfTq7KswLRnJanDPhw2a6902rCZ1V4Ll9Jr6PlJMg/NrEhTcZheVWaA\n6XqUIcN6ELKZw3Tdaxbm4Ud3HwpfmiRvSJIbYGG3v3UTRy7SmQDTq7KcQtR1in4ucTh5GKCf/eHV\nyqs3ZU92dNHL9Jo+Cy/hraS0T7GOI78VuAgd6LJ5/TJQh38W2xOppFiPOo2bFZj2JdY6LC2/0o0b\n6Q6OGrTzbpWugemlIQNM+zV6fVg0w0o9my2sjew//K5IGl9hkR1mmjADTGMw3NMRRyl2f9jcTVgb\n2UDv5sLCM720jQGmsYlmOrGezRbWRla/5vNska5suN+16IeWpssA0zjFkWHNrN5qZAHWhFcXuMHV\nZuGlyzHANGYRLIk1s9ttDLdvOH9YVIU+kJQ3A0zjV/QlsTYMoJFd/lWH2ZxbpF2ZXroyA0yTUtDp\nxG5gBhiZP2xc5u4ayo/RpWthgGmCCphhg2atNjPYvIlUqgk3FqV5wvTSNTPANFmjd00swpDUpdme\n3djcN+8CWPEU51BR8RlgmobilGLrzDTb3VptwACAQL9Vr+fbJgFGl/bOANOUFKQUW2O2Tr/V7NQY\nAD0aK8wfyas1yuR+YChGBpimKt9S7BO85jhzgaRDq8YgIfRonOeQAZYjo0vXzQDTtOVYiq0xu8xC\nQmjSDSQJoUPrbr40zTZoyPtraJ8MMOUjl1KsQ2uF+YRQpx9IBtTWN/fUa9osvLR/BphyM/1SbIX5\nNhs9GsMAW2VuCo+rURZeGhcDTDmb2q2n3sl7n8eBJt0OrTp9oEvzPIcm94jaycJLY2SAKX+jt55i\nYqPbBu0lDgSSNht1+glhg/bT3DiJx9JORpfGzgBTUUz6DoorzF/g4IDaDOvpDo4V5t/Bj439gbSN\nc4aaEANMxTK5hbE1Zi+y2KfRohNIetSXODjGn6+djC5NlAGmIhr7wtidfOUeDjfobdBu0gW6NM/i\nBWCTYnRpCgwwFdR4F8YucvAch2sM0jtxpBvoH+O28bRVI4wuTY0BpkIbV4ytMfssNwyozbJWYzCg\ntsSBP+PrxtnWyjO6NGUGmCKwzxgLfwL3zZ1tHunQatGp0+/SPMfhibS1kowu5cIAUzSGg+Oeh8s1\nOMMzx56zHmYa9IA1ZtefmuHmyTS0Sowu5cgAU3z2XJA9C88ALM0doAl9WCYxvfbH6FLuDDDF6hpj\nLHwrvBFOQRdmoA4deAbLr+tjbqk4DDDF7erziqtwGlqwwuade9fhy/C6aTazDIwuFY0BppK4bEF2\nBp4A4AK0oA8XSN6dVzPjY26psAwwlcqOguxTnHglj0MP5rL5w1M5NjAaw9zC6FJRGWAqp6wg+wOe\nhBm4yOYHp6zC/fk2rdDMLUXEAFO5PcngER6+g7MwAz04S9LPu1HFY24pRgaYyu0peJBBh3ML0IB1\neDiE1w7/uuKDtbmlqE31A92l6Qvhp+E4zAJwMUm+c+Svwui/rEhfMLRUGgaYBGUPs3I/O1WWASZt\nt224J8IR38RSFRhg0lUUP892tpDiNVIaOwNM2rNdAyM10Q6V1+NKxWSASWNzhYAZFzusNGSASZKi\nVMu7AZIkXQ8DTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIU\nJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUD\nTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0yS\nFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQl\nA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUJQNM\nkhQlA0ySFCUDTJIUJQNMkhQlA0ySFCUDTJIUpf8Pn/vvnmKIatYAAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 31, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo4()"]}, {"collapsed": false, "outputs": [], "prompt_number": 32, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Bibliography\n", "------------\n", "<html><a name=\"biblio\"></a></html>\n", "\n", "\n", "* [Schmidt] R. Schmidt, _Multiple emitter location and signal parameter estimation_, IEEE Trans. on Antennas and Propagation, vol. 34, no. 3, pp. 276-280, 1986.\n", "* [LiaoFannjiang14] W. Liao, A. Fannjiang, _MUSIC for Single-Snapshot Spectral Estimation: Stability and Super-resolution_, 2014.\n", "* [Prony] G. R. B. de Prony, _Essai Experimentale et Analytique_, J. de L'Ecole Polytechnique 2, pp. 24-76, 1795.\n", "* [KrimViberg] H. Krim and M. Viberg, _Two decades of array signal processing research: the parametric approach_, IEEE Signal Processing Magazine 13(4), pp.67-94, 1996.\n", "* [RoyKailath] R. Roy and T. Kailath, _ESPRIT ? Estimation of signal parameters via rotation invariance  techniques_, IEEE Trans. Acoust., Speech, Signal Proc., vol. 17, no. 7, July 1989\n", "* [HuaSarkar] Y. Hua and T. K. Sarkar, _Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise_, IEEE Transactions on Acoustics, Speech and Signal Processing 38(5), pp.814-824, 1990.\n", "* [DemNeedNg] L. Demanet, D. Needell and N. Nguyen, _Super-resolution via superset selection and pruning_, Proceedings of the 10th International Conference on Sampling Theory and Applications, 2013."], "metadata": {}, "cell_type": "markdown"}], "metadata": {}}], "nbformat": 3, "metadata": {"kernelspec": {"name": "matlab_kernel", "language": "matlab", "display_name": "Matlab"}, "language_info": {"mimetype": "text/x-matlab", "name": "matlab", "file_extension": ".m", "help_links": [{"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md", "text": "MetaKernel Magics"}]}}}