626 lines
53 KiB
Plaintext
626 lines
53 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
|
||
"grade": false,
|
||
"grade_id": "cell-98f724ece1aacb67",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"# CDS: Numerical Methods Assignments\n",
|
||
"\n",
|
||
"- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
|
||
"\n",
|
||
"- Solutions must be submitted via the Jupyter Hub.\n",
|
||
"\n",
|
||
"- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
|
||
"\n",
|
||
"## Submission\n",
|
||
"\n",
|
||
"1. Name all team members in the the cell below\n",
|
||
"2. make sure everything runs as expected\n",
|
||
"3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
|
||
"4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
|
||
"5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
|
||
"6. submit your solutions **in time (before the deadline)**"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "raw",
|
||
"metadata": {},
|
||
"source": [
|
||
"team_members = \"Koen Vendrig, Kees van Kempen\""
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "2d40509e944f45e5c23cc3e732547aad",
|
||
"grade": false,
|
||
"grade_id": "cell-2a0ab25436430f05",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"## Discrete and Fast Fourier Transforms (DFT and FFT)\n",
|
||
"\n",
|
||
"In the following we will implement a DFT algorithm and, based on that, a FFT algorithm. Our aim is to experience the drastic improvement of computational time in the FFT case."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "90044d2a3233721361765db06afba03b",
|
||
"grade": true,
|
||
"grade_id": "cell-abee6fbaf30772f2",
|
||
"locked": false,
|
||
"points": 0,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"import numpy as np\n",
|
||
"from matplotlib import pyplot as plt\n",
|
||
"import timeit"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "bbe05b6dd0bc5b479cf66199114d7e4d",
|
||
"grade": false,
|
||
"grade_id": "cell-a1c3327dc1cad0db",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"### Task 1\n",
|
||
"\n",
|
||
"Implement a Python function $\\text{DFT(yk)}$ which returns the Fourier transform defined by\n",
|
||
"\n",
|
||
"\\begin{equation}\n",
|
||
"\\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij x_k}\n",
|
||
"\\end{equation}\n",
|
||
"\n",
|
||
"with $x_k = \\frac{2\\pi k}{N}$ and $j = 0, 1, ..., N-1$. The $\\text{yk}$ should represent the array corresponding to $y_k = f(x_k)$. Please note that this definition is slightly different to the one we introduced in the lecture. Here we follow the notation of Numpy and Scipy.\n",
|
||
"\n",
|
||
"Hint: try to write the sum as a matrix-vector product and use $\\text{numpy.dot()}$ to evaluate it."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "2fded00d23ce12e99ba32c09a6370b3c",
|
||
"grade": true,
|
||
"grade_id": "cell-5f6638846212c9d1",
|
||
"locked": false,
|
||
"points": 3,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"def DFT(yk):\n",
|
||
" \"\"\"\n",
|
||
" Return discrete fourier transform (DFT) of yk for N discrete frequency intervals.\n",
|
||
" \"\"\"\n",
|
||
" \n",
|
||
" N = len(yk)\n",
|
||
" xk = 2*np.pi*np.arange(N)/N\n",
|
||
" beta = np.dot(yk, np.exp(np.outer(-np.arange(N), xk*1j)))\n",
|
||
" return beta"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "15164465870c44b9b4e6328f56eaed20",
|
||
"grade": false,
|
||
"grade_id": "cell-74e9ce917ff9d690",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"### Task 2 \n",
|
||
"\n",
|
||
"Make sure your function $\\text{DFT(yk)}$ and Numpy’s FFT function $\\text{numpy.fft.fft(yk)}$ return\n",
|
||
"the same data by plotting $|\\beta_j|$ vs. $j$ for\n",
|
||
"\n",
|
||
"\\begin{equation}\n",
|
||
" y_k = f(x_k) = e^{20i x_k} + e^{40 i x_k}\n",
|
||
"\\end{equation}\n",
|
||
"and\n",
|
||
"\\begin{equation}\n",
|
||
" y_k = f(x_k) = e^{i 5 x_k^2}\n",
|
||
"\\end{equation}\n",
|
||
"\n",
|
||
"using $N = 128$ for both routines."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "14cef25098dec15c058916f4fc58c133",
|
||
"grade": true,
|
||
"grade_id": "cell-7cc28776346ee714",
|
||
"locked": false,
|
||
"points": 1,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAGICAYAAAC5jEaHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAAklEQVR4nO3de5SlZXXv+998VwGtiMilIYbWNAYGcovdpuPG42WAGgMGheyNCBElXsLhDN0KbqMQo8Zztjl6dBjjIGKIuIGEKAbdigx1S1S2aEBshKEIEi5CbIPSoIgXbr3eef54L+tdq1ZVra56q9/1zPX9jNGjq1atqnoeVlOz5vM8cz7m7gIAAAAAoAtZ1wMAAAAAAMwuklIAAAAAQGdISgEAAAAAnSEpBQAAAAB0hqQUAAAAANAZklIAAAAAQGdISgEAAAAAnSEpBQAAAAB0hqQU2IHM7BlmdrWZfc3MPm5mO5WPv9fMrjKzf2g89nozO7DbERfM7GQz2zry2FSPGQCAaVLFSDPbamZXln/WLvL8F5rZW83s76o4C0RFUgrsWD+U9Dx3f66kOyUdZ2ZPk7Sfuz9H0vclnSBJ7n6Ou9/a2UhLZtaT9FIVY68em+oxAwAwbdz9HEmPSvrf7n5k+WfrIs//kru/V9KvJO28o8YJdIGkFNiB3P1ud3+wfPcRSbmk/0PSl8rHvijpWZJkZleWf+9sZtea2Roz+yMzO3cHD/tkSf9cjrWy4JinYLwAAEydKq5LelZ50uivrLBg3DSz0yV9yd1/1cmggR2EpBTogJn9lqQXSvqcpD0kPVB+6OeS9jSzvSXdI0nu/oik8yV9SNIrJP3XHTjOnqQTJV0y8qEFx9zleAEAmEaNuH63pAMkPVfSPpL+80Jxs0xIXyDpqWa2RxfjBnaUua4HAMwaM3u8pH+Q9Cfu/qiZ3S/p8eWHd5f0U0m/I+m7jU/7V0kfVnFkdtt2fK+/lCR3/8tlDvcUSZ9099zMmo8vNeZljRcAgKB+R9J33f1hSQ9Lkpl9WtIRkj6lMXHT3T8i6SPdDBfYsdgpRXhm9jgz65vZExuPHWZmd5vZbqv0Pf/UzG4ys5+b2RfMbJ/y8TlJn5D0Lne/pXz6v6pYCZWkP5D0DTUSvHJ19YOS3q4iSZSZ7Wpmf74DxnyIpFea2RclHWhmH1pqzDtqvAAALGaa4r8GMbL5fZ8j6bZxcROYNSSlCM/df6miGc/TGw+/R9JfufsvFvtcM7vczO5f4M/lC3zOn0s6XdJLJK2V9CNJ/7388MmS/pOkt5f1ly9z9xsk/cTMrpJ0qIoV08MlfcfMdpF0kaQzJL1P0svMbFdJvyfp383sbWa2fnv/m0w6Znd/q7u/0N2PlnSru79BkhYZ8807YrwAACxlyuL/4ZK+I+nZZnZdGT/3U1EeMy9uspiLWWPu3vUYgFVnZhdIut3d/x8ze66kCyQ91d0fMbO/lnSRu1/fwvfZR9IPJG10938rH3umpA+7+8aVfv3G93mLiiM/Jy52PHaS47s7YsyTjhcAgDalGv/N7EhJ6yT9lqSL3f3OlY4RmGbUlGJWfEuD46b/n6S3l40FJOlgFTt8bXi+pF0kXduowTRJKw54I/aQ9FZJZ2mwClt8s2IF99nlu2vKx84o3/+6ux/bwZgXHC8AAKso1fj/DLGYixlCUopZ8S1JbzGz/6IiUfunxsce4+4PmdlZknaV9A4vjxCY2RdU1HyMc5W7HzPy2J6S/qe7v7Td4Q9z97PLN+cleM2kc8JGR6s+5sXGCwDAKko1/rOYi5nC8V3MhLI28wEVrdj/L3f/Qvn4PpI+JmmrpCvc/Z8W/ioTfZ9nqrjm5YXu/u2y0+5Rki7zDv5nm/D47lSNGQCAtsxq/AdSQ6MjzISyBft3Jd1ZBaTS4SoaINyx0oBUfp+rJf3fkj5lZr+UdJOko6c5IKU4ZgAAJkH8B9LATilmgpntLOk2FbUZ1zQeP0PSFkmvknS2u3+nmxECAIC2Ef+BNLBTilnxTknfaAak0uGSvi3pTyT9jZntuaMHBgAAVg3xH0gAO6UIzcyeLumrKu4G+yN3v7fjIQEAgFVG/AfSQlIKAAAAAOgMx3cBAAAAAJ0hKQUAAAAAdIakFAAAAADQmbmuByBJe++9t69fv77rYQAAgrjuuuvudfe1XY8jZcRmAECbFovNU5GUrl+/Xps3b+56GACAIMzsrq7HkDpiMwCgTYvFZo7vAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDNTUVMKACl69NFHtWXLFj300ENdD2VmrVmzRuvWrdNOO+3U9VAAAFOA2Ny95cRmklIAWKYtW7Zot9120/r162VmXQ9n5ri77rvvPm3ZskX7779/18MBAEwBYnO3lhubOb4LAMv00EMPaa+99iLodcTMtNdee7EaDgCoEZu7tdzYTFIKACtA0OsW//0BAKOIDd1azn9/ju8CQMJ6vZ4OP/xwPfroo5qbm9MrX/lKnXnmmcqyTFdeeaWOO+64+vjM3nvvrYMOOkjf+MY39Mgjj+gHP/iBDjroIEnSX/zFX+iEE07ocioAAIRAbN5+JKUAkLDHPOYxuuGGGyRJ99xzj/74j/9YDzzwgN71rndJkp7znOfo8ssvn/d5d955p4499tj6cwEAQDuIzduP47sAEMQ+++yj8847T+ecc47cvevhAAAw84jNk2GnFABa8K7PfU83/ccDrX7NQ37z8Xrniw/drs95ylOeon6/r3vuuUeSdNVVV2nDhg2SpJe+9KV629ve1uoYAQCYVsTmdJCUJug/fvB97fkbT9Kax+za9VBW1U+23K5dH7+nHvf4PboeCpCshY4IAZFtue1G7fvkA7XTzrt0PRQAmIfYPB9JaWLyfl+7XXCkbjjodTrij9/e9XBW1SPnH6s7fuN5eub/+bddDwVY0vaumq6WO+64Q71eT/vss49uvvnmrocD7HC/fOBnWvsPR+qGje/S7x3/uq6HA6BDxOZ0UFOamDzvazd7UP7g/V0PZdU9zn+p3sM/73oYQDK2bt2q008/Xa9//etph4+Z9fCDv9Iu9qj6v/5Z10MBAGLzhNgpTUy/v6140fJ+10NZdZn6MzFPYCUefPBBbdiwoW47/4pXvEJvetObuh4W0Bnv96s3uh0IgJlFbN5+JKWJ8Tyv3uh2IDuAyWWiSxmwmH5/4V+8jzzySB155JFjP7Z+/XrdeOONqzQqoDu5z06cBDCdiM3bj+O7ien3txVvzMAKcM/zmZgnAKA9eRkn60VcAMDUIylNTF4GWZuBe45MLmOlGwCwHfL6RBGLmgCQCpLSxOQzVCvTU05SCgDYLp6XJ4roSQAAySApTYxXQXYGkjVTPhPzBAC0Z7BTGv9EEQBEQVKamKqm1GZgBZidUgDA9pqlE0UAEAVJaWK8TtLirwD3zIvdUgAAJjVDJ4oAIAqS0sTMygrwYJ78UgGk5Pvf/742bNigjRs36vbbb9eHPvQhHXzwwXr5y1+uz3zmM7rpppvGft7DDz+sF7zgBdqwYYMuueQSXXXVVTr00EO1YcMGXX311fr85z+/g2eCVOXUlALAkBRiM0lpYvIyyEY/1jor8wSi+cxnPqMTTjhB119/vX77t39bH/7wh3XFFVfo4osvXjTwXX/99ZKkG264QS972ct08cUX6+yzz9YNN9ygW265haQUE8vz8iQRNaUAICmN2ExSmph6BzH4CnBdOxt8RxhYqTvvvFMHH3yw/vRP/1SHHnqoXvjCF+rBBx/UkUceqc2bN0uS7r33Xq1fv16SdMEFF+j444/X7//+72v9+vU655xz9IEPfEAbN27UEUccoZ/+9KeSisu93/jGN2rDhg067LDDdO211yrPcx144IHaunWrpKKhzAEHHFC///nPf14f/OAHde655+qoo47S6aefrjvuuEPHHHOM3v3ud+uyyy7Tn/3Zn2nDhg26/fbb6zncc889OuWUU/Stb31LGzZs0N/93d/pk5/8pN7+9rfr5JNP1jve8Q5dcskl9UotsJiq+y7xA0BXiM3bb27FXwE7VHUZuAWvKZ2VeSKQL5wl/fi77X7N3zhcOuY9Sz7t1ltv1cc//nH9/d//vU488UR96lOfWvT5N954o66//no99NBDOuCAA/Te975X119/vc4880xddNFFOuOMMyRJv/71r3XDDTfoa1/7ml796lfrxhtv1CmnnKKLL75YZ5xxhv7lX/5FT3va07R27VpJ0ote9CKdfvrpetzjHqc3v/nNkqQvfvGL+upXv6q9995bt956q4499lidcMIJQ+PZZ5999NGPflTvf//7dfnll0uSrr766vq5F1xwgTZv3qxzzjlne/8LYgZ53X2XkzbAzCM2JxOb2SlNTF0rE3wFmJ1SYHL777+/NmzYIEn63d/9Xd15552LPv+oo47SbrvtprVr12r33XfXi1/8YknS4YcfPvS5J598siTpuc99rh544AHdf//9evWrX62LLrpIkvSxj31Mr3rVq1qfD7ASeX824iSA6UZs3j5L7pSa2cckHSvpHnc/rHzsfZJeLOkRSbdLepW7319+7GxJr5HUl/QGd/9fqzP02eQzUmuZs1OK1Eywarpadtlll/rtXq+nBx98UHNzc/X/Rw899NCCz8+yrH4/yzJt27at/piZDX2emelJT3qS9t13X33lK1/Rtddeq4svvrj1+QArUXWpjx4nAUyA2JyMSXZKL5B09MhjV0g6zN1/R9K/STpbkszsEEknSTq0/JwPm1mvtdFC3p+RpLSeJyvdwHKsX79e1113nSTp0ksvXdbXqGpEvv71r2v33XfX7rvvLkl67Wtfq1NOOUUvfelL1etN/iN+t9120y9+8YvtHsdyPw+zybkSBsCUIjYvbMmk1N2/JumnI499yd2rlP0aSevKt4+T9Al3f9jdfyDpNknPaG20aBzfjR1sq+NXWfB5AqvlzW9+s84991xt3LhR995777K+xpo1a7Rx40adfvrpOv/88+vHX/KSl+iXv/xlfTzota99bd24YTEnnXSS3ve+99Ut6T/ykY/oIx/5yJKfd9RRR+mmm26i0REmMqgpZVETwHQhNi/MfIKW6Wa2XtLl1fHdkY99TtIl7v6PZnaOpGvc/R/Lj50v6QvuvuhSwKZNm3yS/2iQbv/uNfrtT/2Bvr3rc/X0P/tc18NZNff9ZIv2OvdQfX+nQ/TUt13d9XCAsW6++WYdfPDBXQ9jVRx55JF6//vfr02bNs372ObNm3XmmWfqqquu6mBk8417HczsOnefP3hMLNXYfOM3PqfDrjhF1+75Yj3jDf/Y9XAA7GDE5jRj84q675rZ2yRtk7TdB5fN7DRJp0nSk5/85JUMY6bUre4VewfROb4LTKX3vOc9Ovfcc5OsV8HSQsTmaqc0+NVpAFCJEJuX3X3XzP5ERQOkl/tgu/VHkp7UeNq68rF53P08d9/k7puqlsVY2qy0us+dRkdAl6688sqxK7FnnXWW7rrrLj372c/uYFRYbRFic90QkPgBIJjIsXlZSamZHS3pLZJe4u6/bnzoMkknmdkuZra/pAMlXbvyYaKS11elBE9K65pSVroBAJOjphQA0jPJlTAfl3SkpL3NbIukd6rotruLpCvKtsTXuPvp7v49M/ukpJtUHOt9nTtRoU2DVvex/7NyJQxS4e7z2rNjx5mkLwJmS9WHMfriLYCFEZu7tZzYvGRS6u4nj3n4/DGPVc9/t6R3b/dIMJFZOZZU1c7SfRfTbM2aNbrvvvu01157Efw64O667777tGbNmq6HginieREfSUqB2URs7tZyY/OKGh1hx6uOJc3OTim/VGB6rVu3Tlu2bNHWrVu7HsrMWrNmjdatW7f0EzEzuKcUmG3E5u4tJzaTlCYmr7vvxt4pzavuuySlmGI77bST9t9//66HAaDJq+7txA9gFhGb07Ts7rvoyIzslFat/Dm+CwDYHnmfkzYAkBqS0sT4jCRr+YzcxwoAaFfV6IjjuwCQDpLSxPiMdKXNy0YVWfB5AgBaVp8oIikFgFSQlCZm0Oo+9vFduu8CAJaj7lIfPE4CQCQkpYmpW90H30F0uu8CAJbDZ+NEEQBEQlKamEFNaewV4Lxf7pSSlAIAtgc7pQCQHJLS1FSt7oOvADsr3QCAZSB+AEB6SEoTU7W6j15rWe0I99gpBQBsB2pKASA9JKWpqXdKYydr1JQCAJal2il1dkoBIBUkpYmpV4CDJ2vVPaXUlAIAtkvVe4H4AQDJIClNjc/G8d3qnrmMlW4AwHaoa0qjx0kACISkNDGzslPqrHQDAJZjRuIkAERCUpqYagU4C95VsKopJSkFAGwXdkoBIDkkpamp7ymNHWzdqSkFACwDJ20AIDkkpYmZla60nhc7wdF3hAEA7aKmFADSQ1KaGp+NFeCqpnTOYs8TANCyGYmTABAJSWlqqhXg6DuIjUvPq91hAACWVHZtDx8nASAQktLEVDuIveArwHl/ML9+f1uHIwEAJKXuvdBf4okAgGlBUpoan42a0uZOac5OKQBgQjYrJ4oAIBCS0tTMSFfBakdYknJ2SgEAE3JqSgEgOSSlqanuKfXgK8CNrol5zhEsAMCE6L4LAMkhKU2Mz+BOab9PUgoAmIzNSJwEgEhIShNjXt3fGTvY+tBOaey5AgBaRE0pACSHpDQxM1Mr09gpdWpKAQCTqspcosdJAAiEpDQ1dbCNvQLcvJuUmlIAwMRISgEgOSSlialqZeYseLAduhKGpBQAMBmblRNFABAISWlimrWWHrnWclbmCQBo14ycKAKASEhKU+PNrrRxay2Hu+/GnScAoGX11WksaAJAKkhKUzMrXWmHdko5vgsAmIxRUwoAySEpTYw1ErQ88g7i0Dz5xQIAMCGSUgBIDklpcgY1MqEbALFTCgBYBqOmFACSQ1KamqFay7jJWjMRzfPAO8IAgFbRfRcA0kNSmhibkZpS88EKN913AQCTK+IHO6UAkA6S0tQ0uu964JpS555SAMAy1Dul5ixqAkAiSEqTM4s1pfxSAQCYzKycKAKASEhKEzPUfTdwUjozXYYBAO1qJKXccw0AaVgyKTWzj5nZPWZ2Y+OxPc3sCjO7tfx7j/JxM7MPmdltZvYdM3v6ag5+Js3IDqI358kF6ACACQ3vlMZdvAWASCbZKb1A0tEjj50l6cvufqCkL5fvS9Ixkg4s/5wm6dx2homKzcoKcLN2ll8qAAATskbXXU7aAEAalkxK3f1rkn468vBxki4s375Q0vGNxy/ywjWSnmBmT2xprJAkzcj9nax0AwCWgZpSAEjPcmtK93X3u8u3fyxp3/Lt/ST9sPG8LeVj85jZaWa22cw2b926dZnDmD1DwbYfN9g2a0oV+D5WAJgmEWLz8Iki4gcApGDFjY7c3aXtvwzM3c9z903uvmnt2rUrHcbMMJ+RndJml2FqSgFgh4gRmxsxI3ScBIA4lpuU/qQ6llv+fU/5+I8kPanxvHXlY2iJDd3fGbhWJm/ex8ovFQCAyWSNOBm69wIABLLcpPQySaeWb58q6bONx19ZduE9QtLPG8d80QYf7CBG7r47tCPsJKUAgMmYN+/zjhsnASCSuaWeYGYfl3SkpL3NbIukd0p6j6RPmtlrJN0l6cTy6Z+X9CJJt0n6taRXrcKYZ9rwTmngZI3uuwCAZbBZaQgIAIEsmZS6+8kLfOj5Y57rkl630kFhMbOxUzo78wQAtIl7SgEgPStudIQdq1krE/n+tWb3XVa6AQCTGu5SHzdOAkAkJKWpGaq1DLyDOCvzBAC0avj47nZfDgAA6ABJaWJmpVZmqNERK90AgAllYqcUAFJDUpqYbEYaHTUbOonuuwCACQ3XlJKUAkAKSEoT02x1r9D3d9LoCACw/YZPFBE/ACAFJKWJaQbbPHCtpXElDABgGbKhe0qJHwCQApLSxJjnyt0kSR54p3RonoGTbwBAu0z9QfxgpxQAkkBSmhjzXNvUkyR55FrLxjzFSjcAYEImr+MHjY4AIA0kpYkxNZLSwMnaUPLNSjcAYELZ0OIt8QMAUkBSmphMubbZnKTYyZo15xl5RxgA0Kqh+BF48RYAIiEpTczwDmLcYGtDx3fjJt8AgHZljRNFNDoCgDSQlCamOL5b7SDGTdaKpJSVbgDA9snc6/gR++o0AIiDpDQxmbv61U5p6AYOeT1PBU6+AQDtytSv40fkq9MAIBKS0sSY+urXtZa+xLPTlXlf22yn4h1qSgEAEzL5IE6yUwoASSApTUwz2CqPu1Nq7uob3XcBANsnU95YvCUpBYAUkJQmJvPBsdbQNaWNXyrYKQUATCpTo8yFngQAkASS0sQ0k7XIwdY8V14npXGTbwBAu7JmmQsnbQAgCSSliSmOJZW1loGDrXmuftU9MXDyDQBoV1HmUsTJyIu3ABAJSWliMnfl2QzslCofzDNwQycAQLt63owfcRdvASASktLEZOort/g1pZnyep4WOPkGALTL5IM4GfrqNACIg6Q0MUWwrY7vxk3WiprS8vgVjY4AABPqqRk/OGkDACkgKU1M1jjWGrkBkM3IPAEA7RqKH4GvTgOASEhKE5PJ6660kWtKM3d5efyKpBQAMKlip5SaUgBICUlpYjL15TOwg2jqy5Wp7xb6mDIAoF098zpORl68BYBISEoTY3L5DOyUFvPM1FcWOvkGALQn7xdxsYqTka9OA4BISEoTU7S6LxsdBU7WMs8ly+TKJBodAQAmkJeLtVWcjLx4CwCRkJQmxuSN47txg60pr3dKje6JAIAJ9MsrYJx7SgEgKSSliekpl2agpjQrk1KXhU6+AQDt8eq4bt19l/gBACkgKU2MKZfXx5ICJ6XukvXUN2pKAQCTGeyUxi9zAYBISEoTM7xTGncFOFNfLpPLZKx0AwAmkI/slFJTCgBpIClNTM9c3itXgAMHW5NLWU+5epKoKQUALK3uvttjpxQAUkJSmpAq2Ko+lhQ3WatqSnNqSgEAE6p3Rum+CwBJISlNSF4H217xd+BkLVNRU5ork7HSDQCYQFVTql78hoAAEAlJaUKqYGvWU98tdLDN1JfKmtLIx5QBAO2proCx+kQR8QMAUkBSmpCq265nPfWVhU7WivtYi5pSo6YUADCBusyFmlIASApJaUIGO6UmVxa6prTnuWSZcqOmFAAwmarMxXrxr04DgEhIShMyaHVf7JRa4GTNqCkFAGyn+TulceMkAESyoqTUzM40s++Z2Y1m9nEzW2Nm+5vZN83sNjO7xMx2bmuws64OtpYVtZaBk7Wecnm5I0xSCgCYRLUzms3A1WkAEMmyk1Iz20/SGyRtcvfDJPUknSTpvZL+2t0PkPQzSa9pY6BotLa3nvqWhV4BNuXFTqlloZNvAEB78rzqvhv/6jQAiGSlx3fnJD3GzOYkPVbS3ZKeJ+nS8uMXSjp+hd8DpbqmNCt2SiPvIPaU1zvCkecJAGiP1zWl1ZUwcRdvASCSZSel7v4jSe+X9O8qktGfS7pO0v3uXi5Vaouk/VY6SBSqVvdFrWUv9A5iz1zKenL1il1TAACW4GWZi2Xxr04DgEhWcnx3D0nHSdpf0m9K2lXS0dvx+aeZ2WYz27x169blDmOmNGtK88A1pUPztLjzBIBpk3pszqt7Si1THvzqNACIZCXHd18g6QfuvtXdH5X0aUnPkvSE8jivJK2T9KNxn+zu57n7JnfftHbt2hUMY3bUre6zqittzGCb582GTjQ6AoAdJfXYXC9qZnNFUkpNKQAkYSVJ6b9LOsLMHmtmJun5km6S9FVJJ5TPOVXSZ1c2RFQGwTZ2993Bfaw9ucVNvgEA7aprSs2Uy4gfAJCIldSUflNFQ6NvS/pu+bXOk/RWSW8ys9sk7SXp/BbGCQ1a3VtZUxp1B7Gap1c7wmKlGwCwtGajo2KnNGacBIBo5pZ+ysLc/Z2S3jny8B2SnrGSr4vxqlb3lsWutRzslFb3lLLSDQBYmg/VlBrddwEgESu9EgY7UH1Pafia0jLZzsrju+yUAgAmkNdXpxX3XEc9UQQA0ZCUJqRudV82AIrawKHZfZeaUgDAxKoylyzj+C4AJISkNCF1q/t6BThmslbvCFuvPL4bM/kGALRrtEs9SSkApIGkNCF5fSl40X036rGkuqY0y+RmyhRzngCAdtW9F6xXxsmYi7cAEA1JaUKGdxB7sqDJWtWoQtZTbnG7DAMAWlYd3+2xUwoAKSEpTYg3jyUF7r7brCmVsrDJNwCgXVWczCwrGwISPwAgBSSlCfFGTakHDrbNmiCneyIAYELVPdfqFb0Xoi7eAkA0JKUJqVvdB+9KW++UZsU8qSkFAEzCvYiTWUZNKQCkhKQ0JSO1MlHv76xWus1i7wgDANrleREXi+67vbBXpwFANCSlCamPtZb3lEZdAa67J7JTCgDYDt6Mk8ZOKQCkgqQ0IfUOYjZXHN8Nu1NaHd/tSYHnCQBoWZmEZr258kQRi5oAkAKS0oT4yA5i1BVg7zdXujNlQecJAGhX3q8WbzPKPwAgISSlCam672Z1992YO4h5s8uw9dgpBQBMptopzei+CwApISlNiA91pbWwtZZV913LMklx5wkAaNfgPm92SgEgJSSlCalb3ffmlFsvbLCtjinLip3SLOg8AQAtq08UzYUucwGAaEhKE1K3urdMCtzAYdDQqVc2dIo5TwBAuwY7peU9pZR/AEASSEoTMhRsLe6xpLxfNXQquu9yfBcAMIm690IvU64eO6UAkAiS0pT4cFIaNVmrfqmouu+y0g0AmMjQ4i07pQCQCpLShFSt7rNeL3QDh3pHuFfslPaCJt8AgHZV5R+93lzRvZ2dUgBIAklpShqt7kPvlDZqSmU9akoBAJPx0e677JQCQApIShPSbHWvwMda87L7rlns5BsA0LLmfd6Br04DgGhISlMy0uo+i3osqdop7fWkrKeMlW4AwASqxdss64W+Og0AoiEpTchw991e2J3S+pcKyySx0g0AmFDVKK83p8hXpwFANCSlCWm2uo+crFU1per15Fkv7DwBAC2rd0qz0FenAUA0JKUpGd0pDRps3Yua0ox7SgEA28Or7rv0JACAlJCUJmS41X3cYOt5cSy56r6bBT2mDABoV7OmNPLVaQAQDUlpSnw2uu/WtbOWySzTnPFLBQBgadUVMFnwxVsAiIakNCVDre4DB9vqPtbenDzrFQ/lQecKAGiN+6CmNPLiLQBEQ1KakOaxpMi1lnm/7J5Y/lIhSf3+ti6HBABIgY+UuUS9Og0AgiEpTUmz1X3gpLTeKc16Mit2SnN2SgEAS7B8UOYS+eo0AIiGpDQlQ63u4zYA8qFfKkySlLNTCgBYgnuubV79ahP36jQAiIakNCWNVveyTFnUroJ17exc0YFXUp5zBAsAsATvK1exmBn56jQAiIakNCGzUlPqjftYBzWlJKUAgCV4rrz81SZ0Q0AACIakNCHNVvehk9Jqp7SXSRk1pQCAyVjer5NSuu8CQDpIShPSbHXvWdyaUo3ZKXVqSgEAS3J2SgEgQSSlKWm0ug+9U5oP5klNKQBgYnlfedkgL3KcBIBoSEoT0mx1b4G771ZXwjTvKSUpBQAsxRo1pSSlAJCOFSWlZvYEM7vUzL5vZjeb2TPNbE8zu8LMbi3/3qOtwc66Zqt7t0yZeb2rGErdfbcnlfeUhpwnAKBd3pfXx3cDL94CQDAr3Sn9G0lfdPenSnqapJslnSXpy+5+oKQvl++jDY1W95EbADW7DFtWdd+lphQAsBRXv7lTypUwAJCEZSelZra7pOdKOl+S3P0Rd79f0nGSLiyfdqGk41c2RNSGjiUVyWnIZK38JcJ6c3Xy7RzfBQAswfK+XNSUAkBqVrJTur+krZL+h5ldb2YfNbNdJe3r7neXz/mxpH1XOkgUhlvdB24AlA+6DFtVU9rnFwsAwBKoKQWAJK0kKZ2T9HRJ57r7Rkm/0shRXXd3aXxBh5mdZmabzWzz1q1bVzCMWTJodT9I1uLulPZ6vbr7LjulALD6Uo/N5vmgpjTy1WkAEMxKktItkra4+zfL9y9VkaT+xMyeKEnl3/eM+2R3P8/dN7n7prVr165gGDOk2ep+RmpKB913AybfADBl0o/NuXJjpxQAUrPspNTdfyzph2Z2UPnQ8yXdJOkySaeWj50q6bMrGiFqo63uJanfj7eDaF6sbGeNe0rpvgsAWEoRJ4vF29BXpwFAMHMr/Pz/KuliM9tZ0h2SXqUi0f2kmb1G0l2STlzh90Cl0eq+2ilVwGOt7oOa0tC1swCAVg0d321cnVZ1cgcATKcVJaXufoOkTWM+9PyVfF0sZNDq3izwVSl1Temcsl75ywU7pQCAJZj3B8d3G2UuPZJSAJhq/JROyFCr+8A1pVbuilpzpzRi8g0AaJf74ERR5KvTACAYktKUjKkpjdiV1j3XNi93hKuaUi5ABwAswZplLpR/AEAySEoT0qyVsSxwsPX+oFFFVvwdMfkGALTN6+O7oa9OA4BgSEqTks9GsG3sCJsVZc8hk28AQKsyn40yFwCIhqQ0Ic1W96qvSonX7t7y/iAprZpTBLz6BgDQMs+Vl8d2I1+dBgDRkJQmpHl8V5F3SuWNpLRc6aamFACwBFM+b6c04tVpABANSWlCmq3uBw2AAgbbvK/cql8qyoZOrHQDAJZgnsvLndLQV6cBQDAkpSlptrrPqp3SeMmaNWpKs6yoKQ2ZfAMAWpX5/J1SakoBYPqRlCak2eq+agAU8qqU5jyzuFffAADaNmgIGPnqNACIhqQ0KY1W9+VVKRF3SiVXf6Sm1FnpBgAsIfNcmoWr0wAgGJLShDRb3Vt1rDVgsLW8MU9WugEAEzLuKQWAJJGUpqTR6r4OthGTtWZNaS/wMWUAQKvM+/Lq+G7gq9MAIBqS0oQ0W91br9pBjJesNa++qWtKWekGACzB5DNydRoAxEJSmpBmq3tZVSsTMdjmjdrZaqeUlW4AwOIy79fJaOir0wAgGJLShDRb3Wf1peAxd0rzep5la/+QyTcAoE1FTWkZHwNfnQYA0ZCUJqXZ6j5uV8Hm8d3BPaXxkm8AQLuK7rtVozziBwCkgqQ0IUOt7kPXlPYHyXfdqCJe8g0AaJdpUOYS++o0AIiFpDQhw63uq2Qt4LFWHzSqyHpxjykDANplyiWLf3UaAERDUpqQZqv7KlmL2OrevN84vss9pQCAyWSNhoChr04DgGBIShPSbHVfBduQO6XNHWFqSgEAEyqO78YvcwGAaEhKEzK+1X28YJt5v9FluPwnyko3AGAJmbxOSmNfnQYAsZCUJqTZ6t4iNwDyvJ5n1purHwMAYDHNhoCRr04DgGhIShPSbHUfuQGQqXEfq1FTCgCYTPP4buSr0wAgGpLShAy1ug/cwMGajSqq5JudUgDAEjLlgzIXakoBIBkkpQlptrqvd0o9XlKaeWOnNPIxZQBAq0wuz6rF28BXpwFAMCSlCWm2us8C379myuvuu4Oa0njzBAC0q6e8PrY7uDqN+AEA046kNCFDtTKBu+82j+9mHN8FAEzINOi9MGgISPwAgGlHUpqQrJGU1slaP94KsA0d3+VKGADAZDLljeO7NMoDgFSQlCYkcx8cSyqTtZA7pY2GTr3y+K67dzkkAEACijg5XP4RMU4CQDQkpQnJ1K93EC1wTWnWuBLGyuTbAs4TANCurFFTWsUP79PoCACmHUlpQkxe15JWO6UWcAW4WVM62CklKQUALK4ocxlZvOWkDQBMPZLShAzXlMbdKW02dKquhKHREQBgKc2d0iwrklNxJQwATD2S0oRkatSUWtya0uGaIJJSAMBkijhZLWpSUwoAqSApTUimvupW91WyFnKntC9v/NPsu4WcJwCgXXOW12Uu9dVpxA8AmHokpQkxed3qPvKxVpMP7mOV1FcWcp4AgPbU95GOnrThnlIAmHokpQnpeT6/1X3AFeCsMU9Jxa4pjY4AAIvoV1126+O73FMKAKkgKU2INWpKq1b3EXcQm42OpGKn1OieCABYRF7uiFodJ6kpBYBUkJQmpNdodd+rGwDFWwHORpJSl4WcJwCgPXm5U+ojV6fRkwAApt+Kk1Iz65nZ9WZ2efn+/mb2TTO7zcwuMbOdVz5MSMUO4qDVfdxamaL7bq9+v2/UlAIAFpeXyaeVi7dVmQvxAwCmXxs7pW+UdHPj/fdK+mt3P0DSzyS9poXvARU7paMNHDzgDmKmfrE7WnKZjJVuAMAi+v0yTmQjV6cRPwBg6q0oKTWzdZL+UNJHy/dN0vMkXVo+5UJJx6/ke2CgZz4ItsG779Yt/SXl6kmiphQAsLB8pPuucc81ACRjpTulH5T0FknVT/y9JN3v7mULPG2RtN+4TzSz08xss5lt3rp16wqHEV9erQCXwbZXHUuKeHx3pKY0p6YUAHaIlGOz1913hxdv2SkFgOm37KTUzI6VdI+7X7ecz3f389x9k7tvWrt27XKHMTOqWpl6BbismYm4ApxpuKY0VyYLOE8AmDYpx+a6pjQbvjotYpwEgGjmVvC5z5L0EjN7kaQ1kh4v6W8kPcHM5srd0nWSfrTyYaLf36Y5NVvdZ+p7zB3ETH1ppKaU7okAgMUMFm9Hr04jfgDAtFv2Tqm7n+3u69x9vaSTJH3F3V8u6auSTiifdqqkz654lJCXx3Q9G95BjLgCbPKRefaKOlMAABbgIzWlPWpKASAZq3FP6VslvcnMblNRY3r+KnyPmdMva2XqY7sqk9KANaU9H3QZlqTcYu4IAwDaU8fJbLSmNF6cBIBoVnJ8t+buV0q6snz7DknPaOPrYqDuKji0g2iygMmaUVMKANhOdUOjbPjqNBY1AWD6rcZOKVbBaPddKe7x3Z5yuTVrSklKAQCLy/tFnLCR7rv0JACA6UdSmggfaeAgVVelxEvWTPnwPC1m8g0AaI+PdN+tr05zehIAwLQjKU3EoFamWWuZhTy+29NwTanL2CkFACwqz8t7Sssd0sHVafHiJABEQ1KaCK+SMovffbdnPlQ76+oVu6cAACzAyzIXq+7zrq9OI34AwLQjKU3ErNSUjp2n8UsFAGBxeRknbPTqNGpKAWDqkZQmIq9rZZo7iPG67w4uP28e36XREQBgcdWi5lCZizJqSgEgASSliah3EOcF21jJ2uA+1kbyHbR2FgDQnoUaAhI/AGD6kZQmorr824Lf31nN00eOX5lY6QYALMzHnCiKuHgLABGRlCai6io42n03WrAd7JSO3lPKSjcAYGE+tqbU6L4LAAkgKU1EfSwpfE1pmWRnI8d32SkFACwirxc1R69Oi7V4CwARkZQmYrTVvSTl6oVr4DCu+y41pQCAJVVlLj2O7wJAakhKEzGu1b1bvJ3ScY0qiuO7sZJvAEC76i71wa9OA4CISEoTsVCre1OsYFvXlGbNnVJTFmyeAIB21Q0Bs7nBYwHLXAAgIpLSRCy8gxgrWasaVQy19LdeuHkCANrl4xoCslMKAEkgKU3E2Fb3AbvvjqspVcAdYQBAu6pFzWz0SrFgcRIAIiIpTcS4VvcRd0rzMcm30z0RALCEqiGggl+dBgARkZQmYlyr++KqlFjBNh/zS4VbRk0pAGBR9U5pj5pSAEgNSWkqxrS6L4JtrGStblQRvHYWANAuH9t9N97VaQAQEUlpIsa3uo/XACgf06iCnVIAwFLG9V6IeHUaAEREUpqIsa3uzcId3627DDd+qZBlMrHSDQBYhI9pCEijPABIAklpIsa1uveAV6VUjSpGa2czVroBAIuoFm+zHuUfAJAaktJEjGt1H7GmNB/XZdh67JQCABZVnbTJjO67AJAaktJEjGt1H7n7bnNHWDJqSgEAi6qvTmOnFACSQ1KaiPGt7jNlwYJtXVNqIzulweYJAGjZQvdcs6gJAFOPpDQR41rdF8daYwXb8d0T6b4LAFhcXVOajd5TSvwAgGlHUpqIhVvdxwq2PqamlO67AICluM8v/4h4dRoARERSmooxre4jNgDK+2WX4dHuu+yUAgAWM677bsCr0wAgIpLSRIxvdW/xrkrJ5zeqEEkpAGAJdffd0cVbdkoBYOqRlCZiXKv7kDulY2pnSUoBAEvy8Yu3JKUAMP1IShMxrtW9LN5VKdWOsDUbVVhPWbDkGwDQMh/XEJDuuwCQApLSVNTHkoaTtWjB1vOyprSZfGe9cFffAADaVS1q9hpXp+VG/ACAFJCUJmKwgzh8KXi0YFvfx9o8vqt4O8IAgJaNaQgosVMKACkgKU1E1eo+6w3XWkYLtt4vGzc1a4KyHkkpAGBx42pKLaOmFAASQFKainE7pRZ4pzQb7b5LTSkAYBF1mctoTSnxAwCmHUlpIqruu81amYjdd6t5Dh2/MnZKAQCLcy/i4VCcVBbv6jQACIikNBV1rUzz+G68WksfcyWMWaY5izVPAEC7LJ8fJ50rxQAgCSSlqRhzrNUj7iDWtbONle5yzlWzJwAARlW9F5o7peL4LgAkYdlJqZk9ycy+amY3mdn3zOyN5eN7mtkVZnZr+fce7Q13dnldKzNaaxkrUcv7Ve3scEMnSer3t3UxJABACsYu3nJ8FwBSsJKd0m2S/pu7HyLpCEmvM7NDJJ0l6cvufqCkL5fvY6XKYGvRV4B9fvJtVryds1MKAFjI2O678XovAEBEy05K3f1ud/92+fYvJN0saT9Jx0m6sHzahZKOX+EYIY3tKijL1Au2U+pja4JMkpSzUwoAWEjeV99t5MF4vRcAIKJWakrNbL2kjZK+KWlfd7+7/NCPJe3bxveYeeUKcG90BTjYlTCD41eDHeGqE2+ecwQLALAAz9Uf+bUmZJwEgIBWnJSa2eMkfUrSGe7+QPNjXvRnH3tuxsxOM7PNZrZ569atKx1GeAvVlMbdKR2epyT1+ySlALCako7N3pfPS0rj9V4AgIhWlJSa2U4qEtKL3f3T5cM/MbMnlh9/oqR7xn2uu5/n7pvcfdPatWtXMoyZYOX9a1nwmlKva4Ia/zQzakoBYEdIOTab+7yd0ohxEgAiWkn3XZN0vqSb3f0DjQ9dJunU8u1TJX12+cNDxX1+TalbL9xOqRbZKXVqSgEAC/G+XMM1peyUAkAa5pZ+yoKeJekVkr5rZjeUj/25pPdI+qSZvUbSXZJOXNEIUahrShsvWRZvBbi6i7Q5T2pKAQBL8lx9m79TSlIKANNv2Umpu39d0mibu8rzl/t1MZ6N6Uobsaa0uhJm3D2lJKUAgIVYPn+nlKQUANLQSvddrD73XNt8dAW4p8y83l0MYczl5yrvKQ01TwBAy1y5esOPWE9ZsBNFABARSWkqvK98zAqwFKsB0Lguw9WuaZ+aUgDAQhaIkxlXwgDA1CMpTYXnykdfrizgsdbylwcbqp2tdkoDzRMA0CobFyc5vgsASSApTYTl/bHBVgq2g5jP7zJs1Y5wn18sAAALoKYUAJJFUpoMn5eUWsRay7rLcPP4LjulAIDF2bia0oyaUgBIAUlpKvK+chtZAQ5YazmupnRQOxtnngCAlvmYOMlOKQAkgaQ0EeNrZar7O+MEXPNiRTsbc09pqB1hAECrxsVJo/suACSBpDQV3pcvUFPqkXZKfX5N6SD55vguAGA883xenHTL4l2dBgABkZQmw9UfXQHO4u2UDmpKBzulWa9MviPNEwDQLs+V22iX+oBxEgACIilNhC3QVVCKVVNq5W6ojdspDTRPAEC7ip3S0ThZvB8pTgJARCSlqRh7T2nZDCjQCrB7rm0+fkfYuQAdALAAUy4f6b5L+QcApIGkNBHjamWq+zv7kbrSel/5yEq3ZcX7XAkDAFiQ5/O67w7uuQ4UJwEgIJLSZMyvlQnZlXZs98SivpSVbgDAQsYt3lJTCgBpIClNRNHqfngF2OsV4DjJmuX9+UlpVV8aaJ4AgHaZ9+u4OHiw6r1A/ACAaUZSmoixx3frWstIwdbHJKXlSjc1pQCABdiY+DHovRApTgJAPCSliTDvzz++G7FWJu/PqwlSVt3Hyi8VAIDxbMx93hawSz0ARERSmgr3+bUyvapWxjsY0OqwMTWlWVbUlMbaEQYAtMnk84/vUlMKAEkgKU3EYivAHqz77vxjytU8SUoBAOMtVlNK/ACA6UZSmgxfpPtupGDr6i9UO8tKNwBgATbmRFHdkyBUnASAeEhKE5F5Xz7Sfbe+FDxQraXl8+dprHQDAJaQKZdzTykAJImkNBWeKy+T0EpWd98NtIM4rqa0V9WUBponAKBVNiZOqj5pE6f3AgBERFKaCFM+f6c0C3hP6dirb6ruu6x0AwDGM+Wa92sNO6UAkASS0kSY5/KRFeCYNaX5mNrZaqeUlW4AwHhFnJyF+7wBIB6S0kRkPn+nNGKtZXElzPA8s6x8P1KXYQBAq4qa0tErYeKdKAKAiEhKk7HYDmKcWstxx3ezgPMEALRrbPmHET8AIAUkpYnIfH6tjJU7iJFWgM3785JvhTymDABo07id0ohxEgAiIilNhI29p7RaAQ4UbMfcM5f1ylpa7ikFACzA5HVjo/qxKk6yqAkAU42kNBHm/TErwFVX2jjJmnl/zPHdeLWzAIB2ZePiZNV9l/gBAFONpDQRpjE7iOWxVkXaKV10RzhO8g0AaJfJ53ep71WLmsQPAJhmJKWJyLw/71hS3VUw0Apw5v15XYarnVIFmicAoF2Zcmn0Pu8ySc3p3g4AU42kNBFFTenwCnDVlTbUTqnn8+fZm6s/BgDAOOPu865PFLFTCgBTjaQ0EUX33ZF7SiPWlGr+faxZwPtYAQDtGntPab1TSvwAgGlGUpoI08IrwJFqLcetdFvVfTfQPAEA7RrbfZeaUgBIAklpIky5ZCM7iL3q/s44tTKZj9kp5Z5SAMASMuXybGRR0+LFSQCIiKQ0Edm4HcQ6WfMuhrRK8nndd6kpBQAsZVyjo8HibaQ4CQDxkJQmwsbUygyuSomzAlzUzo6/jzVUQycAQKsy5XUNacXqngRx4iQARERSmohMPi8pzbJyRThQrYyNuae0R00pAGAJ4+KkBey9AAARrVpSamZHm9ktZnabmZ21Wt9nVmQ+fwW4OtYaqdbSvD/vmHKvnie/VAAAxss8l7L4cRIAIlqVpNSKzgJ/K+kYSYdIOtnMDlmN7zUrxl2VUjVwiLSDaPL58+T4LgBgCeNqSuv4waImAEy11dopfYak29z9Dnd/RNInJB23St9rJmQatwJcXQoeJ1nLFtkpjTRPAEC7xnbfLXsvcE8pAEy3uVX6uvtJ+mHj/S2S/tMqfS9J0k/v+ZG2nvefV/NbdOop/sCYWpni/f2++7e65aZ/6mJYrXvSti26ZZd9hx6rroR5yl3/rFv++5UdjArAjrLtqHfo0Gf9YdfDQEvu/Y+7dN/HTtwh3+sge3TePaVV74W11/yVbvnWOTtkHAAQTf8F79IhRxy9qt9jtZLSJZnZaZJOk6QnP/nJLXy9TA/P7brirzOtbp7bqMc+7Y+GHtv7N35L39r9aO3y8NaORtW+2+YOU37YCUOPZb2ertnnRD32gds7GhWAHWWnuc7CEtR+bFa242Lzd+Z+T3tseMnQY098ymG6brfnaadHf75DxgAAEe3cW/3YbO7t391lZs+U9Jfu/gfl+2dLkrv/v+Oev2nTJt+8eXPr4wAAzCYzu87dN3U9jpQRmwEAbVosNq9WTem3JB1oZvub2c6STpJ02Sp9LwAAAABAolZlL9bdt5nZ6yX9L0k9SR9z9++txvcCAAAAAKRr1Q4Iu/vnJX1+tb4+AAAAACB9q3V8FwAAAACAJZGUAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDMkpQAAAACAzpCUAgAAAAA6Q1IKAAAAAOgMSSkAAAAAoDMkpQAAAACAzpCUAgAAAAA6Y+7e9RhkZlsl3dXSl9tb0r0tfa1pxjxjYZ6xMM/u/Za7r+16ECkjNi8L84yFecbCPLu3YGyeiqS0TWa22d03dT2O1cY8Y2GesTBPYNis/FthnrEwz1iY53Tj+C4AAAAAoDMkpQAAAACAzkRMSs/regA7CPOMhXnGwjyBYbPyb4V5xsI8Y2GeUyxcTSkAAAAAIB0Rd0oBAAAAAIkIk5Sa2dFmdouZ3WZmZ3U9nraY2ZPM7KtmdpOZfc/M3lg+vqeZXWFmt5Z/79H1WNtgZj0zu97MLi/f39/Mvlm+rpeY2c5dj3GlzOwJZnapmX3fzG42s2dGfD3N7Mzy3+yNZvZxM1sT5fU0s4+Z2T1mdmPjsbGvoRU+VM75O2b29O5Gvn0WmOf7yn+73zGz/2lmT2h87OxynreY2R90MmhMFWJz+j/LJWJzpNeT2Exs7mTQEwiRlJpZT9LfSjpG0iGSTjazQ7odVWu2Sfpv7n6IpCMkva6c21mSvuzuB0r6cvl+BG+UdHPj/fdK+mt3P0DSzyS9ppNRtetvJH3R3Z8q6Wkq5hvq9TSz/SS9QdImdz9MUk/SSYrzel4g6eiRxxZ6DY+RdGD55zRJ5+6gMbbhAs2f5xWSDnP335H0b5LOlqTy59JJkg4tP+fD5c9mzChic/o/yxuIzQFeT2IzsVlTHJtDJKWSniHpNne/w90fkfQJScd1PKZWuPvd7v7t8u1fqPghuZ+K+V1YPu1CScd3MsAWmdk6SX8o6aPl+ybpeZIuLZ+S/DzNbHdJz5V0viS5+yPufr8Cvp6S5iQ9xszmJD1W0t0K8nq6+9ck/XTk4YVew+MkXeSFayQ9wcyeuEMGukLj5unuX3L3beW710haV759nKRPuPvD7v4DSbep+NmM2UVsTvRnXBOxOdbrKWLz8Y3Hic1TJEpSup+kHzbe31I+FoqZrZe0UdI3Je3r7neXH/qxpH27GleLPijpLZLy8v29JN3f+J8swuu6v6Stkv5HeRTqo2a2q4K9nu7+I0nvl/TvKgLezyVdp3ivZ9NCr2Hkn0+vlvSF8u3I88TyzMS/CWJziNeV2Bzr9WwiNicyzyhJaXhm9jhJn5J0hrs/0PyYFy2Uk26jbGbHSrrH3a/reiyrbE7S0yWd6+4bJf1KI8eBgryee6hYndtf0m9K2lXzj5qEFeE1XIqZvU3FEcaLux4L0BVicxjE5hkQ4TVcSsqxOUpS+iNJT2q8v658LAQz20lF0LvY3T9dPvyT6phB+fc9XY2vJc+S9BIzu1PFEa/nqajveEJ5xESK8bpukbTF3b9Zvn+pikAY7fV8gaQfuPtWd39U0qdVvMbRXs+mhV7DcD+fzOxPJB0r6eU+uFcs3DyxYqH/TRCbQ/0sJzbHej2biM2JzDNKUvotSQeW3cN2VlHQe1nHY2pFWbtxvqSb3f0DjQ9dJunU8u1TJX12R4+tTe5+truvc/f1Kl6/r7j7yyV9VdIJ5dMizPPHkn5oZgeVDz1f0k0K9nqqOBp0hJk9tvw3XM0z1Os5YqHX8DJJryw7/R0h6eeNo0TJMbOjVRzle4m7/7rxocsknWRmu5jZ/iqaR1zbxRgxNYjNif+MIzbHej1FbCY2T3NsdvcQfyS9SEW3qdslva3r8bQ4r2erOGrwHUk3lH9epKKm48uSbpX0L5L27HqsLc75SEmXl28/RcX/PLdJ+mdJu3Q9vhbmt0HS5vI1/YykPSK+npLeJen7km6U9A+Sdonyekr6uIp6nEdVrLC/ZqHXUJKp6EB6u6Tvquh62PkcVjDP21TUp1Q/jz7SeP7bynneIumYrsfPn+7/EJvT/1nemDOxOcDrSWwmNnc9/oX+WDlYAAAAAAB2uCjHdwEAAAAACSIpBQAAAAB0hqQUAAAAANAZklIAAAAAQGdISgEAAAAAnSEpBRJjZv/a9RgAAMAAsRlYGa6EAQAAAAB0hp1SIDFm9suuxwAAAAaIzcDKkJQCAAAAADpDUgoAAAAA6AxJKQAAAACgMySlAAAAAIDOkJQC6aFlNgAA04XYDKwASSmQEDPbS9JPux4HAAAoEJuBlSMpBRJhZr8p6WpJ7+96LAAAgNgMtMXcOW0AAAAAAOgGO6UAAAAAgM6QlAIAAAAAOkNSCgAAAADoDEkpAAAAAKAzJKUAAAAAgM6QlAIAAAAAOvP/AxVAbdTJLxFzAAAAAElFTkSuQmCC\n",
|
||
"text/plain": [
|
||
"<Figure size 1152x432 with 2 Axes>"
|
||
]
|
||
},
|
||
"metadata": {
|
||
"needs_background": "light"
|
||
},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"N = 128\n",
|
||
"\n",
|
||
"xk = 2*np.pi*np.arange(N)/N\n",
|
||
"yk0 = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
"yk1 = np.exp(5j*xk*2)\n",
|
||
"\n",
|
||
"fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(16,6))\n",
|
||
"\n",
|
||
"ax[0].set_xlabel(\"j\")\n",
|
||
"ax[1].set_xlabel(\"j\")\n",
|
||
"\n",
|
||
"ax[0].set_title(\"$y_k = e^{20ix_k} + e^{40ix_k}$\")\n",
|
||
"ax[0].plot(np.abs(DFT(yk0)), label=\"DFT\")\n",
|
||
"ax[0].plot(np.abs(np.fft.fft(yk0)), label=\"numpy.fft.fft\")\n",
|
||
"ax[0].legend(loc=\"upper right\")\n",
|
||
"\n",
|
||
"ax[1].set_title(\"$y_k = e^{i5x^2_k}$\")\n",
|
||
"ax[1].plot(np.abs(DFT(yk1)), label=\"DFT\")\n",
|
||
"ax[1].plot(np.abs(np.fft.fft(yk1)), label=\"numpy.fft.fft\")\n",
|
||
"ax[1].legend(loc=\"upper right\")\n",
|
||
"\n",
|
||
"# TODO: So the graphs overlap completely. Is this good enough?\n",
|
||
"# To make it more clear, we could mirror one of the graphs (multiply by -1),\n",
|
||
"# like what is often done in spectroscopy, or we could add the difference.\n",
|
||
"\n",
|
||
"fig.show()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "1674c10abab84f4d7c15bc7e5fea53b2",
|
||
"grade": false,
|
||
"grade_id": "cell-e04e5bc7ed412f64",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"### Task 3\n",
|
||
"\n",
|
||
"Analyze the evaluation-time scaling of your $\\text{DFT(yk)}$ function with the help of the timeit\n",
|
||
"module. Base your code on the following example:\n",
|
||
"\n",
|
||
"```python\n",
|
||
"import timeit\n",
|
||
"\n",
|
||
"tOut = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n",
|
||
"tMean = np.mean(tOut)\n",
|
||
"```\n",
|
||
"This example evaluates $\\text{DFT(yk)}$ 5 × 10 times and stores the resulting 5 evaluation times in tOut. Afterwards we calculate the mean value of these 5 repetitions. \n",
|
||
"Use this example to calculate and plot the evaluation time of your $\\text{DFT(yk)}$ function for $N = 2^2, 2^3, ..., 2^M$. Depending on your implementation you might be able to go up to $M = 10$. Be careful and increase M just step by step!"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "0e92c1c2d548a1c9b1476dd295415c2e",
|
||
"grade": true,
|
||
"grade_id": "cell-0ab81532ab86e322",
|
||
"locked": false,
|
||
"points": 4,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"M = 2 gives\n",
|
||
"tOut = [0.00035531632602214813, 0.00026446394622325897, 0.0002669980749487877, 0.000263310968875885, 0.0002625705674290657]\n",
|
||
"\n",
|
||
"M = 3 gives\n",
|
||
"tOut = [0.0002719983458518982, 0.0002478230744600296, 0.0002496456727385521, 0.00025005731731653214, 0.0008301911875605583]\n",
|
||
"\n",
|
||
"M = 4 gives\n",
|
||
"tOut = [0.000346149317920208, 0.0003256509080529213, 0.00032440759241580963, 0.00031386781483888626, 0.0003223838284611702]\n",
|
||
"\n",
|
||
"M = 5 gives\n",
|
||
"tOut = [0.0007905261591076851, 0.0005042999982833862, 0.0004924368113279343, 0.0004901541396975517, 0.0005300091579556465]\n",
|
||
"\n",
|
||
"M = 6 gives\n",
|
||
"tOut = [0.002741251140832901, 0.002402886748313904, 0.0024533523246645927, 0.0024026362225413322, 0.002384801395237446]\n",
|
||
"\n",
|
||
"M = 7 gives\n",
|
||
"tOut = [0.009145548567175865, 0.008851788938045502, 0.00881863571703434, 0.008819866925477982, 0.008821901865303516]\n",
|
||
"\n",
|
||
"M = 8 gives\n",
|
||
"tOut = [0.03434724546968937, 0.03383493982255459, 0.03385954722762108, 0.03491049725562334, 0.03388229105621576]\n",
|
||
"\n",
|
||
"M = 9 gives\n",
|
||
"tOut = [0.14674979075789452, 0.14636401552706957, 0.1461899448186159, 0.14625835418701172, 0.14628244005143642]\n",
|
||
"\n",
|
||
"M = 10 gives\n",
|
||
"tOut = [0.5750202471390367, 0.5727888783439994, 0.5727971633896232, 0.5725235631689429, 0.5941787445917726]\n",
|
||
"\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"for M in range(2, 10+1):\n",
|
||
" N = 2**M\n",
|
||
" xk = 2*np.pi*np.arange(N)/N\n",
|
||
" # Using the first equation for yk from the previous exercise.\n",
|
||
" yk = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
" tOut = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n",
|
||
" tMean = np.mean(tOut)\n",
|
||
" print(\"M =\", M, \"gives\")\n",
|
||
" print(\"tOut =\", tOut)\n",
|
||
" print()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "3063d5b408e133b1930c56fa45fed54e",
|
||
"grade": false,
|
||
"grade_id": "cell-73ca4de972164356",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"### Task 4\n",
|
||
"\n",
|
||
"A very simple FFT algorithm can be derived by the following separation of the sum from\n",
|
||
"above:\n",
|
||
"\n",
|
||
"\\begin{align}\n",
|
||
" \\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij \\frac{2\\pi k}{N}} \n",
|
||
" &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi 2k}{N}} \n",
|
||
" + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi (2k+1)}{N}}\\\\ \n",
|
||
" &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi k}{N/2}}\n",
|
||
" + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi k}{N/2}} e^{-ij \\frac{2\\pi}{N}}\\\\\n",
|
||
" &= \\beta^{\\text{even}}_j + \\beta^{\\text{odd}}_j e^{-ij \\frac{2\\pi}{N}}\n",
|
||
"\\end{align}\n",
|
||
"\n",
|
||
"where $\\beta^{\\text{even}}_j$ is the Fourier transform based on only even $k$ (or $x_k$) and $\\beta^{\\text{odd}}_j$ the Fourier transform based on only odd $k$. In case $N = 2^M$ this even/odd separation can be done again and again in a recursive way. \n",
|
||
"\n",
|
||
"Use the template below to implement a $\\text{FFT(yk)}$ function, making use of your $\\text{DFT(yk)}$ function from above. Make sure that you get the same results as before by comparing the results from $\\text{DFT(yk)}$\n",
|
||
"and $\\text{FFT(yk)}$ for both functions defined in task 2.\n",
|
||
"\n",
|
||
"```python\n",
|
||
"def FFT(yk):\n",
|
||
" \"\"\"Don't forget to write a docstring ...\n",
|
||
" \"\"\"\n",
|
||
" N = # ... get the length of yk\n",
|
||
" \n",
|
||
" assert # ... check if N is a power of 2. Hint: use the % (modulo) operator\n",
|
||
" \n",
|
||
" if(N <= 2):\n",
|
||
" return # ... call DFT with all yk points\n",
|
||
" \n",
|
||
" else:\n",
|
||
" betaEven = # ... call FFT but using just even yk points\n",
|
||
" betaOdd = # ... call FFT but using just odd yk points\n",
|
||
" \n",
|
||
" expTerms = np.exp(-1j * 2.0 * np.pi * np.arange(N) / N)\n",
|
||
" \n",
|
||
" # Remember : beta_j is periodic in j !\n",
|
||
" betaEvenFull = np.concatenate([betaEven, betaEven])\n",
|
||
" betaOddFull = np.concatenate([betaOdd, betaOdd])\n",
|
||
" \n",
|
||
" return betaEvenFull + expTerms * betaOddFull\n",
|
||
"```"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "fdad5c990a077e2bc2c58d4c4da46973",
|
||
"grade": true,
|
||
"grade_id": "cell-ce8233802d8ccb83",
|
||
"locked": false,
|
||
"points": 3,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"def FFT(yk):\n",
|
||
" \"\"\"\n",
|
||
" Return the fast fourier transform (FFT) of array yk by considering odd\n",
|
||
" and even points and making use of discrete fourier transforms (DFTs).\n",
|
||
" \"\"\"\n",
|
||
" \n",
|
||
" N = len(yk)\n",
|
||
"\n",
|
||
" # N should be a power of two\n",
|
||
" assert np.log2(N).is_integer()\n",
|
||
"\n",
|
||
" if(N <= 2):\n",
|
||
" return DFT(yk)\n",
|
||
"\n",
|
||
" else:\n",
|
||
" betaEven = FFT(yk[::2])\n",
|
||
" betaOdd = FFT(yk[1::2])\n",
|
||
"\n",
|
||
" expTerms = np.exp(-1j * 2.0 * np.pi * np.arange(N) / N)\n",
|
||
"\n",
|
||
" # Remember : beta_j is periodic in j !\n",
|
||
" betaEvenFull = np.concatenate([betaEven, betaEven])\n",
|
||
" betaOddFull = np.concatenate([betaOdd, betaOdd])\n",
|
||
"\n",
|
||
" return betaEvenFull + expTerms * betaOddFull"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "33b415b478341202e57888621063fc35",
|
||
"grade": true,
|
||
"grade_id": "cell-a27cf0cb147b31ae",
|
||
"locked": false,
|
||
"points": 1,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEaCAYAAADqqhd6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhjUlEQVR4nO3df5RkZZ3f8fenqgcGUeTXwAKjO+PKwaCocEbEGF0inhVc4rA5oKABdEk4niPrD/YcBdEgJpujceOvZDUhooCHRV38AUtEJSxEooI7KEF+iIwwQrPAtCAgKjJd95s/7r3V1dXdzExX9zz3ufV5nTNnum5VdX2nuudzn/o+z71XEYGZmbVLJ3UBZma29BzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQg53M7MWcrhb60k6XNIPJH1X0qWSVlTbPyrpeklfHNh2hqQD01ZcknSSpKmhbY2u2ZrD4W7j4D7gNRHxamATsF7SS4ADIuJVwE+B4wEi4r9FxF3JKq1I6gInUNZeb2t0zdYsDndrvYh4ICJ+V918CiiAfw58p9r2LeCVAJKuq/7eSdIPJa2U9GeSPruDyz4J+Luq1tqCNTegXmsYh7uNDUl/CPwJ8PfAHsDj1V2PAXtK2hvYDBARTwEXAJ8GTgb+YgfW2QXeCHx56K4Fa05ZrzXTROoCzHYESbsBXwTeGhFbJD0K7Fbd/WzgEeDFwE8GnvZ94DOUrZDp7XitDwFExIcWWe6/Ab4SEYWkwe1bq3lR9Vo7eeRuy0LSMyX1JO03sO1Fkh6Q9Kxles1/J+l2SY9JukrSPtX2CeBLwHkRcWf18O8Dr62+fh3wPQaCshoRfxL4IGXYImlXSe/fATUfDJwi6VvAgZI+vbWad1S9lg+Huy2LiHiCctLvsIHNHwH+U0T8+umeK+lKSY8u8OfKBZ7zfuDtwBuAVcD9wH+s7j4JeDnwwao//aaIuBl4SNL1wAuBrwKHALdI2hm4GHg38DHgTZJ2BV4G3CvpHElrtvc92daaI+J9EfEnEXE0cFdEvBPgaWq+Y0fUa3mRT/lry0XShcDPI+I/SHo1cCHwgoh4StIngIsj4sdL8Dr7APcAh0bEz6ptrwA+ExGHjvr9B17nvcARwBufru2xLW2ZHVHzttZr7eSRuy2nf2Rm5P6fgQ9WE38A/4xyxLkUjgJ2Bn5Yj/ApV5M8tkTfv7YH8D7grOE7Bj9tVPeftZVPGzui5gXrtfbzyN2WjaTDKZfznUnZCz40ql84Sf8nIv5Y0lnArsC/H7jvKuBVC3zb6yPimKHXeQdwZEScsEz/lO2yjSP3RtVs7eORuy2n/wf8AfBfgLMHwnsf4NeSvgDcGxEfjIFRRkQcExHPXODPMfO8zo+AfynpsOr77yZpvYaWmjRMjjVbRhzutmwi4veUq082RcRVA3cdQtmuuTsi/nYJXucHwIeBr0p6ArgdODoa/LE0x5otL27L2LKRtBOwkXJC74aB7e8GJoG3UY7ob0lToVl7eeRuy+lc4HuDwV45hLIt8VbgU5L23NGFmbWdR+625Ko+8rXALcCfRcQvE5dkNnYc7mZmLeS2jJlZCznczcxaqBFnhdx7771jzZo1qcswM8vKTTfd9MuIWDXffY0I9zVr1rBhw4bUZZiZZUXSLxa6z20ZM7MWcribmbWQw93MrIUa0XM3M1sKW7ZsYXJykieffDJ1KUtq5cqVrF69mhUrVmzzcxzuZtYak5OTPOtZz2LNmjW05QSbEcHDDz/M5OQka9eu3ebnuS1jZq3x5JNPstdee7Um2AEksddee233pxGHu5m1SpuCvbaYf5PDPaHHfvVLNvyv/5m6jEV7+KFJfvTtL6Yuw8zm4Z57QndeczGH33oeD687hr32XZ26nO1217f/B4f//L/y5Kv/NSt32TV1OWaN0O12OeSQQ/q3v/GNb7Bp0ybWr1/f75nvvffeHHTQQXzve9/jqaee4p577uGggw4C4AMf+ADHH3/8yHU43BOK6d8DsOWpPGf2o/cUHQW96S2pSzFrjF122YWbb7551rZNmzbxqle9iiuvnHut9E2bNnHsscfOec6o3JZJKIoeAEWvl7iSRarq7+Vav1mLeeSeUpShGMV04kIWqa6/l2n91mrn/f1t3P5Pjy/p9zx4/90491+98Gkf87vf/Y6XvvSlAKxdu5avf/3rAFx//fX97SeccALnnHPOktY2zOGeUn/knmk49kfumdZvtgzma8sAC7ZllovDPaHs2zLVyL0oMq3fWm1rI+y2c889pczDUbnvnMxazOGeUlGUf+W62iSq+t2WMWsct2USUjWRmu/Ivay/N+1wN6s98cQTc7YdeeSRHHnkkfM+fs2aNdx6661LXsdWR+6SPi9ps6RbB7Z9TNJPJd0i6euSdh+472xJGyXdKel1S15xi0R/tUye4T7TVnK4mzXNtrRlLgSOHtp2NfCiiHgx8DPgbABJBwMnAi+snvMZSd0lq7ZtiszbGv22TKY7J7MW22q4R8R3gUeGtn0nIupEugGoj51fD3wpIn4fEfcAG4HDl7DeVlHmI3flvk7frMWWYkL1z4Grqq8PAO4buG+y2mbzibxXm8gjd7PGGincJZ0DTAOXLOK5p0vaIGnD1NTUKGXkq8h85BuZH4Rl1mKLDndJbwWOBd4SEVFtvh94zsDDVlfb5oiI8yNiXUSsW7Vq1WLLyFr+bZly5J5r/WZttqilkJKOBt4L/HFE/HbgriuAv5X0cWB/4EDghyNX2VZ1OGba1pBH7mZz1Kf83bJlCxMTE5xyyim85z3vodPpcN111+2wU/9uNdwlXQocCewtaRI4l3J1zM7A1dUVQm6IiLdHxG2SvgLcTtmueUfU6/1sjv7IPfIMx9w/eZgth8Fzy2zevJk3v/nNPP7445x33nnAwueYWepT/2413CPipHk2X/A0j/8r4K9GKWps1IfvZ3oQkHI/8ZnZMttnn304//zzednLXsaHPvShHfraPkI1oXrkS6YfbkS1WibXCWFrt6vOggd/srTf8w8OgWM+sl1Ped7znkev12Pz5s3Ajjv1r8M9oZmedZG4ksXpt2Uyrd8shR116l+He0L91SbZ9tzr1TJ51m8tt50j7OVy99130+122Weffbjjjjt22Ov6rJAJ9dsyma+W8YSq2fympqZ4+9vfzhlnnEG1+GSH8cg9pczXiXcyr99sOdSX2auXQp588smceeaZO7wOh3tCM0sh8wzHekLV4W424+kuGL8jT/3rtkxCdc+aTMPRJw4zay6He0K596zdljFrLod7QnVbg0xHvm7LWBPNnOqqPRbzb3K4J9TJfuRer3PPc+dk7bNy5UoefvjhVgV8RPDwww+zcuXK7XqeJ1QT6i+FzD3cM63f2mf16tVMTk7SttOIr1y5ktWrV2/9gQMc7gnVPetcTz/QIe8JYWufFStW9M+4OO7clklo5gjPPA/fd8/drLkc7gl1KENRmYZj7p88zNrM4Z7QzLll8gzHeufktoxZ8zjcE8q9Z+117mbN5XBPKPe2Rn/nlGn9Zm3mcE9I/bZGnuvEc//kYdZmDveE6nDMdkLVq2XMGsvhnlCnf5m9PJdCdjP/5GHWZg73hHLvWc/Un+fOyazNHO4J1ROqubZlupnXb9ZmWw13SZ+XtFnSrQPb9pR0taS7qr/3qLZL0qclbZR0i6TDlrP43LVl5J7rOn2zNtuWkfuFwNFD284CromIA4FrqtsAxwAHVn9OBz67NGW2U+5tjW5/tUye9Zu12VbDPSK+CzwytHk9cFH19UXAcQPbL47SDcDukvZbolpbpw5HZTry7WRev1mbLbbnvm9EPFB9/SCwb/X1AcB9A4+brLbNIel0SRskbWjb6Tm3Vf/cMpmO3CeUd1vJrM1GnlCN8qz4231m/Ig4PyLWRcS6VatWjVpGlmYmVPNbSlgMXgTYE6pmjbPYcH+obrdUf2+utt8PPGfgcaurbTaPbsY9997A1ZfcljFrnsWG+xXAqdXXpwKXD2w/pVo1cwTw2ED7xobk3LMeDPccd05mbbfVKzFJuhQ4Ethb0iRwLvAR4CuSTgN+Abyxevg3gdcDG4HfAm9bhppbY2ZCNb9w7E1v6X+dY1vJrO22Gu4RcdICdx01z2MDeMeoRY2LnFfL9AZ77hnunMzazkeoJhJFQUflPHR9ubqchHvuZo3mcE8k9wnJ3Os3azuHeyKzwzG/kXsxsPwxx/rN2s7hnsjgOvEcR76512/Wdg73RAZH7p0MR76Fl0KaNZrDPZFe5iNfj9zNms3hnsjgapNOhqtlisznDMzazuGeyPTgQUAZjnxnt5Xyq9+s7RzuiQxeVDrHkXsUXgpp1mQO90SyXwo52HPPcOdk1nYO90QGwzHHtkbhg5jMGs3hnshgWyPPtszgzim/+s3azuGeyKyRb4bh7tUyZs3mcE+kbssUoSxHvvXIvQj1LxdoZs3hcE+kPjfLFiboZhiOg/V75G7WPA73ROqR+xYm8hy5D9afYVvJrO0c7onUE6pbtCLLtsas+jPcOZm1ncM9kXpCMteR7+z689s5mbWdwz2ROhyntSLLcK8nVKe1Ist1+mZt53BPZCYc8xy5D9af41JOs7ZzuCdST6j2yHvknmv9Zm3ncE+knpCc1gTdDCckB+v3hKpZ84wU7pLeI+k2SbdKulTSSklrJd0oaaOkL0vaaamKbZP+yDfznnuu9Zu13aLDXdIBwDuBdRHxIqALnAh8FPhERDwf+BVw2lIU2jb1OvFeJ8+DmPrh3slzzsCs7UZty0wAu0iaAJ4BPAC8Brisuv8i4LgRX6OVIsq2RqEVdHMMx/r0A7nWb9Zyiw73iLgf+GvgXspQfwy4CXg06uSCSeCA+Z4v6XRJGyRtmJqaWmwZ2YpeGYhFZwVdBVHkFZD9c8t08jwIy6ztRmnL7AGsB9YC+wO7Akdv6/Mj4vyIWBcR61atWrXYMrJV9NsaO1W38wp3Bur3hKpZ84zSlnktcE9ETEXEFuBrwCuB3as2DcBq4P4Ra2yl+gLZRRXug1dmykExUL/bMmbNM0q43wscIekZkgQcBdwOXAscXz3mVODy0Upsp7pzFZ0VwOzzo2ehOio1Ou65mzXRKD33GyknTn8E/KT6XucD7wPOlLQR2Au4YAnqbJ9qtUx0V1Q3Mwv3YqZ+r5Yxa56JrT9kYRFxLnDu0Oa7gcNH+b7joJ6QrEfuvV5ek5KD9XvkbtY8PkI1kajbGt2y5x65jdwH6u9kuNrHrO0c7qn0R755TqhmX79ZyzncE6nbGqp67sV0XuE4XL/D3axZHO6pFEMTqkVm4Riz6y8ymzMwazuHeyL9ke/EzkCGSyGrHntdv0fuZs3icE+lvnpRNaFa9DKbkOzXn+dqH7O2c7inUrVh+iP33NoyxTTT0QF1q5tbEhdkZoMc7onUSwc7E/XIPbORb1FQ0EGdMtzdljFrFod7Iur33Kt17pmN3BU9enSgCvd6DsHMmsHhnkh9EJNyHblHzyN3swZzuKfSXydeh3tm4Vj06Glm5J7dzsms5RzuqUSPIkSnW57eJ7dwVz1yV92Wyat+s7ZzuKdSlD3ruq2R3WqZmD2hmtvOyaztHO6pVCPfeuROZm2NckK1C123ZcyayOGeiOaM3DMLx6Juy1RtpdzqN2s5h3sqUYd7HY55tTX6Pfd6zsAHMZk1isM9ERU9CnVQ1daIzE4/oCjr7+T6ycOs5RzuqUSPgm4/HHNbbaIoKOiibvkr5IOYzJrF4Z5KFLN67rmFYz1y7/fcvVrGrFEc7oloaLVMbuFOFAQDbaXc6jdrOYd7Isp8nbiGTj/gpZBmzeJwT6WekKzXuefYc1e3v9ontzkDs7YbKdwl7S7pMkk/lXSHpFdI2lPS1ZLuqv7eY6mKbZN65N7NtC1T9ty7dNyWMWukUUfunwK+FREvAF4C3AGcBVwTEQcC11S3bYiil3XPWlQ993pC2G0Zs0ZZdLhLejbwauACgIh4KiIeBdYDF1UPuwg4brQS26k/8u3kOXLvDLWVItyWMWuSUUbua4Ep4AuSfizpc5J2BfaNiAeqxzwI7DvfkyWdLmmDpA1TU1MjlJEnFdU694zbMqGZ+otph7tZk4wS7hPAYcBnI+JQ4DcMtWAiIoCY78kRcX5ErIuIdatWrRqhjDyJohr5lm0NMgv3Tr0Usvrk0b9gtpk1wijhPglMRsSN1e3LKMP+IUn7AVR/bx6txHbKfZ17uVpmZudUZHb6BLO2W3S4R8SDwH2SDqo2HQXcDlwBnFptOxW4fKQKW0pREIMj98xGvh2qtkw9oeqeu1mjTIz4/L8ALpG0E3A38DbKHcZXJJ0G/AJ444iv0UoaPrdMdgcxDe2cvFrGrFFGCveIuBlYN89dR43yfcdBpwrHep17fiP3YtaEam5tJbO28xGqifQnVCdWlBsyC8dy5zRwhGpmOyeztnO4J1IvJexmulpG9Ga3ZTKr36ztHO6J1EsJZw4Cyisc67ZMrqdPMGs7h3siGgrH3Ea+M20Zj9zNmsjhnkinOv3AzIRqXuvEO1VbZqKaM/BZIc2axeGeSCdm96yV2ci3EwWo279AtkfuZs3icE+k7lkDTEcn0557Z2ZCOLP6zdrO4Z6IKEDl21/QgSK3tszwhGpe9Zu1ncM9kXpCEigvlF1sSVzR9ulWO6dc20pmbedwT6SekIRq5J7hhCqDI3e3ZcwaxeGeyKyRu7oos3DsREF0JrJdymnWdg73RAYnVHsZjtzrtow6HYqQJ1TNGsbhnkhneEI1s3DsUBCdgZ2TR+5mjeJwT2QwHAs62U1Iljunmfpz2zmZtZ3DPZFuNSEJeU6oTlBAZ3C1j8PdrEkc7okM9twLOllNqEZR0FH0d049j9zNGsfhnkg3Zka+hfIK91591ajB+j1yN2sUh3sis3vWeS2FrMNdA/V75G7WLA73RLqDbZnMRu5Fdb3UWatlMpszMGs7h3siHQro1EshuyijcOyP3Dt5zhmYjQOHeyITKmbaGspr5NurRu7MCvd86jcbBw73BIbbGkGHTkYj36gnVGftnPKp32wcjBzukrqSfizpyur2Wkk3Stoo6cuSdhq9zHbprzapj1BVG9oy+dRvNg6WYuT+LuCOgdsfBT4REc8HfgWctgSv0Sq96er0vgMj95x61sXQUsigg3yZPbNGGSncJa0G/hT4XHVbwGuAy6qHXAQcN8prtNHMyLc8o2KhTnnxjkwMj9x7mX3yMBsHo47cPwm8F/rJtBfwaETUw7hJ4ID5nijpdEkbJG2YmpoasYy8DE9IRman/C2G689snb7ZOFh0uEs6FtgcETct5vkRcX5ErIuIdatWrVpsGVkanpDMbkK1GDqIKbNPHmbjYGKE574SeIOk1wMrgd2ATwG7S5qoRu+rgftHL7Nd5kxIqsNERtcgLbzO3azxFj1yj4izI2J1RKwBTgT+ISLeAlwLHF897FTg8pGrbJmiGG7LTGQ18q3bMqovsSevljFrmuVY5/4+4ExJGyl78Bcsw2tkrd+zrtsyyqstU++cZkbu7rmbNc0obZm+iLgOuK76+m7g8KX4vm01p62hbpYjdzozI/eOR+5mjeIjVBMYXm1CZuFYT6h2Op5QNWsqh3sCwyP3UKc8kVgmZnru5a9Pbks5zcaBwz2B4dUyoW5ePfdeeYStVLVlMlvKaTYOHO4J9NeJ91eb5NVzj3pCtTszZ5DTJw+zceBwT6Df1si9LTPrCNt86jcbBw73BIqhc8ugLt2M2hr9Tx6zVsvkU7/ZOHC4JxBF3iP3uv5Od2DOIKP6zcaBwz2BOatlOhNZhWMMHaGKl0KaNY7DPYHhkTuZjXzrk352OjMTwjmt0zcbBw73BGYO358Z+XZzCvdeWetMW6ZDB/fczZrE4Z5AzHMQUE4j3+GdU271m40Dh3sCM4fvz4zcc2rLUNfvCVWzxnK4J1BPqHb6E6rdrNoyxdC5Zcq2ktsyZk3icE9g5gjPmXXuWY18e0NLITNb7WM2DhzuCcysE6/CvdNlhfIZ+c6pP7e2ktkYcLgnMN9SSBg4FXDDRcwzoepwN2sUh3sC/QnVgZE7zJwtsvGqnVN3Yqb+rlfLmDWKwz2B4cP3cwv3uv5uJ9M5A7Mx4HBPIIbWiavflskj3Bk65W8or9U+ZuPA4Z7CgiP3THru9ci931byhKpZ0zjcE+i3ZTqze+65TKgSs5dySl0m5HA3axKHewr9kXv19mfalul2Zw7CKjdnsnMyGwOLDndJz5F0raTbJd0m6V3V9j0lXS3prurvPZau3HboH6Faj3xznVAdWOcOMD29JVVJZjZklJH7NPCXEXEwcATwDkkHA2cB10TEgcA11W0bFPVSwhXl7U5eI3fNmTMoQz6X+s3GwaLDPSIeiIgfVV//GrgDOABYD1xUPewi4LgRa2yffs+9DEdlFu71QUzdTD95mI2DJem5S1oDHArcCOwbEQ9Udz0I7LvAc06XtEHShqmpqaUoIxtzV5vUI99MetZz6s9rtY/ZOBg53CU9E/gq8O6IeHzwvogIIOZ7XkScHxHrImLdqlWrRi0jL8OrTTIbuRM9eiHUmT0hHLnUbzYGRgp3SSsog/2SiPhatfkhSftV9+8HbB6txBYaWm3SD/cik3AsevQGfnXcljFrnlFWywi4ALgjIj4+cNcVwKnV16cCly++vHYabsv0L5SdS1sjehSDvzr9nVMm9ZuNgYkRnvtK4GTgJ5Jurra9H/gI8BVJpwG/AN44UoVtFLNPmZvbyF3FULhndlZLs3Gw6HCPiP8LaIG7j1rs9x0LRXk0Z7YTqlHM25bJZs7AbAz4CNUE1L9MXaf6O6+2hoppCs1ty/SmHe5mTeFwTyCGVpuoOg1BkcsRntGjoNu/mVtbyWwcONxTKAp6g+Go6opGmYzc57Rlupm1lczGgMM9AcXQUsJuZm2ZodUy/dU+HrmbNYbDPYXhcKxGvrkshVQU87dlMqnfbBw43FMoevQGJiTr87pnM/KN3qwJVfnEYWaN43BPYLit0em3NfIY+ZYj9/naMnnUbzYOHO4pDIUjOfbc5XXuZk3mcE9g7sg9r9UyC0+o5lG/2ThwuKcwdPh+fdGLXM6qqCiIWROqVc/d4W7WGA73BDTnIKDqNASRRzjOacvUB2H1MjkIy2wMONwTGA7HeuSeS896wbZSr0hVkpkNcbgnoJh9hGqOq2Vi1s4ps6WcZmPA4Z7AnJF7faHsbMK9R6GZE4p6QtWseRzuKURBZLxapjO0lHNm5J5H/WbjwOGewJyedTUhmc3IndltGY/czZrH4Z6AoqDQQM89s5Gvoueeu1nDOdwTGA7Hbo5tmVk7J4/czZrG4Z6AGDo3SxWO2axzZ/acgTLbOZmNA4d7Ap3oEQMj325mbZnh+uuRey5zBmbjwOGewJxzs9QXys6kZ90Z7rn7lL9mjeNwT6AzdBBQN7ORb4diaOSe12ofs3GwbOEu6WhJd0raKOms5XqdHA2vlpmoDmLKpS2joXDvTuwE5FO/2ThYlnCX1AX+BjgGOBg4SdLBy/FaOeow3LPO68Rh5SePeXrumdRvNg4mtv6QRTkc2BgRdwNI+hKwHrh9KV/kluu+ym7fPXcpv+UO8dzeg9y54tD+7XpC9aBNl7Dpw1emKmub7V88xEN6Yf92fW6cP7r9b9j04YtTlWWWpQf/6ASOeMvS59hyhfsBwH0DtyeBlw8+QNLpwOkAz33ucxf1Ijvt+mweecbaRZaYziOsRS9+U/92d2KCG1afxk6P3pWwqm33CGtZ+bKT+7f32Hs/btj3RHb6zT8lrMosTxPP2ndZvq8iYum/qXQ8cHRE/Nvq9snAyyPijPkev27dutiwYcOS12Fm1maSboqIdfPdt1wTqvcDzxm4vbraZmZmO8Byhfs/AgdKWitpJ+BE4Iplei0zMxuyLD33iJiWdAbwbaALfD4ibluO1zIzs7mWa0KViPgm8M3l+v5mZrYwH6FqZtZCDnczsxZyuJuZtZDD3cyshZblIKbtLkKaAn6xyKfvDfxyCcvZ0Vx/Wq4/Ldc/mj+MiFXz3dGIcB+FpA0LHaGVA9eflutPy/UvH7dlzMxayOFuZtZCbQj381MXMCLXn5brT8v1L5Pse+5mZjZXG0buZmY2xOFuZtZCWYd7bhfhlvQcSddKul3SbZLeVW3fU9LVku6q/t4jda0LkdSV9GNJV1a310q6sfoZfLk6xXNjSdpd0mWSfirpDkmvyOz9f0/1u3OrpEslrWzyz0DS5yVtlnTrwLZ532+VPl39O26RdFi6yvu1zlf/x6rfn1skfV3S7gP3nV3Vf6ek1yUpupJtuGd6Ee5p4C8j4mDgCOAdVc1nAddExIHANdXtpnoXcMfA7Y8Cn4iI5wO/Ak5LUtW2+xTwrYh4AfASyn9LFu+/pAOAdwLrIuJFlKfTPpFm/wwuBI4e2rbQ+30McGD153TgszuoxqdzIXPrvxp4UUS8GPgZcDZA9X/5ROCF1XM+U+VUEtmGOwMX4Y6Ip4D6ItyNFREPRMSPqq9/TRksB1DWfVH1sIuA45IUuBWSVgN/Cnyuui3gNcBl1UMaWzuApGcDrwYuAIiIpyLiUTJ5/ysTwC6SJoBnAA/Q4J9BRHwXeGRo80Lv93rg4ijdAOwuab8dUugC5qs/Ir4TEdPVzRsorzQHZf1fiojfR8Q9wEbKnEoi53Cf7yLcBySqZbtJWgMcCtwI7BsRD1R3PQgszxVzR/dJ4L1AUd3eC3h04Be96T+DtcAU8IWqtfQ5SbuSyfsfEfcDfw3cSxnqjwE3kdfPABZ+v3P8P/3nwFXV142qP+dwz5akZwJfBd4dEY8P3hfl2tTGrU+VdCywOSJuSl3LCCaAw4DPRsShwG8YasE09f0HqHrT6yl3UvsDuzK3ZZCVJr/fWyPpHMpW6yWpa5lPzuGe5UW4Ja2gDPZLIuJr1eaH6o+f1d+bU9X3NF4JvEHSJsoW2Gso+9e7Vy0CaP7PYBKYjIgbq9uXUYZ9Du8/wGuBeyJiKiK2AF+j/Lnk9DOAhd/vbP5PS3orcCzwlpg5WKhR9ecc7tldhLvqUV8A3BERHx+46wrg1OrrU4HLd3RtWxMRZ0fE6ohYQ/le/0NEvAW4Fji+elgja69FxIPAfZIOqjYdBdxOBu9/5V7gCEnPqH6X6vqz+RlUFnq/rwBOqVbNHAE8NtC+aQxJR1O2J98QEb8duOsK4ERJO0taSzkx/MMUNQIQEdn+AV5POVv9c+Cc1PVsQ73/gvIj6C3AzdWf11P2rq8B7gL+N7Bn6lq38u84Eriy+vp5lL/AG4G/A3ZOXd9Wan8psKH6GXwD2COn9x84D/gpcCvwRWDnJv8MgEsp5we2UH5yOm2h9xsQ5Qq4nwM/oVwV1MT6N1L21uv/w/994PHnVPXfCRyTsnaffsDMrIVybsuYmdkCHO5mZi3kcDczayGHu5lZCznczcxayOFuthWSvp+6BrPt5aWQZmYt5JG72VZIeiJ1DWbby+FuZtZCDnczsxZyuJuZtZDD3cyshRzuZlvnJWWWHYe72dOQtBdzrwFq1ngOd7MFSNof+AHldUvNsuKDmMzMWsgjdzOzFnK4m5m1kMPdzKyFHO5mZi3kcDcza6H/Dw49JMHj2uG3AAAAAElFTkSuQmCC\n",
|
||
"text/plain": [
|
||
"<Figure size 432x288 with 1 Axes>"
|
||
]
|
||
},
|
||
"metadata": {
|
||
"needs_background": "light"
|
||
},
|
||
"output_type": "display_data"
|
||
}
|
||
],
|
||
"source": [
|
||
"# N needs to be a power of two for FFT to work.\n",
|
||
"N = 2**7\n",
|
||
"\n",
|
||
"xk = 2*np.pi*np.arange(N)/N\n",
|
||
"yk = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
"\n",
|
||
"fig, ax = plt.subplots()\n",
|
||
"\n",
|
||
"ax.set_xlabel(\"j\")\n",
|
||
"\n",
|
||
"ax.set_title(\"$y_k = e^{20ix_k} + e^{40ix_k}$\")\n",
|
||
"ax.plot(np.abs(FFT(yk)), label=\"FFT\")\n",
|
||
"ax.plot(np.abs(DFT(yk)), label=\"DFT\")\n",
|
||
"ax.legend(loc=\"upper right\")\n",
|
||
"\n",
|
||
"fig.show()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"deletable": false,
|
||
"editable": false,
|
||
"nbgrader": {
|
||
"cell_type": "markdown",
|
||
"checksum": "c1d2dd0c2543e228d48d84e94720a3a8",
|
||
"grade": false,
|
||
"grade_id": "cell-c61537808f7cce68",
|
||
"locked": true,
|
||
"schema_version": 3,
|
||
"solution": false,
|
||
"task": false
|
||
}
|
||
},
|
||
"source": [
|
||
"### Task 5\n",
|
||
"\n",
|
||
"Analyze the evaluation-time scaling of your $\\text{FFT(yk)}$ function with the help of the timeit module and compare it to the scaling of the $\\text{DFT(yk)}$ function."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {
|
||
"deletable": false,
|
||
"nbgrader": {
|
||
"cell_type": "code",
|
||
"checksum": "bb4f1eee977aad469245f60fca596938",
|
||
"grade": true,
|
||
"grade_id": "cell-aaf90559928426bf",
|
||
"locked": false,
|
||
"points": 1,
|
||
"schema_version": 3,
|
||
"solution": true,
|
||
"task": false
|
||
}
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"M = 2 gives\n",
|
||
"tOutDFT = [0.00021187309175729752, 0.00018186774104833603, 0.0001812456175684929, 0.00018110498785972595, 0.00018095411360263824]\n",
|
||
"tOutFFT = [0.0005802353844046593, 0.0005557788535952568, 0.0006145015358924866, 0.0005566608160734177, 0.0005704071372747421]\n",
|
||
"\n",
|
||
"M = 3 gives\n",
|
||
"tOutDFT = [0.00022258423268795013, 0.00020575150847434998, 0.0002390751615166664, 0.00020509026944637299, 0.00020475033670663834]\n",
|
||
"tOutFFT = [0.001284896396100521, 0.0013337591663002968, 0.0012675542384386063, 0.0013139024376869202, 0.0012842761352658272]\n",
|
||
"\n",
|
||
"M = 4 gives\n",
|
||
"tOutDFT = [0.0003250986337661743, 0.00035554729402065277, 0.0003026863560080528, 0.00030192453414201736, 0.0003027757629752159]\n",
|
||
"tOutFFT = [0.00277732964605093, 0.002760577015578747, 0.0021515777334570885, 0.002049895003437996, 0.0028353100642561913]\n",
|
||
"\n",
|
||
"M = 5 gives\n",
|
||
"tOutDFT = [0.0008856169879436493, 0.0004905443638563156, 0.0004913564771413803, 0.0005037290975451469, 0.0004903236404061317]\n",
|
||
"tOutFFT = [0.004227722063660622, 0.004218194633722305, 0.004212621599435806, 0.006331468001008034, 0.006671415641903877]\n",
|
||
"\n",
|
||
"M = 6 gives\n",
|
||
"tOutDFT = [0.0037633972242474556, 0.002388077788054943, 0.0024485038593411446, 0.0023846719413995743, 0.0024057207629084587]\n",
|
||
"tOutFFT = [0.01122717373073101, 0.011149754747748375, 0.011180805042386055, 0.01120598241686821, 0.011219387874007225]\n",
|
||
"\n",
|
||
"M = 7 gives\n",
|
||
"tOutDFT = [0.00920125376433134, 0.00882750190794468, 0.008815147913992405, 0.00881881546229124, 0.008816582150757313]\n",
|
||
"tOutFFT = [0.022676337510347366, 0.023197690956294537, 0.02264594007283449, 0.022664004936814308, 0.022718658670783043]\n",
|
||
"\n",
|
||
"M = 8 gives\n",
|
||
"tOutDFT = [0.034523812122642994, 0.033904923126101494, 0.033997087739408016, 0.03390498273074627, 0.03390084486454725]\n",
|
||
"tOutFFT = [0.046751600690186024, 0.04573535453528166, 0.04478597640991211, 0.03446220513433218, 0.034463477320969105]\n",
|
||
"\n",
|
||
"M = 9 gives\n",
|
||
"tOutDFT = [0.12952119670808315, 0.1336375381797552, 0.13360701967030764, 0.1335998559370637, 0.13358618132770061]\n",
|
||
"tOutFFT = [0.09165043104439974, 0.07926731836050749, 0.06876837275922298, 0.06870093382894993, 0.06888525560498238]\n",
|
||
"\n",
|
||
"M = 10 gives\n",
|
||
"tOutDFT = [0.5503582386299968, 0.5721757095307112, 0.572254091501236, 0.572665935382247, 0.5729674575850368]\n",
|
||
"tOutFFT = [0.17185171879827976, 0.1377342501655221, 0.13842287193983793, 0.1409124033525586, 0.13889379799365997]\n",
|
||
"\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"for M in range(2, 10+1):\n",
|
||
" N = 2**M\n",
|
||
" xk = 2*np.pi*np.arange(N)/N\n",
|
||
" # Using the first equation for yk from the second exercise.\n",
|
||
" yk = np.exp(20j*xk) + np.exp(40j*xk)\n",
|
||
" tOutDFT = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n",
|
||
" tOutFFT = timeit.repeat(stmt=lambda: FFT(yk), number=10, repeat=5)\n",
|
||
" tMean = np.mean(tOut)\n",
|
||
" print(\"M =\", M, \"gives\")\n",
|
||
" print(\"tOutDFT =\", tOutDFT)\n",
|
||
" print(\"tOutFFT =\", tOutFFT)\n",
|
||
" print()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"For small $M$, `DFT` is faster, but as $M$ increases, `FFT` gets a lot more efficient."
|
||
]
|
||
}
|
||
],
|
||
"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.8.10"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 4
|
||
}
|