path: root/circuitpython/extmod/ulab/docs/numpy-fft.ipynb

{
 "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
}