From 4f0de19fd9617494d8d9711ecd3b8cd6bb398435 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Tue, 15 Feb 2022 14:58:19 +0100 Subject: [PATCH] 05: Finish all tasks except for some TODOs --- ...Discrete and Fast Fourier Transforms.ipynb | 149 +++++++++++++++--- 1 file changed, 124 insertions(+), 25 deletions(-) diff --git a/Week 2/5 Discrete and Fast Fourier Transforms.ipynb b/Week 2/5 Discrete and Fast Fourier Transforms.ipynb index d8945f4..24de6c1 100644 --- a/Week 2/5 Discrete and Fast Fourier Transforms.ipynb +++ b/Week 2/5 Discrete and Fast Fourier Transforms.ipynb @@ -137,6 +137,7 @@ }, "outputs": [], "source": [ + "# TODO: Write docstrings, ugh.\n", "def DFT(yk):\n", " N = len(yk)\n", " xk = 2*np.pi*np.arange(N)/N\n", @@ -292,31 +293,31 @@ "output_type": "stream", "text": [ "M = 2 gives\n", - "tOut = [0.0010804971680045128, 0.0005500782281160355, 0.0005198698490858078, 0.0005232971161603928, 0.0006225360557436943]\n", + "tOut = [0.0003551868721842766, 0.00026587769389152527, 0.0002683410421013832, 0.0002652248367667198, 0.00029958970844745636]\n", "\n", "M = 3 gives\n", - "tOut = [0.0005145492032170296, 0.00047505367547273636, 0.000535009428858757, 0.0003488045185804367, 0.00028009340167045593]\n", + "tOut = [0.00037261005491018295, 0.0003149900585412979, 0.0002723895013332367, 0.00024654995650053024, 0.000258052721619606]\n", "\n", "M = 4 gives\n", - "tOut = [0.0002853143960237503, 0.00024075806140899658, 0.00021690316498279572, 0.00021761376410722733, 0.00023322459310293198]\n", + "tOut = [0.00034819450229406357, 0.00033616088330745697, 0.0003195982426404953, 0.00034148991107940674, 0.00035313330590724945]\n", "\n", "M = 5 gives\n", - "tOut = [0.000910954549908638, 0.0005088187754154205, 0.0005268435925245285, 0.0005177054554224014, 0.0005334652960300446]\n", + "tOut = [0.0007255515083670616, 0.0005129771307110786, 0.0004988498985767365, 0.0004972768947482109, 0.0005336767062544823]\n", "\n", "M = 6 gives\n", - "tOut = [0.002634277567267418, 0.0023930883035063744, 0.0024030981585383415, 0.002396835945546627, 0.0023911641910672188]\n", + "tOut = [0.0025200899690389633, 0.0023925071582198143, 0.0023851143196225166, 0.0023938296362757683, 0.0023842724040150642]\n", "\n", "M = 7 gives\n", - "tOut = [0.0092502785846591, 0.008834738284349442, 0.00884521659463644, 0.008864735253155231, 0.008844294585287571]\n", + "tOut = [0.008896775543689728, 0.00880778580904007, 0.008781216107308865, 0.008825629949569702, 0.008792318403720856]\n", "\n", "M = 8 gives\n", - "tOut = [0.034397597424685955, 0.03395726904273033, 0.03396071586757898, 0.03401851560920477, 0.03401240427047014]\n", + "tOut = [0.03435589000582695, 0.033995422534644604, 0.034102585166692734, 0.03416615817695856, 0.03535666409879923]\n", "\n", "M = 9 gives\n", - "tOut = [0.14702877961099148, 0.148259031586349, 0.14664556831121445, 0.14639647398144007, 0.14646969363093376]\n", + "tOut = [0.1477726474404335, 0.14786271005868912, 0.14967191498726606, 0.1493966607376933, 0.14906080160290003]\n", "\n", "M = 10 gives\n", - "tOut = [0.5762987565249205, 0.5738314474001527, 0.5741532389074564, 0.5737893972545862, 0.574077476747334]\n", + "tOut = [0.5873688040301204, 0.5887598115950823, 0.5843716822564602, 0.5818474553525448, 0.5803691502660513]\n", "\n" ] } @@ -397,7 +398,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": { "deletable": false, "nbgrader": { @@ -415,13 +416,33 @@ "outputs": [], "source": [ "def FFT(yk):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()" + " # TODO: Write a docstring.\n", + " \"\"\"Don't forget to write a docstring ...\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": null, + "execution_count": 6, "metadata": { "deletable": false, "nbgrader": { @@ -436,12 +457,37 @@ "task": false } }, - "outputs": [], + "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": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "# Do your plotting here ...\n", + "# N needs to be a power of two for FFT to work.\n", + "N = 2**7\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "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()" ] }, { @@ -468,7 +514,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": { "deletable": false, "nbgrader": { @@ -483,18 +529,71 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "M = 2 gives\n", + "tOutDFT = [0.00026126671582460403, 0.00021027959883213043, 0.0001835385337471962, 0.00023766234517097473, 0.00022200308740139008]\n", + "tOutFFT = [0.0006446670740842819, 0.0007413318380713463, 0.0006805472075939178, 0.0006410712376236916, 0.0006419029086828232]\n", + "\n", + "M = 3 gives\n", + "tOutDFT = [0.00022369623184204102, 0.00022925715893507004, 0.0002404283732175827, 0.00021736416965723038, 0.0002731001004576683]\n", + "tOutFFT = [0.0013865511864423752, 0.0013830652460455894, 0.0013337815180420876, 0.0012838244438171387, 0.0012847669422626495]\n", + "\n", + "M = 4 gives\n", + "tOutDFT = [0.0003327140584588051, 0.0003035273402929306, 0.0003030272200703621, 0.00030216481536626816, 0.00030276738107204437]\n", + "tOutFFT = [0.002759365364909172, 0.002286495640873909, 0.0020805923268198967, 0.0020598340779542923, 0.0020735692232847214]\n", + "\n", + "M = 5 gives\n", + "tOutDFT = [0.0008888524025678635, 0.000489652156829834, 0.00048649683594703674, 0.00048589520156383514, 0.0004848325625061989]\n", + "tOutFFT = [0.004292245022952557, 0.00426078587770462, 0.004242740571498871, 0.005865932442247868, 0.006726261228322983]\n", + "\n", + "M = 6 gives\n", + "tOutDFT = [0.0036992961540818214, 0.002382858656346798, 0.002396935597062111, 0.002390684559941292, 0.002375013194978237]\n", + "tOutFFT = [0.011373533867299557, 0.011346021667122841, 0.01132622268050909, 0.011333126574754715, 0.011340481229126453]\n", + "\n", + "M = 7 gives\n", + "tOutDFT = [0.009215572848916054, 0.008818225935101509, 0.008802506141364574, 0.008820069953799248, 0.00881765503436327]\n", + "tOutFFT = [0.022970101796090603, 0.023467930033802986, 0.02285583410412073, 0.02295252773910761, 0.022949432022869587]\n", + "\n", + "M = 8 gives\n", + "tOutDFT = [0.03451558295637369, 0.033903918229043484, 0.03388609457761049, 0.033958752639591694, 0.033889370039105415]\n", + "tOutFFT = [0.047724733129143715, 0.04677246883511543, 0.04511415120214224, 0.035544090904295444, 0.03561900556087494]\n", + "\n", + "M = 9 gives\n", + "tOutDFT = [0.12958059646189213, 0.13343970850110054, 0.13339493237435818, 0.1336062252521515, 0.13355298340320587]\n", + "tOutFFT = [0.09326057601720095, 0.08018506038933992, 0.07027726992964745, 0.07104555331170559, 0.0710984542965889]\n", + "\n", + "M = 10 gives\n", + "tOutDFT = [0.5566168483346701, 0.5795272067189217, 0.5755985928699374, 0.5750947641208768, 0.5745749343186617]\n", + "tOutFFT = [0.1753099588677287, 0.14073420222848654, 0.14395484700798988, 0.14109977800399065, 0.14102207031100988]\n", + "\n" + ] + } + ], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "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": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "For small $M$, `DFT` is faster, but as $M$ increases, `FFT` gets a lot more efficient." + ] } ], "metadata": {