-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path5_DFT-Implementation.py
361 lines (361 loc) · 98.4 KB
/
5_DFT-Implementation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The DFT: Numerical Aspects\n",
"\n",
"In this notebook we will look at some numerical issues associated to the DFT; in particular we we will look at the differences in precision between the \"naive\" way of computing the DFT and the FFT algorithm.\n",
" \n",
"As a quick reminder, the definitions of the direct and inverse DFT for a length-$N$ signal are:\n",
"\n",
"\\begin{align*}\n",
" X[k] &= \\sum_{n=0}^{N-1} x[n]\\, e^{-j\\frac{2\\pi}{N}nk}, \\quad k=0, \\ldots, N-1 \\\\\n",
" x[n] &= \\frac{1}{N}\\sum_{k=0}^{N-1} X[k]\\, e^{j\\frac{2\\pi}{N}nk}, \\quad n=0, \\ldots, N-1\n",
"\\end{align*}\n",
"\n",
"The DFT produces a complex-valued vector that we can represent either via its real and imaginary parts or via its magnitude $|X[k]|$ and phase $\\angle X[k] = \\arctan \\frac{\\text{Im}\\{X[k]\\}}{\\text{Re}\\{X[k]\\}}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Direct Implementation\n",
"\n",
"### Numerical errors in real and imaginary parts\n",
"\n",
"The DFT can be easily implemented using the change of basis matrix ${W}_N$. This is an $N\\times N$ complex-valued matrix whose elements are \n",
"\n",
"$$\n",
" {W}_N(n,k)=e^{-j\\frac{2\\pi}{N}nk}\n",
"$$\n",
"\n",
"so that the DFT of a vector $\\mathbf{x}$ is simply $\\mathbf{X} = W_N\\mathbf{x}$. Note that the inverse DFT can be obtained by simply conjugating ${W}_N$ so that $\\mathbf{x} = W_N^*\\mathbf{X}$.\n",
"\n",
"We can easily generate the matrix ${W}_N$ in Python like so:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# first our usual bookkeeping\n",
"%matplotlib inline\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"plt.rcParams[\"figure.figsize\"] = (14,4)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def dft_matrix(N):\n",
" # create a 1xN matrix containing indices 0 to N-1\n",
" a = np.expand_dims(np.arange(N), 0)\n",
" \n",
" # take advantage of numpy broadcasting to create the matrix\n",
" W = np.exp(-2j * (np.pi / N) * (a.T * a))\n",
" \n",
" return W"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Let's try out on a short signal and verify the invertibility of the transform"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-8.88178420e-16+2.66453526e-15j 1.77635684e-15+1.48029737e-16j\n",
" 0.00000000e+00-1.70234197e-15j]\n"
]
}
],
"source": [
"x = np.array([5, 7, 9])\n",
"\n",
"# DFT matrix\n",
"N = len(x)\n",
"W = dft_matrix(N);\n",
"\n",
"# DFT\n",
"X = np.dot(W, x)\n",
"# inverse DFT\n",
"x_hat = np.dot(W.T.conjugate(), X) / N\n",
"\n",
"print(x-x_hat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see, the difference between the original vector and the \"reconstructed\" vector is not exactly zero. This is due to the small numerical errors that accumulate in the $N^2$ multiplications and additions needed by the direct and inverse transforms. \n",
"\n",
"While minor in this case, this numerical imprecision can be very annoying if we switch to a magnitude/phase representation, as we will see now. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Numerical errors in phase\n",
"\n",
"Let's first define a more interesting signal such as a length-128 step signal: \n",
"\n",
"$$\n",
" x[n] = \\begin{cases}\n",
" 1 & \\mbox{for $0 \\leq n < 64$} \\\\\n",
" 0 & \\mbox{for $64 \\leq n < 128$}\n",
" \\end{cases}\n",
"$$\n",
"\n",
"Conveniently, we can compute its DFT analytically (it's just a geometric series) and we have \n",
"\n",
"$$\n",
" X[k] = \\begin{cases}\n",
" 64 & \\mbox{for $k=0$} \\\\\n",
" 0 & \\mbox{for $k \\neq 0$, $k$ even} \\\\\n",
" \\frac{(-1)^{(k-1)/2}\\,e^{-j\\pi\\frac{63}{128}k}}{\\sin(\\frac{\\pi}{128}k)} & \\mbox{for $k$ odd}\n",
" \\end{cases}\n",
"$$\n",
"\n",
"From this it's easy to compute the phase; we will set its value to zero whenever the magnitude is zero (i.e. for even-indexed values) and we have\n",
"\n",
"$$\n",
" \\angle X[k] = \\begin{cases}\n",
" 0 & \\mbox{for $k$ even} \\\\\n",
" -\\pi + \\frac{\\pi}{128}k & \\mbox{for $k$ odd}\n",
" \\end{cases}\n",
"$$\n",
"\n",
"However, let's see what happens if we compute all of this numerically:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAD8CAYAAABO44oNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFO1JREFUeJzt3X+wpXV9H/D3x2Uxq9EQy6bRXXTJhJhsrA3pHdDaSYkxZbEGnNZWGNMaa0IzDSZtLCnUDtm10zGJmaZ2JKmMMSatFQlSsuNgSUbNtOOIZQmNCHSbLUZZILIaIU0l8sNP/7gHvF7O3XsunOWer/t6zTCc53s+PJ/P7jPPvffN85znVncHAABgVE/b7AEAAACeDKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADC0Ezar8cknn9y7du3arPYAAMCCu+mmm77Q3dvXq9u0ULNr164cOHBgs9oDAAALrqo+O0ud288AAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABjauk8/q6r3JHlVknu7+0VT3q8k70jyyiRfTvJj3f0H8x70WLn25rvy9usP5u77HsjzTtqWi89+YZI8bu3Vp+940rXz2Mci91vk2fQ7fmbTb/x+ALBR1d1HL6j6gSR/nuS31gg1r0zypiyHmjOTvKO7z1yv8dLSUm/2I52vvfmuXHrNLXngoUceW9v6tEoqeeiRr/29bNu6JX/3r+3IB2+66wnXzmMfi9xvkWfTz7HWb5x+b/s7f0WwAeAxVXVTdy+tW7deqJnsbFeSD60Rat6V5Pe7+/2T7YNJzurue462z0UINS/7hY/mrvsemKl2S1UemeHv6mi189jHIvdb5Nn0O35m02/sfjtO2paPX/LymfYBwDe+WUPNPD5TsyPJnSu2D0/Wpg11YVUdqKoDR44cmUPrJ+fuGQNNkpm/UR+tdh77WOR+G6nVb+x+G6nVT7+N1G7k6zIAPGoeoaamrE39btXdV3T3Uncvbd++fQ6tn5znnbRt5totNe2PubHaeexjkfttpFa/sfttpFY//TZSu5GvywDwqHmEmsNJTlmxvTPJ3XPY7zF38dkvzLatW75ubevTKlu3fP03221bt+SCM095UrXz2Mci91vk2fQ7fmbTb/x+jz5YAAA2Yt2nn81gf5KLqurKLD8o4P71Pk+zKB79MOrPXf2pPPjIV7NjxZN6Vq+9+vQdWXrBc55U7Tz2scj9Fnk2/Y6f2fQbvx8AbNQsTz97f5Kzkpyc5PNJfj7J1iTp7v8weaTzO5PsyfIjnd/Q3es+AWARHhTwqNe+6xNJkg/845cedW0etd/o/RZ5Nv2On9n0G78fACSzPyhg3Ss13X3BOu93kp/awGwAAABzM4/P1AAAAGwaoQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGhCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaEINAAAwNKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKHNFGqqak9VHayqQ1V1yZT3n19VH6uqm6vqU1X1yvmPCgAA8Hjrhpqq2pLk8iTnJNmd5IKq2r2q7F8luaq7T09yfpJfnfegAAAA08xypeaMJIe6+47ufjDJlUnOW1XTSZ49ef0tSe6e34gAAABrO2GGmh1J7lyxfTjJmatq9ib53ap6U5JnJnnFXKYDAABYxyxXamrKWq/aviDJe7t7Z5JXJvmPVfW4fVfVhVV1oKoOHDlyZOPTAgAArDJLqDmc5JQV2zvz+NvL3pjkqiTp7k8k+aYkJ6/eUXdf0d1L3b20ffv2JzYxAADACrOEmhuTnFZVp1bViVl+EMD+VTWfS/JDSVJV35PlUONSDAAAcMytG2q6++EkFyW5PsntWX7K2a1V9daqOndS9uYkP1FVf5jk/Ul+rLtX36IGAAAwd7M8KCDdfV2S61atXbbi9W1JXjbf0QAAANY30y/fBAAAWFRCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaEINAAAwNKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGhCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoc0UaqpqT1UdrKpDVXXJGjV/v6puq6pbq+o/z3dMAACA6U5Yr6CqtiS5PMkPJzmc5Maq2t/dt62oOS3JpUle1t1fqqpvO1YDAwAArDTLlZozkhzq7ju6+8EkVyY5b1XNTyS5vLu/lCTdfe98xwQAAJhullCzI8mdK7YPT9ZW+q4k31VVH6+qG6pqz7QdVdWFVXWgqg4cOXLkiU0MAACwwiyhpqas9artE5KcluSsJBckeXdVnfS4/6j7iu5e6u6l7du3b3RWAACAx5kl1BxOcsqK7Z1J7p5S8zvd/VB3fybJwSyHHAAAgGNqllBzY5LTqurUqjoxyflJ9q+quTbJDyZJVZ2c5dvR7pjnoAAAANOsG2q6++EkFyW5PsntSa7q7lur6q1Vde6k7PokX6yq25J8LMnF3f3FYzU0AADAo9Z9pHOSdPd1Sa5btXbZited5Gcn/wAAADxlZvrlmwAAAItKqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGhCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaEINAAAwNKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGgzhZqq2lNVB6vqUFVdcpS611RVV9XS/EYEAABY27qhpqq2JLk8yTlJdie5oKp2T6l7VpKfTvLJeQ8JAACwllmu1JyR5FB339HdDya5Msl5U+r+dZJfSvIXc5wPAADgqGYJNTuS3Lli+/Bk7TFVdXqSU7r7Q3OcDQAAYF2zhJqastaPvVn1tCS/kuTN6+6o6sKqOlBVB44cOTL7lAAAAGuYJdQcTnLKiu2dSe5esf2sJC9K8vtV9cdJXpJk/7SHBXT3Fd291N1L27dvf+JTAwAATMwSam5MclpVnVpVJyY5P8n+R9/s7vu7++Tu3tXdu5LckOTc7j5wTCYGAABYYd1Q090PJ7koyfVJbk9yVXffWlVvrapzj/WAAAAAR3PCLEXdfV2S61atXbZG7VlPfiwAAIDZzPTLNwEAABaVUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGhCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaEINAAAwNKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGgzhZqq2lNVB6vqUFVdMuX9n62q26rqU1X1kap6wfxHBQAAeLx1Q01VbUlyeZJzkuxOckFV7V5VdnOSpe5+cZKrk/zSvAcFAACYZpYrNWckOdTdd3T3g0muTHLeyoLu/lh3f3myeUOSnfMdEwAAYLpZQs2OJHeu2D48WVvLG5N8eNobVXVhVR2oqgNHjhyZfUoAAIA1zBJqaspaTy2s+tEkS0nePu397r6iu5e6e2n79u2zTwkAALCGE2aoOZzklBXbO5Pcvbqoql6R5C1J/mZ3f2U+4wEAABzdLFdqbkxyWlWdWlUnJjk/yf6VBVV1epJ3JTm3u++d/5gAAADTrRtquvvhJBcluT7J7Umu6u5bq+qtVXXupOztSb45yW9X1f+sqv1r7A4AAGCuZrn9LN19XZLrVq1dtuL1K+Y8FwAAwExm+uWbAAAAi0qoAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaEINAAAwNKEGAAAYmlADAAAMTagBAACGJtQAAABDE2oAAIChCTUAAMDQhBoAAGBoQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwNCEGgAAYGhCDQAAMDShBgAAGJpQAwAADE2oAQAAhibUAAAAQxNqAACAoQk1AADA0IQaAABgaDOFmqraU1UHq+pQVV0y5f2nV9UHJu9/sqp2zXtQAACAaU5Yr6CqtiS5PMkPJzmc5Maq2t/dt60oe2OSL3X3d1bV+Ul+Mclrj8XAAHzjuvbmu/L26w/m7vseyPNO2paLz35hXn36jg2tJzkmtfodP7PpN3a/RZ5txH6vPn3HE/yK/tRaN9QkOSPJoe6+I0mq6sok5yVZGWrOS7J38vrqJO+squrunuOsAHwDu/bmu3LpNbfkgYceSZLcdd8DufSaW3Lgs3+aD95010zrF//2HyaVPPRIz7VWv/n2W+TZ9Bu73yLPNmK/S6+5JUmGCDa1Xu6oqtck2dPdPz7Z/gdJzuzui1bUfHpSc3iy/X8mNV9Ya79LS0t94MCBOfwRnrzfuOBN+fYjd2b3c5/92Npt9/xZknzd2lrrG6mdxz4Wud8iz6bf8TObfmP2u/lz9+UrDz+S1aoq075XrbU+zTxq9Tt+ZtNv7H6LPNso/e74lh1514vPS5LsOGlbPn7Jy2fax7FQVTd199J6dbNcqakpa6v/dmapSVVdmOTCJHn+858/Q+unxnOe+fQ84/4tX7f2jBO3TK2dtr6R2nnsY5H7LfJs+h0/s+k3Zr9pgSbJmt+QN3IzwDxq9du8Wv30O1a1+q1fe/d9D8y8j800y5WalybZ291nT7YvTZLuftuKmusnNZ+oqhOS/EmS7Ue7/WyRrtQAsPle9gsfzV1TvnluqcojU76drLU+zTxq9Tt+ZtNv7H6LPNuI/Ua5UjPL089uTHJaVZ1aVScmOT/J/lU1+5O8fvL6NUk+6vM0AGzExWe/MNu2fv0VnW1bt+SCM0+ZeX3r0ypbt9Tca/U7fmbTb+x+izzbiP22bd3y2IMFFt2WvXv3HrVg7969X923b98fJXlfkjcl+U/d/cGqeuu+ffuetXfv3oP79u27Jcnr9u3b97Yk35fkJ/fu3fulo+33iiuu2HvhhRfO508BwPC++7nPzs5v3ZZb7ro/f/4XD2fHSdty2Y/szj/5we+ceX3vud+bv7X72+deq9/xM5t+Y/db5NlG7HfZj+ze9IcE7Nu37569e/desV7durefHStuPwMAAI5mnrefAQAALCyhBgAAGJpQAwAADE2oAQAAhibUAAAAQ9u0p59V1ZEkn92U5tOdnOQLmz0ET4hjNzbHb2yO39gcv7E5fmNz/Gbzgu7evl7RpoWaRVNVB2Z5XByLx7Ebm+M3NsdvbI7f2By/sTl+8+X2MwAAYGhCDQAAMDSh5muu2OwBeMIcu7E5fmNz/Mbm+I3N8Rub4zdHPlMDAAAMzZUaAABgaMd9qKmqPVV1sKoOVdUlmz0PR1dVp1TVx6rq9qq6tap+ZrL+nKr6var6o8m/v3WzZ2W6qtpSVTdX1Ycm26dW1Scnx+4DVXXiZs/I2qrqpKq6uqr+1+Q8fKnzbwxV9c8mXzc/XVXvr6pvcv4trqp6T1XdW1WfXrE29VyrZf9+8rPMp6rq+zdvcpI1j9/bJ187P1VV/6WqTlrx3qWT43ewqs7enKnHdlyHmqrakuTyJOck2Z3kgqravblTsY6Hk7y5u78nyUuS/NTkmF2S5CPdfVqSj0y2WUw/k+T2Fdu/mORXJsfuS0neuClTMat3JPmv3f3dSf5qlo+l82/BVdWOJD+dZKm7X5RkS5Lz4/xbZO9NsmfV2lrn2jlJTpv8c2GSX3uKZmRt783jj9/vJXlRd784yf9OcmmSTH6OOT/J907+m1+d/IzKBhzXoSbJGUkOdfcd3f1gkiuTnLfJM3EU3X1Pd//B5PX/zfIPVDuyfNx+c1L2m0levTkTcjRVtTPJ307y7sl2JXl5kqsnJY7dAquqZyf5gSS/niTd/WB33xfn3yhOSLKtqk5I8owk98T5t7C6+78l+dNVy2uda+cl+a1edkOSk6rquU/NpEwz7fh19+9298OTzRuS7Jy8Pi/Jld39le7+TJJDWf4ZlQ043kPNjiR3rtg+PFljAFW1K8npST6Z5C939z3JcvBJ8m2bNxlH8e+S/FySr062/1KS+1Z8kXcOLrbvSHIkyW9MbiF8d1U9M86/hdfddyX55SSfy3KYuT/JTXH+jWatc83PM+P5R0k+PHnt+M3B8R5qasqax8ENoKq+OckHk/zT7v6zzZ6H9VXVq5Lc2903rVyeUuocXFwnJPn+JL/W3acn+X9xq9kQJp+9OC/JqUmel+SZWb5laTXn35h8LR1IVb0ly7fTv+/RpSlljt8GHe+h5nCSU1Zs70xy9ybNwoyqamuWA837uvuayfLnH73UPvn3vZs1H2t6WZJzq+qPs3yr58uzfOXmpMntMIlzcNEdTnK4uz852b46yyHH+bf4XpHkM919pLsfSnJNkr8e599o1jrX/DwziKp6fZJXJXldf+33qjh+c3C8h5obk5w2efrLiVn+kNb+TZ6Jo5h8BuPXk9ze3f92xVv7k7x+8vr1SX7nqZ6No+vuS7t7Z3fvyvK59tHufl2SjyV5zaTMsVtg3f0nSe6sqhdOln4oyW1x/o3gc0leUlXPmHwdffTYOf/Gsta5tj/JP5w8Be0lSe5/9DY1FkdV7UnyL5Kc291fXvHW/iTnV9XTq+rULD/w4X9sxowjO+5/+WZVvTLL/7d4S5L3dPe/2eSROIqq+htJ/nuSW/K1z2X8yyx/ruaqJM/P8jfvv9fdqz9gyYKoqrOS/PPuflVVfUeWr9w8J8nNSX60u7+ymfOxtqr6viw/6OHEJHckeUOW/weZ82/BVdW+JK/N8m0vNyf58Szft+/8W0BV9f4kZyU5Ocnnk/x8kmsz5VybBNV3ZvnJWV9O8obuPrAZc7NsjeN3aZKnJ/nipOyG7v7JSf1bsvw5m4ezfGv9h1fvk6M77kMNAAAwtuP99jMAAGBwQg0AADA0oQYAABiaUAMAAAxNqAEAAIYm1AAAAEMTagAAgKEJNQAAwND+PwIq4Vko3fsbAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x22336f98390>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N = 128\n",
"x = np.zeros(N)\n",
"x[0:64] = 1\n",
"\n",
"plt.stem(x);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAD8CAYAAACsP5F0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGc9JREFUeJzt3X2MZfdZH/Dv49m1M3EaJsYbY4+T2pGsDSlusrAKcV0hcGjXQBqvTFKCKF2lRk4lSkNLN3hBKrbUykFb8dZSVCsBtmrq2HKMbVGBsZwgWlScrFkah5itQwDHuyYeIMvrNvHu/vrH3LX3ZWbn3tm5c+9v5/ORVnPPmefe+5w5c879fee8bLXWAgAA0JOLJt0AAADAqAQZAACgO4IMAADQHUEGAADojiADAAB0R5ABAAC6I8gAAADdEWQAAIDuCDIAAEB3Nq3nm11++eXtmmuuWc+3BAAAOvLkk0/+SWtty0p16xpkrrnmmuzfv3893xIAAOhIVf3RMHVOLQMAALojyAAAAN0RZAAAgO4IMgAAQHcEGQAAoDvreteyafDQgUPZ++jBHD5yNFfNzWb3jq3ZuW1+0m0BAAAj2FBB5qEDh7Lnwady9MXjSZJDR45mz4NPJYkwAwAAHdlQp5btffTgSyHmpKMvHs/eRw9OqCMAAGA1NlSQOXzk6EjzAQCA6bShgsxVc7MjzQcAAKbThgoyu3dszezmmdPmzW6eye4dWyfUEQAAsBob6mL/kxf0f+CBT+crx09k3l3LAACgSxsqyCSLYebeTz6bJLnvfTdMuBsAAGA1NtSpZQAAwIVBkAEAALojyAAAAN0RZAAAgO4MFWSqaq6qHqiq36uqp6vqhqq6rKoeq6pnBl9fM+5mAQAAkuGPyPx0kl9trb0xyZuTPJ3kjiSPt9auS/L4YBoAAGDsVgwyVfXqJN+U5MNJ0lr7SmvtSJJbkuwblO1LsnNcTQIAAJxqmCMyb0iykOQXqupAVX2oqi5NckVr7fkkGXx97Rj7BAAAeMkwQWZTkq9P8nOttW1J/jojnEZWVbdX1f6q2r+wsLDKNgEAAF42TJB5LslzrbUnBtMPZDHYfLGqrkySwdcXlnpya+2e1tr21tr2LVu2rEXPAADABrdikGmt/XGSL1TV1sGstyf5bJJHkuwazNuV5OGxdAgAAHCGTUPW/UCSj1TVxUk+n+S9WQxB91fVbUmeTfLu8bQIAABwuqGCTGvtd5JsX+Jbb1/bdgAAAFY27P8jAwAAMDUEGQAAoDuCDAAA0B1BBgAA6I4gAwAAdEeQAQAAuiPIAAAA3RFkAACA7ggyAABAdwQZAACgO4IMAADQHUEGAADojiADAAB0R5ABAAC6I8gAAADdEWQAAIDuCDIAAEB3BBkAAKA7ggwAANAdQQYAAOiOIAMAAHRHkAEAALqzaZiiqvrDJH+Z5HiSY6217VV1WZL7klyT5A+T/OPW2pfG0yYAAMDLRjki8y2ttbe01rYPpu9I8nhr7bokjw+mAQAAxu58Ti27Jcm+weN9SXaefzsAAAArGzbItCS/VlVPVtXtg3lXtNaeT5LB19cu9cSqur2q9lfV/oWFhfPvGAAA2PCGukYmyY2ttcNV9dokj1XV7w37Bq21e5LckyTbt29vq+gRAADgNEMdkWmtHR58fSHJLyV5a5IvVtWVSTL4+sK4mgQAADjVikGmqi6tqr918nGSf5jkM0keSbJrULYrycPjahIAAOBUw5xadkWSX6qqk/X/vbX2q1X1qST3V9VtSZ5N8u7xtQkAAPCyFYNMa+3zSd68xPw/TfL2cTQFAABwLudz+2UAAICJEGQAAIDuCDIAAEB3BBkAAKA7ggwAANAdQQYAAOiOIAMAAHRHkAEAALojyAAAAN0RZAAAgO4IMgAAQHcEGQAAoDuCDAAA0B1BBgAA6I4gAwAAdEeQAQAAuiPIAAAA3RFkAACA7ggyAABAdwQZAACgO5sm3cC0eOjAoex99GAOHzmaq+Zms3vH1uzcNj/ptgAAYGx6HgMLMllcgXsefCpHXzyeJDl05Gj2PPhUknSzIgEAYBS9j4GHPrWsqmaq6kBV/fJg+tqqeqKqnqmq+6rq4vG1OV57Hz340go86eiLx7P30YMT6ggAAMar9zHwKNfIvD/J06dM/3iSn2ytXZfkS0luW8vG1tPhI0dHmg8AAL3rfQw8VJCpqquTfEeSDw2mK8lNSR4YlOxLsnMcDa6Hq+ZmR5oPAAC9630MPOwRmZ9K8oEkJwbTX53kSGvt2GD6uSRLnkhXVbdX1f6q2r+wsHBezY7L7h1bM7t55rR5s5tnsnvH1gl1BAAA49X7GHjFIFNV70jyQmvtyVNnL1Halnp+a+2e1tr21tr2LVu2rLLN8dq5bT5333p9Lp5Z/HHMz83m7luv7+IiJwAAWI3ex8DD3LXsxiTvrKpvT/KKJK/O4hGauaraNDgqc3WSw+Nrc/x2bpvPvZ98Nkly3/tumHA3AAAwfj2PgVc8ItNa29Nau7q1dk2S9yT5eGvte5J8Ism7BmW7kjw8ti4BAABOMcpdy870w0n+dVV9LovXzHx4bVoCAAA4t5H+Q8zW2q8n+fXB488neevatwQAAHBu53NEBgAAYCIEGQAAoDuCDAAA0B1BBgAA6I4gAwAAdEeQAQAAuiPIAAAA3RFkAACA7ggyAABAdwQZAACgO4IMAADQHUEGAADojiADAAB0R5ABAAC6I8gAAADdEWQAAIDuCDIAAEB3BBkAAKA7ggwAANAdQQYAAOiOIAMAAHRHkAEAALqzYpCpqldU1Ser6v9U1e9W1V2D+ddW1RNV9UxV3VdVF4+/XQAAgOGOyHw5yU2ttTcneUuSm6vqbUl+PMlPttauS/KlJLeNr00AAICXrRhk2qK/GkxuHvxrSW5K8sBg/r4kO8fSIQAAwBmGukamqmaq6neSvJDksSS/n+RIa+3YoOS5JPPLPPf2qtpfVfsXFhbWomcAAGCDGyrItNaOt9bekuTqJG9N8rVLlS3z3Htaa9tba9u3bNmy+k4BAAAGRrprWWvtSJJfT/K2JHNVtWnwrauTHF7b1gAAAJY2zF3LtlTV3ODxbJJvTfJ0kk8kedegbFeSh8fVJAAAwKk2rVySK5Psq6qZLAaf+1trv1xVn03y0ar6d0kOJPnwGPsEAAB4yYpBprX26STblpj/+SxeLwMAALCuRrpGBgAAYBoIMgAAQHcEGQAAoDuCDAAA0B1BBgAA6I4gAwAAdEeQAQAAuiPIAAAA3RFkAACA7ggyAABAdwQZAACgO4IMAADQHUEGAADojiADAAB0R5ABAAC6I8gAAADdEWQAAIDuCDIAAEB3BBkAAKA7ggwAANAdQQYAAOjOikGmql5XVZ+oqqer6ner6v2D+ZdV1WNV9czg62vG3y4AAMBwR2SOJfmh1trXJnlbku+vqjcluSPJ462165I8PpgGAAAYuxWDTGvt+dbabw8e/2WSp5PMJ7klyb5B2b4kO8fVJAAAwKlGukamqq5Jsi3JE0muaK09nyyGnSSvXevmAAAAljJ0kKmqVyX5WJIfbK39xQjPu72q9lfV/oWFhdX0CAAAcJqhgkxVbc5iiPlIa+3BwewvVtWVg+9fmeSFpZ7bWruntba9tbZ9y5Yta9EzAACwwQ1z17JK8uEkT7fWfuKUbz2SZNfg8a4kD699ewAAAGfbNETNjUm+N8lTVfU7g3k/kuSDSe6vqtuSPJvk3eNpEQAA4HQrBpnW2v9KUst8++1r2w4AAMDKRrprGQAAwDQQZAAAgO4IMgAAQHcEGQAAoDuCDAAA0B1BBgAA6I4gAwAAdEeQAQAAuiPIAAAA3RFkAACA7ggyAABAdwQZAACgO5sm3UCPHjpwKHsfPZjDR47mqrnZ7N6xNTu3zU+6LQAANriNNE4VZEb00IFD2fPgUzn64vEkyaEjR7PnwaeS5IL9JQEAYPpttHGqU8tGtPfRgy/9cpx09MXj2fvowQl1BAAAG2+cKsiM6PCRoyPNBwCA9bDRxqmCzIiumpsdaT4AAKyHjTZOFWRGtHvH1sxunjlt3uzmmezesXVCHQEAwMYbp7rYf0QnL5T6wAOfzleOn8j8BX43CAAA+rDRxqmCzCrs3Dafez/5bJLkvvfdMOFuAABg0UYapzq1DAAA6I4gAwAAdGfFIFNVP19VL1TVZ06Zd1lVPVZVzwy+vma8bQIAALxsmCMyv5jk5jPm3ZHk8dbadUkeH0wDAACsixWDTGvtN5L82Rmzb0myb/B4X5Kda9wXAADAslZ7jcwVrbXnk2Tw9bVr1xIAAMC5jf1i/6q6var2V9X+hYWFcb8dAACwAaw2yHyxqq5MksHXF5YrbK3d01rb3lrbvmXLllW+HQAAwMtWG2QeSbJr8HhXkofXph0AAICVbVqpoKruTfLNSS6vqueS/FiSDya5v6puS/JsknePs8mePXTgUPY+ejCHjxzNVXOz2b1ja3Zum590WwAAdMJ4cmkrBpnW2ncv8623r3EvF5yHDhzKngefytEXjydJDh05mj0PPpUkfvkAAFiR8eTyxn6x/0a299GDL/3SnXT0xePZ++jBCXUEAEBPjCeXJ8iM0eEjR0eaDwAApzKeXJ4gM0ZXzc2ONB8AAE5lPLk8QWaMdu/YmtnNM6fNm908k907tk6oIwAAemI8ubwVL/Zn9U5egPWBBz6drxw/kXl3mQAAYATGk8sTZMZs57b53PvJZ5Mk973vhnPWurUeAMDGMMq4b5Tx5EYiyEwJt9YDANgYjPvWhmtkpoRb6wEAbAzGfWtDkJkSbq0HALAxGPetDaeWTYmr5mZzaIlf3nPdWs81NQAA02PYsdlqxn2czRGZKTHqrfVOnlt56MjRtLx8buVDBw6tQ7cAAJxqlLGZWyqvDUdkpsSot9Y717mVSz3H0RsAgNENO4YaZWzmlsprQ5CZIqPcWm+UcyvdGQMAYHSjjKFGve7FLZXPn1PLOrXcOZRLzR/1zhgPHTiUGz/48Vx7x//IjR/8uNPVAIALxijjnFHGUKOMzVgbgkynRjm3cjVHb4a99kboAQAmbdjxyKjjnFHGUK57WX9OLevUKOdWjnJnjFHO7xz1lDXX6QAAwxp23DDKeGTUa4xHGUO57mX9OSLTsZ3b5rPt9XP5xmsvy2/ecdOyG8q4jt6Mcrh1NXdZc7QHAC4s4zhyMsp4ZNTrWEY9yjLs2Iy14YjMBjCuozdrFXrO968rJ+uHPdrjyBAALG/Uz8lJHzkZZTwy6v/f4ijLdBNkNohh74yxe8fW03YyyfJ/eRhX6EnGd4rbNAUk4QtgYxjX/n4ctav5nOwpnIwyzjnJ3cWmlyDDaUb5y8O4Qk8yvqM90xSQegtfPYY6ARDWTo/b9TTst8a1vx9X7ahnUPQWThxhubAIMpxl2L88jCv0JOM72jMNAWmcrz0NH4LTEuqmJQCO87V7q52WPizf6LU9btfTst8a1/5+XLWjnkHRYzhxhOXC4WJ/zsuwF7Xt3Dafu2+9PhfPLP7Kzc/N5u5br1+TGxSMct/2UWrHFZDG+dqjXPA4DbXT0scoF5Wu5hbl43jt3mqnpQ/Lt7raHrfradlvjWt/P67aUf8vlFHmj/LZPuq4wUX2G9N5BZmqurmqDlbV56rqjrVqigvTKDuZUXZgo+wYpyEgjfO1p+FDcFpC3TQEwHG+dm+109KH5VtdbY/b9bTst8a1vx9X7ah36RJOmKRVn1pWVTNJfjbJP0jyXJJPVdUjrbXPrlVz02DSh/OnpXYSfSx16Pdcr3vmIeUkufGDHx+qdue2+bNe+1veuCUfe/LQkoe1z6d2946tyx4y/5Y3bjmr51Fqlzts/1Wzm4euvWpu9qye5165OV/6mxfXvHatej7f5VuqLlkclJxP7e4dW5cd8Bw6cvSsnte79nyXb9p/FpZv7Wp73K6nZb+1Fvv78/3cGfUz6u5brx/qc3I1n8GjfrYvZRrGRRfSeKtn1Vpb3ROrbkhyZ2ttx2B6T5K01u5e7jnbt29v+/fvX9X7raXv+i//O8nZ50WeOf/Mc2CTxQ3/7luvX3IjXO/andvmT+t5XLUbefm+8xvmc/+nnjtr53y+tXffen2S03f8y33QjFL7nd8wf9b8zRdVUsmLx9uKtWvxGj2+XyVZak84N7s5Xz52YtW1s5tn8orNFy05QDrzdSZRe77LN+0/C8u3drU9btfT8n7nu79fq8+dUT+jevu8vtDHI+u5fNOgqp5srW1fse48gsy7ktzcWvu+wfT3JvnG1tq/WO450xJkfuG7fyBfs/CFvOnKV582/7PP/0WSvDT/wLNH8uVjx896/iWbZnLJ5osmXrvt9XOn9TyuWsvX1/K97rLZ/P7CX6e1lks2zeR4azl2/MRQta+7bDZf+LOjS772ppmLcvxEW/Pa8+15rZbvTBddVLmoasnXHqV208xFOdFaTpxYeV+7nrVrtXyj1Fq+taud1PL1tl1Py36rx8+oafg8m5afxUZZvj/e8rq8997/eFbNJAwbZM7nGplaYt5Ze+Squr2q9lfV/oWFhfN4u7Vz2aWX5JUXz5w1/5UXz5w2f7nBzpePHZ+K2jN7HlftWvZs+dam53PVXv6qS3LFqy/J13zVK7Lt9XPLDl6Wqr38VZcs+9rHjp8YS+359rxWy5ckMxct7tYu2TSTN1x+6TkHfsPWHjt+Im+4/NLT6peznrVrtXyj1Fq+/pevt+16WvZbSX+fUdPweXZm7Vr2bPnO7uOySy9ZsmaabchTy4Z14wc/vuQ5sPNzs/nNO27aMLXT0oflG3/ttPRh+cZfOy19WL7V1U5LH9NQOy19WL7x105LHxf68k2D9Tgi86kk11XVtVV1cZL3JHnkPF5v6ozrbli91U5LH5Zv/LXT0oflG3/ttPRh+VZXOy19TEPttPRh+cZfOy19XOjL15OZO++8c1VPvPPOO0/cddddzyT5SJIfSPLfWmsfO9dz7rnnnjtvv/32Vb3fJLzxylfn6tfM5qlDf56/+n/HMj83m3/7j9605B0eLuTaaenD8vlZWD4/C8s3XX1MQ+209GH5/CwulOWbBnfdddfzd9555z0r1a361LLV6O3UMgAAYH2tx6llAAAAEyHIAAAA3RFkAACA7ggyAABAdwQZAACgO+t617KqWkjyR+v2hud2eZI/mXQTrJr11zfrr2/WX9+sv75Zf/2y7ob3t1trW1YqWtcgM02qav8wt3VjOll/fbP++mb99c3665v11y/rbu05tQwAAOiOIAMAAHRnIweZeybdAOfF+uub9dc3669v1l/frL9+WXdrbMNeIwMAAPRrIx+RAQAAOrUhg0xV3VxVB6vqc1V1x6T74dyq6nVV9Ymqerqqfreq3j+Yf1lVPVZVzwy+vmbSvbK0qpqpqgNV9cuD6Wur6onBuruvqi6edI8srarmquqBqvq9wTZ4g22vH1X1rwb7zc9U1b1V9Qrb3/Sqqp+vqheq6jOnzFtye6tFPzMYy3y6qr5+cp2TLLv+9g72n5+uql+qqrlTvrdnsP4OVtWOyXTdtw0XZKpqJsnPJvm2JG9K8t1V9abJdsUKjiX5odba1yZ5W5LvH6yzO5I83lq7Lsnjg2mm0/uTPH3K9I8n+cnBuvtSktsm0hXD+Okkv9pae2OSN2dxPdr2OlBV80n+ZZLtrbWvSzKT5D2x/U2zX0xy8xnzltvevi3JdYN/tyf5uXXqkeX9Ys5ef48l+brW2t9N8n+T7EmSwTjmPUn+zuA5/3kwRmUEGy7IJHlrks+11j7fWvtKko8muWXCPXEOrbXnW2u/PXj8l1kcSM1ncb3tG5TtS7JzMh1yLlV1dZLvSPKhwXQluSnJA4MS625KVdWrk3xTkg8nSWvtK621I7Ht9WRTktmq2pTklUmej+1varXWfiPJn50xe7nt7ZYk/7Ut+q0kc1V15fp0ylKWWn+ttV9rrR0bTP5WkqsHj29J8tHW2pdba3+Q5HNZHKMygo0YZOaTfOGU6ecG8+hAVV2TZFuSJ5Jc0Vp7PlkMO0leO7nOOIefSvKBJCcG01+d5MgpO3bb4PR6Q5KFJL8wODXwQ1V1aWx7XWitHUryH5I8m8UA8+dJnoztrzfLbW/GM/35Z0l+ZfDY+lsDGzHI1BLz3LqtA1X1qiQfS/KDrbW/mHQ/rKyq3pHkhdbak6fOXqLUNjidNiX5+iQ/11rbluSv4zSybgyupbglybVJrkpyaRZPRzqT7a9P9qUdqaofzeKp8h85OWuJMutvRBsxyDyX5HWnTF+d5PCEemFIVbU5iyHmI621Bwezv3jyMPrg6wuT6o9l3ZjknVX1h1k8jfOmLB6hmRuc6pLYBqfZc0mea609MZh+IIvBxrbXh29N8gettYXW2otJHkzy92L7681y25vxTCeqaleSdyT5nvby/3ti/a2BjRhkPpXkusFdWy7O4oVWj0y4J85hcE3Fh5M83Vr7iVO+9UiSXYPHu5I8vN69cW6ttT2ttatba9dkcVv7eGvte5J8Ism7BmXW3ZRqrf1xki9U1dbBrLcn+Wxse714NsnbquqVg/3oyfVn++vLctvbI0n+6eDuZW9L8ucnT0FjelTVzUl+OMk7W2t/c8q3Hknynqq6pKquzeJNGz45iR57tiH/Q8yq+vYs/lV4JsnPt9b+/YRb4hyq6u8n+Z9JnsrL11n8SBavk7k/yeuz+IH97tbamRdJMiWq6puT/JvW2juq6g1ZPEJzWZIDSf5Ja+3Lk+yPpVXVW7J4o4aLk3w+yXuz+Ecw214HququJN+VxVNaDiT5viyeh2/7m0JVdW+Sb05yeZIvJvmxJA9lie1tEE7/UxbvePU3Sd7bWts/ib5ZtMz625PkkiR/Oij7rdbaPx/U/2gWr5s5lsXT5n/lzNfk3DZkkAEAAPq2EU8tAwAAOifIAAAA3RFkAACA7ggyAABAdwQZAACgO4IMAADQHUEGAADojiADAAB05/8DjXdNS8Si49UAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2233939e630>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAD8CAYAAAChIeEzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X2MHdd53/HfoxVFryUEK4VsYy65IY0ItB3R8aYLWyqLImbdUIlcaavEiAM1cZ0EtIGkcYKUyrJCW6d1IAYs0jaJG0SIXaeAwKiQaFqIkjBypCItEUledl1SsszUNWNKS8WiIq9lSwu+LJ/+ce817+7OvXvPzpk7c2a+H0AQ7+xw7lnO6zPnOc8xdxcAAAAApOiashsAAAAAAOtFQAMAAAAgWQQ0AAAAAJJFQAMAAAAgWQQ0AAAAAJJFQAMAAAAgWQQ0AAAAAJJFQAMAAAAgWQQ0AAAAAJJ1bRlfumnTJt++fXsZXw0AAAAgASdOnHjF3TevtV4pAc327ds1OztbxlcDAAAASICZfW2Q9Ug5AwAAAJAsAhoAAAAAySKgAQAAAJAsAhoAAAAAySKgAQAAAJCsUqqcAcCgjs7N69Cx0zq3sKgtY6Pav3enpifHy24WAACoCAKakvGwBvR2dG5eB46c0uKlJUnS/MKiDhw5JUmcJwAaKeS5gWcMNAUBTYl4WAP6O3Ts9HfOj47FS0s6dOw05wiATHV+iA95buAZA03CGJoS9XtYQ39H5+a1++AT2jHzmHYffEJH5+bLbhIKcG5hMWg5gGbrPMTPLyzKdfUhvi73iJDnBp4x0CS5Axoze5OZPWNm/8fMnjOzX4vRsCbgYW196n7DwlVbxkaDlgNotro/xIc8N/CMgSaJ0UNzQdIed/8BSe+SdLuZ3Rphu7XHw9r61P2Ghav2792p0Q0jy5aNbhjR/r07S2oRgCqr+0N8yHMDzxhoktwBjbd8u/1xQ/s/z7vdJuBhbX3qfsPCVdOT47r/7l26bqR1qRofG9X9d+8i/xtApro/xIc8N/CMgSaJUhTAzEYknZD0fZI+6e5Px9hu3XUeyu59+KQuLl3ReM0GLxZly9io5jOCl7rcsELVeQCs1DpPDj9zVpL00EduK7k1AKps/96dywbCS/V6iA95buAZA00SJaBx9yVJ7zKzMUmfNbNb3P3Z7nXMbJ+kfZI0MTER42trgYe1cHW/YYWgig0AXNWEh/iQ5waeMdAUUaucufuCpP8h6faMnz3g7lPuPrV58+aYX4uGIQ3pKsYTAcBy05PjmpwY03t23KTjM3saeW8AmiZ3D42ZbZZ0yd0XzGxU0vsk/UbulgF98NaphfFEAACg6WL00LxF0pNmdlLSFyQ97u5/FGG7ANZQ9wGwAAAAa4lR5eyku0+6+zvd/RZ3/3cxGgZgbVSxAQAATRelKACAcjRhACwAAEA/BDRA4hhPBAAAmixqlTMAAAAAGCYCGgAAAADJIqABAAAAkCzG0AAAACA5R+fmdejYaZ1bWNQWiuI0GgENAAAAknJ0bl4HjpzS4qUlSdL8wqIOHDklSQQ1DUTKGQAAAJJy6Njp7wQzHYuXlnTo2OmSWoQyEdAAAAAgKecWFoOWo94IaAAAAJCULWOjQctRbwQ0AAAASMr+vTs1umFk2bLRDSPav3dnSS1CmSgKAAAAgKR0Bv7f+/BJXVy6onGqnDUaAQ0AAACSMz05rsPPnJUkPfSR20puDcpEyhkAAACAZBHQAAAAAEgWAQ0AAACAZBHQAAAAAEgWAQ0AAACAZFHlDACABjk6N69Dx07r3MKitlDqFkANENAAANAQR+fmdeDIKS1eWpIkzS8s6sCRU5JEUAMgWaScAQDQEIeOnf5OMNOxeGlJh46dLqlFAJAfAQ0AAA1xbmExaDkApICUMwAAGmLL2KjmM4KXLWOjJbQGqWIcFqomdw+NmW0zsyfN7Hkze87MPhajYQAAIK79e3dqdMPIsmWjG0a0f+/OklqE1HTGYc0vLMp1dRzW0bn5spuGBouRcnZZ0q+4+9sl3Srp583sHRG2CwAAIpqeHNf9d+/SdSOt2//42Kjuv3sXb9cxMMZhoYpyp5y5+0uSXmr/+Vtm9rykcUlfyrttAAAQ1/TkuA4/c1aS9NBHbiu5NUgN47BQRVGLApjZdkmTkp7O+Nk+M5s1s9nz58/H/FoAAAAMQa/xVozDQpmiBTRmdoOkRyT9kru/tvLn7v6Au0+5+9TmzZtjfS0AAACGhHFYqKIoVc7MbINawcyD7n4kxjYB1BcVcgAgTZ1r9b0Pn9TFpSsa5xqOCsgd0JiZSfqUpOfd/TfzNwlAnTFTOQCkjXFYqJoYKWe7Jf2UpD1m9sX2fz8aYbsAaogKOQAAIKYYVc7+lySL0BYADUCFHAAAEFPUKmcAsBYq5AAAgJgIaGrq6Ny8dh98QjtmHtPug08wgy8qgwo5AAAgpihVzlAtDLpGlVEhBwAAxERAU0P9Bl3z0BgX5YfXhwo5AKqAazhQDwQ0NcSg6+GgJwwA0sU1vJoIMrEejKGpIQZdDwflh4G4GPuHYeIaXj2dIHN+YVGuq0Em1wKshYCmhhh0PRz0hAHx8CCDYeMaXj0EmeVK+aUSKWc1xKDr4dgyNqr5jBsfPWHpILWhOhj7h2HjGl49BJlxhdzjUk/BpIempqYnxzU5Mab37LhJx2f2JHEwpoaesLTRI1AtPMhg2LiGVw8p8/GE3uNS7x0joAHWaXpyXPffvUvXjbROo/GxUd1/9y6Cx0SkfvGuGx5kMGxcw6uHIHNtg6aFhd7jUn+pRMoZkAPlh9OV+sW7bvbv3bks3UHiQQbF4xpeLU1NmR80NSwkLSz0Hpd6CiYBDYBGSv3iXTdNfZABsFydg8yswEXSwEFKyFjD0Htc6i+VCGgANFLqF+86qvODDIBmWRm8vPdtm/XIiflVgcubNlzTNzWsextZAYqU3esSeo9L/aUSAQ2ARkr94g1guKiKiDypYQ8+dVa+Yr3FS0urgpmOTsDTvQ2TVm1Dyu51Wc89LuWXSgQ0ABor5Ys3gOFJvaQt8gs5BrJSw7ICkX5GzDK3sTKoWavXpSn3OKqcAQAA9EFVxPoqompYSHGZsdENmZXdljw7BHKJynwZ6KEBAADog6qI6ShqMsmQY6DXeJes3pWP3/n9klanhh06djpzG+Njo9p6YyvFrO69LiHooQEA1N6gb2GBLMyTlIYiJ5MMOQZ6zadzz60Tmb0rWZOhMydPGAIaAECthT7kACvxcFmuKkwmGXIM9Jq09RPTu1YFLr0w8WsYUs4AALUWMncDkIWqiOWpymSSocdAjAH5TRrUnxcBDQCg1hj/gBh4uIxr0LEuVZpMkmOgukg5AwDUGuMfgOEYNDUsJA20qLQwibSuOokS0JjZp83sZTN7Nsb2AAD5MRC+hfEPQPFCgpSiBuOvJ0DJGpCP9MTqofmMpNsjbQsAkBMD4a/iLSywfsOcpyVWrwsBSvNEGUPj7n9hZttjbAsAkB8D4Zcj9x1oiTVPy0ox5mmJMRgfzcQYGgCoIQbCA1ip6vO00OuC9RpaQGNm+8xs1sxmz58/P6yvBYBGYiA80Bx1maeFQAXrNbSyze7+gKQHJGlqasqH9b0A0ESh5UgBVMugqWFNnqcF6CDlDABqiDegQLqqUDFMIjUM6YhVtvmwpL+UtNPMXjSzn42xXQDA+vFwAVRLihXDeDGCFMSqcvaTMbYDAACQkrJTw4quGEZqGFJAyhkAAMA6VCE1jLQwgIAGAABgmZRSw0gLA4ZY5QwAAKAMdZ9MkrQwNB09NAAAoLaYTBKoPwIaAACQHCaTBNBByhkAAEgKk0kC6EYPDQAAqIQiel2YTBKoPwIaAABQmEGDlJCxLkwmCaAbKWcAAKAQIalh/XpdVq7LZJIAutFDAwAAgqQ0T4tEWhhQd/TQAADQcLHmackzIL/oXhcA9UUPDQAADcY8LQBSR0ADAEANMU8LgKYg5QwAgEQMmhrGPC0AmoQeGgAAEhCSGsY8LQCahIAGAIASpVgxjNQwAFVCyhkAAJGVnRrGPC0AmoSABgCAiHoFKbNfe1VPfvn8siCnqMkk9+/duawN0tq9LgQoAFJFQAMAwBpC5mnpFaQ8+NRZeftzJ8hZuV5Hr9SwQYMU5mlJR8ixBSAbAQ0AAH2EpIVJvVPDfMXnxUtLGjHTkq/8CRXDmiL02AKQjaIAAIBGKmqell7VwbIsuVMxrMFCjy0A2QhoAAC1MWiQElICOXSelqyqYdajvZ0KYVQMa6bQYwtANgIaAEAtVGWelqyyxvfcOtGzJ4Zel+YKPbYAZIsS0JjZ7WZ22sy+YmYzMbZZtEHf4gEAypXaPC3S6tSwT0zvoicGq6zn2AKwWu6iAGY2IumTkv6xpBclfcHMHnX3L+XddlEYhAcA5ar7PC1ZGKSPlahGB8RhnlFdJWgDZrdJ+ri7721/PiBJ7n5/r78zNTXls7Ozub43j90Hn9D8wqI+cvJzeus3r77p23jtiCYnxvTXf/u6JGn7d1+/7O9lLY+x7pdeek2S9I63fFfU7wvZbpG/X1HbKOrfLXQbWe2o8r9FjO+L8fuFtLnI3y/v/hv2vg79vir+fq98+4K++srrunLl6v3nmmtMb910vV7+1oVl7Z07u6ALl1eXNs66Xq+1bve/Rb82bLphY5TjMGv9XtsI2U8xzp2QtsU4tvKu+8q3L+jMK69r6Ypr47Uj2nbTqDbdsDHK71fUNoo6LmJso+rHchV+v6LWLXJfx7pu/c3mbfrw4d9WFZjZCXefWmu9GGWbxyW90PX5RUnvyWjQPkn7JGliYiLC165fr7d4nRvhGxez5wXIWh5j3TdfN5K5PO/3hWw3xveFrBtjG0X9u4VuI6sdVf63iPF9MX6/kDYX+fvl3X/D3teh3zfM3++Vb1/Q11+7IHfXN16/9J2HzpXrvvDq4rJAQpKuXHG98Oqibrx+w7LlWQFK9/LubW+7aTQzSNl2U6vXpfvfYtMNGyUp8yG51+8Xehxmrd9rGyH7Kca5E9K2GMdWnnVXBp8XLi/pq6+83rfNef/tY2yjqOMixjaqfixX4fcrat0i93Ws69ZN12/M/FmVxeih+YCkve7+c+3PPyXp3e7+L3r9nar00Kw0Pjaq4zN79BO/95eSVqcEZC2PsW4veb8vZLsxvi+0bcP8/Yps26BtiPV9ddnXoW0OWXeY+2/Y+3qYx8XRufmeqTBZ62ZN+nj/3btWpVntmHls1ZwsUqsS2JmDdyxbFnq9jjVJYVHHbFFtqEo7ilh3rWOgqoraH0WqwjFUpLKP5aKl2Oa1DLOH5kVJ27o+b5V0LsJ2CxMy2zIANFEnQLm4dEXS2mMN+w3I33rj8nEpIeNXQq/X05PjjD+oGUobA1hLjCpnX5B0s5ntMLPrJH1Q0qMRtluYTknN8bFRmag2A6A5js7Na+7sgp4+82rUySSLqhrG9RqUNgawltw9NO5+2cx+QdIxSSOSPu3uz+VuWcF4iwegLjpBysWlK9p98Ik1K4YN0usS+lZ8PVXDBk0N43pdnkGPrSKRVQFgLTFSzuTufyzpj2NsCwAwuJAgpV+vy8p1QwIUqf9DZ2cMTTeClOoLTTssSmgADKB5okysCQCIq4jUsCInkyQ1rH5C0w6LND05ruMze3Tm4B06PrOH4wrAMlF6aAAAa4uRGrZS0ZNJhrwVp9elXhiMD9RTFVJJYyOgAYAhCAlSqBiGKghNOwRQfVVJJY2NlDMAWKdeaWFZy6uQGkZaGEKEph0CqL4qpZLGRA8NAKxDr7dcs197VY+cmF+1fOUNpOPcwmKuXhcqhqEoocdWHdNYgLqpayopAQ0AdBn0oazXW67DT7+gJfdVy0fMVi2XwlPDqBiGYRr02KprGgtQN3VNJSXlDGiYQatn1cmgv3Ovh7Ks9Xu9zcoKWjrLSQ1DXdU1jQWom7qmktJDAzRIyMD0uhj2PC29emLG2+k6Wek79LogdXVNYwHqpq7zOhHQAA0SUj2r6vKmhmUFKaGD8bPSwn7s743rkRPzmeliBCmoq7qmsQB1VMd7ESlnQINU/S3qMFPDes3TkqXXYPystLBPTO8iXQyNU9c0FgBpoIcGaJAqv0UddmpYkfO01PHt13pR+aoZ6prGAiAN9NAADVLGW9RBe12Yp6V+QnrSkL7pyXEdn9mjMwfv0PGZPZw3AIaGHhqgQfq9Rc0amJ4l5I17SK9LaGoY87RUX0hPGgAA60VAAzRMngf20Lkmqp4ahmJVfcxWDKTUAUD5SDkDUEhamERqWNOFFFlIESl1AFAN9NAANTXom+Oi0sIkUsOaLrQnLTWk1KFO6G1EyghogArKe2OpQsUwidSwpqt75asmpNShGULTiYGqIaABKqbfjUVSJSeT7BegdNpTxwdarK3OQWqVy6ADIehtROoIaICK6XVj+fijz+nC5StJVQzr/B1uiKijuqfUoTnobUTqCGiAkq1ML8sKLiRpYfHSqmVUDAPKQw8k6oLeRqSOKmeojEErbVVdyO+RlV5mgd9HxTCgPEwmiTooY9JlxFWXZ6j1oocGlVCXAYkx5mlxSdb+f8fohhG9acM1+sYbq3tpqBgGAMiD3sa01eUZKo9cAY2ZfUDSxyW9XdK73X02RqPQPFUfkDho1bHQ36NXfrKr1VPSfWORRGoYAKAQ3DPSVfVnqGHI20PzrKS7Jf1ehLagwcoYkFjleVrGx0Z1fGZP5t/hDRoAAOigqEPOgMbdn5cks9Csf2C5YQ9ILGOelqwAigH5AAAgD4o6UBQAFRFrQOKgg+L6BSkrhc7TkvV7vPdtm3sGUAzIBwAA60VRhwF6aMzs85K+J+NH97n75wb9IjPbJ2mfJE1MTAzcQDRDvwGJMVLDVhr2PC39AigqIwEAgPWiqMMAAY27vy/GF7n7A5IekKSpqSlfY3U0UFY6VazUsK03Lg8+hj1Pyy8/9MXMdZuU3woAAIrR9JR0Us5QiqqmhhU1T0uvPNYm5bcCAAAUIW/Z5n8q6bclbZb0mJl90d33RmkZkjJoWlhn3SKqhsVIDStqnpbQXh6gykLOdwAAipa3ytlnJX02UluQqBiTSa6nathK/YKGw8+cXbX+MLtnyW9FXTCBGwCgakg5Q09FpIVJ1U4NK9L05LiOz+zRmYN3UAgAyQo93wEAKFreiTWRmCpPJlmF1DAA/TGBGwCgaghoGqSMySSzMJkkkC4mcANQd4wTTA8pZzWQUsUwqfqpYbhq0GMLzcEEbgDqrNfLX+5/1UYPTeJSrBjW+TsEMNXG4G9kocAFgDoLyVBBdRDQVNSg3Z1lVAzLQoBSP1zU0Qvne/FIeUkb+y9djBNMEylnQxKSuhPS3dnUimEoHhd1oBykvKSN/Zc2JsJOEwHNEIRe3ELGuoSceKFBCmWGm42LOlAOSmOnjf2XttBxgow1rQYCmhxSm6dFIkjB4Bj8DZSD3tG0sf/SFvLyl9646mAMzQp1n6cFGBTHFpqq7PEPlMZOG/svfYOOE2SsaXUQ0HRJYZ6Ww8+cXbU+A3RRFI4tNE0VqvuFFmVBtbD/moPeuOog5awL87TUE/mtAAZVhfEP3APSxv5rDsaaVgc9NF2Yp6V+qvC2FUA6qvLGlXtA2th/zUBvXHUQ0HRhnpb6Ib8VQAjGPwAYFGNNq4OApktIkMJBnIaqvG0FkAbeuAIIwQvraiCg6RIapHAQVx9vWwGE4GUVAKSHgGYFgpR64W3r8JRd6haIhfsAAKSFKmeoNarNDAeTiwEIQfVJADHRQ4Pa421r8Si+AGBQVJ8EEBs9NAXgzROahuILAAZVhbl+ANQLAU1kpN6giZhcDMCgeAECIDYCmsh484Qm2r93p0Y3jCxbRvEFAFl4AQIgNgKayHjzhCai+AKAQfECBEBsuYoCmNkhSf9E0kVJ/0/Sh919IUbDUsW8J2gqii8Uj9LYqAPm+gEQW94emscl3eLu75T0V5IO5G9S2njzBCDEoEVEGJ+HflIrRjM9Oa7jM3t05uAdOj6zh2AGhUvtHEGYXAGNu/+Zu19uf3xK0tb8TUobqTcABhUSpDA+D70Q7AL9cY7UX8x5aH5G0kMRt5csUm8ADCJk/h7G56EX5oEC+uMcqb81Axoz+7yk78n40X3u/rn2OvdJuizpwT7b2SdpnyRNTEysq7EAUCchQQrj89ALwS7QH+dI/a2Zcubu73P3WzL+6wQzH5L0fkn3uLv32c4D7j7l7lObN2+O9xsAQKJCytcyPg+9UAYZ6I9zpP5yjaExs9sl/aqkO939jThNAoBmCAlSGJ+HXgh2gf44R+ov7xia35G0UdLjZiZJT7n7R3O3CgAaILR8LePzkIUyyEB/nCP1lyugcffvi9UQAGgighTEwHEE9Mc5Um9556EBAAAAgNIQ0AArMPkWAABxcE/FMBDQAF2YfCt93DwBoBq4p2JYCGiALszGnjZungBQHdxTMSwENEAXJt9KGzdPAKgO7qkYFgIaoAuTb6WNm+f6kaoHIDbuqRgWAhqgC5NvpY2b5/qQqgegCNxTMSwENEAXZmNPGzfP9SFVD0ARuKdiWHJNrIlm6qSmXFy6ot0Hn6jdbLtMvpUuZoNeH1L1ABSFeyqGgYAGQXqlpkjigoVK4OYZbsvYqOYzghdS9QAAKSDlDEFITQHqh1Q9AEDKCGgQJDQ1hcpJQPWlmufO9QUAIJFyhkAhqSmkpwHpSC1Vj+sLAKCDHhoECUlNIT0NQFG4vgAAOuihQZCQKlJUTgJQFK4vAIAOApoB1b1UcYhBU1OonIQ64RpQLVxfAAAdpJwNoCqzaKc2AJbKSaiLqlwDcBXXFwBABwHNAKqQq53iA1WqlZOAlapwDcByXF8AAB2knA2gCrna/R6oqnwDT61yEpClCtcArMb1BQAg0UMzkF452cPM1eaBCihPFa4BAAAgGwHNAKqQq80DFVCeKlwDAABANgKaAVQhV5sHKqA8VbgGAACAbLnG0JjZv5d0l6Qrkl6W9M/d/VyMhlVN2bnaIfO/AIiv7GsAAADIlrcowCF3/9eSZGa/KOnfSPpo7lYhEw9UAAAAwHK5Us7c/bWuj9dL8nzNAQAAAIDB5S7bbGa/LumnJX1T0ntztwgAAAAABrRmD42Zfd7Mns347y5Jcvf73H2bpAcl/UKf7ewzs1kzmz1//ny83wAAAABAY5l7nCwxM/teSY+5+y1rrTs1NeWzs7NRvhdxHJ2b170Pn9TFpSsap+AAAAAASmZmJ9x9aq31co2hMbObuz7eKenLebaHchydm9eBI6d0cemKJGl+YVEHjpzS0bn5klsGAAAA9Jd3HpqD7fSzk5J+WNLHIrQJQ3bo2GktXlpatmzx0pIOHTtdUosAAACAweQqCuDuPxarISjPuYXFoOUAAABAVeTtoUENbBkbDVoOAAAAVAUBDbR/706NbhhZtmx0w4j2791ZUosAAACAweSehwbp61QzO3TstM4tLGoLVc4AAACQCAIaSGoFNQQwAAAASA0pZwAAAACSRUADAAAAIFkENAAAAACSRUADAAAAIFkENAAAAACSZe4+/C81Oy/pa0P/4mybJL1SdiOwbuy/tLH/0sb+Sxv7L23sv3Sx7wb3ve6+ea2VSgloqsTMZt19qux2YH3Yf2lj/6WN/Zc29l/a2H/pYt/FR8oZAAAAgGQR0AAAAABIFgGN9EDZDUAu7L+0sf/Sxv5LG/svbey/dLHvImv8GBoAAAAA6aKHBgAAAECyGh3QmNntZnbazL5iZjNltwf9mdk2M3vSzJ43s+fM7GPt5TeZ2eNm9n/b/7+x7LYim5mNmNmcmf1R+/MOM3u6ve8eMrPrym4jspnZmJk9bGZfbp+Dt3HupcPMfrl93XzWzA6b2Zs4/6rLzD5tZi+b2bNdyzLPN2v5rfazzEkz+8HyWg6p5/471L5+njSzz5rZWNfPDrT332kz21tOq9PW2IDGzEYkfVLSj0h6h6SfNLN3lNsqrOGypF9x97dLulXSz7f32YykP3f3myX9efszquljkp7v+vwbkv5je999Q9LPltIqDOI/S/pTd3+bpB9Qaz9y7iXAzMYl/aKkKXe/RdKIpA+K86/KPiPp9hXLep1vPyLp5vZ/+yT97pDaiN4+o9X773FJt7j7OyX9laQDktR+jvmgpO9v/53/0n5GRYDGBjSS3i3pK+7+VXe/KOkPJd1VcpvQh7u/5O7/u/3nb6n1QDWu1n77g/ZqfyBpupwWoh8z2yrpDkm/3/5skvZIeri9CvuuoszsuyT9Q0mfkiR3v+juC+LcS8m1kkbN7FpJb5b0kjj/Ksvd/0LSqysW9zrf7pL037zlKUljZvaW4bQUWbL2n7v/mbtfbn98StLW9p/vkvSH7n7B3c9I+opaz6gI0OSAZlzSC12fX2wvQwLMbLukSUlPS/q77v6S1Ap6JP2d8lqGPv6TpHslXWl//m5JC10XeM7B6nqrpPOS/ms7ZfD3zex6ce4lwd3nJf0HSWfVCmS+KemEOP9S0+t843kmPT8j6U/af2b/RdDkgMYyllHyLQFmdoOkRyT9kru/VnZ7sDYze7+kl939RPfijFU5B6vpWkk/KOl33X1S0usivSwZ7bEWd0naIWmLpOvVSlNaifMvTVxLE2Jm96mVQv9gZ1HGauy/QE0OaF6UtK3r81ZJ50pqCwZkZhvUCmYedPcj7cVf73Svt///clntQ0+7Jd1pZn+tVnrnHrV6bMbaKTAS52CVvSjpRXd/uv35YbUCHM69NLxP0hl3P+/ulyQdkfT3xfmXml7nG88ziTCzD0l6v6R7/Oq8Key/CJoc0HxB0s3tKi/XqTUg69GS24Q+2mMuPiXpeXf/za4fPSrpQ+0/f0jS54bdNvTn7gfcfau7b1frXHvC3e+R9KSkH2+vxr6rKHf/G0kvmNnO9qJ/JOlL4txLxVlJt5rZm9vX0c7+4/xLS6/z7VFJP92udnarpG92UtNQHWZ2u6RflXSnu7/R9aNHJX3QzDaa2Q61ijs8U0YbU9boiTXN7EfVeks8IunT7v7rJTcJfZjZP5D0PyWd0tVxGP9KrXE0/13ShFo37g+4+8rBlKgIM/shSf/S3d8vCqtTAAAAyklEQVRvZm9Vq8fmJklzkv6Zu18os33IZmbvUqugw3WSvirpw2q9FOPcS4CZ/Zqkn1Ar1WVO0s+plafP+VdBZnZY0g9J2iTp65L+raSjyjjf2kHq76hVIesNSR9299ky2o2WHvvvgKSNkv62vdpT7v7R9vr3qTWu5rJa6fR/snKb6K/RAQ0AAACAtDU55QwAAABA4ghoAAAAACSLgAYAAABAsghoAAAAACSLgAYAAABAsghoAAAAACSLgAYAAABAsghoAAAAACTr/wMpTsipNsjI1gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2233938d160>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"W = dft_matrix(N);\n",
"\n",
"# DFT\n",
"X = np.dot(W, x)\n",
"plt.stem(abs(X));\n",
"plt.show();\n",
"\n",
"plt.stem(np.angle(X));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly we have a problem with the phase, although the magnitude looks nice. This is inherent to the fact that the phase is computed by taking the arctangent of a ratio. When the computed DFT values are close to zero, the denominator of the ratio will be also close to zero and any numerical error in its value will lead to large errors in the phase. As we will see in the next section, this problem can be alleviated by using smarter algorithms than the direct naive method.\n",
"\n",
"Let's still verify the inverse DFT:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAEDCAYAAADjkkFLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X2MHdd53/Hfw+WKWkl1VrK2ibkUSzl2qchWrY0JWwqLwmYMrJy4FitXsF3XNVwbUoCmdYSEqtgUsFM0EAsWcVLHCSLYrh1AYGhLDCPYTVjXUmFXiKSQWVmULLFRpZjWUrFWltd6W/Bl+fSPvVdcLmfu3rn3zMw5M98PIIg7vJx79s7MOfOcec5zzd0FAAAAAClaU3cDAAAAAGBQBDQAAAAAkkVAAwAAACBZBDQAAAAAkkVAAwAAACBZBDQAAAAAklVbQGNmXzKz58zs0UD7+wszmzezr+f8/efM7OUQ7wUAAAAgDnU+ofmypOsC7m+3pI9m/YWZbZE0HvC9AAAAAESgtoDG3b8t6YXl28zsZztPWg6Z2XfM7IoC+/uWpJdWbjezES0FO7cO22YAAAAAcVlbdwNWuEPSr7j735jZOyX9gaRtQ+7zVyXd4+7PmtnQDQQAAAAQj2gCGjO7SNIvSPrassBjXefvbpD0nzL+2ay7T/fY53pJN0p6V9DGAgAAAIhCNAGNltLf5t396pV/4e77JO0bYJ9Tkt4k6clOkHSBmT3p7m8aqqUAAAAAohBN2WZ3f1HS02Z2oyTZkrcNuc9vuPvPuPsmd98k6VWCGQAAAKA56izbvEfSX0rabGbPmNknJH1E0ifM7LuSHpN0fYH9fUfS1yT9Ymd/ualoAAAAAJrB3L3uNgAAAADAQKJJOQMAAACAomopCnDppZf6pk2b6nhrAAAAAAk4dOjQ8+4+sdrragloNm3apIMHD9bx1gAAAAASYGbf7+d1pJwBAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASFYtVc4AAAAAxGP/zKx2HziiY/MLWj8+ph3Tm7V9arLuZvWFgAYAAABosf0zs9q577AWTi5KkmbnF7Rz32FJSiKoIeUMAAAAaLHdB468Fsx0LZxc1O4DR2pqUTEENAAAAECLHZtfKLQ9NkMHNGZ2vpk9ZGbfNbPHzOy3QjQMAAAAQPnWj48V2h6bEE9ojkva5u5vk3S1pOvM7JoA+wUAAABQsh3TmzU2OnLWtrHREe2Y3lxTi4oZuiiAu7uklzs/jnb+82H3CwAAAKB83YX/t971iE4sntZkG6ucmdmIpEOS3iTp8+7+YMZrbpJ0kyRt3LgxxNsCAAAACGD71KT2PHRUkrT35mtrbk0xQYoCuPuiu18taYOkd5jZWzNec4e7b3H3LRMTEyHeFgAAAEDLBa1y5u7zkv63pOtC7hcAAAAAsoSocjZhZuOdP49Jeo+kJ4bdLwAAAACsJsQamjdI+kpnHc0aSV91968H2C8AAAAA9BSiytkjkqYCtAUAAAAACgm6hgYAAAAAqkRAAwAAACBZBDQAAAAAkkVAAwAAACBZIaqcAQAAAIjI/plZ7T5wRMfmF7R+fEw7pjdr+9Rk3c0qBQENAAAA0CD7Z2a1c99hLZxclCTNzi9o577DktTIoIaUMwAAAKBBdh848low07VwclG7DxypqUXlIqABAAAAGuTY/EKh7akjoAEAAAAaZP34WKHtqSOgAQAAABpkx/RmjY2OnLVtbHREO6Y319SiclEUAAAAAGiQ7sL/W+96RCcWT2uSKmcAAAAAUrJ9alJ7HjoqSdp787U1t6ZcpJwBAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkEdAAAAAASBYBDQAAAIBkra27AQAAAECT7J+Z1e4DR3RsfkHrx8e0Y3qztk9N1t2sxhr6CY2ZXWZm95nZ42b2mJl9KkTDAAAAgNTsn5nVzn2HNTu/IJc0O7+gnfsOa//MbN1Na6wQT2hOSfp1d/9rM/t7kg6Z2Tfd/XsB9g0AAAAkY/eBI1o4uXjWtoWTi9p94MjQT2l48pNt6IDG3Z+V9Gznzy+Z2eOSJiUR0AAAAKBVjs0vFNrer+6Tn26w1H3yI6n1QU3QogBmtknSlKQHM/7uJjM7aGYH5+bmQr4tAAAAEIX142OFtver15OftgsW0JjZRZLulvRr7v7iyr939zvcfYu7b5mYmAj1tgAAAEA0dkxv1tjoyFnbxkZHtGN681D7LevJTxMECWjMbFRLwcyd7r4vxD4BAACA1GyfmtTtN1yl80aWbrMnx8d0+w1XDZ0WVtaTnyYIUeXMJH1R0uPu/jvDNwkAAABI1/apSU1tHNc7L79E99+2Lcgal7Ke/DRBiCc0WyV9VNI2M3u4898vBdgvAAAAAJX35KcJQlQ5+z+SLEBbAAAAgFYpUop5+9Sk9jx0VJK09+Zrq2xm1EJ8Dw0AAACAgijFHEbQss0AAAAA+kMp5jAIaAAAAIAaUIo5DAIaAAAAoAaUYg6DgAYAAACoAaWYw6AoAAAAAFCD7sL/W+96RCcWT2tylSpnyEZAAwAAANSEUszDI+UMAAAAQLIIaAAAAAAki4AGAAAAQLIIaAAAAAAki4AGAAAAQLIIaAAAAAAki4AGAAAAQLIIaAAAAAAki4AGAAAAQLLW1t0AAADQLvtnZrX7wBEdm1/Q+vEx7ZjerO1Tk3U3C0CiCGgAAEBl9s/Maue+w1o4uShJmp1f0M59hyWJoAbAQEg5AwAAldl94MhrwUzXwslF7T5wpKYWAUgdAQ0AAKjMsfmFQtsBYDUENAAAoDLrx8cKbQeA1RDQAACAyuyY3qyx0ZGzto2NjmjH9OaaWgQgdRQFAAAAleku/L/1rkd0YvG0JqlyBmBIBDQAAKBS26cmteeho5KkvTdfW3NrAKSOlDMAAAAAySKgAQAAAJAsAhoAAAAAySKgAQAAAJAsAhoAAAAAySKgAQAAAJAsAhoAAAAAySKgAQAAAJAsvlgTAACg5fbPzGr3gSM6Nr+g9eNj2jG9WdunJutuFtAXAhoAAALj5hBlKePc2j8zq537Dmvh5KIkaXZ+QTv3HZYkzlskIUhAY2ZfkvQ+Sc+5+1tD7BMAgBRxc4iylHVu7T5w5LV9di2cXNTuA0cy95tiwF6kzSn+fm0X6gnNlyX9vqQ/DrQ/AOiJAQcxyDoPi94cAv0q69w6Nr/Q9/YUA/YibU7x90OgogDu/m1JL4TYFwCspjvgzM4vyHVmwNk/M1t309AieefhbIGbQ6CIIoFHEevHx/re3iuoilWRNqf4+6HCKmdmdpOZHTSzg3Nzc1W9LYAGYsBBDPLOwxGzzNfn3TQC/SoSeBSxY3qzxkZHzto2NjqiHdObz3ltWUFVmYq0OcXfDxUGNO5+h7tvcfctExMTVb0tgAZiwEEM8s63Rfe+bw6BIooEHkVsn5rU7TdcpfNGlm4LJ8fHdPsNV2WmWJUVVJWpSJtT/P3A99AASBADDmKQd751bwb7uTkEiigSeAyy76mN43rn5Zfo/tu25e6zrKCqTEXanOLvB8o2A0jQjunNZy3alBhw0FsZRSR6nYfbpya156GjkqS9N1871PsAy9V9bnWvm1vvekQnFk9rMoGiLEXaXObvRzGb8oQq27xH0rskXWpmz0j6tLt/McS+AWClFAdU1CevatHB77+g+56YG/jmgvMQbVV3UDWIIm0u4/ejelq5ggQ07v7hEPsBgH6lOKCiHnmL9+984Ki88/NqNxd5M6uch4gZTwTiEVM59yaeF6yhAQA0Wt7ifV/xc16lPMqEI0Wct3GJpZhNU88LAhoAQKMVKRaRdXNBmXCkiPM2LrEUs2nqeUFAAwBotKyqRdnfFJN9cxHLzGrs9s/Mauuue3X5bd/Q1l33Jj/jmzrO27jEUj2tqecFAQ0AoNGySt1+5JqNfd9cxDKzGrOmprEsl1rAxnkblzJLbhfR1POCgAYA0Hgrv2PjP2+/qu+bi1hmVmPW1DSWrhQDNs7b+PT7XT9laup5wffQAABaqd8KZZRnXl1T01i6YqpQ1S/O23bpt3JZU88LAhoAAFZBeebe1o+PaTYjeEk9jaUr1YCN87Yden3HTZYmnheknAEAsExqayVi0NQ0lq6mrjtAMzQ95bMfBDQAAHSkuFYiBrEseC5L0wM2pC3VJ4ghkXIGAEBHr5nODRczG99LE9NYupq67gDN0PSUz37whAYA0BjDposx04k8MVSoArLwBJGABkgKuf1AvhDpYqyVAJCapqd89oOABkhEW3P7CeLQrxALY5npBJCitj9BZA0NULN+a8en+D0Iw+pVirKpvzMGFyJdrNdaie76ECBl/Y45QEoIaIAaFblhb2NufxuDOAwu1MLYJi9ujx032+VikghNRcoZUKMiKTJtzO1vYxCHwZEulra2ptVWie8rQVMR0AA1KnLD3sabtTYGcSmIdV0TC2PTxs12+ZgkQlORcgbUqEiKTBu/B2HH9Oaz0iOk5gdxsYs9ZYV0sfj0m0bGzXb5yv6+ElIGURee0AA1KvrUpW1VTJhxjw+z6CiiSBoZT2TLV+aTflIGUScCGqBG3LCvrm1BXOyYRUcRRQLgNqbVVq3MMSfUZEesKa2IGylnQM1IkUFKyk5ZQRpCpJFtuPjsc6aNabV1KGvMCTHZEXtKK+JFQAMA6FtM65rI169Hr5vOlYoGwEzwpCvEZAel+jEoUs4AAH2LJU2SfP36kEaGLCGONSmtGBRPaAAAheTNolf5xKTXTfXKVCaEFSqNrHsOoRlCpAyS0opBEdAAAIZWde57kZtqhEUaGfIMe6xjSmlFWkg5AwAMrepyzpT4rQ9pZChLLCmtSA8BDQBgaFXnvnNTXR9uOtul6jLKRUr1U+IZXaScobWokASEU3XuO2sz6kUaWTsUqWhXNUo8YzkCGrQSHSEQVl7u+7uvmNDWXfeWMnHATTVi1ZQJs5iLb1DiGcuRcoZWqjrfH2i6rDSkD7x9UncfmqW0MlqlSSXFYy6jHHPbUD0CGrQSHSEQ3src9/uemGPiIHGsUSiuSRNmMRffiLltqB4BDVqJjhBY3bA3s0wcpK1JTxqq1KTzPubiGzG3DdUjoEEr0RECvYW4me01ccDMf/ya9KShSk2aMIu5ol3MbUP1KApQkaYsEExRr89+mG80BposxILbXoUCKMoRvyY9aahS074cMubiGzG3DdUioKkAFbXqs9pnT0fYPKlNHsTa3hA3s3kTB1QnSkPVpbhjMsx1WXTCLNY+AP0p8/hxbvQvSEBjZtdJ+j1JI5K+4O67Quy3KRi868Nn3y6pTR7E3N5QN7NZEwe37H0487XM/MclpicNVd7YhfjulX4nzGLuA7C6Mo9fzN8BFKOh19CY2Yikz0t6r6QrJX3YzK4cdr9NwmP7+vDZt0tqOf8xt7fMdWZNWmPQZLGsUai6OEGV12XMfQBWV+bx49woxtx9uB2YXSvpM+4+3fl5pyS5++15/2bLli1+8ODBod43hD/75G9o/NjT2vT6C8/a/rc/ekWS+trez2tnjs7r+KmzT0pJWrd2RFMbx4O/3yCvff7l43r6+Ve0eNq1bu2ILrtkTJdetK6096vq91vts//esy9Kkq58w+sqb9tyWe0I8doYzq0Q+8g7Tiu3P/DUj5TnZ37q/FLaNsxn0au917zx9bV/9s+/fFz/b+4VuZ/dL0j9H5Osbc+/fFxPPf+KTp8+M/6sWWN646UX6uXjpwr9fv28Xx2vzdse83VW5DPO217Ga3v14xdfODrU75c19j353MvnvFfX68ZGM9tcpN8atg/I23dZ10iIfcQyJhYZ80Mcv2E+t177zjsPQx3rv5u4TB/f87nc96+SmR1y9y2rvS5EytmkpB8s+/kZSe/MaNBNkm6SpI0bNwZ42+G98MpxnXfi3E7y1Yxtedv7ee1ll4xlDt6XXTJWyvsVfe3Km4vjpxb11POvlNq2EPsI8dlfcN7IOf++yPs9//Jx/fDF43J3/fiVk6/d8BX9LFa2I2+/eW2O8bMPuY+847Ry+7q1I7k3PlV+Fr2O3/LX92pviLYNu49LL1qXG2D0e0yytnU/i6xJlOdeOn7Ovy96PQzTtlCvzdse83XWb//Ua3sZr826RrrbXz1xbqJJv/1y3ti3dmSNTi2ePme/69aOFBoz+nntIH1A3r6LvDbEuVzWeV90H8O+Nu/1IY7fMJ9br32X3W9dcuG6zL+LWYgnNDdKmnb3T3Z+/qikd7j7v837N7E8ofngH/2lpHPzW7O275+ZzVzg1+8+euX/FmlHGa/duuvezFz5yfExbbh4rJS2Vfn7DZJ73c/7rcxvlZZScm6/4arc3Om8Ni/Xa7957Y71sw+5j36EOibD/B6rHb/lry/y2kHaFuL3q1rRzzM1MV9nMRz/LEXGqCJ9QN5+x8dGdfzU6SB98GqvHaQP6Ge/sQnR55T5WZR1/IbRtL5vUFU+oXlG0mXLft4g6ViA/Uaje1Kd6MzYDLpAMNYTsNc6k+5gkbKyPvte+a3DfG6hChnsn5nVzNF5nVg8ra277m1FdZTu75cVwHZvZspW5Pj1am8bZZ2zFPZAr+IEK6/rIv1y3tj3k4WT+uwHr67kugzVB7Sxv49BmX0440MxIQKav5L0ZjO7XNKspA9J+hcB9huNsm5ciyqrw2pzac5h5A2Gs/MLmnvp+MDHKUQhg15BeNM7w7onD4oev7rb21X3DVHeObuy7+2isEdYdR//XopMVBSZoOs19lV5XQ77Xm3u72NQ5rkSy/iQgqGrnLn7KUm/KumApMclfdXdHxt2vzGpo1JWd3B58OkXtHXXvfqP+w9ndlghqryUWc0ohJWfxf6Z2cxtVcsL+Ewa6jiFqABFdZT69Dp+MZy3WfJuiKpsX945O2KW+XomXMKJ4fivZvvUpO6/bZue3vXLuv+2bbk3eUX6zzrGvjL6APp7IEBAI0nu/j/c/R+6+8+6+2+H2GddsjqbqkuMZg0udz5wtLQOq1uac3J8TKb6SnNmyfosdnztu9px13drH3yzBkOTtHJVWtHjFGKQpVx1ffKO37uvmIj2pjGGG6K8c3PRPeoJlyaI4fiHUqT/rHrsKytwpL8HAn2xZlPkdTYfePuk7j4021f+bghZg0te6YZQHVasjzWzPouTp8/9NOrIqc9Kg8hKX5AG+3b1YfJmSSOsT97xi3ktSAw3RHnn7OSyz4888nLEcPxDKbqOrsqxr6w+gP4eIKA5S15nc98Tc7r9hqsqW2hcZBBpeodV5LOoIgVwZX75ysEwr2rOIN+uPswAF9M3fLdR1vG7Ze/Dma8tujaqqevoep2zsU64NEUMx38Q/fbLsSgrcKS/BwKlnDVFr86m3/zdEHqtzViu22EVycmNNYc/T5EBtcoUwLw0gVjWI8WcRthWw6aulrnOIYbzlnO2PjEc/6JSWPezUlnp61w7AE9ozhLLLFXebMsH3j6p+56YO+spkaS+S0qnWAkl67MYXWOSSScXz6SelTn4plqGN9ZZyrYadha1zJS1WM5bztl6xHL8i4g5hTNPmU9SYr52Yq6gh+YgoFkmlse2RQaXrbvu7bukdIoDQN5nkbWtrN8h1TK8RTDglG/Ym8ay1zmkeN4inNSOf4rrflIMHIeV4kQq0kRAs0xMnU2/g0uRmvspDgBS/mdR1XGJ5cldWRhwqjPMTWPTz0OgiFSvh9QCx2GlOJGKNLGGZoUQa2WqXKdSJCe36vLTTZFifnkRTSrZ2iQr+5F3XzHR6PMQKKLp/XJTpDqRivS0NqApK+ioeqFikU6dAaA/K88NSUEWXMZakIEBJz5Z/cjdh2b1gbdPsvAXEAvhU8FEKqrSypSzMlNsqn68WqTmfkwpdbHKOzduv+Eq3X/btuD7lYY/54Zd/9IrdSPE2hrW5xTXq4T8MOch0CRF0rfoh+oRy9pkNF8rA5oyg446ZruzOvWq6/M3ZbAo69woa78hAqW8ASfvm+1DtQ/5yuxHmnKtAv1inWB9mEhFVVoZ0JR5sxDDQsWqO+8mDRZlnRtl7TdEoDTIN9uvLDgxSPv63UcbldWPNOlaBfrFwvR6hZhIZSIGq2nlGpoyczpjWKdS9SLvJi0qL+vcKGu/oQKlrGIYRfedtUaoSetzqlwDVVY/0qRrFehXk/qhNkrxS1RRvVYGNGUGHTEsVKy6827SYFHWuVHWfssMzovsO2/AGb9gtLT2lSUrcKl6QC2rH2nStQr0i4XpaWMiBv1oZcpZqJzOqtep9KvqtLcY0uxCKSvft6z9lrngste+VxacyBtw1q1do7HRkb72IdWfVpAXuJw/uqby1Lky+pEmXasYXN3XWdVYmJ42JmLQj1YGNNLwNwsx56JX3Xk3bbAoKyAtY79lLrgsUkEvb2D5ycJJffaDV/e1jxiuqbzAbOW2rqwvsI1Z065VFBfDdVY1FqanjYkY9KO1Ac2wQi0yLGOmrOrOm8GiXmU+Eex3370GnH73EcPC3aIzfrEMqP32I1yriOE6q0OZ/WTbnnhVjYkY9IOAZkAhHoGWOVNWddpb3vvR0bdDiAEnhrSCvMBsfGxUx0+d7jt1rkpF+5G6U2JRrxiusyZp4xOvqjERg360sihACCEWGTZ9oRuVSeJURrWuEIvYY1i4m1e84TPvf0vtxT7yNL0fQVgxXGdNwvVXjaxKnHWostoliuEJzYCaMiNdpramNsQs5qeCMaQVrDYTGON52/R+BGHFcJ01Cddfe/A0Lm4ENAMK8Qi06Qvd6OjjE3OQGUtaQWopWb36EVI+sVIs11lTNH0cxxkxj58goBlKE2aky0RHH5/Yg8zUgokY5PUj775igtnERFQdeHKdhdP0cRxnxD5+th1raGoUw5dw9jJsrmiZX2CKwZA/3zx5/ch9T8yR258A1hqmLfZxHOEwfsaNJzQ1i3WmrFeuqCRKxCaK2cSwYknpyupHbtn7cOZrmU2MS680lpS+46jNYh3HERbjZ9wIaJApb5D9zD2P6fip05SITRRBZjixLxAl5TMNvdJYCGiAeDB+xo2ABpnyBtn5hZPnbGNRXFoIMsOIfYEos4lpyAs8f2psNIqnfwDOYPyMF2tokKnoLC5pLGib2BeIktufhqy1hqNrTK+cODX0uhq+MwNAW/CEBpnyZnfPH12jH7967lMa0ljQNimkdDGbGL+sNJZXT5w6p58t+vQv9pRIrC6WNXpACnhCg0x5s7uf/qdvCVK5jJlDpI4qfghl5begz2dMGknFnv7xDfZpo/odUAxPaJCr1+zuMIviVqugNgxmtFAVFoiiLCGe/sWeEoneYl+jB8SGgAaFDZvGUlaZUlIsUDVSulBEvxMuIQo6pJASiXwEpEAxpJyhcmV11KRYAIhVkRSiEAUdSIlMG1/iCBTDExpUrqyZQ2a0AMSqaArRsE//SIlMG2XXgWIIaFC5Xh31noeODrxfUiwAxKqOCRdSItNFQAoUQ0CDyvXqqIcJaJjRAhArJlxQVNMDUor4ICQCGtSijI6aGa1qMAgBxTHhApxBER+ENlRAY2Y3SvqMpJ+T9A53PxiiUcCgmj6jVTcGIWAwZU+4MNGAlFCWGqEN+4TmUUk3SPqjAG0BEDkGIWBwZU24MNGA1FDEB6ENVbbZ3R93d2riIojuDOODT7+grbvu5RuRI8QgBMSHkvVIDWWpEVpl30NjZjeZ2UEzOzg3N1fV2yIRRb6jAfVhEALiw0QDUsP3JCG0VQMaM/tfZvZoxn/XF3kjd7/D3be4+5aJiYnBW4xGYoYxDQxCQHyYaEBqQnx5LLDcqmto3P09VTQE7cYMYxqoJAfEhwpqSBFFfBASZZsRBb6jIR0MQkBcmGgABkeFwGYYtmzzP5P0OUkTkr5hZg+7+3SQlqFVmGEEgMEx0QAUR4XA5hi2ytmfuvsGd1/n7j9NMINBkU8LAKujGiQQDut3m4OUM0SDGUYAyMdsMhAW63ebo7KyzQAAYHDMJgNhUSGwOQhoADQeaTpoAmaTgbD4KoLmIKAB0Gh8aSuagtlkYHVFJrBYv9scrKEB0Gi90nQYtJASqkECvQ2yzoz1u83AExoAjUaaDpqC2WSgN9aZtRdPaAA0Wq8vbeUL1ZAaZpOBfExgtRdPaAA0Wt6iz3dfMcHaGgBoENaZtRcBDYBGy0vTue+JOVITAKBBqFrWXqScAWi8rDSdW/Y+nPlaUhMAIE3dfn73gSM6Nr+g9eNjpBK3BAENgMrFsHal19oaAECaWGfWTqScAahULN8LQ2oCAADNQEADoFKxlNWkBC4AAM1AyhmASsVUVpPUBAAA0scTGgCVoqwmAAAIiYAGQKVYuwIAAEIi5QxApSirCQAAQiKgAVA51q4AAIBQSDkDAAAAkCwCGgAAAADJIqABAAAAkCwCGgAAAADJIqABAAAAkCxz9+rf1GxO0vcrf+Nsl0p6vu5GYGAcv7Rx/NLG8Usbxy9tHL90cez69w/cfWK1F9US0MTEzA66+5a624HBcPzSxvFLG8cvbRy/tHH80sWxC4+UMwAAAADJIqABAAAAkCwCGumOuhuAoXD80sbxSxvHL20cv7Rx/NLFsQus9WtoAAAAAKSLJzQAAAAAkkVAAwAAACBZrQ5ozOw6MztiZk+a2W11twe9mdllZnafmT1uZo+Z2ac62y8xs2+a2d90/n9x3W1FNjMbMbMZM/t65+fLzezBzrHba2bn1d1GZDOzcTO7y8ye6FyD13LtpcPMbun0m4+a2R4zO5/rL15m9iUze87MHl22LfN6syX/rXMv84iZ/Xx9LYeUe/x2d/rPR8zsT81sfNnf7ewcvyNmNl1Pq9PW2oDGzEYkfV7SeyVdKenDZnZlva3CKk5J+nV3/zlJ10j6N51jdpukb7n7myV9q/Mz4vQpSY8v+/m/SPps59j9WNInamkV+vF7kv7C3a+Q9DYtHUeqd87OAAADyUlEQVSuvQSY2aSkfydpi7u/VdKIpA+J6y9mX5Z03YptedfbeyW9ufPfTZL+sKI2It+Xde7x+6akt7r7P5L0fyXtlKTOfcyHJL2l82/+oHOPigJaG9BIeoekJ939KXc/IelPJF1fc5vQg7s/6+5/3fnzS1q6oZrU0nH7SudlX5G0vZ4Wohcz2yDplyV9ofOzSdom6a7OSzh2kTKz10n6J5K+KEnufsLd58W1l5K1ksbMbK2kCyQ9K66/aLn7tyW9sGJz3vV2vaQ/9iUPSBo3szdU01JkyTp+7v4/3f1U58cHJG3o/Pl6SX/i7sfd/WlJT2rpHhUFtDmgmZT0g2U/P9PZhgSY2SZJU5IelPTT7v6stBT0SPr79bUMPfyupFslne78/HpJ88s6eK7BeL1R0pyk/95JGfyCmV0orr0kuPuspP8q6aiWApmfSDokrr/U5F1v3M+k519L+vPOnzl+AbQ5oLGMbdSwToCZXSTpbkm/5u4v1t0erM7M3ifpOXc/tHxzxku5BuO0VtLPS/pDd5+S9IpIL0tGZ63F9ZIul7Re0oVaSlNaiesvTfSlCTGz39RSCv2d3U0ZL+P4FdTmgOYZSZct+3mDpGM1tQV9MrNRLQUzd7r7vs7mH3Yfr3f+/1xd7UOurZLeb2Z/q6X0zm1aemIz3kmBkbgGY/aMpGfc/cHOz3dpKcDh2kvDeyQ97e5z7n5S0j5JvyCuv9TkXW/czyTCzD4m6X2SPuJnvgiS4xdAmwOav5L05k6Vl/O0tCDrnprbhB46ay6+KOlxd/+dZX91j6SPdf78MUl/VnXb0Ju773T3De6+SUvX2r3u/hFJ90n6552Xcewi5e5/J+kHZra5s+kXJX1PXHupOCrpGjO7oNOPdo8f119a8q63eyT9q061s2sk/aSbmoZ4mNl1kv69pPe7+6vL/uoeSR8ys3VmdrmWijs8VEcbU2ZnAsT2MbNf0tIs8YikL7n7b9fcJPRgZv9Y0nckHdaZdRj/QUvraL4qaaOWBu4b3X3lYkpEwszeJek33P19ZvZGLT2xuUTSjKR/6e7H62wfspnZ1Voq6HCepKckfVxLk2Jcewkws9+S9EEtpbrMSPqklvL0uf4iZGZ7JL1L0qWSfijp05L2K+N66wSpv6+lClmvSvq4ux+so91YknP8dkpaJ+lHnZc94O6/0nn9b2ppXc0pLaXT//nKfaK3Vgc0AAAAANLW5pQzAAAAAIkjoAEAAACQLAIaAAAAAMkioAEAAACQLAIaAAAAAMkioAEAAACQLAIaAAAAAMn6/8JRmKN4xaQiAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x22339e8b0f0>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAEDCAYAAAD0sNGCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+0XlV95/HP1xDwirWXkKjkkpioTCw21WvvAimzZiHaCaCLpKlMoY5FR1faWTq2nRYm1FnquNqVdDGrra1Wy6ijruVCLELMFDRVcRaWETQY5YeYMYJCblACGBDIEJJ854/nXHJz85znec49v/be5/1aKyv3ee7Jc3bOec45+7v3d+9t7i4AAAAASNlz2i4AAAAAANSNwAcAAABA8gh8AAAAACSPwAcAAABA8gh8AAAAACSPwAcAAABA8oIPfMzsk2b2kJndVdHnfdnM9pnZP+X8/u/M7Ikq9gUAAAAgDMEHPpI+Jem8Cj/vSklv7fcLM5uSNF7hvgAAAAAEIPjAx91vlvTo7PfM7GVZz83tZvYNM3tFgc/7mqRfzH3fzBaoFxRdXrbMAAAAAMJyXNsFmKerJP2Bu//QzM6U9PeSzi35me+WtNXdHzSz0gUEAAAAEI7oAh8ze76k35D0j7MClBOy362X9ME+/2za3dcM+Mylki6SdE6lhQUAAAAQhOgCH/XS8/a5+6vn/sLdr5N03Tw+c1LSyyXtyoKp55nZLnd/eamSAgAAAAhC8GN85nL3xyXdZ2YXSZL1vKrkZ97g7i929xXuvkLSUwQ9AAAAQDqCD3zM7GpJ35S0ysx2m9k7JL1F0jvM7HuS7pa0tsDnfUPSP0p6ffZ5uSlwAAAAANJg7t52GQAAAACgVsH3+AAAAABAWUFPbrB48WJfsWJF28UAAAAAEKjbb7/9YXdfMmy7oAOfFStWaPv27W0XAwAAAECgzOwno2xHqhsAAACA5BH4AAAAAEgegQ8AAACA5BH4AAAAAEgegQ8AAACA5JWe1c3Mlkn6jKQXSzos6Sp3/9CcbUzShyRdIOkpSW9z9++U3TcAAEDdtuyY1pXbdmrPvv1aOj6my9as0rrJibaLBaCgKqazPijpT9z9O2b2S5JuN7OvuPv3Z21zvqTTsj9nSvpo9jcAAECwtuyY1hXX3an9zxySJE3v268rrrtTkgh+gMiUTnVz9wdnem/c/ReS7pE0906wVtJnvOdWSeNmdkrZfQMAANTpym07nw16Zux/5pCu3LazpRIBmK9Kx/iY2QpJk5Jum/OrCUkPzHq9W8cGRzOfscHMtpvZ9r1791ZZPAAAgEL27Ntf6H0A4aos8DGz50v6gqQ/cvfH5/66zz/xfp/j7le5+5S7Ty1ZsqSq4gEAABS2dHys0PsAwlVJ4GNmC9ULej7r7tf12WS3pGWzXp8qaU8V+wYAAKjLZWtWaWzhgqPeG1u4QJetWdVSiQDMV+nAJ5ux7ROS7nH3v8rZbKuk37Oe10p6zN0fLLtvAACAOq2bnNCm9at1/IJelWlifEyb1q9mYgMgQlXM6na2pLdKutPMvpu992eSlkuSu39M0o3qTWW9S73prN9ewX4BAABqt25yQld/635J0jW/f1bLpQEwX6UDH3f/F/UfwzN7G5f0rrL7AgAAKIt1eYBuqqLHBwAAIAqsywN0V6XTWQMAAISMdXmA7iLwAQAAncG6PEB3EfgAAIDOYF0eoLsIfAAAQGewLg/QXUxuAAAAOmNmAoPLr71DBw4d1gSzugGdQeADAACiV2SKatblAbqJwAcAAESNKarnjzWN0CUEPgAAIGqDpqhuuhIfUyBBwIiuYXIDAAAQtVCmqJ4JJKb37ZfrSCCxZcd0o+UYFWsaoWsIfAAAQNRCmaI6tkAilIARaAqBDwAAiFooU1THFkiEEjACTSHwAQAAUVs3OaFN61fr+AW9as3E+Jg2rV/d+DiV2AKJUAJGoCkEPgAAIHrrJic0uXxcZ65cpFs2ntvK4PzYAolQAkagKczqBgAAUIEYF0dlTSN0CYEPAABARQgkgHCR6gYAAAAgefT4AAAABC6mhVGBUBH4AAAABGxmYdSZNYJmFkaVRPADFECqGwAAQMBiWxgVCFUlgY+ZfdLMHjKzu3J+f46ZPWZm383+vK+K/QIAAKQutoVRgVBV1ePzKUnnDdnmG+7+6uzPByvaLwAAQNJiWxgVCFUlgY+73yzp0So+CwAAAEfEtjAqEKomx/icZWbfM7MvmdkrG9wvAABAtNZNTmjT+tU6fkGv2jYxPqZN61czsQFQUFOzun1H0kvc/Qkzu0DSFkmn9dvQzDZI2iBJy5cvb6h4AAAAzSoyRTULowLlNdLj4+6Pu/sT2c83SlpoZotztr3K3afcfWrJkiVNFA8AAKBRM1NUT+/bL9eRKaq37Jhuu2hAshoJfMzsxWZm2c9nZPt9pIl9A+iGLTumdfbmm7Ry4w06e/NNVB4ABI0pqoHmVZLqZmZXSzpH0mIz2y3p/ZIWSpK7f0zSmyX9RzM7KGm/pIvd3avYNwCwuB+A2DBFNdC8SgIfd79kyO8/LOnDVewLAOYa1HJK4AMgREvHxzTdJ8hhimqgPk3O6gYAtaDlFEBsmKIaaB6BD4DosbgfgNgwRTXQPAIfANGj5RRAjNZNTmhy+bjOXLlIt2w8l6AHqFlT6/gAQG1mKguXX3uHDhw6rIkh62EAAIDuIfABkAQW9wMAAIOQ6gYAAAAgeQQ+AAAAAJJHqhsAAACisWXHtK7ctlN79u3XUsZ0ogACHwAAAERhy45pXXHdnc8uWj29b7+uuO5OSSL4wVCkugEAACAKV27b+WzQM2P/M4d05badLZUIMSHwAQAAQBT27Ntf6H1gNgIfAAAARGHp+Fih94HZGOMDAACQkJQH/1+2ZtVRY3wkaWzhAl22ZlXf7VM+FiiOwAcAACARqQ/+n/k/XH7tHTpw6LAmBgQzqR8LFEfgAyBptPYB6JJBg/9Tufetm5zQ1d+6X5J0ze+flbtdF44FiiHwAZAsWvsAdA2D/4/gWGAuAh8AyaK1DziC3s9uWDo+puk+FfsuDv6fz7HgOkkbs7oBSBatfUDPTO/n9L79ch3p/dyyY7rtoqFil61ZpbGFC456b9Dg/5QVPRZcJ+kj8AGQLKY9BXpY9LE71k1OaNP61Tp+Qa+KNzE+pk3rV3ey16LoseA6SR+pbgCSVXTaUyBV9H52y6iD/7ugyLHgOkkfgQ+AZBWZ9hRIWdGxDoxzQBcxPip9laS6mdknzewhM7sr5/dmZn9rZrvM7A4ze00V+wWAYdZNTmhy+bjOXLlIt2w8l8obOqnIWAfGOaCrGB+VvqrG+HxK0nkDfn++pNOyPxskfbSi/QIAgCGKjHVgnAO6ivFR6ask1c3dbzazFQM2WSvpM+7ukm41s3EzO8XdH6xi/wAAYLBRxzowzgFdxviotDU1xmdC0gOzXu/O3jsm8DGzDer1Cmn58uWNFA4AAPQwzgEIU11j77o0pq+p6aytz3veb0N3v8rdp9x9asmSJTUXCwAAzMY4ByA8dY2969qYvqYCn92Sls16faqkPQ3tGwAAjIhxDkB46hp717UxfU2lum2V9G4z+5ykMyU9xvgeAIhLl9Ihuo5xDsjDfaAddY2969qYvkoCHzO7WtI5khab2W5J75e0UJLc/WOSbpR0gaRdkp6S9PYq9gsAaMZMOsRMy+BMOoQkKj1AR3AfaE9dY++6NqavklQ3d7/E3U9x94Xufqq7f8LdP5YFPfKed7n7y9x9tbtvr2K/AIBmdC0dAsCxuA+0p66xd10b09dUqhsAIGJdS4dAf6Q5dRv3gfbMXGeXX3uHDhw6rImKrr+6PjdUBD7oJB7eQDFdS4eIRZP3MtKcwH2gXUXG3hW5N3RpTF9Ts7oBweja1I1AFbqWDhGDpu9lpDmB+0AcqOfkI/BB5/DwBopjiuPwVHUv27JjWmdvvkkrN96gszfflFs5Is0Jdd4HRv0eYjjqOflIdUPnhP7wJg0PoepSOkQMqriXDUpfm4s0J0j13AdIo6xW6PWcNtHjg87Je0iH8PCmexrAqKq4lxVpGS6a5kQLPkbV5R6KOq6TkOs5bSPwQeeEnKPc5Zv/qKhMAT1V3MuKtAwXSXOiEQdFdLWHoq7rJOR6TtsIfNA5IY9V6OrNf1RUpoAjqriXFW0ZXjc5ocnl4zpz5SLdsvHc3H3RiIMiutpDUdd1EnI9p20EPuikUR/eTevqzX9UVKaAo+Xdy0btGa2rZZhGHBTR1R6KOq+TUOs5bSPwAQLS1Zv/qKhMAcMV6Rmtq2WYRhwU0dUeCq6T5hH4AAHp6s1/VDwkgOGK9ozW0TJMIw6K6mIPBddJ8wh8gMB08eY/Kh4SwHAh9IzSiAMMx3XSPNbxARCNmYfB5dfeoQOHDmuCdY6AY7Sx3k7e+mOs+wQMVvY6Ye2/YujxARAVesSAwZruGWW2RaAdXHvFEfgAAJCQptNnqphtkfW5gOKY6bQ4Ut0AAEhMk2lmg8YUnXrS8PS6mVbrmQrcTKu1JHp0gQFCGM8XG3p8AAC1oBW/G8rOtkirNTA/zHRaHIEPAKBybeSeE2i1o+yYIlqtgflhptPiCHyQBCo8QFiabsVnkG97yo4potUamB+mwy6OMT6IHvnhQHiabsUfFGhxH6hfmTFFl61ZddQ9XKLVGhhVkWuPqa8JfJAAKjzN4aaJUTW9lgzpUvFifS6gfjQS91SS6mZm55nZTjPbZWYb+/z+bWa218y+m/15ZxX7BSQqPE0hlQhFNJ17TrpU3FifK02koYejqvTj2M9p6cDHzBZI+oik8yWdLukSMzu9z6bXuPursz8fL7tfYAYVnmYw8xKKaDr3nEG+QFhoLAtLFY3EKZzTKnp8zpC0y93vdfcDkj4naW0FnwuMhApPM+hZQ1FNtuIzyBd5Ym+hjhWNZWGpopE4hXNaReAzIemBWa93Z+/N9dtmdoeZXWtmy/I+zMw2mNl2M9u+d+/eCoqH1FHhaQY9awgd6VKYK4UW6ljRWBaWKhqJUzinVQQ+1uc9n/P6f0la4e6/Jumrkj6d92HufpW7T7n71JIlSyooHrqACk/96FlDSugF6IYUWqhjRWNZWKpoJE7hnFYR+OyWNLsH51RJe2Zv4O6PuPvT2cv/IenXK9gvgAbRs4ZUxNgLQKA2Pym0UMeKxrLwlG0kTuGcVjGd9bclnWZmKyVNS7pY0u/O3sDMTnH3B7OXF0q6p4L9AmhYmbU6ACmMKdHzegE+sPXu1svWD9PQzl/T06rjCKYpT08K57R0j4+7H5T0bknb1AtoPu/ud5vZB83swmyz95jZ3Wb2PUnvkfS2svsFAMQllJ6WvNb+ffufab1s/ZCuNX8ptFDHjDT09MR+TitZx8fdb3T3f+XuL3P3v8jee5+7b81+vsLdX+nur3L317n7D6rYLwAgHqFU4Edt7Q8luCBda/5I0QUwWyWBDwAAw4RSge/XC5AnhOAihQHFbYq9hRpAdQh8gIgx4BkxCaUC368X4KTnLey7bQjBBelaAFCNKiY3ANACBjzHL4SB/k26bM2qo76zUnsV+LkTdcy9ntos21wpDCiuWteuHQDVIPABIjVovAQVgPDlBa7bf/Kovv6DvUlW6EKuwIdcNokZFWcb1OgDAIMQ+ACRCmW8BOYnL3D97K33P7sCdIq9eCFX4EMuG44Y1Ohz6kntpyYCCBdjfIBIhTJegnFG85MXoPqc16HMLAaEgkYfAPNF4ANEKoQBz6GsyxKjIgEqFbpuoTFhsFAafQDEh8AHiFQI61OEsi5LjPoFrpazbdEKHRXneNGYMFwIjT5AUdyXw8AYHyBibY9JIOVk/voNpn/dK5boC7dPl5pZrI3Z/phhqzqMXxlu0EQUM/dDICTMwhoOAh8Axxi1Irt0fEzTfYKcLqScVFHZ7xe4Tr1kUamZxZqe7Y8HerUGNSYQ+BzRdqMPUASzsIaDVDcARymSatPVlJM605HKrjLfdC8c6Y7VYvwKkB6yI8JB4APULLa83iIV2RDGGbWhjcr+qN+jpivOPNCPVvZ672pjApAyGjTCQeCD6MQUSMQ4ULloRbZsD0WMmq7sh9wLxwP9iCqu9642JgApo0EjHAQ+iEpsgUSMaUBUZIdr+hiF3AvHA/2Iqq73LjYmACmjQSMcBD6oRV29MrEFEjGmAVGRHa7pYxRyLxwP9CNivN4BNIMGjTAwqxsqV+csT0VnPGp7mt2QZj0b9VgMmioWPU0fo5C+R/3UNcNW29dvUaGfJwDoOgIf9FWmwlHntI1FKhYhTLN72ZpVR5VBaqf3pOixYKrY4Zo8Rm18j9oOOkK4fosK5XpHM9q+RgAUR6objlF2HE2d6R5FUoxCSIsLJQ0ohGOB+Wv6exTCWLoYv7OhXO+oXwjXCIDi6PGpUCqtP2V7bAb1ypQ9RkVW7A4l3z6E3pNQjgXmr8nv0aB7QFOLaMb6nQ3hesf8jfqMYkFKIE4EPhUpmpYRcpBUtsKRl+7xulcsqSR1ZdSKRRv59qGeV8YeoIiiY+nqwHcWTSvyHI81MAe6rpJUNzM7z8x2mtkuM9vY5/cnmNk12e9vM7MVVew3JEXSMkLvIi87VW9eusfXf7C30dSVpmfeCvm8MlMbighhSnO+s2haked4CNcIgOJKBz5mtkDSRySdL+l0SZeY2elzNnuHpJ+7+8sl/bWkvyy739AUaf0JPXe9igpHv2kbm24hazrfPuTzGuvYg5gWq01JCEFHrN9ZxKvIMyqEawTt4xkVH3P3ch9gdpakD7j7muz1FZLk7ptmbbMt2+abZnacpJ9KWuJDdj41NeXbt28vVb4qfPGdf6rxPfdpxcknHvX+jx95UpK04uQTteP+fXr64KFj/u0Jxy3Q5PLxo7a99d5Hcvf12peefNS2efsb9F4V2z78xNP60d4n5e464bgFWrZoTIuff0Khz/j+g49Lkk4/5QWSVOgYDSvz3M/Oe6/o+2WO56Dz+oKxhX3LUOT/V9e2edsXOW5FPjfv/bnvPfzE07r34Sd1+PCR28RznmN66eIT9dAvnh65bHnvV3HtlP3cvM+o4v9R5JwUuQdUcZ1V8T2s63qo4rtc9tgXLVsVx77p/dV1rsuUbdgzau72eddI3v7qvC83uW0V/7+Qj8WoZR70jHri6YOF/n913bequK7zjudPlyzT26/+O4XCzG5396lh21UxxmdC0gOzXu+WdGbeNu5+0Mwek3SypIfnfpiZbZC0QZKWL19eQfHKe/TJp3X8gWNvhk/Nem/ZorG+F8CyRWPHbHvCcQtyb65zt83b36D3qth28fNP6HvhFvmM5x1/dGtYkWM0bH9zPzvvvaLvlzmeg85rXhmK/P/q2jZv+yLHrej/b5Tj+cCj+4/6rkjS4cOuBx7dr5NOXDhy2fLer+LaKfu5eZ9Rxf+jyDkpcg+o4jqr4ntY1/VQxXe57LEvWrYqjn3T+6vrXJcp27Bn1Nzt866RvP0V+X88/MTT+tnjT8vd9fMnn3k2qKrreKb0jGrq/zfoGXXCwv4JVXXd2+u8rvOO56ITT+i7feiq6PG5SNIad39n9vqtks5w9/80a5u7s212Z69/lG2T30SucHp8fucfvinp2IH0c98fNLB99rZzB1BKvS7ymTSOUfdXpGxFt81T9jNGPUbzKVsVyhzPYee1yP6a3HY+24+qzPFcufEG9bs7maT7Nr+x1bJV9blVqOJcN32d5anjO75lx3QlC83Wdf8M5dg3rc77WRl1TlBT5Dsb47MklGdUU9sOekadsXLRyJ9bRdnqFEo5hmmyx2e3pGWzXp8qaU/ONruzVLdflvRoBfsOyrrJiZFukDPbhDj7V91GPUYx6vJ5rQOzeqGsmQrkgUOHJcWxCCraFcIziqmy49CFZ9SWHdPacf8+HTh0WGdvvimJOk0Vs7p9W9JpZrbSzI6XdLGkrXO22Srp0uznN0u6adj4ntStm5zQLRvP1X2b3/js4P9BZr58t933KAPoAlb0vNYhle8Kg4ebkcr3pZ+qJhxJ+RghPEyVHYfUn1F5DUex3/9KBz7uflDSuyVtk3SPpM+7+91m9kEzuzDb7BOSTjazXZL+s6RjprxGvlS/fLNRsahGSt+VmVm9JsbHZGJWrzqk9H3pp4oKZOrHCOFhquw4pP6MCnmm2jIqWcDU3W+UdOOc99436+f/J+miKvbVRSGsol6nvIrF9p88mlwXa91SS5EIIe0kxa7+Gal9X+aqIhUl9WOE8OQtAp5KT0JKQnhG1SXVnsdKFjBFvVL98s3Iq1h89tb7aWUtKPXvStNSb+1P/ftSRSrKoGNETzXqkHpPQii4fgdLteeRwCcCRb98sV3MeRWLuYPAUuhirVuqN6q2pNrVPyP170sVFci8Y/HLYwuTDorRrhDGi6Ys9UatKqQ6honAJwJFvnwxXsxFKlmhtESHGlzGeqMK9XjSIxK/shXIvGNkpqSDYiBlqTdqVSHVnsdKxvh0UZN5/3nTJEs6pgwxjgfql89sOrbHRwqjJTrkKXJjnFI75OOZ+nSlMX5fmpZ3jP74mu/23T6VoBjVSnmsYIxSb9SqSopjmAh85qGNitrcL19eGeYGPTP27NsfbODTr2Lxulcs0Rdunw5ycGcbg52LPDRju1GFPHi8C4OMY/u+tKHfMbpy286kg2JUJ+TGna5KvVEL+Uh1m4cQukjzyrDArO/2oV/Mc9NR/nzd6mC7WJtuKYoxfbGIkFveUu3qR3ldSBNENUKoM+BoXL/dRY/PPIRQUcvb1yF3jS1c0LeF+upv3d9U8SoRakt00y1FIfeIVCH0lrdQv4ddFkLaUNE0wRDKjHaEUGfA0Ujz7S4Cn3kIoaKWV4aJ7OLtdzHHFviEqun0p9Qfml1IJ0N1QkobGjUoHlRmpC+EOgOORaNWN5HqNg8hdJEOKgPTYFZr7oxjkhpNf2LKYeCIGNOGYiwzqhNCnQHlhDrzKIqjx2ceQugiDaEMXZDXUrtp/WrdsvHcRsrQhR4RWt4wqhh7QAeVOdRJZ1AdntdxC6mXGeUR+MxTCBW1ImXIyy8n73ywEMbX8NAEjogxbSjGMqNaIdQZMD8h1ANQHVLdOiCvteK/brkz6dnCqhBK6zLpi3HIS4cgTaI6MaYNxVhmIEZ13GtDqQegGgQ+Q1RxEbVd6clrrbj6tgfIOx8i9fE1qA4NDM2IcUxYXpklERADFalr6QfqAWkh1W2AKvI6Q8gNHTT1dZHtu6gL42vq0rU0ykENDHOvtarSJLp2jGfEmDY06iLUM9sCKKaulDTqAWmhx2eAKmbiCWE2n7xWiVgXO21SjK3LIRjU8tZ2D2hdmm5gSH1h29SF8GwAUlJXShr1gLTQ4zNAFRdRCLmhea0Vv/3rE/rC7dO0YgxRV+tyyq31eZW6D2y9W08fPJxkK3feAPYFZn2Dn7INDAy4jVsIzwYgJXVOIhJjLzP6o8dngCryOkPIDc1rrfjzdatpxWhJ6q31eZW3ffufSbaVO28A+yVnLqtlYDsV57iF8GwAUsIkIhgFPT4DVJHXGUpuaF5rBa0Y7Qiptb6Onqe8lrc8KVTWB007PvWSRZVPR84UyXEL5dkAVCGEDIaUln4I4XimisBngCouopQuRFQnlNb6ugZY51XqnrvwOfr5U88cs30qlfUmGxioOMeNZwNSEdJEHSk05oZ0PFNE4DNEFRdRChciqhVKa31dPU95lTpJVNYrklrFuYstnDwbkIKQMhhSwPGsF4EP0IJQWuvr7HkaVKlLpbLetlQqzrRwAvEKJYMhFRzPepUKfMxskaRrJK2Q9GNJ/87df95nu0OS7sxe3u/uF5bZLxC7UFrr2+h5SqWyjurQwgnEK5QMhlRwPOtVdla3jZK+5u6nSfpa9rqf/e7+6uwPQc8AMa5xEmOZQ7BuckK3bDxX921+o27ZeG4rFTxmwUEIaOEE4sVzpFocz3qVDXzWSvp09vOnJa0r+XmdFuMUxzGWGUewMBtCwNTOwPy13fjIc6RaHM96lR3j8yJ3f1CS3P1BM3thznbPNbPtkg5K2uzuW/I+0Mw2SNogScuXLy9ZvLjEmO4RY5lxNFLP0LZQxrwBsQllfBzPkWpxPOsztMfHzL5qZnf1+bO2wH6Wu/uUpN+V9Ddm9rK8Dd39KnefcvepJUuWFNhF/GJM94ixzADCQgsnMD+DGh/LarsnCajD0B4fd39D3u/M7GdmdkrW23OKpIdyPmNP9ve9Zva/JU1K+tH8ipyuGAe0xVhmAOGhhRMorq7Gx1B6koCqlR3js1XSpdnPl0r64twNzOwkMzsh+3mxpLMlfb/kfpMU44C2GMsM5KGFE0BM6hofV2dPEtCmsoHPZkm/aWY/lPSb2WuZ2ZSZfTzb5lckbTez70n6unpjfAh8+ogx3SPGMgP9MFEHgNjU1fhIGjtSVWpyA3d/RNLr+7y/XdI7s5//j6TVZfbTJTGme8RYZmCuLkzUMdOjdeDQYZ29+SYWjwUiV9eacKSxI1VlZ3UDgCSk3sJJzj6QpjoaH5lpEakqm+oGAElIfS0ZcvYBjIo0dqSKHh8gQxpQt6Xewpl6jxaAapHGjhTR4wOIge1Iv4Uz9R4tAACGoccHUDcGtmO4lFs4U+/RAgBgGHp8AJEGhPSl3qMFIE2sr4Yq0eMDiKk70Q0p92gBSA+zUaJq9PgAqm8ROAAAMD/MRomq0eMDqL5F4AAAwPyQho6qEfgAGdKAAAAIB2noqBqpbgAAAAgOaeioGj0+AAAACA5p6KgagQ8AAACCRBo6qkSqGwAAAIDkEfgAAAAASB6BD4DOYSVwAAC6h8AHQKfkrQRO8AMAQNoIfAB0CiuBAwDQTQQ+ADqFlcABAOgmAh8AnZK34jcrgQMAkLZSgY+ZXWRmd5vZYTObGrDdeWa208x2mdnGMvsEgDJYCRwAgG4q2+Nzl6T1km7O28DMFkj6iKTzJZ0u6RIzO73kfgFgXtZNTmjT+tWaGB+TSZoYH9Om9atZIA8AgMQdV+Yfu/s9kmRmgzY7Q9Iud7832/ZzktZK+n6ZfaN5M1MAHzh0WGdvvkmXrVlFZRFRYiVwAAC6p4kxPhOSHphG9WJpAAAIfklEQVT1enf2HiLCFMAAAACI2dDAx8y+amZ39fmzdsR99OsO8gH722Bm281s+969e0fcBerGFMAAAACDsUB22Iamurn7G0ruY7ekZbNenyppz4D9XSXpKkmamprKDZDQLKYABgAAyJeXHSOJ9OpANJHq9m1Jp5nZSjM7XtLFkrY2sF9UiCmAAQAA8pEdE76y01n/lpntlnSWpBvMbFv2/lIzu1GS3P2gpHdL2ibpHkmfd/e7yxUbTWMKYAAAgHxkx4Sv7Kxu10u6vs/7eyRdMOv1jZJuLLMvtGumi/bKbTu1Z99+LR0fY1Y3AABQSMozxC4dH9N0nyCH7JhwlAp80C1MAQwAAOYr9TEwl61ZpSuuu/OodDeyY8LSxBgfAAAAdFzqY2BYIDt89PgAAAB0QNtpZl0YA0N2TNjo8QEAAEhcCAuRM0Ms2kbgAwAAkLgQ0syYIRZtI9UNAAAgcSGkmTFDLNpG4AMAAJC4UKZaZgwM2kSqGwAAQOJIMwPo8QEAAEgeaWYAgQ8AAEAnkGaGriPVDQAAAEDyCHwAAECQZhbcvO2+R3X25psaXXMGQHoIfAAAQGNGDWZCWHATQFoIfAAAQCOKBDMhLLgJIC0EPgAAoBFFgpkQFtwEkBYCHwAA0IgiwUzewppNL7gJIB0EPgAAoBFFghkW3ARQNQIfAADQiCLBzLrJCW1av1oT42MySRPjY9q0fjXr0ACYNxYwBQAAjZgJWq7ctlN79u3X0vExXbZmVW4ww4KbAKpE4AMAABpDMAOgLaS6AQAAAEheqcDHzC4ys7vN7LCZTQ3Y7sdmdqeZfdfMtpfZJwAAAAAUVTbV7S5J6yX9wwjbvs7dHy65PwAAAAAorFTg4+73SJKZVVMaAAAAAKhBU2N8XNI/m9ntZrZh0IZmtsHMtpvZ9r179zZUPAAAAAApG9rjY2ZflfTiPr96r7t/ccT9nO3ue8zshZK+YmY/cPeb+23o7ldJuirb914z+8mI+6jbYkmk6sWL8xc3zl/cOH/x4tzFjfMXN87f6F4yykZDAx93f0PZkrj7nuzvh8zseklnSOob+Mz5d0vK7rsqZrbd3XMncEDYOH9x4/zFjfMXL85d3Dh/ceP8Va/2VDczO9HMfmnmZ0n/Vr1JEQAAAACgEWWns/4tM9st6SxJN5jZtuz9pWZ2Y7bZiyT9i5l9T9K3JN3g7l8us18AAAAAKKLsrG7XS7q+z/t7JF2Q/XyvpFeV2U8grmq7ACiF8xc3zl/cOH/x4tzFjfMXN85fxczd2y4DAAAAANSqqemsAQAAAKA1BD4AAAAAkkfgMwIzO8/MdprZLjPb2HZ5kM/MlpnZ183sHjO728z+MHt/kZl9xcx+mP19UttlRT4zW2BmO8zsn7LXK83stuz8XWNmx7ddRvRnZuNmdq2Z/SC7Ds/i+ouHmf1xdu+8y8yuNrPncv2Fy8w+aWYPmdlds97re71Zz99mdZk7zOw17ZUcUu75uzK7f95hZteb2fis312Rnb+dZramnVLHjcBnCDNbIOkjks6XdLqkS8zs9HZLhQEOSvoTd/8VSa+V9K7sfG2U9DV3P03S17LXCNcfSrpn1uu/lPTX2fn7uaR3tFIqjOJDkr7s7q9Qb2Kbe8T1FwUzm5D0HklT7v6rkhZIulhcfyH7lKTz5ryXd72dL+m07M8GSR9tqIzI9ykde/6+IulX3f3XJP1fSVdIUlaXuVjSK7N/8/dZHRUFEPgMd4akXe5+r7sfkPQ5SWtbLhNyuPuD7v6d7OdfqFfpmlDvnH062+zTkta1U0IMY2anSnqjpI9nr03SuZKuzTbh/AXKzF4g6d9I+oQkufsBd98nrr+YHCdpzMyOk/Q8SQ+K6y9Y7n6zpEfnvJ13va2V9BnvuVXSuJmd0kxJ0U+/8+fu/+zuB7OXt0o6Nft5raTPufvT7n6fpF3q1VFRAIHPcBOSHpj1enf2HgJnZiskTUq6TdKL3P1BqRccSXpheyXDEH8j6XJJh7PXJ0vaN+tBwDUYrpdK2ivpf2apih/PFq7m+ouAu09L+u+S7lcv4HlM0u3i+otN3vVGfSY+/0HSl7KfOX8VIPAZzvq8xxzggTOz50v6gqQ/cvfH2y4PRmNmb5L0kLvfPvvtPptyDYbpOEmvkfRRd5+U9KRIa4tGNhZkraSVkpZKOlG99Ki5uP7ixL00Imb2XvXS9z8781afzTh/BRH4DLdb0rJZr0+VtKelsmAEZrZQvaDns+5+Xfb2z2a69LO/H2qrfBjobEkXmtmP1UsrPVe9HqDxLPVG4hoM2W5Ju939tuz1teoFQlx/cXiDpPvcfa+7PyPpOkm/Ia6/2ORdb9RnImFml0p6k6S3+JEFNzl/FSDwGe7bkk7LZrU5Xr2BZVtbLhNyZONBPiHpHnf/q1m/2irp0uznSyV9semyYTh3v8LdT3X3Fepdaze5+1skfV3Sm7PNOH+BcvefSnrAzFZlb71e0vfF9ReL+yW91syel91LZ84f119c8q63rZJ+L5vd7bWSHptJiUM4zOw8Sf9F0oXu/tSsX22VdLGZnWBmK9WbpOJbbZQxZnYkkEQeM7tAvVbnBZI+6e5/0XKRkMPM/rWkb0i6U0fGiPyZeuN8Pi9puXoP94vcfe6AUATEzM6R9Kfu/iYze6l6PUCLJO2Q9O/d/ek2y4f+zOzV6k1McbykeyW9Xb1GNq6/CJjZf5P0O+ql2OyQ9E71xhFw/QXIzK6WdI6kxZJ+Jun9kraoz/WWBbMfVm9GsKckvd3dt7dRbvTknL8rJJ0g6ZFss1vd/Q+y7d+r3rifg+ql8n9p7mdiMAIfAAAAAMkj1Q0AAABA8gh8AAAAACSPwAcAAABA8gh8AAAAACSPwAcAAABA8gh8AAAAACSPwAcAAABA8v4/1DT4jchd7L0AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x22339f44ef0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x_hat = np.dot(W.T.conjugate(), X) / N\n",
"\n",
"plt.stem(np.real(x_hat - x));\n",
"plt.show();\n",
"plt.stem(np.imag(x_hat));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, the error is very small but clearly not zero."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The FFT Algorithm\n",
"\n",
"The FFT algorithm computes the DFT recursively by successively splitting the data vector into smaller pieces and recombining the results. The most well-known version of the FFT operates on data lengths that are a power of two but efficient algorithms exist for all lengths that are factorizable into powers of small primes. \n",
"\n",
"The FFT algorithm is not only much faster than the direct method but it's also better conditioned numerically. This is because in the FFT implementation great care is applied to minimizing the number of trigonometric factors.\n",
"\n",
"As we can see in the examples below, the phase is now accurate and the reconstruction error is almost two orders of magnitude smaller, basically equal to the numerical precision of floating point variables. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAD8CAYAAAC2A3H6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAG/ZJREFUeJzt3X2wnFd9H/DvD/kFhTYjiE2CZCt2UtcEUBMnd3ipOhkgUBvKYEclxYS2hIZR0gmTNNM6tctMk2YmI6XutEkLTeoSCnQYQ4aAUYsTJUF0SNvgIEcJb8aJYxNbkhMbjKCJNTaST//YvfG1tHvvXe377uczo9E+zx7tc5599OjuT+d7zlZrLQAAAIvsadPuAAAAwLgpfAAAgIWn8AEAABaewgcAAFh4Ch8AAGDhKXwAAICFp/ABAAAWnsIHAABYeAofAABg4Z037Q6s56KLLmqXXXbZtLsBAADMqDvvvPNLrbWLN2o304XPZZddlsOHD0+7GwAAwIyqqj/dTDtRNwAAYOEpfAAAgIWn8AEAABaewgcAAFh4Ch8AAGDhzfSqbgAAwGy47cix3Hzw7hw/cTLbt23NDVdfmeuu2jHtbm2awgcAAFjXbUeO5aYPfSYnv346SXLsxMnc9KHPJMncFD+ibgAAwLpuPnj3XxU9q05+/XRuPnj3lHo0uJGM+FTVu5K8JslDrbUX9Hj+pUk+kuS+7q4PtdZ+dhTHBgAAzs1m42vHT5zs+ef77Z9Fo4q6vTvJ25O8d502v9Nae82IjgcAAAxhkPja9m1bc6xHkbN929bxd3RERhJ1a619Iskjo3gtAABg/AaJr91w9ZXZev6Wp+zbev6W3HD1lWPt4yhNco7PS6rqD6vq16vq+f0aVdXeqjpcVYcffvjhCXYPAACWxyDxteuu2pF9e3blgi2d8mHHtq3Zt2fX3CxskExuVbffT/KtrbW/qKpXJ7ktyRW9GrbWbklyS5KsrKy0CfUPAAAWwmbn7QwaX7vuqh259ffuT5J84EdeMtpOT8BERnxaa19rrf1F9/HtSc6vqosmcWwAAFgWq/N2jp04mZYn5+3cduTYWW0XIb42iIkUPlX1LVVV3ccv7B73y5M4NgAALItB5u0sQnxtEKNazvrWJC9NclFVHU3y00nOT5LW2i8neV2Sf1pVp5KcTHJ9a02MDQAANmFcy07Pe3xtECMpfFprb9jg+bens9w1AAAwgGVbdnpcJrmqGwAAMKBlW3Z6XCa1qhsAALDGOOJrq3/+pz746Tx++onsWOd1l43CBwAAJmyc8bVlmrczCFE3AACYMPG1yTPiAwAAI7DZ6FoivjYNCh8AABjSING1RHxtGkTdAABgSINE1xLxtWkw4gMAAH2M84tDE/G1SVL4AABAD+P+4lDxtckSdQMAgB6svLZYjPgAALBUfHHoclL4AACwNHxx6PISdQMAYGmIry0vIz4AAMw98TU2ovABAGCuia+xGaJuAADMNfE1NsOIDwAAM0l8jVFS+AAAMHPE1xg1UTcAAGaO+BqjZsQHAICZI77GqCl8AACYiM3O2UnE1xg9UTcAAMZudc7OsRMn0/LknJ3bjhzr2V58jVFT+AAAMHaDzNlJOiM4+/bsygVbOh9Xd2zbmn17domvcc5GEnWrqncleU2Sh1prL+jxfCX5xSSvTvJokh9qrf3+KI4NAMD0jGPJ6VXia4zSqEZ83p3kmnWef1WSK7q/9ib5pREdFwCAKRkkvtZvbk6//TBqIyl8WmufSPLIOk2uTfLe1vHJJNuq6jmjODYAANNhyWnmyaRWdduR5IE120e7+x48s2FV7U1nVCg7d+6cSOcAAHjSOOJrlpxm2iZV+FSPfa1Xw9baLUluSZKVlZWebQAAGI/V+NrqSM5qfC3JWUWKJaeZJ5Na1e1okkvXbF+S5PiEjg0AwCaJr7GoJjXicyDJW6vq/UlelOSrrbWzYm4AAIyH+BrLblTLWd+a5KVJLqqqo0l+Osn5SdJa++Ukt6ezlPU96Sxn/eZRHBcAgI2Jr8GICp/W2hs2eL4l+bFRHAsAgMGsF187s/C54eorn1IkJeJrLIZJRd0AABgx8TXYPIUPAMAcEl+DwUxqVTcAAEbI6mswGCM+AAAzRHwNxkPhAwAwI8TXYHxE3QAAZoT4GoyPER8AgDHabHQtEV+DcVL4AACMySDRtUR8DcZJ1A0AYEwGia4l4mswTkZ8AADGZJDoWiK+BuOk8AEAGNBm5+0MGl1LxNdgXETdAAAGsDpv59iJk2l5ct7ObUeOndVWdA1mh8IHAGAAg8zbue6qHdm3Z1cu2NL5yLVj29bs27NLdA2mQNQNACCbj6+dy7wd0TWYPiM+AMDSGyS+1m9+znrzdoDpU/gAAEtvkPiaeTswn0TdAICFNY74miWnYT4pfACAhbQaX1sdyVmNryU5q0gZdNlp83Zg/oi6AQALSXwNWMuIDwAwV8TXgHOh8AEA5ob4GnCuRN0AgLkhvgacq5EUPlV1TVXdXVX3VNWNPZ7/oap6uKr+oPvrLaM4LgAw/247ciy79x/K5Td+NLv3H+r53TmrBo2v7duzKxds6Xzc2bFta/bt2SW+Bktq6KhbVW1J8o4kr0xyNMmnqupAa+3zZzT9QGvtrcMeDwBYHINE1xLxNeDcjWLE54VJ7mmt3dtaezzJ+5NcO4LXBQAW3CDRtUR8DTh3o1jcYEeSB9ZsH03yoh7t/n5VfW+SP0ryk621B3q0AQAWwDhWXkusvgacu1EUPtVjXztj+38kubW19lhV/WiS9yR5ec8Xq9qbZG+S7Ny5cwTdAwAmaZwrr62+hvgaMKhRRN2OJrl0zfYlSY6vbdBa+3Jr7bHu5n9N8j39Xqy1dktrbaW1tnLxxRePoHsAwCRZeQ2YRaMY8flUkiuq6vIkx5Jcn+QH1zaoque01h7sbr42yV0jOC4AMEG+OBSYZ0MXPq21U1X11iQHk2xJ8q7W2ueq6meTHG6tHUjy41X12iSnkjyS5IeGPS4AMDm+OBSYdyP5Hp/W2u2ttb/ZWvv21trPdff9627Rk9baTa2157fWvrO19rLW2hdGcVwAYDLE14B5N4qoGwCw4MTXgHmn8AGAJbbZeTvia8C8G0nUDQCYP6vzdo6dOJmWJ+ft3Hbk2FltxdeAeafwAYAlNci8neuu2pF9e3blgi2djw47tm3Nvj27xNeAuSHqBgALZhzLTifia8B8M+IDAAtkkPhav/k5/fYDzDOFDwAsEMtOA/Qm6gYAM26z0bXEstMA/Sh8AGCGrUbXVkdxVqNrSSw7DTAAUTcAmGGDRNcS8TWAfoz4AMAUjHPltUR8DeBMCh8AmLBB4muDRtdWX0N8DeCpRN0AYMKsvAYweUZ8AGBExhFfE10DGA2FDwCMwDjja6JrAMMTdQOAERBfA5htRnwAYB3iawCLQeEDAH2IrwEsDlE3AOhDfA1gcRjxAWDpiK8BLB+FDwBLRXwNYDmJugGwVMTXAJaTER8Alor4GsByGknhU1XXJPnFJFuSvLO1tv+M5y9M8t4k35Pky0le31r74iiODQDJ5uftiK8BLKehC5+q2pLkHUlemeRokk9V1YHW2ufXNPvhJF9prf2Nqro+yc8nef2wx56Ezf4gHWfbWemH8/NeOD/vxayeX5Ke83YO/+kj+fgXHn5K2xuuvvIpbZNOfO1lz704u/cfmsnzm+W2s9IP5+e9cH7Tey/mxShGfF6Y5J7W2r1JUlXvT3JtkrWFz7VJfqb7+INJ3l5V1VprIzj+2AwyAXZcbWelH87Pe+H8vBezfH5PP/9pPeftvO+T92f1B81q2317dmXfnl1Pia+97LkX59fuPDaS85uFtrNyrZ3fbJ+f92J5zm+c78U8qWFrj6p6XZJrWmtv6W7/oyQvaq29dU2bz3bbHO1u/0m3zZfWe+2VlZV2+PDhofo3jN37D+XYiZP5kU9/JN/21WN/tf/C87bkqp3b8sUv/2WS5LJvekaO3H8ij506fdZrDNs2ybrtn/mM88fS1vl5L5yf92La53felqeltZbTT7RceN6WXPqsrXngkZM92w5itc+ff/BrSZLnPecbNzy/tW03Or8Lz3/a1NsOc37z8F44P++F85vue/FnF1+aN9/6n856flqq6s7W2spG7UYx4lM99p1ZTW2mTadh1d4ke5Nk586dw/VsSP0mwK7+RXj08dNn7Rt1243aP/r408bSdpg+Oz/vxTKen/di9Od36vQTT2lz75f+Mk88MXxQYPV433DBlrP2babtRu1XC7tptj2zz+NqO8o+O7/R9Nl74fwm8V486xkX9nx+1o1ixOclSX6mtXZ1d/umJGmt7VvT5mC3ze9W1XlJ/izJxRtF3WZlxOdMO7Ztzf+58eV5/X/53SSdya7jartRPy555taxtHV+3gvn572Y9vn1sqUqp3v86Ni29fw8duqJp8TdKr3/h221z2ttdH5nGqT9IredlX44v/G3nZV+OL/xtz2X9tO22RGfUXyPz6eSXFFVl1fVBUmuT3LgjDYHkryp+/h1SQ7N+vyeZLDvbxhX21nph/Mbf9tZ6YfzG3/bWenHNM7vZc+9OEfuP5E77nsku/cfym1HjvVs28/p1nq+7s+89vnZt2dXdmzbmkrnh/MbX7zT9Rtz21nph/Mbf9tZ6YfzG3/bc2k/L4aOurXWTlXVW5McTGc563e11j5XVT+b5HBr7UCSX0ny36vqniSPpFMczbzVyV6bWdFiXG03ar+6xOqo2zo/74Xz816M+vxWFxB4vBthO3OxgbVtH338VL7y6NfPet3V79Hp148z+7Pyrc+a+PVbpraz0g/n571wftN7L+bJ0FG3cZp21G0jZ0Y7FqntrPTD+Y2/7az0w/mNv+20+zFIdOLMFYWSzv827tuza+5/8AIwWpOMugGwxG47cuys+Fov/RaM6bX/uqt2nBVdU/QAMIxRrOoGwJJaHZk5M76WnB07275ta88Rn+3btvZ87euu2qHQAWBkjPgAcM5uPnh3zy8Ovfng3We1XdTJsgDMB4UPAGcRXwNg0Yi6AfAU4msALCIjPgA8hfgaAItI4QOwJMTXAFhmom4AS0B8DYBlZ8QHYAmIrwGw7BQ+AHNqs9G1RHwNAETdAObQING1RHwNAIz4AMyhQaJrifgaACh8AGbIOFZeS8TXAEDUDWBGjHPltdXXUOgAsKyM+ADMCCuvAcD4GPEBmBGDrryWdIql4ydOZvu2rbnh6iuN6ABAHwofgDFbnbfz+Oknsnv/ob4FipXXAGB8RN0AxqjfvJ1eixaIrwHA+Ch8AMZokHk7Vl4DgPERdQM4B5uNr53LstMKHQAYPSM+AAMaJL7Wb37OestOAwCjp/ABGJBlpwFg/ih8ALpW42t33PdIdu8/1HMEJxl82WnzdgBg+oaa41NVz0rygSSXJflikn/QWvtKj3ank3ymu3l/a+21wxwXYNT6xdeSnFWkWHYaAObPsCM+Nyb5WGvtiiQf6273crK19l3dX4oeYOaIrwHAYhu28Lk2yXu6j9+T5LohXw9gpMTXAIBk+OWsv7m19mCStNYerKpn92n39Ko6nORUkv2ttduGPC7AhsTXAIBVG474VNVvV9Vne/y6doDj7GytrST5wSS/UFXfvs7x9lbV4ao6/PDDDw9wCICnEl8DAFZtOOLTWntFv+eq6s+r6jnd0Z7nJHmoz2sc7/5+b1X9ryRXJfmTPm1vSXJLkqysrLQNzwBYKpv94tBk8Pha0imWjp84me3btq772gDAfBk26nYgyZuS7O/+/pEzG1TVM5M82lp7rKouSrI7yb8d8rjAEhokupaIrwEATxp2cYP9SV5ZVX+c5JXd7VTVSlW9s9vmO5Icrqo/TPLxdOb4fH7I4wJLaJDoWiK+BgA8aagRn9bal5N8X4/9h5O8pfv4/ybZNcxxgMW22fjaING1RHwNAHjSsFE3gKGMc+W11ddQ6AAAw0bdAIZi5TUAYBIUPsBY+OJQAGCWiLoBI+eLQwGAWWPEBxg58TUAYNYY8QFGzheHAgCzRuEDbNpml50WXwMAZo2oG7Ap/ebt9Fq0QHwNAJg1Ch9gUwaZt2P1NQBg1oi6wZLbbHxtkHk7ifgaADBbjPjAEhskvtZvfk6//QAAs0ThA0vMstMAwLJQ+MCCWY2u3XHfI9m9/1DP0ZtVgy47bd4OADCvzPGBBdIvupbEstMAwFIz4gMLZJDoWiK+BgAsD4UPzIHNxtfOZeU18TUAYBmIusGMGyS+Nmh0bfU1FDoAwKIz4gMzzsprAADDU/jAlIwjvia6BgDQm6gbTME442uiawAAZzPiA1MgvgYAMFkKHxgh8TUAgNkk6gYjIr4GADC7hhrxqaofqKrPVdUTVbWyTrtrquruqrqnqm4c5pgwq8TXAABm17BRt88m2ZPkE/0aVNWWJO9I8qokz0vyhqp63pDHhYkRXwMAmH9DRd1aa3clSVWt1+yFSe5prd3bbfv+JNcm+fwwx4ZJEF8DAFgMk1jcYEeSB9ZsH+3ug5knvgYAsBg2HPGpqt9O8i09nnpba+0jmzhGr+Ggts7x9ibZmyQ7d+7cxMvD+AwaX0s6xdLxEyezfdvW3HD1lUZ1AABmwIaFT2vtFUMe42iSS9dsX5Lk+DrHuyXJLUmysrLSt0CCYazO23n89BPZvf9Q3wJFfA0AYDFMIur2qSRXVNXlVXVBkuuTHJjAcaGnfvN2ei1aIL4GALAYhl3O+vur6miSlyT5aFUd7O7fXlW3J0lr7VSStyY5mOSuJL/aWvvccN2GczfIvB2rrwEALIZhV3X7cJIP99h/PMmr12zfnuT2YY4F69lsdC0ZbN5OIr4GALAIJhF1g7EaJLqW9J+f028/AADzT+HD3BskupaYtwMAsIwUPsys1fjaHfc9kt37D/UdwTmX6Jp5OwAAy2WoOT4wLv3ia0nOKlAGXXJ69TUUOgAAy8OIDzNpkPia6BoAABtR+DBR44ivia4BALARUTcmZpzxNdE1AADWY8SHiRFfAwBgWhQ+DE18DQCAWSfqxlDE1wAAmAdGfBiK+BoAAPNA4UNP4msAACwSUTfOIr4GAMCiMeLDWcTXAABYNAqfJbHZ6FoivgYAwOIRdVsCg0TXEvE1AAAWjxGfJTBIdC0RXwMAYPEY8VkCg0TXkidHgW4+eHeOnziZ7du25oarrzSqAwDA3FL4zLHVeTuPn34iu/cf6lucDBpdS8TXAABYLKJuc6rfvJ1eixaIrgEAsOwUPnNqkHk7Vl4DAGDZibrNmM3G185l3o5CBwCAZWXEZ4YMEl/rNz9nvXk7AACwrIYqfKrqB6rqc1X1RFWtrNPui1X1mar6g6o6PMwxF9kg8TXzdgAAYPOGHfH5bJI9ST6xibYva619V2utb4G0qFbja3fc90h27z/UcwQnGSy+Zt4OAABs3lBzfFprdyVJVY2mNwuoX3wtyVlFyqDLTpu3AwAAmzOpOT4tyW9W1Z1VtXdCx5wJ4msAADB9G474VNVvJ/mWHk+9rbX2kU0eZ3dr7XhVPTvJb1XVF1prPeNx3cJob5Ls3Llzky8/eeNYfW31z9988O4cP3Ey27dt7fu6AADA5m1Y+LTWXjHsQVprx7u/P1RVH07ywvSZF9RauyXJLUmysrLShj32OIivAQDAfBl71K2qnlFVf331cZK/m86iCHNLfA0AAObLsMtZf39VHU3ykiQfraqD3f3bq+r2brNvTvK/q+oPk/xeko+21n5jmONOm9XXAABgvgy7qtuHk3y4x/7jSV7dfXxvku8c5jizRnwNAADmy6RWdVso4msAADBfhhrxWVZWXwMAgPmi8DlH4msAADA/RN0AAICFp/ABAAAWnsIHAABYeAofAABg4Sl8AACAhVettWn3oa+qejjJn067H10XJfnStDvBOXP95pvrN99cv/nl2s0312++uX6b962ttYs3ajTThc8sqarDrbWVafeDc+P6zTfXb765fvPLtZtvrt98c/1GT9QNAABYeAofAABg4Sl8Nu+WaXeAobh+8831m2+u3/xy7eab6zffXL8RM8cHAABYeEZ8AACAhafw2YSquqaq7q6qe6rqxmn3h/6q6tKq+nhV3VVVn6uqn+juf1ZV/VZV/XH392dOu6/0V1VbqupIVf3P7vblVXVH9/p9oKoumHYf6a2qtlXVB6vqC9378CXuv/lRVT/Z/bfzs1V1a1U93f03u6rqXVX1UFV9ds2+nvdbdfzH7meZT1fVd0+v5yR9r9/N3X8/P11VH66qbWueu6l7/e6uqqun0+v5pvDZQFVtSfKOJK9K8rwkb6iq5023V6zjVJJ/3lr7jiQvTvJj3et1Y5KPtdauSPKx7jaz6yeS3LVm++eT/Ifu9ftKkh+eSq/YjF9M8huttecm+c50rqP7bw5U1Y4kP55kpbX2giRbklwf998se3eSa87Y1+9+e1WSK7q/9ib5pQn1kf7enbOv328leUFr7W8l+aMkNyVJ97PM9Ume3/0z/7n7GZUBKHw29sIk97TW7m2tPZ7k/UmunXKf6KO19mBr7fe7j/9fOh+6dqRzzd7TbfaeJNdNp4dspKouSfL3kryzu11JXp7kg90mrt+MqqpvTPK9SX4lSVprj7fWTsT9N0/OS7K1qs5L8g1JHoz7b2a11j6R5JEzdve7365N8t7W8ckk26rqOZPpKb30un6ttd9srZ3qbn4yySXdx9cmeX9r7bHW2n1J7knnMyoDUPhsbEeSB9ZsH+3uY8ZV1WVJrkpyR5Jvbq09mHSKoyTPnl7P2MAvJPmpJE90t78pyYk1Pwjcg7Pr25I8nOS/daOK76yqZ8T9Nxdaa8eS/Lsk96dT8Hw1yZ1x/82bfvebzzPz558k+fXuY9dvBBQ+G6se+yyFN+Oq6q8l+bUk/6y19rVp94fNqarXJHmotXbn2t09mroHZ9N5Sb47yS+11q5K8pcRa5sb3bkg1ya5PMn2JM9IJx51JvfffPJv6RypqrelE99/3+quHs1cvwEpfDZ2NMmla7YvSXJ8Sn1hE6rq/HSKnve11j7U3f3nq0P63d8fmlb/WNfuJK+tqi+mEyt9eTojQNu60ZvEPTjLjiY52lq7o7v9wXQKIffffHhFkvtaaw+31r6e5ENJ/nbcf/Om3/3m88ycqKo3JXlNkje2J793xvUbAYXPxj6V5IruqjYXpDOx7MCU+0Qf3fkgv5Lkrtbav1/z1IEkb+o+flOSj0y6b2ystXZTa+2S1tpl6dxrh1prb0zy8SSv6zZz/WZUa+3PkjxQVVd2d31fks/H/Tcv7k/y4qr6hu6/pavXz/03X/rdbweS/OPu6m4vTvLV1Ugcs6OqrknyL5O8trX26JqnDiS5vqourKrL01mk4vem0cd55gtMN6GqXp3O/zpvSfKu1trPTblL9FFVfyfJ7yT5TJ6cI/Kv0pnn86tJdqbzw/0HWmtnTghlhlTVS5P8i9baa6rq29IZAXpWkiNJ/mFr7bFp9o/equq70lmY4oIk9yZ5czr/yeb+mwNV9W+SvD6diM2RJG9JZx6B+28GVdWtSV6a5KIkf57kp5Pclh73W7eYfXs6K4I9muTNrbXD0+g3HX2u301JLkzy5W6zT7bWfrTb/m3pzPs5lU6U/9fPfE3Wp/ABAAAWnqgbAACw8BQ+AADAwlP4AAAAC0/hAwAALDyFDwAAsPAUPgAAwMJT+AAAAAtP4QMAACy8/w+MGufkLoRbpQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2233a65f2b0>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAEDCAYAAADjkkFLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XFwZWd53/Hf47uyK0wZ2bEprHbttVNHQDFBVMEQNylxSOUkHlvjCYNdoC5JxsAkhLRBDor/SJsJY3eUSUInCe2OMaRTz5rYK4QnDSiASUM7ZY0WYS9gRBnM7vou4IWtgMKtV7779A9d7Urac67ue897dM977/cz42H16tVznnPe855zH+695zV3FwAAAACk6IJeJwAAAAAA3aKgAQAAAJAsChoAAAAAyaKgAQAAAJAsChoAAAAAyaKgAQAAAJCsnhU0Zna/mT1jZl+MFO/jZrZiZn+9pd3M7L1m9lUze9LMfivG9gAAAAD0Xi/fofmQpBsjxpuV9JaM9n8taa+kl7j7SyU9GHGbAAAAAHqoZwWNu/+9pFMb28zsx1vvtBw2s8+Y2UsC4n1K0g8yfvUOSX/g7mda/Z4pkjcAAACA6qjad2j2S3qnu/9TSe+W9BcRYv64pDea2aKZfczMrokQEwAAAEAF7Op1AuvM7PmSflrSQ2a23nxR63e3SvqDjD+ru/vkNqEvkvT/3H2iFed+ST8TJ2sAAAAAvVSZgkZr7xatuPsrt/7C3eckzXUZ92lJB1v//oikD3YZBwAAAEDFVOYjZ+7+fUlPmdkbpLNPJ/vJCKHnJd3Q+vc/l/TVCDEBAAAAVIC5e282bHZA0uskXSbp25J+X9Kjkt4v6cWShiQ96O5ZHzXLivcZSS+R9HxJ35X0a+6+YGYjkh6QdIWk/yvp7e7+eNy9AQAAANALPStoAAAAAKCoynzkDAAAAABC9eShAJdddpnv27evF5sGAAAAkIDDhw9/x90v365fTwqaffv2aXFxsRebBgAAAJAAMzvaST8+cgYAAAAgWRQ0AAAAAJJFQQMAAAAgWRQ0AAAAAJJFQQMAAAAgWT15yhkAoP/ML9U1u7CsEysN7R4Z1vTkmKbGR3udVq7U8gUAZKOgAQAUNr9U18zcETVWm5Kk+kpDM3NHJKmSRUJq+QIA8vGRMwBAYbMLy2eLg3WN1aZmF5Z7lFF7qeULAMhHQQMAKOzESiOovddSyxcAkI+CBgBQ2O6R4aD2XkstXwBAPgoaAEBh05NjGh6qbWobHqppenKsRxm1l1q+AIB8PBQAAFDY+hfp73r4CZ1untFoxZ8allq+AIB8FDQAgCimxkd14LFjkqQPv+21Pc5me6nlCwDIRkEDAKg01osBALRDQQMAqCzWiwEAbCfKQwHMbMTMHjazr5jZk2bGe/cAgMJYLwYAsJ1Y79C8T9LH3f1XzOxCSc+LFBcAMMBYLwYAsJ3C79CY2Qsk/aykD0iSu59295WicQEAYL0YAMB2Ynzk7GpJJyV90MyWzOw+M7t4ayczu9PMFs1s8eTJkxE2CwDod6wXAwDYToyCZpekV0l6v7uPS/qhpPds7eTu+919wt0nLr/88gibBQD0u6nxUd1z67W6sLZ2uxodGdY9t17LAwEAAGfF+A7N05KedvdDrZ8fVkZBAwBAN1gvBgDQTuF3aNz9W5KOm9n6+/8/L+nLReMCAAAAwHZiPeXsnZIeaD3h7OuS3hopLgAAQM+xwCtQXVEKGnf/gqSJGLEAAACqhAVegWqLsrAmAABAv2KBV6DaKGgAAADaYIFXoNooaAAAANpggVeg2ihoAAAA2mCBV6DaYj3lDAAAoC+tf/H/roef0OnmGY3ylDOgUihoAAAAtsECr0B18ZEzAAAAAMmioAEAAACQLAoaAAAAAMmioAEAAACQLAoaAAAAAMmioAEAAACQLAoaAAAAAMmioAEAAACQLBbW7ND8Ul2zC8s6sdLQ7tYKwZLOa6vKqsFZ+U6Nj+a2A6iWkLma4rxOLefU8g2V4v6lmDOAclDQdGB+qa6ZuSNqrDYlSfWVhqYfelwyabXpZ9tm5o5IUs8vqFn5zswd0eLRUzp4uH5eu9T7nAGckzeHpfPnakjfqkgt59TyDZXi/qWYM4Dy8JGzDswuLJ+9aK5bPeNni5l1jdWmZheWdzK1TFn5NlabOnDoeGZ7FXIGcE7eHM6aqyF9qyK1nFPLN1SK+5dizgDKQ0HTgRMrjVL6liUvh6Z7ZnsVcgZwTt6czGoP6VsVqeWcWr6hUty/FHMGUB4Kmg7sHhkupW9Z8nKomQX1B9AbeXMyqz2kb1WklnNq+YZKcf9SzBlAeShoOjA9OabhodqmtqELTEO1zQXC8FDt7MMCeikr3+Ghmm6/bm9mexVyBnBO3hzOmqshfasitZxTyzdUivuXYs4AysNDATqw/gXDux5+QqebZzS64SlnW9uq8GXEvHynxkc1ceWllcwZwDnt5nCRvlWRWs6p5Rsqxf1LMWcA5YlW0JhZTdKipLq73xQrblVMjY/qwGPHJEkffttrz7ZntVVBXr557QCqJWSupjivU8s5tXxDpbh/KeYMoBwx36F5l6QnJb0gYkwAAJCwFNeLSTFnYJBFKWjMbI+kX5b0Xkn/NkZMAACQthTXi2mXM4BqivVQgD+VdJekM3kdzOxOM1s0s8WTJ09G2iwAAKiqFNeLSTFnYNAVLmjM7CZJz7j74Xb93H2/u0+4+8Tll19edLMAAKDiUlwvJsWcgUEX4x2a6yXdbGbfkPSgpBvM7L9GiAsAABKW4noxKeYMDLrCBY27z7j7HnffJ+k2SY+6+5sLZwYAAJKW4noxKeYMDDoW1gQAAKWYGh/VPbdeqwtray83RkeGdc+t11b2gQBSmjkDgy7qwpru/neS/i5mTAAAkK4U14tJMWdgkPEODQAAAIBkRX2HBgCAKmKhxPKFHOO8vowTgG5Q0AAA+hoLJZYvZAHNvL6LR0/p4OF6UotwAqgGPnIGAOhrLJRYvpBjnNf3wKHjjBOArlDQAAD6Ggslli/kGOf1bboHxQaAdRQ0AIC+xkKJ5Qs5xnl9a2ZBsQFgHQUNAKCvsVBi+UKOcV7f26/byzgB6AoFDQCgr7FQYvlCjnFe3z+cupZxAtAVnnIGAOh7LJRYvpBjnNeXcQLQDQoaAAA2YC2UtMVYDwdAWihoAABoYc2atMVYDyerL4Bq4zs0AAC0sGZN2mKsh8NYA+mhoAEAoIU1a9IWYz0cxhpIDwUNAAAtrFmTthjr4TDWQHooaAAAaGHNmrTFWA+HsQbSQ0EDAEALa9akLcZ6OIw1kB6ecgYAwAashZK2GOvhAEgLBU1CynxeftHYPMsfeVI7N/LyzWqXFLRvMWKUtX87HQPYSVW+fwJ5OLc6R0GTiDKfl180Ns/yR57Uzo28fBePntLBw/VN7dMPPS6ZtNr0TX2l7H3Lih0ao6z9C9leamMKVPn+CeTh3ArDd2gSUebz8ovG5ln+yJPauZGX74FDx89rXz3jZwuRjX3z9i0rdmiMomKMR2pjClT5/gnk4dwKU7igMbO9ZvZpM3vSzL5kZu+KkRg2K/N5+UVj8yx/5Ent3MjLq+me2R4SI2Sfyzo+ZeZW1TEFqnz/BPJwboWJ8Q7Nc5J+x91fKuk1kn7DzF4WIS42KPN5+UVj8yx/5Ent3MjLq2ZWOEbIPpd1fMrMrapjClT5/gnk4dwKU7igcfdvuvvnW//+gaQnJfHhvsjKfF5+0dg8yx95Ujs38vK9/bq957UPXWAaqtl5ffP2LSt2aIyiYoxHamMKVPn+CeTh3AoT9aEAZrZP0rikQxm/u1PSnZJ0xRVXxNzsQFj/AthdDz+h080zGo34tIuiscvMDWlL7dxol+/ElZee157XNyR2SIwy928nYwA7qcr3TyAP51aYaAWNmT1f0kFJv+3u39/6e3ffL2m/JE1MTHT+gXScVebz8ovG5ln+yJPauZGXb157yL7FiFFUjPFIbUyBKt8/gTycW52LUtCY2ZDWipkH3H0uRkwAAACJ9TgAtFe4oDEzk/QBSU+6+x8XTwkAAGBNu/U4AECK85Sz6yW9RdINZvaF1n+/FCEuAAAYcKzHAWA7hd+hcff/IanzZ5oCAAB0qN16HHsu4RG2AOK8QwMAAFAK1uMAsB0KGgAAUFmsxwFgOxQ0AACgsqbGR3XPrdfqwtraS5bRkWHdc+u1POUMwFlRF9YEAACIjfU4ALTDOzQAAAAAksU7NAAAJCLGApMhMfL6stBl78QYv53MoUqx0b8oaAAASEC7BSY7fcEXskhlXt/Fo6d08HCdhS57IOQciHG+FM2hSrHR3/jIGQAACYixwGRIjLy+Bw4dZ6HLHokxfkXHqcyFTllEFd2ioAEAIAHtFpgsI0Ze36Z74TzQnRjjV3Scyopbdmz0NwoaAAASEGOByZAYeX1rZoXzQHdijF/RcSpzoVMWUUW3KGgAAEhAjAUmQ2Lk9b39ur0sdNkjMcav6DiVudApi6iiWxQ0AAAkIMYCkyEx8vr+4dS1LHTZIzHGr+g4lbnQKYuools85QwAgETEWGAyJEZeXxa67J0Y47eTOVQpNvoXBQ0AAMAOyFpjRRLrruwA1rfpbxQ0AAAAJctaY2X6occlk1abfraNdVfiY32b/sd3aAAAAEqWtcbK6hk/W8ysY92V+Fjfpv9R0AAAAJQsxnpB6A7r2/Q/ChoAAICSxVgvCN1hfZv+R0EDAABQsqw1VoYuMA3VNi9Uyror8bG+Tf/joQAAAAAlW//y+V0PP6HTzTMa3fCUs61tfFE9rrxjz3HuHxQ0AAAAOyBvjRXWXSkf69v0tygFjZndKOl9kmqS7nP3e2PELVM/PY98p/eln45dP8gbj5B2qfN1EGKMf2jORWOE7F8MZc2R0LhF8yhzrMvaXlVwXR5szJ3BkeI4VSWPmAoXNGZWk/Tnkn5B0tOSPmdmj7j7l4vGLks/PY+83b70w/bQXt54LB49pYOH6x21h6yDEGPuhOYckkfR/YuhrOtLaNyic7XMsc6KMSjX5ZTOOXSHuTM4UhynquQRm7n79r3aBTB7raR/5+6TrZ9nJMnd78n7m4mJCV9cXCy03SKuv/dR1VcaetsTH9XV36ufbb9oV03jV4zoG9/9oSRp349dvOnvvvzN70uSXvbiF7Rtk5QZIy9uSPvW7S0dW9Gzz21+tvr6vlw0dEFmbnk5Z7VvzaGb7YXsXxX6prS9vPEwM2XN7bz2LFnzod34dzp3QnO+aFdNl1w8tCluaIxO9y8r37z9yGvb7hiFXEc6OW55cYvO1W72o0jO3W4vZO50e+y36xtyndx6LsfYj262V/RYxIjRr9uLcQ0oe+6U8XqkXd/Q7RXNOSRukZxD7okx7p8hueW1dZLHyu6rdMt9f3Ren14ws8PuPrFdvxgfORuVdHzDz09Lui4joTsl3SlJV1xxRYTNdi/vuePrA/yj0+cPtCQ978JaR215MfLihrRv3V7WSbnevn4j2y5Gu/atOXSzvZD9q0LflLaXNx55L+pD/g+MrPnQbvy39t1o47kVmvOzzzX1o9MXnNcWEiMvrtTZPGvXHjInO42R1R4at+hc7WY/iuTc7fZC5k63x367viHXya3ncl6MvLasPLrZXtFjESNGv24vxjWg7LlTxuuRdn1Dt1c055C4oTG6vSfGuH+G5JbX1kkep374bGafKovxDs0bJE26+6+3fn6LpFe7+zvz/qYq79BsNToyrP/5nhv0xv/8vyQV+9JYVoy8uKHtG223L0VtzaGb7YXsXxX6prS9vPGomamZMbfz2rNkzYcYcyc059GRYe25ZHhT3NAYWWLO943KmpOhcYvO1Rj7ERKj2+2FzJ2yhFwnt57LeTHy2rJ0sz2UJ4W5U8brkXZ9y4obur3YOYfcE2O+9ixyLHbiNXBMnb5DE2Mdmqcl7d3w8x5JJyLELU0/PY98p/eln45dP8gbj9uv29txe8g6CDHGPzTnkDyK7l8MZc2R0LhF8yhzrMs6t6qC6/JgY+4MjhTHqSp5xBbjI2efk3SNmV0lqS7pNkn/MkLc0qx/6akfnvCw0/vST8euH7Qbj4krL+24PS9GyPbKzHn9UZvdxiiac4iy5kho3KJ5lD3WZWyvKtrty9ZzOcXtoT3mzuBIcZyqkkdshQsad3/OzH5T0oLWHtt8v7t/qXBmJZsaH01+8Nbt9L7007HrB3nj0U17ke2FCM0tRox+mCOhcYvmUeZYl7W9quC6PNiYO4MjxXGqSh4xRVmHxt3/RtLfxIgFAAAAAJ2K8R0aAAAAAOgJChoAAAAAyaKgAQAAAJAsChoAAAAAyaKgAQAAAJAsChoAAAAAyaKgAQAAAJAsChoAAAAAyaKgQa75pbqWjq3o0FOndP29j2p+qV5qjLy+MfJAd2KM307mUKXYAJCSGPfgrL7c2+PiuGXb1esEUE3zS3XNzB3R6eYZSVJ9paGZuSOSpKnx0egx8vouHj2lg4frmTFQrhjjl9W3rByqFBsAUtLNPbiT+8D0Q49LJq02veu4OIf7Vj7eoUGm2YVlNVabm9oaq03NLiyXEiOv74FDxwvnge7EGL+i41RW3LJjA0BKYtyDs2KsnvGzxUy3cXEO9618FDTIdGKlEdReNEZe36Z7ZntIHuhOjPErOk5lxS07NgCkJMY9OOTayb29O9y38lHQINPukeGg9qIx8vrWzArnge7EGL+i41RW3LJjA0BKYtyDQ66d3Nu7w30rHwUNMk1Pjml4qLapbXiopunJsVJi5PW9/bq9hfNAd2KMX9FxKitu2bEBICUx7sFZMYYuMA3VNhcv3Nu7x30rHw8FQKb1L5fNLizrxEpDu0eGNT05FvSls5AY7fpOXHlpZvuBx45F2FPkiTV+O5VDlWIDQEq6uQd3GqNoXJzDfSsfBQ1yTY2PRnlB2mmMvL4x8kB3YozfTuZQpdgAkJIY9+B2MYrExTkct2x85AwAAABAsihoAAAAACSLggYAAABAsihoAAAAACSLggYAAABAsgoVNGY2a2ZfMbMnzOwjZjYSKzEAAAAA2E7Rd2g+Ienl7v4KSV+VNFM8JQAAAADoTKGCxt3/1t2fa/34WUl7iqcEAAAAAJ2J+R2aX5X0sbxfmtmdZrZoZosnT56MuNm0zS/VtXRsRYeeOqXr731U80v1Xqc08ELGJKtv3t/HGGvOF5SFc6s7ZR03rhcoS4rnRch9NXT/YsQoa/92OkbKdm3Xwcw+KelFGb+6290/2upzt6TnJD2QF8fd90vaL0kTExPeVbZ9Zn6prpm5IzrdPCNJqq80NDN3RFL2qrooX7sx6aTv9EOPSyatNn3T3y8ePaWDh+uFxprzBWXh3OpOWcctRlzGFFlSPC/ycs66r+bdg6Xs/Qu5j+fFKGv/QrYX8tqlX237Do27v97dX57x33oxc4ekmyS9yd0pVALMLiyrsdrc1NZYbWp2YblHGSFkTLL6rp7xsxfBjX9/4NDxwmPN+YKycG51p6zjFiMuY4osKZ4XeTln3Vfz7sF5+xdyHy/rGDHf4yj6lLMbJf2upJvd/UdxUhocJ1YaQe0oX8iYhIxTM6fWD4nB+YKycG51p6zjFiMuY4osKZ4Xebnl3VdDYsSYU0Ux3+Mo+h2aP5P0DyV9wsy+YGb/KUJOA2P3yHBQO8oXMiYh41QzC9peSF/OFxTFudWdso5bjLiMKbKkeF7k5ZZ3Xw2JEWNOFcV8j6PoU87+sbvvdfdXtv57e6zEBsH05JiGh2qb2oaHapqeHOtRRggZk6y+QxeYhmqbL7LDQzXdft3ewmPN+YKycG51p6zjFiMuY4osKZ4XeTln3Vfz7sF5+xdyHy/rGDHf49j2oQAoz/qXvWYXlnVipaHdI8Oanhyr7BfzBkG7MTnw2LGO+ub9/cSVlxYaa84XlIVzqztlHbcYcRlTZEnxvGiXc9Z9Na9vSOyQGGXuX4wYW1+79CsKmh6bGh+t9IVkEIWMSV7fvLYYL3Q4X1AGzq3ulHXcuF6gLCmeF+3utZ3eg7uJvVOY78XFXIcGAAAAAHYUBQ0AAACAZFHQAAAAAEgWBQ0AAACAZFHQAAAAAEgWBQ0AAACAZFHQAAAAAEgWBQ0AAACAZFHQlGB+qa6lYys69NQpXX/vo5pfqvc6pb4Tcozz+jJOQP9hXgPA4NnV6wT6zfxSXTNzR3S6eUaSVF9paGbuSI+z6i/tjvHWVXLz+i4ePaWDh+sdxQCQhpBrAwCgf/AOTWSzC8tqrDY3tTVWm5pdWO5RRv0n5Bjn9T1w6DjjBPQZrr8AMJgoaCI7sdIIake4kGOc17fpHhQbQPVx/QWAwURBE9nukeGgdoQLOcZ5fWtmQbEBVB/XXwAYTBQ0kU1Pjml4qLapbXiopunJsR5l1H9CjnFe39uv28s4AX2G6y8ADCYeChDZ+hdPZxeWdWKlod0jw5qeHNPU+KgOPHasx9n1h3bHOKTvxJWXdhQDQBpCrg0AgP5BQVOCqfFRbqAlCznGeX0ZJ6D/MK8BYPDwkTMAAAAAyaKgAQAAAJCsKAWNmb3bzNzMLosRDwAAAAA6UbigMbO9kn5BEt94BwAAALCjYrxD8yeS7pKUvVIhAAAAAJSkUEFjZjdLqrv74x30vdPMFs1s8eTJk0U2CwAAAACSOnhss5l9UtKLMn51t6Tfk/QvOtmQu++XtF+SJiYmeDcHAAAAQGHbFjTu/vqsdjO7VtJVkh43M0naI+nzZvZqd/9W1CwBAAAAIEPXC2u6+xFJL1z/2cy+IWnC3b8TIS8AAAAA2Bbr0FTU/FJdS8dWdOipU7r+3kc1v1TvdUpRlbV/qcUtOzaA6sib61wDgP6T4rxOMed1Xb9Ds5W774sVa9DNL9U1M3dEp5tnJEn1lYZm5o5IkqbGR3uZWhRl7V9qcbeLDaB/5M31xaOndPBwvW+v98AgSvF1XIo5b8Q7NBU0u7CsxmpzU1tjtanZheUeZRRXWfuXWtyyYwOojry5fuDQca4BQJ9J8d6eYs4bUdBU0ImVRlB7asrav9Tilh0bQHXkzemmZz/0k2sAkK4U7+0p5rwRBU0F7R4ZDmpPTVn7l1rcsmMDqI68OV1be0pox/0BVF+K9/YUc96IgqaCpifHNDxU29Q2PFTT9ORYjzKKq6z9Sy1u2bEBVEfeXL/9ur1cA4A+k+K9PcWcN4r2UADEs/7lq9mFZZ1YaWj3yLCmJ8eS+FJWJ8rav9Tibhf7wGPHCscHUA3t5vrElZf27fUeGEQpvo5LMeeNKGgqamp8NJmTqBtl7V9qccuODaA68uY61wCg/6Q4r1PMeR0fOQMAAACQLAoaAAAAAMmioAEAAACQLAoaAAAAAMmioAEAAACQLAoaAAAAAMmioAEAAACQLAoaAAAAAMmioAEqaH6prqVjKzr01Cldf++jml+q9zoloHKYJwB2Etec6qKgASpmfqmumbkjOt08I0mqrzQ0M3eECyewAfMEwE7imlNtFDRAxcwuLKux2tzU1lhtanZhuUcZAdXDPAGwk7jmVBsFDVAxJ1YaQe3AIGKeANhJXHOqjYIGqJjdI8NB7cAgYp4A2Elcc6qNggaomOnJMQ0P1Ta1DQ/VND051qOMgOphngDYSVxzqm1XrxMAsNnU+Kiktc/rnlhpaPfIsKYnx862A2CeANhZXHOqrXBBY2bvlPSbkp6T9N/c/a7CWQEDbmp8lIsksA3mCYCdxDWnugoVNGb2c5JukfQKd3/WzF4YJy0AAAAA2F7R79C8Q9K97v6sJLn7M8VTAgAAAIDOFC1ofkLSz5jZITP772b2U3kdzexOM1s0s8WTJ08W3CwAAAAAdPCRMzP7pKQXZfzq7tbfXyLpNZJ+StJfmdnV7u5bO7v7fkn7JWliYuK83wMAAABAqG0LGnd/fd7vzOwdkuZaBcxjZnZG0mWSeAsGAAAAQOmKfuRsXtINkmRmPyHpQknfKZoUAAAAAHSi6GOb75d0v5l9UdJpSXdkfdwMAAAAAMpQqKBx99OS3hwpFwAAAAAIUvQjZwAAAADQMxQ0W8wv1bV0bEWHnjql6+99VPNL9UrHBQAAALYapNeeFDQbzC/VNTN3RKebZyRJ9ZWGZuaOFD4ByooLAAAAbDVorz0paDaYXVhWY7W5qa2x2tTswnIl4wIAAABbDdprTwqaDU6sNILaex0XAAAA2GrQXntS0Gywe2Q4qL3XcQEAAICtBu21JwXNBtOTYxoeqm1qGx6qaXpyrJJxAQAAgK0G7bVn0YU1+8rU+Kiktc8dnlhpaPfIsKYnx862Vy0uAAAAsNWgvfY0d9/xjU5MTPji4uKObxcAAABAGszssLtPbNePj5wBAAAASBYFDQAAAIBkUdAAAAAASBYFDQAAAIBkUdAAAAAASFZPnnJmZiclHd3xDWe7TNJ3ep0Eusb4pY3xSxvjlzbGL22MX7oYu85d6e6Xb9epJwVNlZjZYiePg0M1MX5pY/zSxviljfFLG+OXLsYuPj5yBgAAACBZFDQAAAAAkkVBI+3vdQIohPFLG+OXNsYvbYxf2hi/dDF2kQ38d2gAAAAApIt3aAAAAAAki4IGAAAAQLIGuqAxsxvNbNnMvmZm7+l1PmjPzPaa2afN7Ekz+5KZvavVfqmZfcLM/nfrfy/pda7IZmY1M1sys79u/XyVmR1qjd2HzezCXueIbGY2YmYPm9lXWnPwtcy9dJjZv2ldN79oZgfM7B8w/6rLzO43s2fM7Isb2jLnm635j63XMk+Y2at6lzmk3PGbbV0/nzCzj5jZyIbfzbTGb9nMJnuTddoGtqAxs5qkP5f0i5JeJul2M3tZb7PCNp6T9Dvu/lJJr5H0G60xe4+kT7n7NZI+1foZ1fQuSU9u+Pk/SPqT1tj9H0m/1pOs0In3Sfq4u79E0k9qbRyZewkws1FJvyVpwt1fLqkm6TYx/6rsQ5Ju3NKWN99+UdI1rf/ulPT+HcoR+T6k88fvE5Je7u6vkPRVSTOS1Hodc5ukf9L6m79ovUZFgIEtaCS9WtLX3P3r7n5a0oOSbulxTmjD3b/p7p9v/fsHWntBNaq1cfvLVre/lDTVmwzRjpntkfTLku5r/WySbpBygzuzAAADPklEQVT0cKsLY1dRZvYCST8r6QOS5O6n3X1FzL2U7JI0bGa7JD1P0jfF/Kssd/97Sae2NOfNt1sk/Rdf81lJI2b24p3JFFmyxs/d/9bdn2v9+FlJe1r/vkXSg+7+rLs/JelrWnuNigCDXNCMSjq+4eenW21IgJntkzQu6ZCkf+Tu35TWih5JL+xdZmjjTyXdJelM6+cfk7Sy4QLPHKyuqyWdlPTB1kcG7zOzi8XcS4K71yX9kaRjWitkvifpsJh/qcmbb7yeSc+vSvpY69+MXwSDXNBYRhvPsE6AmT1f0kFJv+3u3+91Ptiemd0k6Rl3P7yxOaMrc7Cadkl6laT3u/u4pB+Kj5clo/Vdi1skXSVpt6SLtfYxpa2Yf2niWpoQM7tbax+hf2C9KaMb4xdokAuapyXt3fDzHkknepQLOmRmQ1orZh5w97lW87fX315v/e8zvcoPua6XdLOZfUNrH++8QWvv2Iy0PgIjMQer7GlJT7v7odbPD2utwGHupeH1kp5y95PuvippTtJPi/mXmrz5xuuZRJjZHZJukvQmP7cQJOMXwSAXNJ+TdE3rKS8Xau0LWY/0OCe00frOxQckPenuf7zhV49IuqP17zskfXSnc0N77j7j7nvcfZ/W5tqj7v4mSZ+W9CutboxdRbn7tyQdN7OxVtPPS/qymHupOCbpNWb2vNZ1dH38mH9pyZtvj0j6V62nnb1G0vfWP5qG6jCzGyX9rqSb3f1HG371iKTbzOwiM7tKaw93eKwXOabMzhWIg8fMfklr/y9xTdL97v7eHqeENszsn0n6jKQjOvc9jN/T2vdo/krSFVq7cb/B3bd+mRIVYWavk/Rud7/JzK7W2js2l0pakvRmd3+2l/khm5m9UmsPdLhQ0tclvVVr/6cYcy8BZvbvJb1Rax91WZL061r7nD7zr4LM7ICk10m6TNK3Jf2+pHllzLdWkfpnWntC1o8kvdXdF3uRN9bkjN+MpIskfbfV7bPu/vZW/7u19r2a57T2cfqPbY2J9ga6oAEAAACQtkH+yBkAAACAxFHQAAAAAEgWBQ0AAACAZFHQAAAAAEgWBQ0AAACAZFHQAAAAAEgWBQ0AAACAZP1/+Zw0blAD2NoAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2233a6caef0>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAEDCAYAAAD0sNGCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+0XWV95/HPh5DgrdQGJCq5IQ22kZZKNe0dRsu0uhQG7HQRhqUV+ws72rTTsT+mFSVjl+3Y5YKazvTHlLqaQQt2OqAyNGY62BTBTru6hHKZKAhMCgUL+SGJYLRKhkDynT/OPuTmcM6595y9z9nPs/f7tVbWPWff5+793Xs/+9nnyfPdz3FECAAAAACa7IS6AwAAAACASaPjAwAAAKDx6PgAAAAAaDw6PgAAAAAaj44PAAAAgMaj4wMAAACg8ZLv+Nj+qO39tr9Y0fr+wvZB23/es/w624/Y/nzx79VVbA8AAABA/ZLv+Ei6TtJFFa5vi6SfHPC7KyLi1cW/z1e4TQAAAAA1Sr7jExF/LenJhctsf0cxcnO37b+x/V0jrO82Sf9UdZwAAAAA0pV8x2eArZJ+ISK+X9K7Jf1hRev9oO17bP+O7ZMqWicAAACAmp1YdwCjsn2ypB+Q9Enb3cUnFb+7VNIH+vzZnoi4cJFVb5b0ZUkr1OlYvXfAugAAAABkJruOjzqjVAcj4nmTD0TEzZJuHmelEbGvePm07T9WZyQJAAAAQANkl+oWEV+X9Ijtt0iSO15Vdr22T++uT9IlkiqZRQ4AAABA/RwRdccwlO0bJL1e0mmSHpf065Jul/RhSadLWi7pxohYUlqa7b+R9F2STpb0hKR3RMQO27dLWiXJkj4v6eci4hvV7g0AAACAOiTf8QEAAACAsrJLdQMAAACAUSU7ucFpp50W69atqzsMAAAAAAm7++67vxIRqxYrl2zHZ926dZqfn687DAAAAAAJs/2PSylHqhsAAACAxqPjAwAAAKDx6PgAAAAAaDw6PgAAAAAaj44PAAAAgMZLdlY3AAAAdGzbuUdbduzS3oOHtHrljK648CxdsmG27rCArFQy4mP7Itu7bD9k+8o+v19r+7O2d9q+x/YPV7FdAACAptu2c48233yv9hw8pJC05+Ahbb75Xm3buafu0ICslO742F4m6RpJb5J0tqS32T67p9ivSfpERGyQdJmkPyy7XQAAgDbYsmOXDj1z5Lhlh545oi07dtUUEZCnKkZ8zpX0UEQ8HBGHJd0oaWNPmZD0ouL1t0naW8F2AQAAGm/vwUMjLQfQXxUdn1lJjy14v7tYttBvSPoJ27sl3SLpF/qtyPYm2/O25w8cOFBBaAAAAHlbvXJmpOUA+qui4+M+y6Ln/dskXRcRayT9sKQ/sf28bUfE1oiYi4i5VatWVRAaAABA3q648CzNLF923LKZ5ct0xYVn1RQRkKcqZnXbLemMBe/X6PmpbO+QdJEkRcTnbL9A0mmS9lewfQAAgMbqzt72npvu0eEjRzXLrG7AWKro+Nwlab3tMyXtUWfygh/rKfOopDdKus72d0t6gSRy2QAAAJbgkg2zuuHvHpUkffxnX1tzNECeSqe6RcSzkt4laYekB9SZve0+2x+wfXFR7Fcl/YztL0i6QdLbI6I3HQ4AAAAAJqKSLzCNiFvUmbRg4bL3L3h9v6TzqtgWAAAAAIyqko4PAAAAsNC2nXu0Zccu7T14SKt5LgkJoOMDAACASm3buUebb773uS9e3XPwkDbffK8k0flBbaqYzhoAAAB4zpYdu57r9HQdeuaItuzYVVNEAB0fAAAAVGzvwUMjLQemgY4PAAAAKrV65cxIy4FpoOMDAACASl1x4VmaWb7suGUzy5fpigvPqikigMkNAAAAULHuBAbvuekeHT5yVLPM6oYE0PFB6+Q4vWaOMQMA2u2SDbO64e8elSR9/GdfK4n7GepFxwetkuP0mjnGDABAL+5nqBvP+KBVcpxeM8eYAQDoxf0MdaPjg1bJcXrNHGMGAKAX9zPUjY4PWiXH6TVzjBkAgF7cz1A3Oj5olRyn18wxZgAAenE/Q92Y3ACtkuP0mjnGDABAL+5nqBsdH7ROv+k1U5djzACA5ik7HfUk72dMlY3FVJLqZvsi27tsP2T7ygFlftT2/bbvs/3fq9guAAAApqM7HfWeg4cUOjYd9bade+oOLenYkI7SHR/byyRdI+lNks6W9DbbZ/eUWS9ps6TzIuJ7JP1y2e0CAABgelKejjrl2JCOKkZ8zpX0UEQ8HBGHJd0oaWNPmZ+RdE1EfFWSImJ/BdsFAADAlKQ8HXXKsSEdVXR8ZiU9tuD97mLZQq+Q9Arbf2v7DtsX9VuR7U22523PHzhwoILQAAAAUIWUp6NOOTako4qOj/ssi573J0paL+n1kt4m6VrbK5/3RxFbI2IuIuZWrVpVQWgAAACoQsrTUaccG9JRxaxuuyWdseD9Gkl7+5S5IyKekfSI7V3qdITuqmD7AAAAmLCUp6NOOTako4qOz12S1ts+U9IeSZdJ+rGeMtvUGem5zvZp6qS+PVzBtgEAAFpr2lM4p/z1CmVjYzrs5ivd8YmIZ22/S9IOScskfTQi7rP9AUnzEbG9+N2/tH2/pCOSroiIJ8puGwAAoK26Uzh3ZzPrTuEsiQ/sI+JYtkMl3+MTEbdExCsi4jsi4oPFsvcXnR5Fx69ExNkRcU5E3FjFdgEAANqKKZyrw7Fsh0o6PgAAAJgupnCuDseyHej4AAAAZIgpnKvDsWwHOj4AAAAZYgrn6nAs26GKWd0AAAAwZUzhXB2OZTvQ8QEAAMhUytNLjyKFqaRzPJYpHLec0PEBAABAbZhKejwct9HxjA8AAABqw1TS4+G4jY4Rn8QwZAkAANqEqaTHw3EbHSM+CekOWe45eEihY0OW23buqTs0AACAiWAq6fFw3EZHxychDFkCAIC2YSrp8XDcRkeqW0IYsgQAAG3DVNLj4biNjo5PQlavnNGePp0chizrxXNXAAAMVsV9cpSppLkvH5PjFNx1ItUtIQxZpofnrgAAGGza90nuyyiDjk9CLtkwq6suPUcrlnVOy+zKGV116Tmt/V+MFPDcFQAAg037Psl9GWWQ6pYYhizTwnNXAKaNNB7kZNr3Se7LS0M70h8jPsAQTBUJYJpI40Fupn2f5L68ONqRwSrp+Ni+yPYu2w/ZvnJIuTfbDttzVWwXmDSeuwIwTaTxIDfTvk9yX14c7chgpVPdbC+TdI2kCyTtlnSX7e0RcX9PuW+V9IuS7iy7TbRDCsO0TBUJYJpI42mXFO5zZWOY9n2S+/LiaEcGq+IZn3MlPRQRD0uS7RslbZR0f0+535T0IUnvrmCbaLjuMG33fyy6w7SSaun88NwVgGngaw3aI4X7XFUxTPs+yX15ONqRwapIdZuV9NiC97uLZc+xvUHSGRHx58NWZHuT7Xnb8wcOHKggNOSKYVoAbUQaT3ukcJ9LIQZUj3ZksCo6Pu6zLJ77pX2CpN+R9KuLrSgitkbEXETMrVq1qoLQkCuGaQG0EV9r0B4p3OdSiAHVox0ZrIpUt92Szljwfo2kvQvef6ukV0r6K9uS9DJJ221fHBHzFWwfDcQwLYC2Io2nHVK4z6UQAyZjUu1ICs+llVHFiM9dktbbPtP2CkmXSdre/WVEfC0iTouIdRGxTtIdkuj0YCiGaQEATZbCfS6FGJCPJkyTXbrjExHPSnqXpB2SHpD0iYi4z/YHbF9cdv1oJ4ZpAQBNlsJ9LoUYkI8mPBNWRaqbIuIWSbf0LHv/gLKvr2KbaL4c0z1yHwIGAExPCulIOd5rUY8mPBNWSccHQBpTkwIA2o17ESalCc+EVfGMDwA1YwgYAJA37kWYlCY8E8aID1CRJgwBAwDyxr0Ik9IdMXzPTffo8JGjms0wpZ+OT02a9CxIk/aljCYMAQMA8sa9KE1N+ayU+zNhpLrVoAnTAXY1aV/KasIQMAAgb9yL0sNnpXTQ8alBk/Jvm7QvZTEtKACgbtyL0sNnpXSQ6laDqvJvUxg2bVIucRXHc9pDwP1illR7vcB0pNAGpI5jhDbKPR2paZr0WSl3dHxqUEX+bSrTVTYllziV4zmKfjFf8ckvSJaeORLPLUt9PzCeHOvstHGMAKSgKZ+VmoBUtxpUkX+byrBpU3KJUzmeo+gX8zNH47lOT1fq+4Hx5Fhnp41jBCAFTfms1ASM+NSgiukAUxk2TX1qw6WmuaRyPEcxSmwp7wfGk2OdnTaOEYAU0l1T/6zUJnR8alI2/zalYdNUc4mHpbn0Sul4LtWgmAeVRbPkWGenjWMEtFtK6a6pflZqG1LdMsWw6eJGSXPJ8Xj2i3n5CdbyZT5uWer7gfHkWGenjWMEtBvprujFiE+mGDZd3LA0lzWnHP8/vjkez0Ex91uW8n5gPDnW2WnjGAHtRroretHxyVhThk0nlX87appLjsdzUMy57QfGk2OdnbZUjlHZdi6F5xSA3JDuejzaEVLdULNJfpsxaS4AUlC2neNb34Hx8DngGNqRjko6PrYvsr3L9kO2r+zz+1+xfb/te2zfZvvbq9gu8jfJ/Fu+vRpACsq2czynAIyHzwHH0I50lE51s71M0jWSLpC0W9JdtrdHxP0Liu2UNBcRT9n+t5I+JOmtZbeN/E06/zaVNBcA7VW2neM5BWB8VXwOaEKK2DjtSBP2u1cVIz7nSnooIh6OiMOSbpS0cWGBiPhsRDxVvL1D0poKtosGGJRn29b8WwDNU7ado50E6tOUFLFR25Gm7HevKjo+s5IeW/B+d7FskHdI+nQF20UDkH8LoOnKtnO0k0B9mpIiNmo70pT97lXFrG7usyz6FrR/QtKcpNcN+P0mSZskae3atRWEhtQx3SyApivbztFOAvVpSqrpqO3IKF8JkpMqOj67JZ2x4P0aSXt7C9k+X9L7JL0uIp7ut6KI2CppqyTNzc317Tylool5j3Uh//Z4TdoXAB1l2zmeV2wP7gFpSWlK7LJ1Y5R2JKX9rlIVqW53SVpv+0zbKyRdJmn7wgK2N0j6I0kXR8T+CrZZq6bmPeaqSeejSfsCABgN94D0pJJqOu26kcp+V610xycinpX0Lkk7JD0g6RMRcZ/tD9i+uCi2RdLJkj5p+/O2tw9YXRaamveYqyadjybtCwBgNNwD0pPKlNjTrhup7HfVqkh1U0TcIumWnmXvX/D6/Cq2k4qm5HuOKtXh9yadj6bm1LZdqtcOgLSkdD9Lod1KIQYpjVTTOj4fpLDfVavkC0zbpo1Ti6Y8/N6k89GkfUFHytcOgLSkcg9Iod1KIYaUpFI3ckfHZwxNzXscJuXh9yadjybtCzpSvnYApCWVe0AK7VYKMaQklbqRu0pS3domx6lFyw4XpzT83ivH8zHIsH3pDjcjLylfOwDSksr9LIV2K4UY6jLsMxufD8qh4zOmnPIeu8PF3f856Q4XS1pyY5r6tIY5nY/FNGlfkP61AyAtKdwDUmi3UoihDot9Zqu7buSOVLcWqGK4mCFWYDxcOwByk0K7lUIMdSDFb7IY8WmBKoaLUxl+B3LDtQMgNym0WynEUIc2p/hNAx2fFqhquJghVmA8qVw7qUwNCyB9KbRbKcQwbW1N8ZsWUt1aoK3DxQCOYWpYAEgfn9kmi45PCzT123cBLB154wCQPj6zTRapbi2RynAxqTZAPVLJG0+lDUglDtSD84+UpfKZrYno+GBqhk3RCGCyUsgbr2Jq/SbFgXpw/oH2ItUNU0OqDVCfFPLGU2kDUokD9eD8A+1FxwdTk0qqDdBGKeSNp9IGpBIH6sH5B9qLVLdFpJIHnEocZaSQaoPpaUKdbZq688ZTaQNSiQP14PwD7cWIzxCpTP+aShxlpZBqg+loSp1FtVJpA1KJA/Xg/APtRcdniFTygFOJo6wUUm0wHU2ps6hWKm1AKnGgHpx/oL0qSXWzfZGk35O0TNK1EXF1z+9PkvQxSd8v6QlJb42IL1Wx7UlKJQ94WBxrTklzaH5QmlMVqTYppFBVEUMK+zEpqVw7TdKU+lJ3ul1qcUzTJOtQCvVzlBgGnf8U9gP5oL7kp3THx/YySddIukDSbkl32d4eEfcvKPYOSV+NiO+0fZmk35L01rLbnrRU8oBTiWOpJjlVaApTYlexf02fTjW3Opu6FOo98lZXuzyt9ox2GdNGfcmTI6LcCuzXSvqNiLiweL9ZkiLiqgVldhRlPmf7RElflrQqhmx8bm4u5ufnS8VWVrdS/9TdN+vlX+s8m3DCCdbLT3uhTjv5JN2/7+uSpLNPf9Fzf/OlJ74pSVr34hcet65+y/v9fb/lX/nG03r4K9/U0aPHDlc3jv3/9HTfdYyyvX5lB+3HUmLe+ehBPf3s8WlOknTSicu0Ye3KvutY6rEYtu6Tlp+wpHWMsr2q9q/3eE7yGA0rO8q5HmV577JhdfYbTz9b+famvX/jlC1zTsap96PEkdOxqHJ70zwWo26v7HXdu3yxNqfM/g1b9ykvXL6kdYyyvVFjGLR/o9xfRtmPUfev6nvUYmXLfB6pKuZJlR11HWWOxTj1pWwbMOr+VdGODFrHl1edoZ++4b8oFbbvjoi5xcpVkeo2K+mxBe93S/rng8pExLO2vybpxZK+srCQ7U2SNknS2rVrKwitnG6P/cHP/5mkTmU+49QZnXbySZKkb1mx7Hl/89Th518Eg5b3+/t+y7vbe+Qr39SRo3FcHN0PkWW216/soP1YSsz9GoKFy/utY6nHYti6uw3NUta91O31Wz7O/vUez0keo2FlRznXoyzvXTasznY761Vur6qyk9xemXMyTr0fJY6cjoXU6Vg//vWnFRH66jefea5ulWm3Fou5d/mgGIatY5Ttlb2ue5cv1uaUOdfD1v3U4f6PElddt8bZv1HuL6Psx6DlVZzrsveBQfV2UtubdtlR11Hms9k49aVsGzBo+aSOxbB1nPrCk/qWT10VIz5vkXRhRLyzeP+Tks6NiF9YUOa+oszu4v0/FGWeGLTeFEZ8xvHWP/qcpOfnjA9aPu04llq2TLznXX173zSn2ZUz+tsr3zDy+qa17knG0Hs869qPUc71KMvL1reqtjft/avi2lnKOqqoc6PGsNTYxik7Shy9y3rTS6TOjFxXXXpOJc/sLGX/hsVwyYbZSo591TEvVofKnOth6+4+hzrpujXO/vWqaj9G3b9pWazettGkPusMqi9Vx1ClVOIY11JHfKqY1W23pDMWvF8jae+gMkWq27dJerKCbSNBk5wqNIVpSKuIIYX9QD6oL8ekMGNgCjGMinZ5OutIWY71NmVNry9NVUXH5y5J622faXuFpMskbe8ps13S5cXrN0u6fdjzPchbd6rQ2ZUzsqqdKnSS655mDCnsB/JBfTkmhRkDU4hhVLTL01lHynKstylren1pqtLP+BTP7LxL0g51prP+aETcZ/sDkuYjYrukj0j6E9sPqTPSc1nZ7eZk28492vnoQR0+clTnXX17K6Y77E5dndu6JxHDoPOfwn4gH9SXjhRmDEwhhnEMqkNV3KNSqJ9VxJDCfkxKrvU2ZU2pL236nFrJF5hGxC0R8YqI+I6I+GCx7P1Fp0cR8f8i4i0R8Z0RcW5EPFzFdnPQzak9fOSoJL7Bvm04/0C1UkgvSSGGqtBGtUeT6i2q07Y2oJKODwYjp7bdOP9AtVJIL0khhqrQRrVHk+otqtO2NqCK6awxBDm17cb5B6qXQnpJCjFUYVgb1Z2ZCs3RlHrbFCmkmLXtcwojPhM2KHeWnNp24PwDSBltFFCPVFLM2tYG0PGZMHJq243zDyBltFFAPVJJMWtbG0Cq24R1hyy37NilvQcPafXKmUbPloHjcf4BpGxYG9X9QlgA1Uslxaxtn1Po+EwBObXtxvlHylLIMW+rVI49bRQwfSlNL96mNoBUNwBoqVRyzNuIYw+0W9tSzFJBxwcAWiqVHPM24tgD7cb04vUg1Q1A7VJJ+UkljmlJJce8jao69m2rs0BVUrh22pRilgo6PgBqNSjlR9JUbwjD4miqlHLM26aKY5/KtQPkhmunvUh1A1CrVFJ+Uoljmsgxr08Vx76NdRaoAtdOezHiA6BWqaRbtfEb7Ns2jWlKqjj2qVw7QG64dtqLjg/QYinkOKeSbpVKHNM2KMc8hboxitzilcrn96deZ3M8J/00ZT9wTOrXDiaHVDegpVKZTjeVdKtU4khBKnVjqXKLtyop19mmnJOm7AeOl/K1g8mi4wO0VCo5zqlM6ZlKHClIpW4sVW7xViXlOtuUc9KU/cDxUr52MFmlUt1snyrp45LWSfqSpB+NiK/2lHm1pA9LepGkI5I+GBEfL7NdAOWllOOcypSeKcSRQlpNHc87ldnvlOrytKVQZ/tpyjlpyn7g+VK9djBZZUd8rpR0W0Ssl3Rb8b7XU5J+KiK+R9JFkn7X9sqS2wVQ0qBcZnKc65NKWs2060bZ/aYup6cp56Qp+wGgo2zHZ6Ok64vX10u6pLdARPx9RDxYvN4rab+kVSW3C6AkcpzTk0pazbTrRtn9pi6npynnpCn7AaCjbMfnpRGxT5KKny8ZVtj2uZJWSPqHktvFEN2UkTsfeVLnXX37yP9bXPbvkQdynNOTSlpNFXVjlHak7H5Tl6tX9j7QlHPSlP0A0LHoMz62PyPpZX1+9b5RNmT7dEl/IunyiDg6oMwmSZskae3ataOsHoWy30bMtxm3CznOaUlpitUydWPUdqSK/aYuV6eq+0BTzklT9gPAEkZ8IuL8iHhln3+fkvR40aHpdmz291uH7RdJ+l+Sfi0i7hiyra0RMRcRc6tWkQ03jrIpI6mk2gBt1JS0mlHbkabsd1NwHwDQVGVT3bZLurx4fbmkT/UWsL1C0p9J+lhEfLLk9rCIsikjqaTaAG3UlLSaUduRpux3U3AfANBUpaazlnS1pE/YfoekRyW9RZJsz0n6uYh4p6QflfRDkl5s++3F3709Ij5fctvoo2zKSEqpNkBuqpiKuglpNeO0I03Y76aY9H0ghSnbAbRTqRGfiHgiIt4YEeuLn08Wy+eLTo8i4r9FxPKIePWCf3R6JqRsyggpJ8B4UpmKOgW0I3mb5PnjOgFQp7KpbkhM2ZQRUk6A8fBcxDG0I3mb5PnjOgFQp7KpblggleH7sikjpJwAo+O5iOO1tR1J5T5Q1qTO37DrZM0pz0+ly/F49otZUnb7gfHkWGfbhI5PRZgGGmg3no8D94HFjXKdDDueqeoX8xWf/IJk6Zkj8dwy6kUz0Qakj1S3ijB8D7Qbz7WA+8DiRrlOcjye/WJ+5mg81+npSn0/MJ4c62zb0PGpCGkuQLvxXAu4DyxulOskx+M5Smwp7wfGk2OdbRtS3SpCmguAVJ5rIcd8PGWPW0r3gZTrwFKvk5SO51INinlQWTRLjnW2bRjxqQhpLgBSwHTB46niuKVyH2hKHUjleI6iX8zLT7CWL/Nxy1LfD4wnxzrbNnR8KkKaC4AUkGM+niqOWyr3gabUgVSO5yj6xbzlLa/Slje/Kqv9wHhyrLNtQ6pbhVJJcwGmKeWUmjYix3w8VR23FO4DTaoDVRzPabdRg2Kuu15gOlJoAzAYIz4AxtaUlJomGZRLTo75cE06bk3al7JoowAsRMcHwNiaklLTJOSYj6dJx61J+1IWbRSAhUh1AzC2JqXUNEU3xWLLjl3ae/CQVq+cyTb9cJopSk06bk3al7JoowAsRMcHwNiYujNNTcgxr+Mb0Jtw3LqatC9l0EYBWIhUNwBjI6UGk0KKEqpAGwVgIUZ8AIyNlBpMCilKqAJtFICF6PgALTGp5yVGSalh6mssFSlKacrxGibtD0BXqVQ326favtX2g8XPU4aUfZHtPbb/oMw2AYwuhSldU4gB+SBFKT1cwwByV/YZnysl3RYR6yXdVrwf5Dcl/e+S2wMwhhSel0ghBuSDb0BPD9cwgNyVTXXbKOn1xevrJf2VpPf2FrL9/ZJeKukvJM2V3CYqlGPaAkaXwvMSKcSAyUghjRKTxzUMIHdlR3xeGhH7JKn4+ZLeArZPkPSfJF2x2Mpsb7I9b3v+wIEDJUPDYkhbaI8Uvsk9hRhQPdqR9uAaBpC7RTs+tj9j+4t9/m1c4jZ+XtItEfHYYgUjYmtEzEXE3KpVq5a4eoyLtIX2SOF5iRRiQPVoR9qDaxhA7hZNdYuI8wf9zvbjtk+PiH22T5e0v0+x10r6Qds/L+lkSStsfyMihj0PhCkgbaE9UpjSNYUYUD3akfbgGgaQu7LP+GyXdLmkq4ufn+otEBE/3n1t++2S5uj0pIHpYtslheclqohh2s+l8RzccLQj7ZJCO8I1CWBcZZ/xuVrSBbYflHRB8V6252xfWzY4TBZpC8jNtJ8n4fmVxdGOYJq4JgGUUarjExFPRMQbI2J98fPJYvl8RLyzT/nrIuJdZbaJ6jBdLHIz7edJeH5lcbQjmCauSQBllE11Q+ZSSFsAlmraz5Pw/MrSNL0dIbUqHVyTAMoom+oGAFMz7el0mb4XpFalhWsSQBl0fABkY9rPk/D8CkitSgvXJIAy6PgAyEZVz5N0U5fufORJnXf17QP/957nV4631OPWJKRWpYVrEkAZPOMDICtlnycZlLrUXXfV22uKUY9bUzBdd3q4JgGMixEfAK1C6tJ42nrcSK0CgOZgxAdAq5C6NJ62HrfuyMKWHbu09+AhrV45w6xuAJApOj4AWoXUpfG0+biRWgUAzUCqG4BWIXVpPBw3AEDuGPEB0CqkLo2H4wYAyB0dHwBL1pRvsE8ldSm345nKcQMAYBykugFYEr7BvlocTwAApouOD4Alaet0xpPC8QQAYLro+ABYkrZOZzwpHE8AAKaLjg+AJRk0bXEbpjOeBI4nAADTRccHwJIwnXG1OJ4AAExXqY6P7VNt32r7weLnKQPKrbX9l7YfsH2/7XVltgtg+i7ZMKurLj1HsytnZEmzK2d01aXnMMvXmDieAABMV9nprK+UdFtEXG37yuL9e/uU+5ikD0bErbZPlnS05HYB1GDa0xmnPN1zFbExPTQAANNTNtVto6Tri9fXS7qkt4DtsyWdGBG3SlJEfCMiniq5XQANl/J0zynHBgAA+iu+HyOBAAAJTUlEQVTb8XlpROyTpOLnS/qUeYWkg7Zvtr3T9hbby/qUk+1Ntudtzx84cKBkaABylvJ0zynHBgAA+ls01c32ZyS9rM+v3jfCNn5Q0gZJj0r6uKS3S/pIb8GI2CppqyTNzc3FEtcPoIFSnu455dgAAEB/i3Z8IuL8Qb+z/bjt0yNin+3TJe3vU2y3pJ0R8XDxN9skvUZ9Oj4A0LV65Yz29OlIpDDdc8qxAQCA/sqmum2XdHnx+nJJn+pT5i5Jp9heVbx/g6T7S24XQMOlPN1zyrEBAID+ynZ8rpZ0ge0HJV1QvJftOdvXSlJEHJH0bkm32b5XkiX915LbBdBwKU/3nHJsAACgv1LTWUfEE5Le2Gf5vKR3Lnh/q6TvLbMtAO1TdrrnSU6HzVTUQHVSnrp+kBxjBtqu7IgPACSJKaeBPOR4reYYMwA6PgAaiimngTzkeK3mGDMAOj4AGoopp4E85Hit5hgzADo+ABpq0NTSTDkNpCXHazXHmAHQ8QHQUEw5DeQhx2s1x5gBlJzVDQBS1Z1dacuOXdp78JBWr5xh1iUgQTleqznGDICOD4CGGDS1LB9EgPTleK3mGDPQdqS6AcgeU8sCAIDF0PEBkD2mlgUAAIuh4wMge0wtCwAAFkPHB0D2mFoWAAAsho4PgOwxtSwAAFgMs7oByB5TywIAgMXQ8QHQCEwtCwAAhiHVDQAAAEDjler42D7V9q22Hyx+njKg3Ids32f7Adu/b9tltgsAAAAAoyg74nOlpNsiYr2k24r3x7H9A5LOk/S9kl4p6Z9Jel3J7QKYoG0792jnowd15yNP6ryrb+eLQAEAQPbKdnw2Srq+eH29pEv6lAlJL5C0QtJJkpZLerzkdgFMyLade7T55nt1+MhRSdKeg4e0+eZ76fwAAICsle34vDQi9klS8fMlvQUi4nOSPitpX/FvR0Q80G9ltjfZnrc9f+DAgZKhARjHlh27dOiZI8ctO/TMEW3ZsaumiAAAAMpbdFY325+R9LI+v3rfUjZg+zslfbekNcWiW23/UET8dW/ZiNgqaaskzc3NxVLWD6Baew8eGmk5AABADhbt+ETE+YN+Z/tx26dHxD7bp0va36fYv5Z0R0R8o/ibT0t6jaTndXwA1G/1yhnt6dPJWb1ypoZoAAAAqlE21W27pMuL15dL+lSfMo9Kep3tE20vV2dig76pbgDqd8WFZ2lm+bLjls0sX6YrLjyrpogAAADKK9vxuVrSBbYflHRB8V6252xfW5S5SdI/SLpX0hckfSEi/mfJ7QKYkEs2zOqqS8/R7MoZWdLsyhlddek5fDkoAADImiPSfJRmbm4u5ufn6w4DAAAAQMJs3x0Rc4uVKzviAwAAAADJo+MDAAAAoPHo+AAAAABoPDo+AAAAABqPjg8AAACAxkt2VjfbByT9Y91xLHCapK/UHQTGxvnLG+cvX5y7vHH+8sb5yxvnb+m+PSJWLVYo2Y5PamzPL2WaPKSJ85c3zl++OHd54/zljfOXN85f9Uh1AwAAANB4dHwAAAAANB4dn6XbWncAKIXzlzfOX744d3nj/OWN85c3zl/FeMYHAAAAQOMx4gMAAACg8ej4AAAAAGg8Oj6LsH2R7V22H7J9Zd3xYDjbZ9j+rO0HbN9n+5eK5afavtX2g8XPU+qOFYPZXmZ7p+0/L96fafvO4vx93PaKumNEf7ZX2r7J9v8trsPXcv3lw/a/L9rOL9q+wfYLuP7SZfujtvfb/uKCZX2vN3f8fvF55h7b31df5JAGnr8tRft5j+0/s71ywe82F+dvl+0L64k6b3R8hrC9TNI1kt4k6WxJb7N9dr1RYRHPSvrViPhuSa+R9O+Kc3alpNsiYr2k24r3SNcvSXpgwfvfkvQ7xfn7qqR31BIVluL3JP1FRHyXpFepcx65/jJge1bSL0qai4hXSlom6TJx/aXsOkkX9SwbdL29SdL64t8mSR+eUowY7Do9//zdKumVEfG9kv5e0mZJKj7LXCbpe4q/+cPicypGQMdnuHMlPRQRD0fEYUk3StpYc0wYIiL2RcT/KV7/kzofumbVOW/XF8Wul3RJPRFiMbbXSPpXkq4t3lvSGyTdVBTh/CXK9osk/ZCkj0hSRByOiIPi+svJiZJmbJ8o6Vsk7RPXX7Ii4q8lPdmzeND1tlHSx6LjDkkrbZ8+nUjRT7/zFxF/GRHPFm/vkLSmeL1R0o0R8XREPCLpIXU+p2IEdHyGm5X02IL3u4tlyIDtdZI2SLpT0ksjYp/U6RxJekl9kWERvyvpPZKOFu9fLOngghsB12G6Xi7pgKQ/LlIVr7X9QnH9ZSEi9kj6bUmPqtPh+Zqku8X1l5tB1xufafLzbyR9unjN+asAHZ/h3GcZ839nwPbJkv6HpF+OiK/XHQ+WxvaPSNofEXcvXNynKNdhmk6U9H2SPhwRGyR9U6S1ZaN4FmSjpDMlrZb0QnXSo3px/eWJtjQjtt+nTvr+n3YX9SnG+RsRHZ/hdks6Y8H7NZL21hQLlsj2cnU6PX8aETcXix/vDukXP/fXFR+GOk/Sxba/pE5q6RvUGQFaWaTeSFyHKdstaXdE3Fm8v0mdjhDXXx7Ol/RIRByIiGck3SzpB8T1l5tB1xufaTJh+3JJPyLpx+PYF25y/ipAx2e4uyStL2a0WaHOQ2Xba44JQxTPg3xE0gMR8Z8X/Gq7pMuL15dL+tS0Y8PiImJzRKyJiHXqXG+3R8SPS/qspDcXxTh/iYqIL0t6zPZZxaI3SrpfXH+5eFTSa2x/S9GWds8f119eBl1v2yX9VDG722skfa2bEod02L5I0nslXRwRTy341XZJl9k+yfaZ6kxS8Xd1xJgzH+tIoh/bP6zO/zgvk/TRiPhgzSFhCNv/QtLfSLpXx54R+Q/qPOfzCUlr1bm5vyUieh8IRUJsv17SuyPiR2y/XJ0RoFMl7ZT0ExHxdJ3xoT/br1ZnYooVkh6W9NPq/Ccb118GbP9HSW9VJ8Vmp6R3qvMcAddfgmzfIOn1kk6T9LikX5e0TX2ut6Iz+wfqzAj2lKSfjoj5OuJGx4Dzt1nSSZKeKIrdERE/V5R/nzrP/TyrTir/p3vXieHo+AAAAABoPFLdAAAAADQeHR8AAAAAjUfHBwAAAEDj0fEBAAAA0Hh0fAAAAAA0Hh0fAAAAAI1HxwcAAABA4/1/pXNuT1dc+1sAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2233b876978>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"X = np.fft.fft(x)\n",
"x_hat = np.fft.ifft(X)\n",
"\n",
"plt.stem(np.angle(X));\n",
"plt.show();\n",
"\n",
"plt.stem(np.real(x_hat - x));\n",
"plt.show();\n",
"\n",
"plt.stem(np.imag(x_hat));"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}