{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import time\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['text.usetex'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MA934\n",
    "\n",
    "## Fourier transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fourier Transform pairs\n",
    "\n",
    "Fourier transform represents a function as a sum of harmonics. It is a linear invertible transformation between the time-domain, $h(t)$, and the frequency domain, $H(f)$, representations of a function:\n",
    "$$\n",
    "\\begin{align*}\n",
    "H(f) &= \\int_{-\\infty}^\\infty h(t)\\,\\mathrm{e}^{2\\,\\pi\\,i\\, f\\, t}\\,dt\\\\\n",
    "h(t) &=  \\int_{-\\infty}^\\infty H(f)\\,\\mathrm{e}^{-2\\,\\pi\\,i\\,f\\, t}\\,df.\n",
    "\\end{align*}\n",
    "$$\n",
    "Different notations used due to 1-parameter family of equivalent normalisations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: Gaussian function\n",
    "\n",
    "$$\n",
    "h(t) = \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}}\\,\\mathrm{e}^{-\\frac{t^2}{2\\,\\sigma^2}}.\n",
    "$$\n",
    "FT can be calculated analytically by \"completing the square\" in the exponent: \n",
    "$$\n",
    "\\begin{align*}\n",
    "H(f) &= \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}}\\,\\int_{-\\infty}^\\infty \\mathrm{e}^{-\\frac{t^2}{2\\,\\sigma^2}} \\,\\mathrm{e}^{2\\,\\pi\\,i\\, f\\, t}\\,dt\\\\\n",
    "&=  \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}}\\,\\int_{-\\infty}^\\infty  \\mathrm{e}^{-\\frac{1}{2\\,\\sigma^2}\\left[t^2 - 4\\pi i \\sigma^2 t f \\right]} d t\\\\\n",
    "&=  \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}}\\,\\int_{-\\infty}^\\infty  \\mathrm{e}^{-\\frac{1}{2\\,\\sigma^2}\\left[t^2 - 2(2\\pi i \\sigma^2 f) t + (2\\pi i \\sigma^2 f)^2 - (2\\pi i \\sigma^2 f)^2 \\right]} d t\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: Gaussian function (continued)\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&=  \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}} \\left[ \\int_{-\\infty}^\\infty  \\mathrm{e}^{-\\frac{1}{2\\,\\sigma^2}(t - 2\\pi i \\sigma^2 f)^2} d t \\right]\\,\\mathrm{e}^{\\frac{1}{2\\,\\sigma^2}\\,(2\\pi i \\sigma^2 f)^2}\\\\\n",
    "&=  \\frac{1}{\\sqrt{2\\,\\pi\\,\\sigma^2}} \\left[ \\sqrt{2\\,\\pi\\,\\sigma^2}\\right]\\,\\mathrm{e}^{-2\\pi^2 \\sigma^2 f^2}\\\\\n",
    "&= \\mathrm{e}^{-2\\pi^2 \\sigma^2 f^2}.\n",
    "\\end{align*}\n",
    "$$\n",
    "where we let have $z = t - 2\\pi i \\sigma^2 f$ and used the fact that\n",
    "$$\n",
    "\\int_{-\\infty}^\\infty \\mathrm{e}^{-\\frac{z^2}{2\\,\\sigma^2}} = \\sqrt{2\\,\\pi\\,\\sigma^2}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fourier tranform of real-valued functions\n",
    "\n",
    "Although $H(f)$ in this example turned out to be real, the Fourier Transform of a real-valued function is generally complex. \n",
    "\n",
    "If $h(t)$ is real then, under complex conjugation, $H(f)$ has the symmetry \n",
    "\n",
    "$$\n",
    "H^*(f) = H(-f).\n",
    "$$\n",
    "\n",
    "To see this:\n",
    "\n",
    "$$\n",
    "H^*(f) = \\int_{-\\infty}^\\infty h(t)^*\\,\\mathrm{e}^{-2\\,\\pi\\,i\\, f\\, t}\\,dt = \\int_{-\\infty}^\\infty h(t)\\,\\mathrm{e}^{2\\,\\pi\\,i\\, (-f)\\, t}\\,dt = H(-f).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Applications of the Fourier Transform\n",
    "\n",
    "* Signal analysis : *power spectrum*\n",
    "$$\n",
    "P(f) = H(f)\\,H^*(f),\n",
    "$$\n",
    "summarises the frequencies present in a signal.\n",
    "* Efficient computation of convolutions,\n",
    "$$\n",
    "(f*g)(t) = \\int_{-\\infty}^\\infty f(\\tau)\\,g(t-\\tau)\\,d\\tau.\n",
    "$$\n",
    "* Differentiation: derivative in time domain is equivalent to multiplication by $f$ in the frequency domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Translation property\n",
    "Let $g(t) = h(t+\\tau).$  \n",
    "Writing both sides in terms of the Fourier transforms we have\n",
    "\n",
    "$$\n",
    " \\int_{-\\infty}^\\infty G(f)\\,\\mathrm{e}^{-2\\,\\pi\\,i\\,f\\, t}\\,df =  \\int_{-\\infty}^\\infty H(f)\\,\\mathrm{e}^{-2\\,\\pi\\,i\\,f\\, (t+\\tau)}\\,df,\n",
    "$$\n",
    "\n",
    "from which we conclude that\n",
    "\n",
    "$$\n",
    "G(f) = \\mathrm{e}^{-2\\pi i \\tau}\\,h(f).\n",
    "$$\n",
    "\n",
    "Translation in the time-domain therefore corresponds to multiplication (by a complex number) in the frequency domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convolution theorem\n",
    "$$\n",
    "\\begin{align*}\n",
    "(f*g)(t) &= \\int_{-\\infty}^\\infty f(\\tau)\\,g(t-\\tau)\\,d\\tau\\\\\n",
    "&=  \\int_{-\\infty}^\\infty \\left[\\int_{-\\infty}^\\infty F(f_1)\\,\\mathrm{e}^{-2\\,\\pi\\,i\\,f_1\\, \\tau}\\,df_1\\right] \\left[ \\int_{-\\infty}^\\infty G(f_2)\\,\\mathrm{e}^{-2\\,\\pi\\,i\\,f_2\\, (t-\\tau)}\\,df_2 \\right] d\\tau\\\\\n",
    "&= \\int_{-\\infty}^\\infty F(f_1)\\, G(f_2)\\, \\mathrm{e}^{-2\\,\\pi\\,i\\,f_2\\, t} \\left[ \\int_{-\\infty}^\\infty \\mathrm{e}^{-2\\,\\pi\\,i\\,(f_1-f_2)\\,\\tau} d\\tau \\right]\\,df_1 df_2\\\\\n",
    "& =  \\int_{-\\infty}^\\infty F(f_1)\\, G(f_1)\\, \\mathrm{e}^{-2\\,\\pi\\,i\\,f_1\\, t} df_1.\n",
    "\\end{align*}\n",
    "$$\n",
    "Convolution in time domain is pointwise product in the frequency domain.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Discrete sampling and Nyquist frequency\n",
    "\n",
    "Sampling: $h(t)$, represented by discrete values, $H = \\left\\{h_n = h(t_n),\\  n\\in \\mathbb{Z} \\right\\}$ measured at discrete points $T = \\left\\{t_n,\\ n \\in \\mathbb{Z} \\right\\}$. \n",
    "\n",
    "Assume samples uniformly spaced in time: $t_n = n\\,\\Delta$. \n",
    "\n",
    "$\\Delta$ is called the sampling interval and its reciprocal, $1/\\Delta$, is the sampling rate. \n",
    "\n",
    "The quantity\n",
    "$$\n",
    "f_c = \\frac{1}{2\\Delta}.\n",
    "$$\n",
    "is known as the Nyquist frequency. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sampling Theorem\n",
    "\n",
    "How frequently must we sample to represent $h(t)$ to a given degree of accuracy?\n",
    "\n",
    "**Shannon Sampling Theorem**: Let $h(t)$ be a function and  $\\left\\{h_n\\right\\}$ be the samples of $h(t)$ obtained with sampling interval $\\Delta$. If $H(f) = 0$ for all $\\left| f \\right| \\geq f_c$ then $h(t)$ is *completely* determined by its samples. \n",
    "\n",
    "An explicit reconstruction is provided by the Whittaker-Shannon interpolation formula:\n",
    "$$\n",
    "h(t) = \\Delta \\sum_{n=-\\infty}^\\infty h_n\\, \\frac{\\sin \\left[2\\pi f_c\\,(t-n\\Delta )\\right]}{\\pi\\,(t-n\\Delta)}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Aliasing\n",
    "\n",
    "From the samples alone, it impossible *in principle* to determine how much power is outside the Nyquist interval.\n",
    "\n",
    "The process of sampling moves frequencies which lie * outside* of the Nyquist interval $[-f_c, f_c]$ into this interval via a process known as aliasing.\n",
    "\n",
    "To see this, suppose we have two harmonics, having frequencies $f_1$ and $f_2$, that agree at all their sample points:\n",
    "\n",
    "$$\n",
    "\\mathrm{e}^{2\\pi i f_1 n \\Delta} = \\mathrm{e}^{2\\pi i f_1 n \\Delta}\\ \\ \\  \\forall n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Aliasing\n",
    "$$\n",
    "\\begin{align*}\n",
    "&\\mathrm{e}^{2\\pi i f_1 n \\Delta} = \\mathrm{e}^{2\\pi i f_2 n \\Delta} &\\forall n\\\\\n",
    "\\Rightarrow &\\mathrm{e}^{2\\pi i (f_1-f_2) n \\Delta} = 1  &\\forall n\\\\\n",
    "\\Rightarrow &\\mathrm{e}^{2\\pi i (f_1-f_2) \\Delta} = 1\\\\\n",
    "\\Rightarrow & (f_1-f_2)\\,\\Delta = m &m \\in \\mathbb{Z}\\\\\n",
    "\\Rightarrow & (f_1-f_2) = \\frac{m}{\\Delta} &m \\in \\mathbb{Z}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The frequencies do not have to be equal! They can differ by an integer multiple of the width of the Nyquist interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Discrete Fourier Transform (DFT)\n",
    "\n",
    "Suppose with have $N$, samples of $h(t)$: $\\left\\{h_n= h(t_n),\\  n=0,\\ldots N-1 \\right\\}$. DFT can be thought of as a discrete version continuous FT obtained by approximating the integral for $H(f)$ using the $N$ values, $h_n$.\n",
    "\n",
    "At what points should we approximate $f$? From the sampling theorem, if the sampling rate is high enough, the Nyquist interval, $[-f_c, f_c]$ contains all the information needed to reconstruct $h(t)$.\n",
    "\n",
    "Therefore we should estimate $H(f)$ at the freqencies\n",
    "\n",
    " $f_n = \\frac{n}{\\Delta N}$ for $n = -\\frac{N}{2},\\ldots +\\frac{N}{2}$,\n",
    " \n",
    " that uniformly discretise the Nyquist interval.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Discrete Fourier Transform \n",
    "\n",
    "Replace integral for $H(f_n)$ by a Riemann sum:\n",
    "$$\n",
    "\\begin{align*}\n",
    "H(f_n) &\\approx \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{2\\pi i f_n t_k} \\Delta \\\\\n",
    "& = \\Delta\\, \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{2\\pi i \\left(\\frac{n}{N\\Delta}\\right) (k\\Delta)}\\\\\n",
    "&= \\Delta\\, \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{\\frac{2\\pi i n k}{N}}\\\\\n",
    "&\\equiv \\Delta\\, H_n.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The array of numbers, $\\left\\{H_n\\right\\}$, defined as\n",
    "$$\n",
    "H_n = \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{\\frac{2\\pi i n k}{N}}\n",
    "$$\n",
    "is called the Discrete Fourier Transform (DFT) of the array $\\left\\{h_n\\right\\}$.\n",
    "\n",
    "The inverse is obtained by changing the sign of the exponent.\n",
    "\n",
    "As written the index of $H_n$ runs over $n = -\\frac{N}{2},\\ldots +\\frac{N}{2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Periodicity of DFT\n",
    "\n",
    "When considered as a function of $n$, the DFT is periodic with period $N$. That is:\n",
    "\n",
    "$$\n",
    "H_{n+N} = H_n.\n",
    "$$\n",
    "\n",
    "Thus only $N$ of the $N+1$ values are independent.\n",
    "\n",
    "We can see this by direct computation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "H_{n+N} &= \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{\\frac{2\\pi i (n+N) k}{N}}\\\\\n",
    "&= \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{\\frac{2\\pi i n k}{N}}\\, \\mathrm{e}^{2\\pi i k}\\\\\n",
    "& = \\sum_{k=0}^{N-1} h_k\\,\\mathrm{e}^{\\frac{2\\pi i n k}{N}}\\\\\n",
    "& = H_n.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### DFT indexing\n",
    "\n",
    "We can use periodicity to shift the index, $n$, by $\\frac{N}{2}$ to get an index running over  $n = 0\\ldots N$ rather than $n = -\\frac{N}{2},\\ldots +\\frac{N}{2}$. \n",
    "\n",
    "This is more convenient for a computer and most implementations of the DFT adopt this convention. \n",
    "\n",
    "One consequence of this convention is that if we want to plot the DFT as a function of the frequencies (for example to compare against an analytic expression for the corresponding continuous Fourier transform) we need to obtain the negative frequencies from $H_{-n} = H_{N-n}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### FFT in Python\n",
    "\n",
    "FFT is made available in Python via the ```pyfftw``` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "cell_style": "split"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.6881789171900392+0.5916479882309668j)\n",
      "(0.9047019733848708-0.4387366153747087j)\n",
      "(0.8277040544279829+0.2793681357976642j)\n",
      "(0.27539319505285476-0.002330991182591074j)\n",
      "(1.1955843835618385+1.6491945185167878j)\n",
      "(0.5326351908397333-1.9874427848958942j)\n",
      "(-0.7180623891102259+1.1617246580474099j)\n",
      "(0.13139467893387619+0.6283549973711337j)\n",
      "===========================================\n",
      "(2.461172169900891+1.8817799065107685j)\n",
      "(-1.955712910010033-1.4271632722049337j)\n",
      "(-2.654439605405103-0.2307995773351923j)\n",
      "(-0.513537101037336-1.4177563328076226j)\n",
      "(-1.2270779065217787+5.482090694674889j)\n",
      "(-3.5765267359932134-3.7794626754431255j)\n",
      "(3.449967207513188+1.8302990031405537j)\n",
      "(-1.4892764559669278+2.394196159312399j)\n",
      "===========================================\n",
      "(2.461172169900891+1.8817799065107685j)\n",
      "(-1.9557129100100332-1.4271632722049339j)\n",
      "(-2.654439605405103-0.2307995773351923j)\n",
      "(-0.513537101037336-1.4177563328076226j)\n",
      "(-1.2270779065217787+5.482090694674889j)\n",
      "(-3.576526735993214-3.7794626754431255j)\n",
      "(3.449967207513188+1.8302990031405537j)\n",
      "(-1.4892764559669278+2.3941961593123984j)\n"
     ]
    }
   ],
   "source": [
    "# Execute \"conda install -c conda-forge pyfftw\" in a terminal first\n",
    "import pyfftw\n",
    "\n",
    "# Source and further guidance: https://hgomersall.github.io/pyFFTW/sphinx/tutorial.html\n",
    "\n",
    "# Create and fill a complex array, a, of length 8 (small so that we can inspect it)\n",
    "\n",
    "# pyfftw.empty_aligned() is a helper function that works like numpy.empty() \n",
    "# but returns the array aligned to a particular number of bytes in memory, in this case 16. \n",
    "\n",
    "# If the alignment is not specified then the library inspects the CPU for an appropriate alignment value. \n",
    "# Having byte aligned arrays allows FFTW to performed vector operations, \n",
    "# potentially speeding up the FFT (a similar pyfftw.byte_align().\n",
    "\n",
    "a = pyfftw.empty_aligned(8, dtype='complex128', n=16)\n",
    "a[:] = np.random.randn(8) + 1j*np.random.randn(8)\n",
    "\n",
    "# Two ways to summon the in-built fft, via either pyfftw \n",
    "b = pyfftw.interfaces.numpy_fft.fft(a)\n",
    "# ... or numpy directly\n",
    "c = np.fft.fft(a)\n",
    "\n",
    "print(*a, sep = \"\\n\")\n",
    "print(\"===========================================\")\n",
    "print(*b, sep = \"\\n\")\n",
    "print(\"===========================================\")\n",
    "print(*c, sep = \"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "cell_style": "split"
   },
   "outputs": [],
   "source": [
    "# Direct implementation\n",
    "def mydft(x):\n",
    "\n",
    "    N = len(x)\n",
    "    n = np.arange(N)\n",
    "    k = n.reshape((N, 1))\n",
    "    e = np.exp(-2j * np.pi * k * n / N)\n",
    "    \n",
    "    X = np.dot(e, x)\n",
    "    \n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computational cost of DFT\n",
    "\n",
    "DFT of an array, $\\left\\{h_n\\right\\}$, of length, $N$, can be thought of as a matrix-vector multiplication\n",
    "$$\n",
    "H_n = \\sum_{k=0}^{N-1} W_{nk}\\, h_k,\n",
    "$$\n",
    "where $W_{pq}$ are the elements of the $N\\times N$ matrix\n",
    "$$\n",
    "W_{pq} = \\mathrm{e}^{\\frac{2\\pi i p q}{N}}.\n",
    "$$\n",
    "\n",
    "Direct computation is $\\mathcal{O}(N^2)$. We shall now show that it is possible to exploit regularities in the structure of the matrix $W$ to perform the calculation in $\\mathcal{O}(N\\log N)$ operations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Danielson-Lanczos lemma\n",
    "\n",
    "DFT of length $N$ is a sum of 2 DFTs of length $N/2$. \n",
    "\n",
    "The sub-transforms of length $N/2$ come from the odd/even elements of the original array:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "H_k = \\sum_{j=0}^{N-1} h_j\\,\\mathrm{e}^{\\frac{2\\pi i\\, k\\, j}{N}} &=  \\left[ \\sum_{j=0}^{\\frac{N}{2}-1} h_{2j}\\,\\mathrm{e}^{\\frac{2\\pi i\\, k\\, (2 j)}{N}} \\right] + \\left[ \\sum_{j=0}^{\\frac{N}{2}-1} h_{2j+1}\\,\\mathrm{e}^{\\frac{2\\pi i\\, k\\, (2 j+1)}{N}} \\right]\\\\\n",
    "&=\\left[ \\sum_{j=0}^{\\frac{N}{2}-1} h^{(e)}_{j}\\,\\mathrm{e}^{\\frac{2\\pi i\\, k\\, j}{N/2}} \\right] + \\mathrm{e}^{\\frac{2 \\pi i\\, k}{N}}\\,\\left[ \\sum_{j=0}^{\\frac{N}{2}-1} h^{(o)}_{j}\\,\\mathrm{e}^{\\frac{2\\pi i\\, k\\, j}{N/2}} \\right].\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Danielson-Lanczos lemma\n",
    "\n",
    "Here $\\left\\{ h^{(o/e)}_j, j=0,\\ldots \\frac{N}{2}-1\\right\\}$ are the arrays of length $\\frac{N}{2}$ constructed from the odd/even elements of the original array respectively. \n",
    "\n",
    "The last line is simply  a linear combination of the DFTs of these two smaller arrays (indices greater than $\\frac{N}{2}-1$ are interpreted via periodicity):\n",
    "\n",
    "$$\n",
    "H_k = H^{(e)}_k + w_k\\, H^{(o)}_k,\n",
    "$$\n",
    "\n",
    "where $w_k = \\mathrm{e}^{\\frac{2 \\pi i\\, k}{N}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fast Fourier Transform (FFT)\n",
    "\n",
    "The Danielson-Lanczos Lemma is the basis for a divide-and-conquer approach to calculating the DFT. \n",
    "\n",
    "Computational complexity satisfies\n",
    "\n",
    "$$\n",
    "F(N) = 2\\,F\\left(\\frac{N}{2} \\right) +2\\,N\n",
    "$$\n",
    "\n",
    "with $F(1)=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split"
   },
   "source": [
    "We know from before that \n",
    "\n",
    "$$\n",
    "F(N) = \\mathcal{O}(N\\,\\log N).\n",
    "$$\n",
    "\n",
    "The divide-and-conquer version of DFT is called the Fast Fourier Transform (FFT).\n",
    "\n",
    "Big difference: FT of 100 MP image ($N=10^8$)) at 1 ns per operation:\n",
    "\n",
    "* $N \\log N = 1.8\\times 10^{9}$ - takes $\\sim 2$  s.\n",
    "* $N^2 = 10^{16}$ - takes $\\sim 115$ days!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recursive implementation of FFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "cell_style": "center"
   },
   "outputs": [],
   "source": [
    "def myfft(x):\n",
    "    # Note that this is restricted to inputs of power of 2 size\n",
    "    N = len(x)\n",
    "    \n",
    "    if N == 1:\n",
    "        return x\n",
    "    else:\n",
    "        X_even = myfft(x[::2])\n",
    "        X_odd = myfft(x[1::2])\n",
    "        factor = \\\n",
    "          np.exp(-2j*np.pi*np.arange(N)/ N)\n",
    "        \n",
    "        X = np.concatenate(\\\n",
    "            [X_even+factor[:int(N/2)]*X_odd,\n",
    "             X_even+factor[int(N/2):]*X_odd])\n",
    "        return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.6881789171900392+0.5916479882309668j)\n",
      "(0.9047019733848708-0.4387366153747087j)\n",
      "(0.8277040544279829+0.2793681357976642j)\n",
      "(0.27539319505285476-0.002330991182591074j)\n",
      "(1.1955843835618385+1.6491945185167878j)\n",
      "(0.5326351908397333-1.9874427848958942j)\n",
      "(-0.7180623891102259+1.1617246580474099j)\n",
      "(0.13139467893387619+0.6283549973711337j)\n",
      "===========================================\n",
      "(2.461172169900891+1.8817799065107685j)\n",
      "(-1.9557129100100332-1.4271632722049343j)\n",
      "(-2.654439605405103-0.23079957733519252j)\n",
      "(-0.5135371010373359-1.4177563328076226j)\n",
      "(-1.227077906521779+5.482090694674889j)\n",
      "(-3.576526735993213-3.779462675443126j)\n",
      "(3.449967207513188+1.8302990031405544j)\n",
      "(-1.4892764559669287+2.394196159312398j)\n",
      "===========================================\n",
      "(2.461172169900891+1.8817799065107685j)\n",
      "(-1.9557129100100332-1.4271632722049343j)\n",
      "(-2.654439605405103-0.23079957733519252j)\n",
      "(-0.5135371010373359-1.4177563328076226j)\n",
      "(-1.227077906521779+5.482090694674889j)\n",
      "(-3.576526735993213-3.779462675443126j)\n",
      "(3.449967207513188+1.8302990031405544j)\n",
      "(-1.4892764559669287+2.394196159312398j)\n"
     ]
    }
   ],
   "source": [
    "print(*a, sep = \"\\n\")\n",
    "print(\"===========================================\")\n",
    "testDFT = mydft(a)\n",
    "testFFT = myfft(a)\n",
    "\n",
    "print(*testFFT, sep = \"\\n\")\n",
    "print(\"===========================================\")\n",
    "print(*testFFT, sep = \"\\n\")"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}