{"cells": [{"cell_type": "markdown", "id": "c760c855", "metadata": {}, "source": ["# ONNX FFTs\n", "\n", "Implementation of a couple of variations of FFT (see [FFT](https://www.tensorflow.org/xla/operation_semantics#fft) in ONNX."]}, {"cell_type": "code", "execution_count": 1, "id": "ddecaddb", "metadata": {}, "outputs": [{"data": {"text/html": ["\n", ""], "text/plain": [""]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "code", "execution_count": 2, "id": "6f17f4f5", "metadata": {}, "outputs": [], "source": ["%matplotlib inline"]}, {"cell_type": "code", "execution_count": 3, "id": "75f2064c", "metadata": {}, "outputs": [], "source": ["%load_ext mlprodict"]}, {"cell_type": "markdown", "id": "05a249a3", "metadata": {}, "source": ["## Signature\n", "\n", "We try to use function [FFT](https://www.tensorflow.org/xla/operation_semantics#fft) or [torch.fft.fftn](https://pytorch.org/docs/stable/generated/torch.fft.fftn.html#torch.fft.fftn)."]}, {"cell_type": "code", "execution_count": 4, "id": "6e2fc017", "metadata": {}, "outputs": [{"data": {"text/plain": ["1302"]}, "execution_count": 5, "metadata": {}, "output_type": "execute_result"}], "source": ["import numpy\n", "from numpy.testing import assert_almost_equal\n", "\n", "def numpy_fftn(x, fft_type, fft_length, axes):\n", " \"\"\"\n", " Implements FFT\n", "\n", " :param x: input\n", " :param fft_type: string (see below)\n", " :param fft_length: length on each axis of axes\n", " :param axes: axes\n", " :return: result\n", " \n", " * `'FFT`': complex-to-complex FFT. Shape is unchanged.\n", " * `'IFFT`': Inverse complex-to-complex FFT. Shape is unchanged.\n", " * `'RFFT`': Forward real-to-complex FFT.\n", " Shape of the innermost axis is reduced to fft_length[-1] // 2 + 1 if fft_length[-1]\n", " is a non-zero value, omitting the reversed conjugate part of \n", " the transformed signal beyond the Nyquist frequency.\n", " * `'IRFFT`': Inverse real-to-complex FFT (ie takes complex, returns real).\n", " Shape of the innermost axis is expanded to fft_length[-1] if fft_length[-1] \n", " is a non-zero value, inferring the part of the transformed signal beyond the Nyquist\n", " frequency from the reverse conjugate of the 1 to fft_length[-1] // 2 + 1 entries.\n", " \"\"\"\n", " if fft_type == 'FFT':\n", " return numpy.fft.fftn(x, fft_length, axes=axes)\n", " raise NotImplementedError(\"Not implemented for fft_type=%r.\" % fft_type)\n", " \n", "\n", "def test_fct(fct1, fct2, fft_type='FFT', decimal=5):\n", " cases = list(range(4, 20))\n", " dims = [[c] for c in cases] + [[4,4,4,4], [4,5,6,7]]\n", " lengths_axes = [([c], [0]) for c in cases] + [\n", " ([2, 2, 2, 2], None), ([2, 6, 7, 2], None), ([2, 3, 4, 5], None),\n", " ([2], [3]), ([3], [2])]\n", " n_test = 0\n", " for ndim in range(1, 5):\n", " for dim in dims:\n", " for length, axes in lengths_axes:\n", " if axes is None:\n", " axes = range(ndim)\n", " di = dim[:ndim]\n", " axes = [min(len(di) - 1, a) for a in axes]\n", " le = length[:ndim]\n", " if len(length) > len(di):\n", " continue\n", " mat = numpy.random.randn(*di).astype(numpy.float32)\n", " try:\n", " v1 = fct1(mat, fft_type, le, axes=axes)\n", " except Exception as e:\n", " raise AssertionError(\n", " \"Unable to run %r mat.shape=%r ndim=%r di=%r fft_type=%r le=%r \"\n", " \"axes=%r exc=%r\" %(\n", " fct1, mat.shape, ndim, di, fft_type, le, axes, e))\n", " v2 = fct2(mat, fft_type, le, axes=axes)\n", " try:\n", " assert_almost_equal(v1, v2, decimal=decimal)\n", " except AssertionError as e:\n", " raise AssertionError(\n", " \"Failure mat.shape=%r, fft_type=%r, fft_length=%r\" % (\n", " mat.shape, fft_type, le)) from e\n", " n_test += 1\n", " return n_test\n", "\n", "\n", "test_fct(numpy_fftn, numpy_fftn)"]}, {"cell_type": "code", "execution_count": 5, "id": "84993aa6", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["1.81 s \u00b1 0 ns per loop (mean \u00b1 std. dev. of 1 run, 1 loop each)\n"]}], "source": ["%timeit -n 1 -r 1 test_fct(numpy_fftn, numpy_fftn)"]}, {"cell_type": "code", "execution_count": 6, "id": "9de9a43f", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["2.07 s \u00b1 0 ns per loop (mean \u00b1 std. dev. of 1 run, 1 loop each)\n"]}], "source": ["import torch\n", "\n", "def torch_fftn(x, fft_type, fft_length, axes):\n", " xt = torch.tensor(x)\n", " if fft_type == 'FFT':\n", " return torch.fft.fftn(xt, fft_length, axes).cpu().detach().numpy()\n", " \n", "%timeit -n 1 -r 1 test_fct(numpy_fftn, torch_fftn)"]}, {"cell_type": "markdown", "id": "e55d6dbf", "metadata": {}, "source": ["## Numpy implementation"]}, {"cell_type": "code", "execution_count": 7, "id": "aa74068d", "metadata": {}, "outputs": [], "source": ["import numpy\n", "\n", "\n", "def _dft_cst(N, fft_length, dtype):\n", " def _arange(dim, dtype, resh):\n", " return numpy.arange(dim).astype(dtype).reshape(resh)\n", "\n", " def _prod(n, k):\n", " return (-2j * numpy.pi * k / fft_length) * n\n", "\n", " def _exp(m):\n", " return numpy.exp(m)\n", " \n", " n = _arange(N, dtype, (-1, 1))\n", " k = _arange(fft_length, dtype, (1, -1))\n", " M = _exp(_prod(n, k))\n", " return M\n", "\n", "\n", "def custom_fft(x, fft_type, length, axis, dft_fct=None):\n", " # https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/fft/_pocketfft.py#L56\n", " if dft_fct is None:\n", " dft_fct = _dft_cst\n", " if fft_type == 'FFT':\n", " if x.shape[axis] > length:\n", " # fft_length > shape on the same axis\n", " # the matrix is shortened\n", " slices = [slice(None)] * len(x.shape)\n", " slices[axis] = slice(0, length)\n", " new_x = x[tuple(slices)]\n", " elif x.shape[axis] == length:\n", " new_x = x\n", " else:\n", " # other, the matrix is completed with zeros\n", " shape = list(x.shape)\n", " shape[axis] = length\n", " slices = [slice(None)] * len(x.shape)\n", " slices[axis] = slice(0, length)\n", " zeros = numpy.zeros(tuple(shape), dtype=x.dtype)\n", " index = [slice(0, i) for i in x.shape]\n", " zeros[tuple(index)] = x\n", " new_x = zeros\n", "\n", " cst = dft_fct(new_x.shape[axis], length, x.dtype)\n", " perm = numpy.arange(len(x.shape)).tolist() \n", " if perm[axis] == perm[-1]:\n", " res = numpy.matmul(new_x, cst).transpose(perm)\n", " else:\n", " perm[axis], perm[-1] = perm[-1], perm[axis] \n", " rest = new_x.transpose(perm)\n", " res = numpy.matmul(rest, cst).transpose(perm)\n", " perm[axis], perm[0] = perm[0], perm[axis]\n", " return res\n", " raise ValueError(\"Unexpected value for fft_type=%r.\" % fft_type)\n", "\n", "\n", "def custom_fftn(x, fft_type, fft_length, axes, dft_fct=None):\n", " if len(axes) != len(fft_length):\n", " raise ValueError(\"Length mismatch axes=%r, fft_length=%r.\" % (\n", " axes, fft_length))\n", " if fft_type == 'FFT':\n", " res = x\n", " for i in range(len(fft_length) - 1, -1, -1):\n", " length = fft_length[i]\n", " axis = axes[i]\n", " res = custom_fft(res, fft_type, length, axis, dft_fct=dft_fct)\n", " return res\n", " raise ValueError(\"Unexpected value for fft_type=%r.\" % fft_type)\n", "\n", " \n", "shape = (4, )\n", "fft_length = [5,]\n", "axes = [0]\n", "rnd = numpy.random.randn(*shape) + numpy.random.randn(*shape) * 1j\n", "custom_fftn(rnd, 'FFT', fft_length, axes), numpy_fftn(rnd, 'FFT', fft_length, axes)\n", "assert_almost_equal(custom_fftn(rnd, 'FFT', fft_length, axes),\n", " numpy_fftn(rnd, 'FFT', fft_length, axes), decimal=5)\n", "\n", "shape = (4, 3)\n", "fft_length = [3, 2]\n", "axes = [0, 1]\n", "rnd = numpy.random.randn(*shape) + numpy.random.randn(*shape) * 1j\n", "custom_fftn(rnd, 'FFT', fft_length, axes), numpy_fftn(rnd, 'FFT', fft_length, axes)\n", "assert_almost_equal(custom_fftn(rnd, 'FFT', fft_length, axes),\n", " numpy_fftn(rnd, 'FFT', fft_length, axes), decimal=5)"]}, {"cell_type": "code", "execution_count": 8, "id": "5c454666", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["2.35 s \u00b1 0 ns per loop (mean \u00b1 std. dev. of 1 run, 1 loop each)\n"]}], "source": ["%timeit -n 1 -r 1 test_fct(numpy_fftn, custom_fftn, decimal=4)"]}, {"cell_type": "markdown", "id": "f27bd70d", "metadata": {}, "source": ["## Benchmark"]}, {"cell_type": "code", "execution_count": 9, "id": "507d8348", "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 24/24 [00:06<00:00, 3.91it/s]\n"]}, {"data": {"text/html": ["\n", "\n", "
\n", " \n", " \n", " | name | \n", " custom_fftn | \n", " numpy_fftn | \n", " torch_fftn | \n", "
\n", " \n", " | length | \n", " | \n", " | \n", " | \n", "
\n", " \n", " \n", " \n", " | 8 | \n", " 0.000585 | \n", " 0.000911 | \n", " 0.003643 | \n", "
\n", " \n", " | 16 | \n", " 0.001669 | \n", " 0.001373 | \n", " 0.004087 | \n", "
\n", " \n", " | 24 | \n", " 0.002682 | \n", " 0.003273 | \n", " 0.005745 | \n", "
\n", " \n", " | 32 | \n", " 0.004288 | \n", " 0.003275 | \n", " 0.004657 | \n", "
\n", " \n", " | 40 | \n", " 0.004818 | \n", " 0.003831 | \n", " 0.005198 | \n", "
\n", " \n", "
\n", "
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["piv.plot(logy=True, logx=True, title=\"FFT benchmark\", figsize=(12, 4));"]}, {"cell_type": "code", "execution_count": 11, "id": "9c29380c", "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:13<00:00, 1.33s/it]\n"]}, {"data": {"text/html": ["\n", "\n", "
\n", " \n", " \n", " | name | \n", " custom_fftn | \n", " numpy_fftn | \n", " torch_fftn | \n", "
\n", " \n", " | length | \n", " | \n", " | \n", " | \n", "
\n", " \n", " \n", " \n", " | 2 | \n", " 0.000434 | \n", " 0.001167 | \n", " 0.023980 | \n", "
\n", " \n", " | 4 | \n", " 0.001117 | \n", " 0.001671 | \n", " 0.022530 | \n", "
\n", " \n", " | 8 | \n", " 0.001428 | \n", " 0.002077 | \n", " 0.022102 | \n", "
\n", " \n", " | 16 | \n", " 0.004654 | \n", " 0.002874 | \n", " 0.019792 | \n", "
\n", " \n", " | 32 | \n", " 0.003172 | \n", " 0.002689 | \n", " 0.017474 | \n", "
\n", " \n", " | 64 | \n", " 0.006966 | \n", " 0.004612 | \n", " 0.018116 | \n", "
\n", " \n", " | 128 | \n", " 0.030904 | \n", " 0.011608 | \n", " 0.023369 | \n", "
\n", " \n", " | 256 | \n", " 0.123821 | \n", " 0.025853 | \n", " 0.023532 | \n", "
\n", " \n", " | 512 | \n", " 0.476802 | \n", " 0.043352 | \n", " 0.033228 | \n", "
\n", " \n", " | 1024 | \n", " 1.527917 | \n", " 0.109868 | \n", " 0.052858 | \n", "
\n", " \n", "
\n", "
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["piv.plot(logy=True, logx=True, title=\"FFT benchmark (power2)\", figsize=(12, 4));"]}, {"cell_type": "markdown", "id": "616099d8", "metadata": {}, "source": ["## Profiling"]}, {"cell_type": "code", "execution_count": 13, "id": "531658bb", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["f -- 1 1 -- 0.01752 0.54515 -- :8:f (f)\n", " custom_fftn -- 100 100 -- 0.00234 0.52763 -- :57:custom_fftn (custom_fftn)\n", " custom_fft -- 100 100 -- 0.19936 0.52516 -- :20:custom_fft (custom_fft)\n", " _dft_cst -- 100 100 -- 0.31917 0.32366 -- :4:_dft_cst (_dft_cst)\n", " _arange -- 200 200 -- 0.00088 0.00449 -- :5:_arange (_arange)\n", " -- 200 200 -- 0.00128 0.00128 -- ~:0: ()\n", " -- 200 200 -- 0.00064 0.00064 -- ~:0: ()\n", " -- 200 200 -- 0.00169 0.00169 -- ~:0: () +++\n", " -- 100 100 -- 0.00011 0.00011 -- ~:0: () +++\n", " -- 100 100 -- 0.00024 0.00024 -- ~:0: ()\n", " -- 100 100 -- 0.00076 0.00076 -- ~:0: ()\n", " -- 100 100 -- 0.00102 0.00102 -- ~:0: () +++\n", " -- 300 300 -- 0.00013 0.00013 -- ~:0: () +++\n", " -- 400 400 -- 0.00024 0.00024 -- ~:0: ()\n", " -- 300 300 -- 0.00271 0.00271 -- ~:0: ()\n"]}], "source": ["from pyquickhelper.pycode.profiling import profile2graph, profile\n", "\n", "shape = [512, 128]\n", "fft_length = [128]\n", "axes = [1]\n", "rnd = numpy.random.randn(*shape) + numpy.random.randn(*shape) * 1j\n", "\n", "def f():\n", " for i in range(100):\n", " custom_fftn(rnd, 'FFT', fft_length, axes)\n", "\n", "stat, text = profile(f)\n", "gr = profile2graph(stat)\n", "print(gr[0].to_text(fct_width=40))"]}, {"cell_type": "markdown", "id": "7690454d", "metadata": {}, "source": ["We can see that function `_dft_cst` is the bottle neck and more precisely the exponential. We need to use the symmetries of the matrix it builds."]}, {"cell_type": "markdown", "id": "4e250ff9", "metadata": {}, "source": ["## Faster _dft_cst\n", "\n", "The function builds the matrix $M_{nk} = \\left( \\exp\\left(\\frac{-2i\\pi nk}{K}\\right) \\right)_{nk}$ where $1 \\leqslant n \\leqslant N$ and $1 \\leqslant k \\leqslant K$. So it computes powers of the unity roots.\n", "\n", "$$\n", "\\exp\\left(\\frac{-2i\\pi nk}{K}\\right) = \\exp\\left(\\frac{-2i\\pi k}{K}\\right)^n = \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{nk}\n", "$$\n", "\n", "We use that expression to reduce the number of exponentiels to compute."]}, {"cell_type": "code", "execution_count": 14, "id": "8b60fd16", "metadata": {}, "outputs": [{"data": {"text/plain": ["((3, 4), dtype('complex64'))"]}, "execution_count": 15, "metadata": {}, "output_type": "execute_result"}], "source": ["import numpy\n", "from numpy.testing import assert_almost_equal\n", "\n", "def _dft_cst(N, fft_length, dtype=numpy.float32):\n", " def _arange(dim, dtype, resh):\n", " return numpy.arange(dim).astype(dtype).reshape(resh)\n", "\n", " n = _arange(N, dtype, (-1, 1))\n", " k = _arange(fft_length, dtype, (1, -1))\n", " M = (-2j * numpy.pi * k / fft_length) * n\n", " numpy.exp(M, out=M)\n", " return M\n", "\n", "\n", "M = _dft_cst(3, 4, numpy.float32)\n", "M.shape, M.dtype"]}, {"cell_type": "code", "execution_count": 15, "id": "0600a293", "metadata": {}, "outputs": [{"data": {"text/plain": ["((4, 3), dtype('complex128'))"]}, "execution_count": 16, "metadata": {}, "output_type": "execute_result"}], "source": ["M = _dft_cst(4, 3, numpy.float64)\n", "M.shape, M.dtype"]}, {"cell_type": "code", "execution_count": 16, "id": "38d760e2", "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[ 1. +0.00000000e+00j, 1. +0.00000000e+00j, 1. +0.00000000e+00j],\n", " [ 1. +0.00000000e+00j, -0.5-8.66025404e-01j, -0.5+8.66025404e-01j],\n", " [ 1. +0.00000000e+00j, -0.5+8.66025404e-01j, -0.5-8.66025404e-01j],\n", " [ 1. +0.00000000e+00j, 1. +2.44929360e-16j, 1. +4.89858720e-16j]])"]}, "execution_count": 17, "metadata": {}, "output_type": "execute_result"}], "source": ["M"]}, {"cell_type": "code", "execution_count": 17, "id": "30466d1b", "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[ 1. +0.00000000e+00j, 1. +0.00000000e+00j, 1. +0.00000000e+00j],\n", " [ 1. +0.00000000e+00j, -0.5-8.66025404e-01j, -0.5+8.66025404e-01j],\n", " [ 1. +0.00000000e+00j, -0.5+8.66025404e-01j, -0.5-8.66025404e-01j],\n", " [ 1. +0.00000000e+00j, 1. +6.10622664e-16j, 1. +1.22124533e-15j]])"]}, "execution_count": 18, "metadata": {}, "output_type": "execute_result"}], "source": ["def _dft_cst_power(N, fft_length, dtype=numpy.float32):\n", " if dtype == numpy.float32:\n", " ctype = numpy.complex64\n", " else:\n", " ctype = numpy.complex128\n", " M = numpy.empty((N, fft_length), dtype=ctype)\n", " M[0, :] = 1\n", " M[1, 0] = 1\n", " root = numpy.exp(numpy.pi / fft_length * (-2j))\n", " current = root\n", " M[1, 1] = root\n", " for i in range(2, M.shape[1]):\n", " current *= root\n", " M[1, i] = current\n", " for i in range(2, M.shape[0]):\n", " numpy.multiply(M[i-1, :], M[1, :], out=M[i, :])\n", " return M\n", "\n", "M_pow = _dft_cst_power(4, 3, numpy.float64)\n", "M_pow"]}, {"cell_type": "code", "execution_count": 18, "id": "965651d7", "metadata": {}, "outputs": [], "source": ["assert_almost_equal(M, M_pow)"]}, {"cell_type": "code", "execution_count": 19, "id": "a7194b92", "metadata": {}, "outputs": [], "source": ["dims = (10, 15)\n", "assert_almost_equal(_dft_cst(*dims, dtype=numpy.float32), \n", " _dft_cst_power(*dims, dtype=numpy.float32),\n", " decimal=5)"]}, {"cell_type": "markdown", "id": "36b1a3d7", "metadata": {}, "source": ["## Benchmark again"]}, {"cell_type": "code", "execution_count": 20, "id": "40a61cb1", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["1.46 s \u00b1 0 ns per loop (mean \u00b1 std. dev. of 1 run, 1 loop each)\n"]}], "source": ["def custom_fftn_power(*args, **kwargs):\n", " return custom_fftn(*args, dft_fct=_dft_cst_power, **kwargs)\n", "\n", "\n", "%timeit -r 1 -n 1 test_fct(numpy_fftn, custom_fftn_power, decimal=4)"]}, {"cell_type": "code", "execution_count": 21, "id": "3a2707eb", "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 24/24 [00:07<00:00, 3.19it/s]\n"]}, {"data": {"text/html": ["\n", "\n", "
\n", " \n", " \n", " | name | \n", " custom_fftn | \n", " custom_fftn_power | \n", " numpy_fftn | \n", " torch_fftn | \n", "
\n", " \n", " | length | \n", " | \n", " | \n", " | \n", " | \n", "
\n", " \n", " \n", " \n", " | 8 | \n", " 0.000991 | \n", " 0.000837 | \n", " 0.001177 | \n", " 0.007033 | \n", "
\n", " \n", " | 16 | \n", " 0.002758 | \n", " 0.002591 | \n", " 0.002069 | \n", " 0.006228 | \n", "
\n", " \n", " | 24 | \n", " 0.003087 | \n", " 0.002816 | \n", " 0.002499 | \n", " 0.005564 | \n", "
\n", " \n", " | 32 | \n", " 0.003767 | \n", " 0.003068 | \n", " 0.003306 | \n", " 0.005985 | \n", "
\n", " \n", " | 40 | \n", " 0.004710 | \n", " 0.003975 | \n", " 0.004044 | \n", " 0.005733 | \n", "
\n", " \n", "
\n", "
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["piv.plot(logy=True, logx=True, title=\"FFT benchmark 2\", figsize=(12, 4));"]}, {"cell_type": "code", "execution_count": 23, "id": "ed6c4d1a", "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["f -- 1 1 -- 0.02624 0.57688 -- :8:f (f)\n", " custom_fftn_power -- 100 100 -- 0.00094 0.55064 -- :1:custom_fftn_power (custom_fftn_power)\n", " custom_fftn -- 100 100 -- 0.00609 0.54970 -- :57:custom_fftn (custom_fftn)\n", " custom_fft -- 100 100 -- 0.46378 0.54342 -- :20:custom_fft (custom_fft)\n", " _dft_cst_power -- 100 100 -- 0.07599 0.07726 -- :1:_dft_cst_power (_dft_cst_power)\n", " -- 100 100 -- 0.00126 0.00126 -- ~:0: ()\n", " -- 100 100 -- 0.00008 0.00008 -- ~:0: () +++\n", " -- 100 100 -- 0.00025 0.00025 -- ~:0: ()\n", " -- 100 100 -- 0.00096 0.00096 -- ~:0: ()\n", " -- 100 100 -- 0.00109 0.00109 -- ~:0: ()\n", " -- 300 300 -- 0.00020 0.00020 -- ~:0: () +++\n", " -- 400 400 -- 0.00027 0.00027 -- ~:0: ()\n"]}], "source": ["from pyquickhelper.pycode.profiling import profile2graph, profile\n", "\n", "shape = [512, 128]\n", "fft_length = [128]\n", "axes = [1]\n", "rnd = numpy.random.randn(*shape) + numpy.random.randn(*shape) * 1j\n", "\n", "def f():\n", " for i in range(100):\n", " custom_fftn_power(rnd, 'FFT', fft_length, axes)\n", "\n", "stat, text = profile(f)\n", "gr = profile2graph(stat)\n", "print(gr[0].to_text(fct_width=40))"]}, {"cell_type": "markdown", "id": "9ef4af25", "metadata": {}, "source": ["## Cooley\u2013Tukey FFT algorithm\n", "\n", "See [Cooley\u2013Tukey FFT algorithm](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm).\n", "\n", "The FFT matrix is defined by the matrix computation $F_{ak} = X_{an} M_{nk}$, then one coefficient is ($1 \\leqslant n, k \\leqslant K$):\n", "\n", "$$\n", "F_{ak} = \\sum_n X_{an} M_{nk} = \\sum_n X_{an} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{nk}\n", "$$\n", "\n", "Let's assume K is even, then $\\exp\\left(\\frac{-2i\\pi k}{K}\\right) = -\\exp\\left(\\frac{-2i\\pi \\left(k + \\frac{K}{2}\\right)}{K}\\right)$."]}, {"cell_type": "code", "execution_count": 24, "id": "803fcc3a", "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAADSCAYAAADAMi7MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6yElEQVR4nO3de3xU1b3w/8+XBCVYEdBDJQQVDiUJkDBKSLmJiIaLQW6lrRBExCiIHH9tH0F8VXgsPr5Aoe052CB9JByOBAVr0dIIoRgSavKggBohgqEBI01IEQPINRiT7++PmYyTZHIfQib5vl+veTH7stZee2ezv3utvWZtUVWMMcYYf9bmahfAGGOMaSwLZsYYY/yeBTNjjDF+z4KZMcYYv2fBzBhjjN+zYGaMMcbvWTAzrZ6IrBaRRVe7HMaYhhP7nZkx3xOREUCSqoY00fZmAvGqOqwptmdMS2U1M2MaSUQCrnYZjGntLJiZFkFEVER6eUyvE5H/4/o+QkTyReR/ichXIlIoIg9XXldErgO2AcEict71CfayrXUi8oqIbBWRC8DdIhIuIukickZEPhOR8R7r3yAir4nISRH5UkSeFZE2IhIOrAYGu7Z1xrX+fSJyUETOiUiBiDx1pY6bMS2FBTPTWtwM3AB0Ax4BEkSkk+cKqnoBGAscV9UfuD7Hq8lvGvACcD3wIfBX4G9AF+A/gA0iEupa92XXtnsCdwEzgIdV9RAwB9jt2lZH1/qJwGxVvR7oB+xs7M4b09JZMDOtRQmwRFVLVHUrcB4IrSVNTf6iqpmqWgY4gB8Ay1T1W1XdCSQDU11NkA8Az6jqOVXNA34LPFhLWfuISAdVPa2qHzeinMa0ChbMTGtRpKrfeUxfxBmAGuqfHt+DgX+6Alu5L3HWAm8C2rqmKy+rzk+A+4AvRWSXiAxuRDmNaRUsmJmW4iLQ3mP65gbmU9fuvZ7rHQe6i4jn/6dbgALga5w1rVu9LPO6PVXdq6oTcDZZvgO8WccyGdNqWTAzLUUWME1EAkRkDM5nUw1xArhRRG6oR5oPcQbTBSLS1tW9/35go6qW4gxGL4jI9SJyK/ArIMljeyEicg2AiFwjInEicoOqlgBngTKMMTWyYGZaiv8PZwA5A8ThrNHUm6p+DrwBHHX1TKzSm9FLmm9d2x6Lsya2CpjhygucHUIuAEeBDOB1YK1r2U7gM+BfIvK1a96DQJ6InMXZQSSuIftiTGtiP5o2xhjj96xmZowxxu9ZMDPGGOP3LJgZY4zxexbMjDHG+D0LZsYYY/xe4NUuQHVuuukmve222652MYy54lSVf/zjH/z7v/87AQEVB+AvKirixIkTqCoBAQHccssttG/fvpqcnPLy8vjhD39IUFBQhflnz54lPz8fVaV9+/bcdtttiIi7DIcOHaJPnz7k5eXxzTffEBgYSN++fd3p8/PzOXPmDG3atOHaa6/l1ltvJTAwkEuXLnHixAns/2vz8NFHH32tqv92tcvR5FS1WX4GDBigxrQGycnJ+otf/MLrsszMTD116pSqqm7dulWjo6MbtI3S0lINCQnRnJwcVVVdtGiRrlmzxr18586dOm/ePFVV3bVrl3700Ufat2/fCnls375dS0pKVFV1wYIFumDBAveye+65R7/88ssGlc34FrBPm8E1vKk/1sxoTAMlJSURHR2Nw+Fg9uzZlJaWsnfvXiIjIykuLubChQv07duX7Oxs0tPTGT58OLGxsYSGhjJnzhzKypwDe2zYsIEJEyZ43caQIUPo1Mk5uP+gQYPIz88HnLWvsLAw4uLiCA8PZ8qUKVy8eBGAESNGsG/fvgr5FBUVcc0119C7d28AYmJi+POf/+xenpKSwtixYwEYPnw4nTt3rlKWUaNGERgYWKUsAPfffz8bN26s/0E0xkcsmBnTAIcOHWLTpk1kZmaSlZVFQEAAGzZsYODAgYwfP55nn32WBQsWMH36dPr16wfAnj17ePnllzl48CBHjhxh8+bNAGRmZjJgwIBat5mYmOgOOAA5OTnMnTuXQ4cO0aFDB1atWlVt2ptuuonvvvvOHeTeeust/vnP78dKTktLY8SIEXXe/7Vr11YoS1RUFO+//36d0xvjaz4JZiKy1vXSw+xqlouIrBSRXBHZLyJ3+GK7xjS1dz4pYOiynQx78r/Yvms3vfs5cDgcpKamcvToUQAWL17Mjh072LdvHwsWLHCnjY6OpmfPngQEBDB16lQyMjIAOHXqFNdff32N201LSyMxMZEXX3zRPa979+4MHToUgOnTp7vz80ZE2LhxI7/85S+Jjo7m+uuvdz+fKygooHPnzrU+iyv3wgsvEBgYSFzc96NsdenShePHq3v1mzFXnq86gKwD/gC8Vs3yscCPXJ8fA6+4/jXGb7zzSQHPbD7ApZJSFAjqezft7n2E5yZHMPH279/oUlRUxPnz5ykpKaG4uJjrrrsOwN3Zolz5dGBgIGVlZbRp04aEhAReffVVALZu3UpwcDD79+8nPj6ebdu2ceONN1ZJX910ZYMHD3bXnv72t79x+PBhwNnEOHr06Dodg3Xr1pGcnExqamqF7RUXF1fpcGJMU/JJzUxV/w6cqmGVCcBrrueTHwAdRaSrL7ZtTFNZvj2HSyWlALS7tT8XczI5f6aI5dtzOHXqFF9+6Xxl2ezZs3n++eeJi4vj6aefdqffs2cPX3zxBWVlZWzatIlhw4YBEBoa6q7VPfHEE2RlZZGVlUVwcDDHjh1j8uTJrF+/3v28q9yxY8fYvXs3AK+//ro7v+p89dVXAFy+fJkXX3yROXPmABWfl9UkJSWFl156iS1btlSpxR0+fNjdnGrM1dBUz8y6UfFlhvl4eTmhiDwmIvtEZN/JkyebqGjG1M3xM5fc36+56RY63vkgJ95cxN7fPUJMTAyFhYW89tprtG3blmnTprFw4UL27t3Lzp07ARg4cCDz5s0jPDycHj16MGnSJABiY2NJT0/3us0lS5ZQVFTE3LlzcTgcREVFuZeFhoaSkJBAeHg4p0+f5vHHH6+x/MuXLyc8PJzIyEjuv/9+Ro4cSWlpKbm5uYSFhbnXmzp1KoMHDyYnJ4eQkBASExMBmDdvHufOnSMmJgaHw+EOhuBsBo2Nja3fATXGh3w2ar6I3AYkq2qV2zMRScb5SvkM13Qq8LSq7qu8brmoqCit3CPLmKtp6LKdFHgEtHLdOgaRuXBkjWnT09NZsWIFycnJVZYVFhYyY8YMduzYUeey5OXlMW7cOLKzvT6mrrOMjAySkpJYvXp1g/O4fPkyd911FxkZGe7ejubqEZGPVDWq9jVblqaqmRUA3T2mQ/j+TbvG+IX5o0MJalvxR81BbQOYPzq0Ufl27dqVRx99lLNnzzYqn4YYNmxYowIZOJs7ly1bZoHMXFVNFcy2ADNcvRoHAd+oamETbds0MVVl5MiRXi/OGzZsIDIykoiICIYMGcKnn35aa37x8fEcPHiwyvzU1FTuuOMOHA4Hw4YNIzc3172ssLCQUaNGATBmzBg6duzIuHHjKqSPi4sjNDSUfv36MWvWLEpKSgBITk5m8eLFVbY38fZuLJ0cQbeOQQjOGtnSSp0/qjNixAivtbJyP/vZz+jQoUOt+ZS77bbbGl0r85Uf/ehH9erW31hNdX794Q9/oFevXogIX3/9tXu+qvLkk0/Sq1cvIiMj+fjjjyukGzt2LPn5+dWmr66M3377LcOHD+e7776r87EwHnzxy2ucb+YtBEpwPg97BOcbcue4lguQABwBDgBRteVpI4D4r6YY0UJV9Uc/+pEePHhQVVUTEhL0oYceci9bu3atrlixQlVV33vvPd2yZYvGxsZWSP/uu+9qWVmZlpWV6QMPPKCrVq1SVdWysjJ1OBx64cKFBpfNXDlNdX59/PHH+sUXX+itt96qJ0+edM9/9913dcyYMVpWVqa7d++usI2LFy/qwIEDa0xfUxmfe+45TUpKanCZVW0EkMYGxKmq2lVV26pqiKomqupqVV3tWq6q+oSq/ruqRmgNz8rM1eNPI1qAsyt6+d35N998Q3BwsHuZZw+9e+65x+vvuO677z5EBBEhOjraXRYRqbUmZerP386v22+/3et4k3/5y1+YMWMGIsKgQYM4c+YMhYXOhqb09HR3LbW69NWVEWDixIls2LChDkfTVGYjgBjA/0a0AFizZg333XcfISEhrF+/noULFwJQWlpKTk4Offr0qdO+l5SUsH79esaMGeOeZyNa+JY/nl/VKSgooHv377sAhISEUFDg7AKwbdu2CudRfcvYr18/9u7d26BytXYWzFo5fx3RAuD3v/89W7duJT8/n4cffphf/epXAHz44Yf8+Md1/03+3LlzGT58OHfeead7no1o4Rv+fH41RGZmZq2/96upjAEBAVxzzTWcO3fO52Vr6az7USvmzyNanDx5kk8//dQdtH7+85+774jrc3f8m9/8hpMnT/LHP/6xwnwb0aLx/Pn8qkm3bt0qjGuZn59Pt27dOHr0KN27d+eaa66pNY/qygjOnzq0a9euQWVrzaxm1or584gWnTp14ptvvnEPybRjxw7Cw8MBZy/He++9t9b9X7NmDdu3b+eNN96gTZuK/xVsRIvG8+fzqybjx4/ntddeQ1X54IMPuOGGG+jatWudb6JqKmNRURE33XQTbdu2bVDZWjMLZq2YP49oERgYyKuvvspPfvIT+vfvz/r161m+fDknT56kXbt2FZqh7rzzTn7605+SmppKSEgI27dvB2DOnDmcOHGCwYMH43A4WLJkiTuNjWjReP58fgGsXLmSkJAQ8vPziYyMJD4+HnB2HOrZsye9evXi0UcfdT97S0lJqRDMqktfUxntvGs4n40A4ms2AsiV1xJHtEhKSiI/P9/dGaQhTpw4wbRp00hNTW1UWVq7lnh+Vefy5csMHTrUa6/I+pg8eTLLli2rUmOrDxsBxLQ6LXFEi+nTpzcqkIGzGei3v/2tj0rUerXE86s61157baMD2bfffsvEiRMbFchaM6uZtXLvfFLA8u05HD9zieCOQcwfHVqnES2MqQs7v5pea62ZWTAzxpgWpLUGM2tmNMYY4/csmBljjPF7FsyMMcb4PQtmxhhj/J4FM2OMMX7Pgpkxxhi/Z8HMGGOM37NgZowxxu9ZMDPGGOP3LJgZY4zxexbMjDHG+D0LZsYYY/yeBTNjjDF+z4KZMcYYv2fBzBhjjN+zYHYVqCojR470+qbcDRs2EBkZSUREBEOGDOHTTz+tNb/4+HgOHjxYZf6dd96Jw+HA4XAQHBzMxIkT3ctKSkq44447AJg1axZdunShX79+FdLPnz+fsLAwIiMjmTRpEmfOnAHgwIEDzJw5s+47bIypM19fH0RkjYj08TJ/nojkioiKyE0e88NEZLeIXBaRp7ykWy0iQ0XkpyLymYiUiUiUx/IYEflIRA64/h3psew9EelUtyNRPxbMroKtW7fSv39/OnToUGVZjx492LVrFwcOHGDRokU89thjtea3Zs0a+vSpcq7y/vvvk5WVRVZWFoMHD2by5MnuZRkZGQwdOhSAmTNnkpKSUiV9TEwM2dnZ7N+/n969e7N06VIAIiIiyM/P59ixY3XeZ2NM3fj6+qCq8apa9W4XMoF7gS8rzT8FPAmsqCbLQcAHQDYwGfh7peVfA/eragTwELDeY9l6YG6thW4AC2b1kJSURHR0NA6Hg9mzZ1NaWsrevXuJjIykuLiYCxcu0LdvX7Kzs0lPT2f48OHExsYSGhrKnDlzKCsrA5x3VxMmTPC6jSFDhtCpk/PGZdCgQeTn5wOQl5dHWFgYcXFxhIeHM2XKFC5evAjAiBEjqOmt3GfPnmXnzp0VamYpKSmMHTsWgOHDh9O5c+cq6UaNGkVgYGCVsgDcf//9bNy4sa6HzpgWr7lcH4C+IvKWiLQHEJF0z5pTOVX9RFXzvMz/SlX3AiWVl4lIOHBYVUtV9ZCq5lST73HX5GdAkIhc65reAkyt4TA2mAWzOjp06BCbNm0iMzOTrKwsAgIC2LBhAwMHDmT8+PE8++yzLFiwgOnTp7ub6/bs2cPLL7/MwYMHOXLkCJs3bwYgMzOTAQMG1LrNxMREd8AByMnJYe7cuRw6dIgOHTqwatWqOpX9nXfe4Z577qlwp5eWlsaIESPqvP9r166tUJaoqCjef//9Oqc3piVrTtcHnAHkLFemBjQWqNqMU72fAB+r6mUAVT0NXCsiN/q6YIG+yERExgD/BQQAa1R1WaXlM4HlQIFr1h9UdY0vtn2lvfNJAcu35/B56puc+3A3vfs5uCGoLZcuXaJLly4ALF68mIEDB9KuXTtWrlzpThsdHU3Pnj0BmDp1KhkZGUyZMoVTp05x/fXX17jdtLQ0EhMTycjIcM/r3r27u2lw+vTprFy5kqeeqtKkXcUbb7xBfHy8e7qgoIDOnTvTvn37Oh2DF154gcDAwPK7PgC6dOnC8ePHa0hlTMvXXK8PQBI1NxU21Gjg4bqsKCJ9gReBUZUWfQUEA0W+LFijg5mIBAAJQAyQD+wVkS1e2mg3qeq8xm6vKb3zSQHPbD7ApZJSFAjqezft7n2E5yZHMPH2bu71ioqKOH/+PCUlJRQXF3PdddcBICIV8iufDgwMpKysjDZt2pCQkMCrr74KONvKg4OD2b9/P/Hx8Wzbto0bb7yxSvrqpr35+uuv2bNnD2+//bZ7XkpKCqNHj67TMVi3bh3JycmkpqZW2F5xcTFBQUF1ysOYlqi5Xx8A9eX+upotO3o0Ida0bgjwNjBDVY9UWtwOuOTLsoFvmhmjgVxVPaqq3wIbAe8Nvn5m+fYcLpWUAtDu1v5czMnk/Jkilm/P4dSpU3z5pfO56ezZs3n++eeJi4vj6aefdqffs2cPX3zxBWVlZWzatIlhw4YBEBoaytGjRwF44okn3J00goODOXbsGJMnT2b9+vX07t27QnmOHTvG7t27AXj99dfd+dXkrbfeYty4cbRr1849z/N5WU1SUlJ46aWX2LJlS5Va3OHDh6v0fjSmNWnO1wdgGpCBb90NpNW2koh0BN4FFqpqZqVlAtwM5Pm4bD4JZt2Af3pM57vmVfYTEdnvejDZ3VtGIvKYiOwTkX0nT570QdEa5/iZ728errnpFjre+SAn3lzE3t89QkxMDIWFhbz22mu0bduWadOmsXDhQvbu3cvOnTsBGDhwIPPmzSM8PJwePXowadIkAGJjY0lPT/e6zSVLllBUVMTcuXNxOBxERX3/3DY0NJSEhATCw8M5ffo0jz/+eK37sHHjRqZO/f55a2lpKbm5uYSFhbnnTZ06lcGDB5OTk0NISAiJiYkAzJs3j3PnzhETE4PD4WDOnDnuNGlpacTGxtbhKBrTMjXX6wPQF+gEvFJT+UXkSRHJB0KA/SKyxjX/Ztf8XwHPiki+iHSg0vMyEZnkWm8w8K6IbHctmgf0AhaLSJbr08W1bADwgap+V8vhrTdRbVxNVESmAGNUNd41/SDwY88mRdfDvvOqellEZgM/V9WR3nN0ioqK0pp66DWFoct2UnCmam24W8cgMhfWWHzS09NZsWIFycnJVZYVFhYyY8YMduzYUeey5OXlMW7cOLKzs+ucxpuMjAySkpJYvXp1g/O4fPkyd911FxkZGe7ejsa0Ns31+iAiH6lqld6LjSUiH+O8tlfp5ViPPP4L2KKqqb4rmZMvamYFgGdNK4TvO3oAoKpF5b1ZgDU4o3OzN390KEFtAyrMC2obwPzRoY3Kt2vXrjz66KNefxR5pQ0bNqxRgQyczRnLli2zQGZatZZ4faiJqt7RmEDmkn0lAhn4pmYWCBwG7sEZxPYC01T1M491uqpqoev7JOBpVR1UU77NoWYG3/dWOn7mEsEdg5g/OrTCw11jTOvVHK8PV6pm1tw1OpgBiMh9wH/i7Jq/VlVfEJElwD5V3SIiS4HxwHc4f13+uKp+XlOezSWYGWOMP7Fg1sxYMDPGmPprrcHMRgAxxhjj9yyYGWOM8XsWzIwxxvg9C2bGGGP8ngUzY4wxfs+CmTHGGL9nwcwYY4zfs2BmjDHG71kwM8YY4/csmBljjPF7FsyMMcb4PQtmxhhj/J4FM2OMMX7Pgpkxxhi/Z8HMGGOM37NgZowxxu9ZMDPGGOP3LJgZY4zxexbMjDHG+D0LZsYYY/yeBTNjjDF+z4KZMcYYv2fBzBhjjN+zYGaMMcbvWTCrA1Vl5MiRnD17tsqyDRs2EBkZSUREBEOGDOHTTz+tNb/4+HgOHjzodTu//vWv6d27N+Hh4axcudK9rKSkhDvuuAOAWbNm0aVLF/r161ch/fz58wkLCyMyMpJJkyZx5swZAA4cOMDMmTPrscfGmJZCnHaKSAcvy+JEZL+IHBCR/yci/euQ3xoR6eNlfg8R+VBEckVkk4hc47Gsq4j8zfU9RUTOiEhypfQbRCRHRLJFZK2ItHXNHyciS2orlwWzOti6dSv9+/enQ4cq5wI9evRg165dHDhwgEWLFvHYY4/Vmt+aNWvo06fKucC6dev45z//yeeff86hQ4d44IEH3MsyMjIYOnQoADNnziQlJaVK+piYGLKzs9m/fz+9e/dm6dKlAERERJCfn8+xY8fqvM/GmBbjPuBTVa16Nw5fAHepagTwPPB/a8tMVeNVterdOLwI/F5VewGngUc8lo0Btru+Lwce9JJ+AxAGRABBQLxr/rvA/SLSvqZytehglpSURHR0NA6Hg9mzZ1NaWsrevXuJjIykuLiYCxcu0LdvX7Kzs0lPT2f48OHExsYSGhrKnDlzKCsrA5y1rwkTJnjdxpAhQ+jUqRMAgwYNIj8/H4C8vDzCwsKIi4sjPDycKVOmcPHiRQBGjBjBvn37quT1yiuvsHjxYtq0cf5ZunTp4l6WkpLC2LFjARg+fDidO3eukn7UqFEEBgZWKQvA/fffz8aNG+t3AI0xV01Dr19APxFZLSLl1/c44C/etqGq/09VT7smPwBCAETkNhH53FVbOiQib5UHExFJF5Eoz3xERICRwFuuWf8DTPRYZQywzbXNVOCcl7JsVRdgT3lZXNPpwLiajleLDWaHDh1i06ZNZGZmkpWVRUBAABs2bGDgwIGMHz+eZ599lgULFjB9+nR3c92ePXt4+eWXOXjwIEeOHGHz5s0AZGZmMmDAgFq3mZiY6A44ADk5OcydO5dDhw7RoUMHVq1aVWP6I0eOsGnTJqKiohg7diz/+Mc/3MvS0tIYMWJEnfd/7dq1FcoSFRXF+++/X+f0xpirpzHXLyAb+Hdgsiu7ocBHddjsI7gCjksosEpVw4GzwNwa0t4InFHV71zT+UA3ABEJAEKrqc1V4WpefBDwbH7aB9xZU7rAumReh42PAf4LCADWqOqySsuvBV4DBgBFwM9VNc8X267snU8KWL49h89T3+Tch7vp3c/BDUFtuXTpkrums3jxYgYOHEi7du0qPJeKjo6mZ8+eAEydOpWMjAymTJnCqVOnuP7662vcblpaGomJiWRkZLjnde/e3d00OH36dFauXMlTTz1VbR6XL1+mXbt27Nu3j82bNzNr1izef/99CgoK6Ny5M+3b11jLdnvhhRcIDAwkLi7OPa9Lly4cP368TumNMVeHr65fwBvAMJw1pc6qWqUm5ElE7sYZzIZ5zP6nqma6vicBTwIrGrBbPwY+rMf6q4C/q6rn3fdXQHBNiRodzFxRNwGIwRmN94rIlkpR+BHgtKr2EpEHcLat/ryx267snU8KeGbzAS6VlKJAUN+7aXfvIzw3OYKJt3dzr1dUVMT58+cpKSmhuLiY6667rnxfKu8bAIGBgZSVldGmTRsSEhJ49dVXAeeztODgYPbv3098fDzbtm3jxhtvrJK+uunKQkJCmDzZeTM1adIkHn74YcDZxDh69Og6HYN169aRnJxMampqhe0VFxcTFBRUpzyMMU3P19cvQF3/ficibVS1TESeAB51zb9PVY+LSCSwBhirqkVe0lc37akI6Cgiga7aWQhQ4Fo2loq1rGqJyP8G/g2YXWlRO+BSTWl90cwYDeSq6lFV/RbYCFR+wDQBZxsqOO8U7pHaruwNsHx7DpdKSgFod2t/LuZkcv5MEcu353Dq1Cm+/PJLAGbPns3zzz9PXFwcTz/9tDv9nj17+OKLLygrK2PTpk0MG+a8SQkNDeXo0aMAPPHEE2RlZZGVlUVwcDDHjh1j8uTJrF+/nt69e1coz7Fjx9i9ezcAr7/+uju/6kycOJG0tDQAdu3a5c7P83lZTVJSUnjppZfYsmVLlVrc4cOHq/R+NMY0H766frn8HChvJsoBegKoaoKqOlyf4yJyC7AZeFBVD1cq0i0iMtj1fZpHflW4nmulAVNcsx7i++d09wDv1bb/IhIPjAamqmpZpcW9cTafVssXwawb8E+PaXdbqbd1XFH7G5xtrBWIyGMisk9E9p08ebLeBTl+5vvAfc1Nt9Dxzgc58eYi9v7uEWJiYigsLOS1116jbdu2TJs2jYULF7J371527twJwMCBA5k3bx7h4eH06NGDSZMmARAbG0t6errXbS5ZsoSioiLmzp2Lw+EgKur756KhoaEkJCQQHh7O6dOnefzxx2ss/8KFC/nzn/9MREQEzzzzDGvWrKG0tJTc3FzCwsLc602dOpXBgweTk5NDSEgIiYmJAMybN49z584RExODw+Fgzpw57jRpaWnlD4eNMc2Qr65fQD+cvRTfdmX3LjCims0uxnktXiUiWSLi2TMtB3hCRA4BnYBXatmFp4FfiUiuK89EEfk3oNizmVNE3gf+hLNSky8i5c1Oq4EfArtdZVnskffdrv2oljgDasOJyBRgjKrGu6YfBH6sqvM81sl2rZPvmj7iWufr6vKNiopSbz3+ajJ02U4KzlStiXbrGETmwpE1pk1PT2fFihUkJydXWVZYWMiMGTPYsWNHncuSl5fHuHHjyM6u8WaiVhkZGSQlJbF69eoG53H58mXuuusuMjIy3L0djTHNi6+uXyLykaq676pFpCvwmqrG1LUsInIbkKyqjWrOEZHpQEjlfhT1zOOHwOuqek9N6/miZlYAdPeY9mwrrbKOiAQCN+BsY/Wp+aNDCWobUGFeUNsA5o8ObVS+Xbt25dFHH/X6o+krbdiwYY0KZOBs7ly2bJkFMmOasSt1/VLVQuBVbz+avtJUNakxgczlFuB/1baSL2pmgcBhnO2iBcBeYJqqfuaxzhNAhKrOcXUAmayqP6sp34bUzOD73kDHz1wiuGMQ80eHVnh4aowxzZUvrl+Va2atRaODGYCI3Af8J86u+WtV9QXX8CP7VHWLiLQD1gO3A6eAB1T1aE15NjSYGWNMa9Zag5lP2p1UdSuwtdK8xR7fi4Gf+mJbxhhjTGUtdgQQY4wxrYcFM2OMMX7Pgpkxxhi/Z8HMGGOM37NgZowxxu9ZMDPGGOP3LJgZY4zxexbMjDHG+D0LZsYYY/yeBTNjjDF+z4KZMcYYv2fBzBhjjN+zYGaMMcbvWTAzxhjj9yyYGWOM8XsWzIwxxvg9C2bGGGP8ngUzY4wxfs+CmTHGGL9nwcwYY4zfs2BmjDHG71kwM8YY4/csmBljjPF7FsyMMcb4PQtmV4GqMnLkSM6ePVtl2YYNG4iMjCQiIoIhQ4bw6aef1ppffHw8Bw8erDJ/5syZ9OjRA4fDgcPhICsry72spKSEO+64A4BZs2bRpUsX+vXrVyH9/PnzCQsLIzIykkmTJnHmzBkADhw4wMyZM+u+w8aYOvP19UFE1ohIHy/z54lIroioiNxUaVlbEfnY9X2tiHwlItmV1lkuIp+LyH4ReVtEOrrmR4jIuvrssy9YMLsKtm7dSv/+/enQoUOVZT169GDXrl0cOHCARYsW8dhjj9Wa35o1a+jTp8q5CsDy5cvJysoiKysLh8Phnp+RkcHQoUMBZ9BLSUmpkjYmJobs7Gz2799P7969Wbp0KQARERHk5+dz7NixuuyuMaYefH19UNV4Va16twuZwL3Al16WDXMtB1gHjPGyzg6gn6pGAoeBZ1zbOwCEiMgttRbOhyyY1UNSUhLR0dE4HA5mz55NaWkpe/fuJTIykuLiYi5cuEDfvn3Jzs4mPT2d4cOHExsbS2hoKHPmzKGsrAxw3l1NmDDB6zaGDBlCp06dABg0aBD5+fkA5OXlERYWRlxcHOHh4UyZMoWLFy8CMGLECPbt21evfUlJSWHs2LEADB8+nM6dO1dZZ9SoUQQGBlYpC8D999/Pxo0b67VNY1qy5nJ9APqKyFsi0h5ARNJFJKpyXqr6iarmVbM7Y4BtrvX+Dpzykv5vqvqda/IDIMRj8V+BB2o6Xr5mwayODh06xKZNm8jMzCQrK4uAgAA2bNjAwIEDGT9+PM8++ywLFixg+vTp7ua6PXv28PLLL3Pw4EGOHDnC5s2bAcjMzGTAgAG1bjMxMdEdcABycnKYO3cuhw4dokOHDqxatarWPH79618TGRnJL3/5Sy5fvuyen5aWxogRI+q8/2vXrq1QlqioKN5///06pzemJWtO1wfgM+AsMLcRu3Q3kF6P9WfhCn4u+4A7G7H9egtsTGIR6QxsAm4D8oCfqeppL+uVAgdck8dUdXxjttuU3vmkgOXbc/g89U3Ofbib3v0c3BDUlkuXLtGlSxcAFi9ezMCBA2nXrh0rV650p42OjqZnz54ATJ06lYyMDKZMmcKpU6e4/vrra9xuWloaiYmJZGRkuOd1797d3TQ4ffp0Vq5cyVNPPVVtHkuXLuXmm2/m22+/5bHHHuPFF19k8eLFFBQU0LlzZ9q3b1+nY/DCCy8QGBhYftcHQJcuXTh+/Hid0hvTUjXX6wOQBDwJrKjvPolIN+CUql6s4/q/Br4DNnjM/goIru+2G6NRwQxYCKSq6jIRWeiaftrLepdU1dHIbTW5dz4p4JnNB7hUUooCQX3vpt29j/Dc5Agm3t7NvV5RURHnz5+npKSE4uJirrvuOgBEpEJ+5dOBgYGUlZXRpk0bEhISePXVVwFnW3lwcDD79+8nPj6ebdu2ceONN1ZJX910ZV27dgXg2muv5eGHH2bFCud5nZKSwujRo+t0DNatW0dycjKpqakVtldcXExQUFCd8jCmJWru1wdAG7hrY4DtdVlRRGYC44B7VNVze+2ASw3cfoM0tplxAvA/ru//A0xsZH7NyvLtOVwqKQWg3a39uZiTyfkzRSzfnsOpU6f48kvnc9PZs2fz/PPPExcXx9NPfx/L9+zZwxdffEFZWRmbNm1i2LBhAISGhnL06FEAnnjiCXcHjeDgYI4dO8bkyZNZv349vXv3rlCeY8eOsXv3bgBef/11d37VKSwsBJy9o9555x1384bn87KapKSk8NJLL7Fly5YqtbjDhw9X6f1oTGvSnK8PwDQgg4ZxPy+riYiMARYA473U4noD2VVTXTmNrZn9UFULXd//BfywmvXaicg+nFXRZar6jreVROQx4DGAW25p0o4wXh0/8/2NxTU33ULHOx/kxJuLOKFKzJ86k5CQwK5du2jbti3Tpk2jtLSUIUOGsHPnTtq0acPAgQOZN28eubm53H333UyaNAmA2NhY0tPT6dWrV5VtLlmyhKKiovK2bwIDA92dO0JDQ0lISGDWrFn06dOHxx9/vMbyx8XFcfLkSVQVh8PB6tWrKS0tJTc3l7CwMPd6U6dOJT09na+//pqQkBB+85vf8MgjjzBv3jwuX75MTEwM4HzgvHr1asDZzBEbG9uIo2uMf2uu1wegL84eiq/UVH4ReRJnMLoZ2C8iW4HZQC9V/dxjvTeAEcBNIpIP/G9VTQT+AFwL7HDVCj9Q1TmuZHcD79bjcDaaVKwZellB5D2cO1vZr4H/UdWOHuueVtVOXvLopqoFItIT2ImzSnqkpu1GRUVpfXvo+drQZTspOFO1ptytYxCZC0fWmDY9PZ0VK1aQnJxcZVlhYSEzZsxgx44ddS5LXl4e48aNIzu7cTc7GRkZJCUluYNSQ1y+fJm77rqLjIwMd29HY1qb5np9EJGPVLVK78W6EJFhwHSPoNSQPK4FdgHDPHo7XnG1NjOq6r2q2s/L5y/ACRHpCuD696tq8ihw/XsUZw+Z2322B1fQ/NGhBLUNqDAvqG0A80eHNirfrl278uijj3r9UeSVNmzYsEYFMnA2ZyxbtswCmWnVWuL1QVUzGhPIXG4BFjZlIIM61MxqTCyyHCjy6ADSWVUXVFqnE3BRVS+7fmW+G5hQzY/43JpDzQy+7610/MwlgjsGMX90aIWHu8aY1qs5Xh8aUzPzZ40NZjcCb+KMxF/i7Jp/yvUDvTmqGi8iQ4A/AmU4a4L/6WpvrVFzCWbGGONPWmswa1Q7kaoWAfd4mb8PiHd9/39ARGO2Y4wxxtTERgAxxhjj9yyYGWOM8XsWzIwxxvg9C2bGGGP8ngUzY4wxfs+CmTHGGL9nwcwYY4zfs2BmjDHG71kwM8YY4/csmBljjPF7FsyMMcb4PQtmxhhj/J4FM2OMMX7Pgpkxxhi/Z8HM+JyqMnLkSK9vyt2wYQORkZFEREQwZMgQPv3001rzi4+P5+DB6t/l+uSTT/KDH/ygwrzCwkJGjRoFwJgxY+jYsSPjxo2rsE5cXByhoaH069ePWbNmUVJSAkBycjKLFy+utVzm6miq86u68wOgpKSEO+64A4BZs2bRpUsX+vXrVyH9/PnzCQsLIzIykkmTJnHmzBkADhw4wMyZM+uxx6YuLJgZn9u6dSv9+/enQ4cOVZb16NGDXbt2ceDAARYtWsRjjz1Wa35r1qyhT58+Xpft27eP06dPV5mfkpLC6NGjAedFZf369VXWiYuL4/PPP+fAgQNcunSJNWvWABAbG8tf//pXLl68WGvZTNNrqvOruvMDICMjg6FDhwIwc+ZMUlJSqqSPiYkhOzub/fv307t3b5YuXQpAREQE+fn5HDt2rM77bGpnwcy4JSUlER0djcPhYPbs2ZSWlrJ3714iIyMpLi7mwoUL9O3bl+zsbNLT0xk+fDixsbGEhoYyZ84cysrKAOfd8YQJE7xuY8iQIXTq1AmAQYMGkZ+fD0BeXh5hYWHExcURHh7OlClT3MFkxIgReHvreGlpKfPnz+ell16qsiwlJYWxY8cCcM8993D99ddXWee+++5DRBARoqOj3WUREUaMGEFycnJ9D6Gpgb+dX9WdH1Dx/Bo+fDidO3eukn7UqFEEBgZWKQvA/fffz8aNG+t9DE31LJgZAA4dOsSmTZvIzMwkKyuLgIAANmzYwMCBAxk/fjzPPvssCxYsYPr06e7mlD179vDyyy9z8OBBjhw5wubNmwHIzMxkwIABtW4zMTHRfUEAyMnJYe7cuRw6dIgOHTqwatWqGtP/4Q9/YPz48XTt2rXC/NLSUnJycqqtzVVWUlLC+vXrGTNmjHteVFQU77//fp3Sm9r54/lVztv5kZaWxogRI+q8/2vXrq1QFju/fC/wahfAXF3vfFLA8u05fJ76Juc+3E3vfg5uCGrLpUuX6NKlCwCLFy9m4MCBtGvXjpUrV7rTRkdH07NnTwCmTp1KRkYGU6ZM4dSpU15rQp7S0tJITEwkIyPDPa979+7uppvp06ezcuVKnnrqKa/pjx8/zp/+9CfS09OrLPvwww/58Y9/XOdjMHfuXIYPH86dd97pntelSxeOHz9e5zyMd/56fnmqfH4UFBTQuXNn2rdvX6dj8MILLxAYGEhcXJx7np1fvmfBrBV755MCntl8gEslpSgQ1Pdu2t37CM9NjmDi7d3c6xUVFXH+/HlKSkooLi7muuuuA5zNcZ7KpwMDAykrK6NNmzYkJCTw6quvAs5nHcHBwezfv5/4+Hi2bdvGjTfeWCV9ddOePvnkE3Jzc+nVqxcAFy9epFevXuTm5rJt27YKd9E1+c1vfsPJkyf54x//WGF+cXExQUFBdcrDeOfP51c5b+eH5/PY2qxbt47k5GRSU1MrbM/OL9+zZsZWbPn2HC6VlALQ7tb+XMzJ5PyZIpZvz+HUqVN8+eWXAMyePZvnn3+euLg4nn76aXf6PXv28MUXX1BWVsamTZsYNmwYAKGhoRw9ehSAJ554gqysLLKysggODubYsWNMnjyZ9evX07t37wrlOXbsGLt37wbg9ddfd+fnTWxsLP/617/Iy8sjLy+P9u3bk5ubC0Bqair33ntvrfu/Zs0atm/fzhtvvEGbNhX/Kxw+fLhK7zRTP/58fkH154fn87KapKSk8NJLL7Fly5YqtTg7v3zPglkrdvzMJff3a266hY53PsiJNxex93ePEBMTQ2FhIa+99hpt27Zl2rRpLFy4kL1797Jz504ABg4cyLx58wgPD6dHjx5MmjQJcAYab81/AEuWLKGoqIi5c+ficDiIiopyLwsNDSUhIYHw8HBOnz7N448/Xu99OnnyJO3atavQDHXnnXfy05/+lNTUVEJCQti+fTsAc+bM4cSJEwwePBiHw8GSJUvcadLS0oiNja339s33/P388nZ+lJaWkpubS1hYmHu9qVOnMnjwYHJycggJCSExMRGAefPmce7cOWJiYnA4HMyZM8edxs4v3xNVvdpl8CoqKkq99TAyvjN02U4KPC445bp1DCJz4cga06anp7NixQqvPf4KCwuZMWMGO3bsqHNZ8vLyGDduHNnZ2XVO401SUhL5+fksXLiwwXmcOHGCadOmkZqa2qiytHYt8fzKyMggKSmJ1atXNziPy5cvc9ddd5GRkeHu7ehLIvKRqkbVvmbLYjWzVmz+6FCC2gZUmBfUNoD5o0MblW/Xrl159NFHvf6o9UqbPn16owIZOJujfvvb3/qoRK1XSzy/hg0b1qhABs7za9myZVckkLVqqtosPwMGDFBz5b39cb4OWZqqtz2drEOWpurbH+df7SK1OmVlZXr33XfrN998U2VZUlKSRkREaL9+/XTw4MGalZVVa36PPPKIfvbZZ1Xmz5o1SyMjIzUiIkJ/8pOf6Llz59zLjh8/rjExMaqqOnr0aL3hhhs0Nja2Qvpp06Zp7969tW/fvvrwww/rt99+q6qqf/3rX3XRokVey2LnV9MD9mkzuIY39eeqF6C6jwUz01okJyfrL37xC6/LMjMz9dSpU6qqunXrVo2Ojm7wdjyD5S9/+UtdunSpe3rt2rW6YsUKVVV97733dMuWLVWC2bvvvqtlZWVaVlamDzzwgK5atUpVncHY4XDohQsXGlw24zutNZhZM6MxDeRvI1qUD/+kqly6dKlCV3EbMcX4OwtmxjSAv45o8fDDD3PzzTfz+eef8x//8R+AjZhiWoZGBTMR+amIfCYiZSJSbe8ZERkjIjkikisijXs6b8xV9M4nBQxdtpNhT/4X23c5R7RwOBykpqa6f/u0ePFiduzYwb59+1iwYIE7bfmIFgEBAe4RLYB6jWjx4osvuudVHtHCc7SL6vz3f/83x48fJzw8nE2bNgE2YoppGRpbM8sGJgN/r24FEQkAEoCxQB9gqojU7RbQmGakfESLgjOXvh/R4me/5bn/fpecnByee+454PsRLc6dO0dxcbE7fW0jWgAkJCTgcDgDZHlwKB/R4i9/+UujR7QACAgI4IEHHuDPf/4zQINGTPnd735XYb6NaGGutkYFM1U9pKo5tawWDeSq6lFV/RbYCHh/QGBMM+bPI1qoqnuEFFVly5Yt7h/+2ogppiVoih86dAP+6TGdD3ht0xCRx4DHAG655ZYrXzJj6qG6ES1OqBLzp84kJCSwa9cu94gWpaWlDBkyhJ07d9KmTRv3iBa5ubncfffdVUa0KB9n0pPniBbgrMWVd+4oH9Fi1qxZ9OnTp8YRLVSVhx56iLNnz6Kq9O/fn1deeaXaEVM+//xzzp8/7x7RYvTo0cyZM4dbb72VwYMHAzB58mT3S0zT0tLc7+sy5mqodQQQEXkPuNnLol+r6l9c66QDT6lqlS5UIjIFGKOq8a7pB4Efq+q8mrZrI4CY5qYljmhhI6a0PK11BJBaa2aqWnv7Q80KgO4e0yGuecb4lfmjQ92jwJfz9YgW3t6efCVNnz690XnYiCmmOfDJ2Iy11MwCgcPAPTiD2F5gmqp+VlOeVjMzzVH5+7mOn7lEcMcg5o8OrfA6E2OuNquZNYCITAJeBv4NeFdEslR1tIgEA2tU9T5V/U5E5gHbgQBgbW2BzJjmauLt3Sx4GdMMNSqYqerbwNte5h8H7vOY3gpsbcy2jDHGmOrYCCDGGGP8ngUzY4wxfq/ZvpxTRE4CXzYii5uAr31UHF+yctWPlat+rFz10xLLdauq/psvC+MPmm0waywR2dcce/RYuerHylU/Vq76sXK1HNbMaIwxxu9ZMDPGGOP3WnIw+79XuwDVsHLVj5Wrfqxc9WPlaiFa7DMzY4wxrUdLrpkZY4xpJfw6mDX2Tdci0kNEPnTN3yQi1/ioXJ1FZIeI/MP1bycv69wtIlken2IRmehatk5EvvBY5miqcrnWK/XY9haP+VfzeDlEZLfr771fRH7uscxnx6u2t6KLyLWufc91HYvbPJY945qfIyKjG1qGBpbrVyJy0HVsUkXkVo9lXv+eTVi2mSJy0qMM8R7LHnL93f8hIg81YZl+71GewyJyxmPZFTteIrJWRL4SEa+vOxCnla5y7xeROzyWXZFj1WKoqt9+gHAgFEgHoqpZJwA4AvQErgE+Bfq4lr0JPOD6vhp43EfleglY6Pq+EHixlvU7A6eA9q7pdcCUK3C86lQu4Hw186/a8QJ6Az9yfQ8GCoGOvjxeNZ0rHuvMBVa7vj8AbHJ97+Na/1qghyufAB8dn7qU626P8+fx8nLV9PdswrLNBP7gJW1n4Kjr306u752aokyV1v8PnGPGNsXxGg7cAWRXs/w+YBsgwCDgwyt5rFrSx69rZtqIN12LiAAjgbdc6/0PMNFHRZvgyq+u+U4BtqnqRR9tvzr1LZfb1T5eqnpYVf/h+n4c+ArnANe+VJe3onuW9S3gHtexmQBsVNXLqvoFkOvKr0nKpappHufPBzhftdQUGvMm+dHADlU9paqngR3AmKtQpqnAGz7Ybq1U9e84b1yrMwF4TZ0+ADqKSFeu3LFqMfw6mNWRtzdddwNuBM6o6neV5vvCD1W10PX9X8APa1n/Aar+Z3rB1czwexG5tonL1U5E9onIB+VNnzSj4yUi0TjvuI94zPbF8aruXPG6jutYfIPz2NQlbUPVN+9HcN7dl/P29/SVupbtJ66/z1siUv5+wyt1zOqcr6s5tgew02P2lTxetamu7Ffy/GoRGjVqflOQOrzp+mqoqVyeE6qqIlJtl1HXXVcEzlfklHsG50X9GpxddJ8GljRhuW5V1QIR6QnsFJEDOC/aDebj47UeeEhVy1yzG3y8WhoRmQ5EAXd5zK7y91TVI95zuCL+CryhqpdFZDbOmm3Nr+ZuOg8Ab6lqqce8q328TAM0+2CmV+5N10U4q/CBrjvser0Bu6ZyicgJEemqqoWui+9XNWT1M+BtVS3xyLu8lnJZRP4beKopy6WqBa5/j4rzxau3A3/mKh8vEekAvIvzRuYDj7wbfLwqqctb0cvXyRfni2dvwHkuXck3qtcpbxG5F+fNwV2qerl8fjV/T19dnGstm6oWeUyuwfmMtDztiEpp05uiTB4eAJ7wnHGFj1dtqiv7lTpWLUZraGbcC/xInD3xrsF58m5RVQXScD6vAngI8FVNb4srv7rkW6W93nVBL39ONRHw2vPpSpRLRDqVN9OJyE3AUODg1T5err/d2zifJ7xVaZmvjpfXc6WGsk4BdrqOzRbgAXH2duwB/AjY08By1LtcInI78EdgvKp+5THf69/TR+Wqa9m6ekyOBw65vm8HRrnK2AkYRcUWiitWJle5wnB2ptjtMe9KH6/abAFmuHo1DgK+cd2sXalj1XJc7R4ojfkAk3C2HV8GTgDbXfODga0e690HHMZ5d/Vrj/k9cV5wcoE/Adf6qFw3AqnAP4D3gM6u+VE438Bdvt5tOO+42lRKvxM4gPOinAT8oKnKBQxxbftT17+PNIfjBUwHSoAsj4/D18fL27mCs8lyvOt7O9e+57qORU+PtL92pcsBxvr4XK+tXO+5/g+UH5sttf09m7BsS4HPXGVIA8I80s5yHctc4OGmKpNr+jlgWaV0V/R44bxxLXSdy/k4n2/OAea4lguQ4Cr3ATx6aV+pY9VSPjYCiDHGGL/XGpoZjTHGtHAWzIwxxvg9C2bGGGP8ngUzY4wxfs+CmTHGGL9nwcwYY4zfs2BmjDHG71kwM8YY4/f+f92M3WGf+jRwAAAAAElFTkSuQmCC\n", "text/plain": [""]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import matplotlib.pyplot as plt\n", "fig, ax = plt.subplots(1, 1, figsize=(6, 3))\n", "a = numpy.arange(0, 12) * (-2 * numpy.pi / 12)\n", "X = numpy.vstack([numpy.cos(a), numpy.sin(a)]).T\n", "ax.plot(X[:, 0], X[:, 1], 'o');\n", "for i in range(0, 12):\n", " ax.text(X[i, 0], X[i, 1], \"exp(-2pi %d/12)\" % i)\n", "ax.set_title('unit roots');"]}, {"cell_type": "markdown", "id": "f257dc6e", "metadata": {}, "source": ["Then:\n", "\n", "$$\n", "\\begin{array}{rcl}\n", "F_{a,k + \\frac{K}{2}} &=& \\sum_{n=1}^{N} X_{an} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{n\\left(k + \\frac{K}{2}\\right)} \\\\\n", "&=&\\sum_{n=1}^{N} X_{an} (-1)^n \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{nk} \\\\\n", "&=&\\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} - \\sum_{m=1}^{\\frac{N}{2}} X_{a,2m-1} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{(2m-1)k} \\\\\n", "&=&\\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} - \\sum_{m=1}^{\\frac{N}{2}} X_{a,2m-1} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{-k}\n", "\\end{array}\n", "$$\n", "\n", "Then:\n", "\n", "$$\n", "\\begin{array}{rcl}\n", "F_{a,k} + F_{a,k+\\frac{K}{2}} &=& 2\\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk}\n", "= 2\\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{\\frac{K}{2}}\\right)^{mk}\n", "\\end{array}\n", "$$\n", "\n", "Finally:\n", "\n", "$$\n", "\\begin{array}{rcl}\n", "F_{a,k} &=& \\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} + \\sum_{m=1}^{\\frac{N}{2}} X_{a,2m-1} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{-k} \\\\\n", "F_{a,k + \\frac{K}{2}} &=&\\sum_{m=1}^{\\frac{N}{2}} X_{a,2m} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} - \\sum_{m=1}^{\\frac{N}{2}} X_{a,2m-1} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{2mk} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{-k}\n", "\\end{array}\n", "$$"]}, {"cell_type": "markdown", "id": "356b585d", "metadata": {}, "source": ["Now, what happen when *K* is odd, fallback to the original computation.\n", "\n", "$$\n", "F_{ak} = \\sum_n X_{an} M_{nk} = \\sum_n X_{an} \\exp\\left(\\frac{-2i\\pi}{K}\\right)^{nk}\n", "$$"]}, {"cell_type": "code", "execution_count": 25, "id": "971be7bd", "metadata": {"scrolled": false}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["1.5 s \u00b1 0 ns per loop (mean \u00b1 std. dev. of 1 run, 1 loop each)\n"]}], "source": ["import functools\n", "\n", "\n", "def cooley_fft_2p(x, fft_length):\n", " cst = _dft_cst_power(x.shape[-1], fft_length, x.dtype)\n", " return numpy.matmul(x, cst)\n", "\n", "\n", "@functools.cache\n", "def _build_fact(p2_2, fft_length, dtype):\n", " first = numpy.exp(-2j * numpy.pi / fft_length)\n", " fact = numpy.ones(p2_2, dtype=dtype)\n", " for k in range(1, p2_2):\n", " fact[k] = fact[k-1] * first\n", " return fact.reshape((1, -1))\n", "\n", "\n", "def build_fact(p2_2, fft_length, dtype):\n", " return _build_fact(p2_2, fft_length, dtype)\n", "\n", "\n", "def cooley_fft_recursive(x, fft_length):\n", " if len(x.shape) != 2:\n", " raise RuntimeError(\n", " \"Unexpected x.shape=%r.\" % (x.shape, ))\n", " dtype = numpy.complex128 if x.dtype == numpy.float64 else numpy.complex64\n", " if fft_length == 1:\n", " return x[:, :1].astype(dtype)\n", "\n", " if fft_length % 2 == 0:\n", " def split(x):\n", " even = x[:, ::2]\n", " odd = x[:, 1::2]\n", " return even, odd\n", "\n", " def tmp1(even, odd, fft_length):\n", " p2_2 = fft_length // 2\n", " fft_even = cooley_fft_recursive(even, p2_2)\n", " fft_odd = cooley_fft_recursive(odd, p2_2)\n", " return fft_even, fft_odd, p2_2\n", "\n", " def tmp2(x, fft_even, fft_odd, p2_2):\n", " fact = build_fact(p2_2, fft_length, fft_even.dtype)\n", "\n", " fact_odd = fft_odd * fact\n", " return numpy.hstack([fft_even + fact_odd, fft_even - fact_odd])\n", "\n", " # inplace\n", " # result = numpy.empty((x.shape[0], fft_length), dtype=fft_even.dtype)\n", " # numpy.multiply(fft_odd, fact, out=result[:, :p2_2])\n", " # numpy.subtract(fft_even, result[:, :p2_2], out=result[:, p2_2:])\n", " # numpy.add(fft_even, result[:, :p2_2], out=result[:, :p2_2])\n", " # return result\n", " \n", " even, odd = split(x)\n", " fft_even, fft_odd, p2_2 = tmp1(even, odd, fft_length)\n", " result = tmp2(x, fft_even, fft_odd, p2_2)\n", " else:\n", " result = cooley_fft_2p(x, fft_length)\n", " \n", " return result\n", "\n", "\n", "\n", "def cooley_fft(x, fft_length):\n", " return cooley_fft_recursive(x, fft_length)\n", "\n", "\n", "def custom_fft_cooley(x, fft_type, length, axis):\n", " # https://github.com/numpy/numpy/blob/4adc87dff15a247e417d50f10cc4def8e1c17a03/numpy/fft/_pocketfft.py#L56\n", " if fft_type == 'FFT':\n", " if x.shape[axis] > length:\n", " # fft_length > shape on the same axis\n", " # the matrix is shortened\n", " slices = [slice(None)] * len(x.shape)\n", " slices[axis] = slice(0, length)\n", " new_x = x[tuple(slices)]\n", " elif x.shape[axis] == length:\n", " new_x = x\n", " else:\n", " # other, the matrix is completed with zeros\n", " shape = list(x.shape)\n", " shape[axis] = length\n", " slices = [slice(None)] * len(x.shape)\n", " slices[axis] = slice(0, length)\n", " zeros = numpy.zeros(tuple(shape), dtype=x.dtype)\n", " index = [slice(0, i) for i in x.shape]\n", " zeros[tuple(index)] = x\n", " new_x = zeros\n", "\n", " if axis == len(new_x.shape) - 1:\n", " if len(new_x.shape) != 2:\n", " xt = new_x.reshape((-1, new_x.shape[-1]))\n", " else:\n", " xt = new_x\n", " res = cooley_fft(xt, length)\n", " if len(new_x.shape) != 2:\n", " res = res.reshape(new_x.shape[:-1] + (-1, ))\n", " else:\n", " perm = numpy.arange(len(x.shape)).tolist() \n", " perm[axis], perm[-1] = perm[-1], perm[axis] \n", " rest = new_x.transpose(perm)\n", " shape = rest.shape[:-1]\n", " rest = rest.reshape((-1, rest.shape[-1]))\n", " res = cooley_fft(rest, length)\n", " res = res.reshape(shape + (-1, )).transpose(perm)\n", " perm[axis], perm[0] = perm[0], perm[axis]\n", " return res\n", " raise ValueError(\"Unexpected value for fft_type=%r.\" % fft_type)\n", "\n", "\n", "def custom_fftn_cooley(x, fft_type, fft_length, axes):\n", " if len(axes) != len(fft_length):\n", " raise ValueError(\"Length mismatch axes=%r, fft_length=%r.\" % (\n", " axes, fft_length))\n", " if fft_type == 'FFT':\n", " res = x\n", " for i in range(len(fft_length) - 1, -1, -1):\n", " length = fft_length[i]\n", " axis = axes[i]\n", " res = custom_fft_cooley(res, fft_type, length, axis)\n", " return res\n", " raise ValueError(\"Unexpected value for fft_type=%r.\" % fft_type)\n", " \n", "\n", "shape = (4, )\n", "fft_length = [3,]\n", "axes = [0]\n", "rnd = numpy.random.randn(*shape) + numpy.random.randn(*shape) * 1j\n", "assert_almost_equal(custom_fftn_cooley(rnd, 'FFT', fft_length, axes),\n", " numpy_fftn(rnd, 'FFT', fft_length, axes),\n", " decimal=5)\n", "%timeit -n 1 -r 1 test_fct(numpy_fftn, custom_fftn_cooley)"]}, {"cell_type": "code", "execution_count": 26, "id": "441d6c41", "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 24/24 [00:10<00:00, 2.35it/s]\n"]}, {"data": {"text/html": ["\n", "\n", "
\n", " \n", " \n", " | name | \n", " custom_fftn_cooley | \n", " custom_fftn_power | \n", " numpy_fftn | \n", " torch_fftn | \n", "
\n", " \n", " | length | \n", " | \n", " | \n", " | \n", " | \n", "
\n", " \n", " \n", " \n", " | 8 | \n", " 0.002873 | \n", " 0.000685 | \n", " 0.001482 | \n", " 0.005463 | \n", "
\n", " \n", " | 16 | \n", " 0.007197 | \n", " 0.002121 | \n", " 0.001922 | \n", " 0.005063 | \n", "
\n", " \n", " | 24 | \n", " 0.009443 | \n", " 0.002903 | \n", " 0.002739 | \n", " 0.005169 | \n", "
\n", " \n", " | 32 | \n", " 0.012783 | \n", " 0.002556 | \n", " 0.002003 | \n", " 0.004076 | \n", "
\n", " \n", " | 40 | \n", " 0.014142 | \n", " 0.003916 | \n", " 0.003937 | \n", " 0.005118 | \n", "
\n", " \n", "
\n", "
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["piv.plot(logy=True, logx=True, title=\"FFT benchmark 3\", figsize=(12, 4));"]}, {"cell_type": "code", "execution_count": 28, "id": "eb4f30d6", "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:11<00:00, 1.15s/it]\n"]}, {"data": {"text/html": ["\n", "\n", "
\n", " \n", " \n", " | name | \n", " custom_fftn_cooley | \n", " custom_fftn_power | \n", " numpy_fftn | \n", " torch_fftn | \n", "
\n", " \n", " | length | \n", " | \n", " | \n", " | \n", " | \n", "
\n", " \n", " \n", " \n", " | 2 | \n", " 0.000575 | \n", " 0.000471 | \n", " 0.000722 | \n", " 0.019371 | \n", "
\n", " \n", " | 4 | \n", " 0.001153 | \n", " 0.000328 | \n", " 0.001130 | \n", " 0.018366 | \n", "
\n", " \n", " | 8 | \n", " 0.003678 | \n", " 0.000624 | \n", " 0.001779 | \n", " 0.019295 | \n", "
\n", " \n", " | 16 | \n", " 0.006843 | \n", " 0.002255 | \n", " 0.002192 | \n", " 0.020169 | \n", "
\n", " \n", " | 32 | \n", " 0.015574 | \n", " 0.003045 | \n", " 0.002736 | \n", " 0.017193 | \n", "
\n", " \n", "
\n", "