{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementing DFT\n",
    "\n",
    "Copyright 2019 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with a known result.  The DFT of an impulse is a constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 4\n",
    "x = [1, 0, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fft.fft(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Literal translation\n",
    "\n",
    "The usual way the DFT is expressed is as a summation.  Following the notation on [Wikipedia](https://en.wikipedia.org/wiki/Discrete_Fourier_transform):\n",
    "\n",
    "$ X_k = \\sum_{n=0}^N x_n \\cdot e^{-2 \\pi i n k / N} $\n",
    "\n",
    "Here's a straightforward translation of that formula into Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi = np.pi\n",
    "exp = np.exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1+0j)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = 0\n",
    "sum(x[n] * exp(-2j * pi * k * n / N) for n in range(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wrapping this code in a function makes the roles of `k` and `n` clearer, where `k` is a free parameter and `n` is the bound variable of the summation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft_k(x, k):\n",
    "    return sum(x[n] * exp(-2j * pi * k * n / N) for n in range(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we usually we compute $X$ all at once, so we can wrap this function in another function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(x):\n",
    "    N = len(x)\n",
    "    X = [dft_k(x, k) for k in range(N)]\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1+0j), (1+0j), (1+0j), (1+0j)]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dft(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the results check out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DFT as a matrix operation\n",
    "\n",
    "It is also common to express the DFT as a [matrix operation](https://en.wikipedia.org/wiki/DFT_matrix):\n",
    "\n",
    "$ X = W x $\n",
    "\n",
    "with \n",
    "\n",
    "$ W_{j, k} = \\omega^{j k} $\n",
    "\n",
    "and\n",
    "\n",
    "$ \\omega = e^{-2 \\pi i / N}$\n",
    "\n",
    "If we recognize the construction of $W$ as an outer product, we can use `np.outer` to compute it.\n",
    "\n",
    "Here's an implementation of DFT using outer product to construct the DFT matrix, and dot product to compute the DFT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(x):\n",
    "    N = len(x)\n",
    "    ks = range(N)\n",
    "    args = -2j * pi * np.outer(ks, ks) / N\n",
    "    W = exp(args)\n",
    "    X = W.dot(x)\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the results check out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dft(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing FFT\n",
    "\n",
    "Finally, we can implement the FFT by translating from math notation:\n",
    "\n",
    "$ X_k = E_k + e^{-2 \\pi i k / N} O_k $\n",
    "\n",
    "Where $E$ is the FFT of the even elements of $x$ and $O$ is the DFT of the odd elements of $x$.\n",
    "\n",
    "Here's what that looks like in code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft(x):\n",
    "    N = len(x)\n",
    "    if N == 1:\n",
    "        return x\n",
    "    \n",
    "    E = fft(x[::2])\n",
    "    O = fft(x[1::2])\n",
    "    \n",
    "    ks = np.arange(N)\n",
    "    args = -2j * pi * ks / N\n",
    "    W = np.exp(args)\n",
    "    \n",
    "    return np.tile(E, 2) + W * np.tile(O, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The length of $E$ and $O$ is half the length of $W$, so I use `np.tile` to double them up.\n",
    "\n",
    "And the results check out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}