fft_explorations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ddfad7bb-e6c0-42cd-be9e-b1af91dd10d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "bc33978c-c010-42cd-a168-f8dd1acf264a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f1(n):\n",
    "    return 5\n",
    "\n",
    "def f2(n):\n",
    "    return 2 * math.sin( 2 * math.pi * n )#- math.pi/2)\n",
    "\n",
    "def f3(n):\n",
    "    return 3 * math.cos(4 * math.pi * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9b635013-906b-4b84-9346-10cb756e0090",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = [b/4 for b in range(4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "000766a9-17e9-41ab-b8fc-924cbbd85138",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[8.0, 4.0, 8.0, 0.0]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[f1(n) + f2(n) + f3(n) for n in n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "1240af4c-5f48-4ee3-88c9-0f15965cc14b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "df8eaa3b-2500-453a-b23c-556e173102f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5, 5, 5, 5]\n",
      "[0, 2, 0, -2]\n",
      "[3, -3, 3, -3]\n"
     ]
    }
   ],
   "source": [
    "funcs = [f1, f2, f3]\n",
    "for f in funcs:\n",
    "    print([int(f(n)) for n in n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c8a46ecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def DFT_slow(x):\n",
    "    \"\"\"compute the DFT of the 1D array x\"\"\"\n",
    "    x = np.asarray(x, dtype=float)\n",
    "    N = x.shape[0]\n",
    "    n = np.arange(N)\n",
    "    k = n.reshape((N,1))\n",
    "    M = np.exp(-2j * np.pi * k * n /N)\n",
    "    return np.dot(M, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "77523440",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102 ms ± 6.56 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "13.4 µs ± 304 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "x = np.random.random(1024)\n",
    "np.allclose(DFT_slow(x), np.fft.fft(x))\n",
    "%timeit DFT_slow(x)\n",
    "%timeit np.fft.fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f0070ae",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18e0e793",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "283cb078",
   "metadata": {},
   "outputs": [],
   "source": [
    "def FFT(x):\n",
    "    \"\"\"A recursive implementation of the 1D Cooley-Tukey FFT\"\"\"\n",
    "    x = np.asarray(x, dtype=float)\n",
    "    N = x.shape[0]\n",
    "\n",
    "    if N % 2 > 0:\n",
    "        raise ValueError(\"size of x must be a power of 2\")\n",
    "    elif N <= 32:\n",
    "        return DFT_slow(x)\n",
    "    else:\n",
    "        X_even = FFT(x[::2])\n",
    "        X_odd = FFT(x[1::2])\n",
    "        factor = np.exp(-2j * np.pi * np.arange(N)/N)\n",
    "        return np.concatenate([X_even + factor[:N // 2] * X_odd,\n",
    "                               X_even + factor[N // 2:] * X_odd])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "8c78ec40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.random.random(1024)\n",
    "np.allclose(FFT(x), np.fft.fft(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "4a0f03cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102 ms ± 11.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "2.39 ms ± 51.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "12.8 µs ± 274 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit DFT_slow(x)\n",
    "%timeit FFT(x)\n",
    "%timeit np.fft.fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "id": "bc5f0d55",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Permute2(x):\n",
    "    N = len(x)\n",
    "    j = 0\n",
    "    m=0\n",
    "    for i in range(0, N):\n",
    "        print(f\" i:{i}\\n j:{j}\\n m:{m} \\n x:{x}\" )\n",
    "        print(i)\n",
    "        if (j > i):\n",
    "            x[i], x[j] = x[j], x[i]\n",
    "        \n",
    "        m = N\n",
    "        \n",
    "        while (m>=2 and j > m):\n",
    "            j -= m\n",
    "            m = m//2\n",
    "        \n",
    "        j += m \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dab32de8",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = list(range(32))\n",
    "t.insert(0,0)\n",
    "t.append(0)\n",
    "t = Permute(t)\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "00c5ed3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 2]\n"
     ]
    }
   ],
   "source": [
    "a = [2,3]\n",
    "a[0], a[1] = a[1], a[0]\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "id": "02771062",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Bit Reversal\"\"\"\n",
    "\n",
    "\n",
    "def Permute(x):\n",
    "    m = 0\n",
    "    N = len(x)-2\n",
    "    i = 1\n",
    "    for n in range(1, N+1):\n",
    "        if n > i:\n",
    "            #print(i, n, x[i], x[n])\n",
    "            x[i], x[n] = x[n], x[i]\n",
    "        print(f\" n:{n}\\n i:{i}\\n m:{m} \\n x:{x}\" )\n",
    "        \n",
    "        m = N//2\n",
    "        while m >=2 and i > m: \n",
    "            i = i - m\n",
    "            m = m //2\n",
    "        \n",
    "        i = i + m \n",
    "    return x\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "id": "10aa2b3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def permute_zerod(x):\n",
    "    m = 0\n",
    "    N = len(x)\n",
    "    i = 0\n",
    "    for n in range(N):\n",
    "        if n > i:\n",
    "            x[i], x[n] = x[n], x[i]\n",
    "\n",
    "        m = N//2\n",
    "        while m >=2 and i >= m:\n",
    "            i -= m\n",
    "            m //= 2\n",
    "        i += m \n",
    "    \n",
    "    return x\n",
    "\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "id": "a2893e1a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15]\n"
     ]
    }
   ],
   "source": [
    "t = permute_zerod(list(range(16)))\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1133fc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nth_roots_of_unity(N):\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d510d004",
   "metadata": {},
   "outputs": [],
   "source": [
    "#DANIELSSONLANCZOS FFT\n",
    "\n",
    "def danielsonLanczos(x, w, N):\n",
    "    M = 1\n",
    "    while N > M:\n",
    "        Istep = M << 1\n",
    "        for m in range(1, M):\n",
    "            for i in range(m, N, Istep):\n",
    "                j = i + M\n",
    "                Temp = Wm * Xj\n",
    "                Xj = xi - Temp\n",
    "                xi = xi + Temp\n",
    "            \n",
    "        M = Istep"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 ('fft_hdl': venv)",
   "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.10.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "44613a307f027ce3c5c0ffaa502fd519eb545912dd13be465f11902b1a6911bd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}