diff options
Diffstat (limited to 'circuitpython/extmod/ulab/docs/numpy-fft.ipynb')
| -rw-r--r-- | circuitpython/extmod/ulab/docs/numpy-fft.ipynb | 546 |
1 files changed, 546 insertions, 0 deletions
diff --git a/circuitpython/extmod/ulab/docs/numpy-fft.ipynb b/circuitpython/extmod/ulab/docs/numpy-fft.ipynb new file mode 100644 index 0000000..803c923 --- /dev/null +++ b/circuitpython/extmod/ulab/docs/numpy-fft.ipynb @@ -0,0 +1,546 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-01T09:27:13.438054Z", + "start_time": "2020-05-01T09:27:13.191491Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Notebook magic" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2022-01-07T18:24:48.499467Z", + "start_time": "2022-01-07T18:24:48.488004Z" + } + }, + "outputs": [], + "source": [ + "from IPython.core.magic import Magics, magics_class, line_cell_magic\n", + "from IPython.core.magic import cell_magic, register_cell_magic, register_line_magic\n", + "from IPython.core.magic_arguments import argument, magic_arguments, parse_argstring\n", + "import subprocess\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2020-07-23T20:31:25.296014Z", + "start_time": "2020-07-23T20:31:25.265937Z" + } + }, + "outputs": [], + "source": [ + "@magics_class\n", + "class PyboardMagic(Magics):\n", + " @cell_magic\n", + " @magic_arguments()\n", + " @argument('-skip')\n", + " @argument('-unix')\n", + " @argument('-pyboard')\n", + " @argument('-file')\n", + " @argument('-data')\n", + " @argument('-time')\n", + " @argument('-memory')\n", + " def micropython(self, line='', cell=None):\n", + " args = parse_argstring(self.micropython, line)\n", + " if args.skip: # doesn't care about the cell's content\n", + " print('skipped execution')\n", + " return None # do not parse the rest\n", + " if args.unix: # tests the code on the unix port. Note that this works on unix only\n", + " with open('/dev/shm/micropython.py', 'w') as fout:\n", + " fout.write(cell)\n", + " proc = subprocess.Popen([\"../../micropython/ports/unix/micropython\", \"/dev/shm/micropython.py\"], \n", + " stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n", + " print(proc.stdout.read().decode(\"utf-8\"))\n", + " print(proc.stderr.read().decode(\"utf-8\"))\n", + " return None\n", + " if args.file: # can be used to copy the cell content onto the pyboard's flash\n", + " spaces = \" \"\n", + " try:\n", + " with open(args.file, 'w') as fout:\n", + " fout.write(cell.replace('\\t', spaces))\n", + " printf('written cell to {}'.format(args.file))\n", + " except:\n", + " print('Failed to write to disc!')\n", + " return None # do not parse the rest\n", + " if args.data: # can be used to load data from the pyboard directly into kernel space\n", + " message = pyb.exec(cell)\n", + " if len(message) == 0:\n", + " print('pyboard >>>')\n", + " else:\n", + " print(message.decode('utf-8'))\n", + " # register new variable in user namespace\n", + " self.shell.user_ns[args.data] = string_to_matrix(message.decode(\"utf-8\"))\n", + " \n", + " if args.time: # measures the time of executions\n", + " pyb.exec('import utime')\n", + " message = pyb.exec('t = utime.ticks_us()\\n' + cell + '\\ndelta = utime.ticks_diff(utime.ticks_us(), t)' + \n", + " \"\\nprint('execution time: {:d} us'.format(delta))\")\n", + " print(message.decode('utf-8'))\n", + " \n", + " if args.memory: # prints out memory information \n", + " message = pyb.exec('from micropython import mem_info\\nprint(mem_info())\\n')\n", + " print(\"memory before execution:\\n========================\\n\", message.decode('utf-8'))\n", + " message = pyb.exec(cell)\n", + " print(\">>> \", message.decode('utf-8'))\n", + " message = pyb.exec('print(mem_info())')\n", + " print(\"memory after execution:\\n========================\\n\", message.decode('utf-8'))\n", + "\n", + " if args.pyboard:\n", + " message = pyb.exec(cell)\n", + " print(message.decode('utf-8'))\n", + "\n", + "ip = get_ipython()\n", + "ip.register_magics(PyboardMagic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## pyboard" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-07T07:35:35.126401Z", + "start_time": "2020-05-07T07:35:35.105824Z" + } + }, + "outputs": [], + "source": [ + "import pyboard\n", + "pyb = pyboard.Pyboard('/dev/ttyACM0')\n", + "pyb.enter_raw_repl()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-19T19:11:18.145548Z", + "start_time": "2020-05-19T19:11:18.137468Z" + } + }, + "outputs": [], + "source": [ + "pyb.exit_raw_repl()\n", + "pyb.close()" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": { + "ExecuteTime": { + "end_time": "2020-05-07T07:35:38.725924Z", + "start_time": "2020-05-07T07:35:38.645488Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "%%micropython -pyboard 1\n", + "\n", + "import utime\n", + "import ulab as np\n", + "\n", + "def timeit(n=1000):\n", + " def wrapper(f, *args, **kwargs):\n", + " func_name = str(f).split(' ')[1]\n", + " def new_func(*args, **kwargs):\n", + " run_times = np.zeros(n, dtype=np.uint16)\n", + " for i in range(n):\n", + " t = utime.ticks_us()\n", + " result = f(*args, **kwargs)\n", + " run_times[i] = utime.ticks_diff(utime.ticks_us(), t)\n", + " print('{}() execution times based on {} cycles'.format(func_name, n, (delta2-delta1)/n))\n", + " print('\\tbest: %d us'%np.min(run_times))\n", + " print('\\tworst: %d us'%np.max(run_times))\n", + " print('\\taverage: %d us'%np.mean(run_times))\n", + " print('\\tdeviation: +/-%.3f us'%np.std(run_times)) \n", + " return result\n", + " return new_func\n", + " return wrapper\n", + "\n", + "def timeit(f, *args, **kwargs):\n", + " func_name = str(f).split(' ')[1]\n", + " def new_func(*args, **kwargs):\n", + " t = utime.ticks_us()\n", + " result = f(*args, **kwargs)\n", + " print('execution time: ', utime.ticks_diff(utime.ticks_us(), t), ' us')\n", + " return result\n", + " return new_func" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__END_OF_DEFS__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# numpy.fft\n", + "\n", + "Functions related to Fourier transforms can be called by prepending them with `numpy.fft.`. The module defines the following two functions:\n", + "\n", + "1. [numpy.fft.fft](#fft)\n", + "1. [numpy.fft.ifft](#ifft)\n", + "\n", + "`numpy`: https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.ifft.html\n", + "\n", + "## fft\n", + "\n", + "Since `ulab`'s `ndarray` does not support complex numbers, the invocation of the Fourier transform differs from that in `numpy`. In `numpy`, you can simply pass an array or iterable to the function, and it will be treated as a complex array:" + ] + }, + { + "cell_type": "code", + "execution_count": 341, + "metadata": { + "ExecuteTime": { + "end_time": "2019-10-17T17:33:38.487729Z", + "start_time": "2019-10-17T17:33:38.473515Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([20.+0.j, 0.+0.j, -4.+4.j, 0.+0.j, -4.+0.j, 0.+0.j, -4.-4.j,\n", + " 0.+0.j])" + ] + }, + "execution_count": 341, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fft.fft([1, 2, 3, 4, 1, 2, 3, 4])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**WARNING:** The array returned is also complex, i.e., the real and imaginary components are cast together. In `ulab`, the real and imaginary parts are treated separately: you have to pass two `ndarray`s to the function, although, the second argument is optional, in which case the imaginary part is assumed to be zero.\n", + "\n", + "**WARNING:** The function, as opposed to `numpy`, returns a 2-tuple, whose elements are two `ndarray`s, holding the real and imaginary parts of the transform separately. " + ] + }, + { + "cell_type": "code", + "execution_count": 114, + "metadata": { + "ExecuteTime": { + "end_time": "2020-02-16T18:38:07.294862Z", + "start_time": "2020-02-16T18:38:07.233842Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "real part:\t array([5119.996, -5.004663, -5.004798, ..., -5.005482, -5.005643, -5.006577], dtype=float)\r\n", + "\r\n", + "imaginary part:\t array([0.0, 1631.333, 815.659, ..., -543.764, -815.6588, -1631.333], dtype=float)\r\n", + "\r\n", + "real part:\t array([5119.996, -5.004663, -5.004798, ..., -5.005482, -5.005643, -5.006577], dtype=float)\r\n", + "\r\n", + "imaginary part:\t array([0.0, 1631.333, 815.659, ..., -543.764, -815.6588, -1631.333], dtype=float)\r\n", + "\n" + ] + } + ], + "source": [ + "%%micropython -pyboard 1\n", + "\n", + "from ulab import numpy as np\n", + "\n", + "x = np.linspace(0, 10, num=1024)\n", + "y = np.sin(x)\n", + "z = np.zeros(len(x))\n", + "\n", + "a, b = np.fft.fft(x)\n", + "print('real part:\\t', a)\n", + "print('\\nimaginary part:\\t', b)\n", + "\n", + "c, d = np.fft.fft(x, z)\n", + "print('\\nreal part:\\t', c)\n", + "print('\\nimaginary part:\\t', d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### ulab with complex support\n", + "\n", + "If the `ULAB_SUPPORTS_COMPLEX`, and `ULAB_FFT_IS_NUMPY_COMPATIBLE` pre-processor constants are set to 1 in [ulab.h](https://github.com/v923z/micropython-ulab/blob/master/code/ulab.h) as \n", + "\n", + "```c\n", + "// Adds support for complex ndarrays\n", + "#ifndef ULAB_SUPPORTS_COMPLEX\n", + "#define ULAB_SUPPORTS_COMPLEX (1)\n", + "#endif\n", + "```\n", + "\n", + "```c\n", + "#ifndef ULAB_FFT_IS_NUMPY_COMPATIBLE\n", + "#define ULAB_FFT_IS_NUMPY_COMPATIBLE (1)\n", + "#endif\n", + "```\n", + "then the FFT routine will behave in a `numpy`-compatible way: the single input array can either be real, in which case the imaginary part is assumed to be zero, or complex. The output is also complex. \n", + "\n", + "While `numpy`-compatibility might be a desired feature, it has one side effect, namely, the FFT routine consumes approx. 50% more RAM. The reason for this lies in the implementation details." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ifft\n", + "\n", + "The above-mentioned rules apply to the inverse Fourier transform. The inverse is also normalised by `N`, the number of elements, as is customary in `numpy`. With the normalisation, we can ascertain that the inverse of the transform is equal to the original array." + ] + }, + { + "cell_type": "code", + "execution_count": 459, + "metadata": { + "ExecuteTime": { + "end_time": "2019-10-19T13:08:17.647416Z", + "start_time": "2019-10-19T13:08:17.597456Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "original vector:\t array([0.0, 0.009775016, 0.0195491, ..., -0.5275068, -0.5357859, -0.5440139], dtype=float)\n", + "\n", + "real part of inverse:\t array([-2.980232e-08, 0.0097754, 0.0195494, ..., -0.5275064, -0.5357857, -0.5440133], dtype=float)\n", + "\n", + "imaginary part of inverse:\t array([-2.980232e-08, -1.451171e-07, 3.693752e-08, ..., 6.44871e-08, 9.34986e-08, 2.18336e-07], dtype=float)\n", + "\n" + ] + } + ], + "source": [ + "%%micropython -pyboard 1\n", + "\n", + "from ulab import numpy as np\n", + "\n", + "x = np.linspace(0, 10, num=1024)\n", + "y = np.sin(x)\n", + "\n", + "a, b = np.fft.fft(y)\n", + "\n", + "print('original vector:\\t', y)\n", + "\n", + "y, z = np.fft.ifft(a, b)\n", + "# the real part should be equal to y\n", + "print('\\nreal part of inverse:\\t', y)\n", + "# the imaginary part should be equal to zero\n", + "print('\\nimaginary part of inverse:\\t', z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that unlike in `numpy`, the length of the array on which the Fourier transform is carried out must be a power of 2. If this is not the case, the function raises a `ValueError` exception." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### ulab with complex support\n", + "\n", + "The `fft.ifft` function can also be made `numpy`-compatible by setting the `ULAB_SUPPORTS_COMPLEX`, and `ULAB_FFT_IS_NUMPY_COMPATIBLE` pre-processor constants to 1." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computation and storage costs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### RAM\n", + "\n", + "The FFT routine of `ulab` calculates the transform in place. This means that beyond reserving space for the two `ndarray`s that will be returned (the computation uses these two as intermediate storage space), only a handful of temporary variables, all floats or 32-bit integers, are required. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Speed of FFTs\n", + "\n", + "A comment on the speed: a 1024-point transform implemented in python would cost around 90 ms, and 13 ms in assembly, if the code runs on the pyboard, v.1.1. You can gain a factor of four by moving to the D series \n", + "https://github.com/peterhinch/micropython-fourier/blob/master/README.md#8-performance. " + ] + }, + { + "cell_type": "code", + "execution_count": 494, + "metadata": { + "ExecuteTime": { + "end_time": "2019-10-19T13:25:40.540913Z", + "start_time": "2019-10-19T13:25:40.509598Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "execution time: 1985 us\n", + "\n" + ] + } + ], + "source": [ + "%%micropython -pyboard 1\n", + "\n", + "from ulab import numpy as np\n", + "\n", + "x = np.linspace(0, 10, num=1024)\n", + "y = np.sin(x)\n", + "\n", + "@timeit\n", + "def np_fft(y):\n", + " return np.fft.fft(y)\n", + "\n", + "a, b = np_fft(y)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The C implementation runs in less than 2 ms on the pyboard (we have just measured that), and has been reported to run in under 0.8 ms on the D series board. That is an improvement of at least a factor of four. " + ] + } + ], + "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.5" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "382.797px" + }, + "toc_section_display": true, + "toc_window_display": true + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} |