fourier_no_ei.ipynb

{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  },
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Fourier without the ei\n",
      "\n",
      "## Introduction\n",
      "\n",
      "The standard equations for discrete Fourier transforms (DFTs) involve\n",
      "exponentials to the power of $i$ - the imaginary unit.  I personally find\n",
      "these difficult to think about, and it turns out, the DFT is fairly easy to\n",
      "recast in terms of $\\sin$ and $\\cos$.  This page goes through this process,\n",
      "and tries to show how thinking in this way can explain some features of the\n",
      "DFT.\n",
      "\n",
      "### How hard is the mathematics?\n",
      "\n",
      "You will not need heavy math to follow this page.  If you don't remember the\n",
      "following concepts you might want to brush up on them.  There are also links to\n",
      "proofs and explanations for these ideas in the page as we go along:\n",
      "\n",
      "* basic trigonometry (SOH CAH TOA, Pythagoras' theorem);\n",
      "* the angle sum rule;\n",
      "* basic algebra with sums (e.g. $\\sum_{i=0}^{N-1} (x_i + 4) = \\sum_{i=0}^{N-1}\n",
      "x_i + N 4$);\n",
      "* vector dot products;\n",
      "* vector projection using the dot product.\n",
      "\n",
      "You will not need to understand complex numbers in any depth.\n",
      "\n",
      "### Loading and configuring code libraries\n",
      "\n",
      "Load and configure the Python libraries we will use in this notebook:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "1"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "2"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "3"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - tell numpy to print numbers to 4 decimal places only\n",
      "np.set_printoptions(precision=4, suppress=True)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "4"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Some actual numbers to start\n",
      "\n",
      "Let us start with a DFT of some data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# An example input vector\n",
      "x = np.array(\n",
      "    [ 0.4967, -0.1383,  0.6477,  1.523 , -0.2342, -0.2341,  1.5792,\n",
      "      0.7674, -0.4695,  0.5426, -0.4634, -0.4657,  0.242 , -1.9133,\n",
      "     -1.7249, -0.5623, -1.0128,  0.3142, -0.908 , -1.4123,  1.4656,\n",
      "     -0.2258,  0.0675, -1.4247, -0.5444,  0.1109, -1.151 ,  0.3757,\n",
      "     -0.6006, -0.2917, -0.6017,  1.8523])\n",
      "N = 32  # the length of the time-series\n",
      "plt.plot(x)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "5"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Now the DFT\n",
      "X = np.fft.fft(x)\n",
      "X"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "6"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that `X` - the output of the forward DFT - is a vector of complex\n",
      "numbers. We will go into this in detail later.\n",
      "\n",
      "When we do the inverse DFT on `X` we return the original values of our input\n",
      "`x`, but as complex numbers with imaginary part 0:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Apply the inverse DFT to the output of the forward DFT \n",
      "x_back = np.fft.ifft(X)\n",
      "x_back"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "7"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Rewriting the DFT without the $e^i$\n",
      "\n",
      "### Some notation\n",
      "\n",
      "We will use $\\triangleq$ to mean \"is defined as\".  For example, read:\n",
      "\n",
      "$$\n",
      "i \\triangleq \\sqrt{-1}\n",
      "$$\n",
      "\n",
      "as \"$i$ is defined as the square root of -1.\"\n",
      "\n",
      "$\\in$ means \"in\" or \"one of\".  For example:\n",
      "\n",
      "$$\n",
      "k \\in 0, 2, 4\n",
      "$$\n",
      "\n",
      "means that the value $k$ can take any of the three values 0, 2 or 4.\n",
      "\n",
      "$\\notin$ means \"not in\" or \"not one of\".  For example:\n",
      "\n",
      "$$\n",
      "k \\notin 1, 3, 5\n",
      "$$\n",
      "\n",
      "means that the value $k$ can *not* take any of the three values 1, 3 or 5, but,\n",
      "all other things being equal, can take any other value.\n",
      "\n",
      "We write a *vector* of numbers $[v_0, v_1, ... v_{N}]$ as $\\vec{v}$.\n",
      "\n",
      "### DFT and FFT\n",
      "\n",
      "The fast fourier transform (FFT) refers to a particular set of - er - fast\n",
      "algorithms for calculating the DFT.  It is common, but confusing, to use \"FFT\"\n",
      "to mean DFT.\n",
      "\n",
      "### Introducing the discrete Fourier transform\n",
      "\n",
      "Let us say we have a vector of $N$ values in time, or space $\\vec{x} = [x_0, x_1\n",
      "... x_{N-1}]$.  We generally index $\\vec{x}$ with subscript $n$, so the sample\n",
      "at index $n$ is $x_n$.\n",
      "\n",
      "The DFT converts $\\vec{x}$ from a vector in time, or space, to a vector\n",
      "$\\vec{X}$ representing temporal or spatial frequency components.\n",
      "\n",
      "We will call our original $\\vec{x}$ the *signal*, meaning, the signal not\n",
      "transformed to frequency.\n",
      "\n",
      "The DFT converts $\\vec{x}$ to $\\vec{X}$ where $\\vec{X} = [X_0, X_1, ...\n",
      "X_{N-1}]$.  We generally index $\\vec{X}$ with subscript $k$, so the sample at\n",
      "index $k$ is $X_k$.\n",
      "\n",
      "Here is the equation for the discrete Fourier transform:\n",
      "\n",
      "$$\n",
      "X_k = \\sum_{n=0}^{N-1} x_n \\; e^{-i 2 \\pi \\frac{k}{N} n}\n",
      "$$\n",
      "\n",
      "This is the transform from signal to frequency.  We will call this the\n",
      "*forward* Fourier transform.\n",
      "\n",
      "Here is the equation for the inverse Fourier transform:\n",
      "\n",
      "$$\n",
      "x_n = \\frac{1}{N} \\sum_{k=0}^{N-1} X_k \\; e^{i 2 \\pi \\frac{k}{N} n}\n",
      "$$\n",
      "\n",
      "The inverse Fourier transform converts from frequency back to signal.\n",
      "\n",
      "### Scrubbing the $e^i$\n",
      "\n",
      "The forward and inverse equations are very similar; both share a term $e^{iz}$,\n",
      "where $z = -2 \\pi \\frac{k}{N} n$ for the forward transform; $z = 2 \\pi\n",
      "\\frac{k}{N} n$ for the inverse transform.\n",
      "\n",
      "Some people are used to looking at the form $e^{iz}$ and thinking \"aha, that's a\n",
      "rotation around a circle\".  Apparently this is an intuition that builds up over\n",
      "time working with these sorts of equations.\n",
      "\n",
      "Unfortunately, some of us find it hard to think in complex exponentials, or in\n",
      "terms of complex numbers.\n",
      "\n",
      "So, in this tutorial, we will express the Fourier transform in terms of $\\sin$\n",
      "and $\\cos$.  We will be using complex numbers, but almost entirely as a pair of\n",
      "numbers to represent two components of the same thing, rather than a single\n",
      "number with a real and imaginary part.\n",
      "\n",
      "Having said that, we will need some very basic properties of complex and\n",
      "imaginary numbers.\n",
      "\n",
      "### Lightning refresher on complex numbers\n",
      "\n",
      "$\\newcommand{R}[1]{\\operatorname{Re}\\,(#1)}\\newcommand{I}[1]{\\operatorname{Im}\\,\n",
      "(#1)}$\n",
      "First we remember that $i$ is the *imaginary unit*.  $i$ is the imaginary number\n",
      "such that:\n",
      "\n",
      "$$\n",
      "i \\triangleq \\sqrt{-1}\n",
      "$$\n",
      "\n",
      "Therefore:\n",
      "\n",
      "$$\n",
      "i^2 = -1\n",
      "$$\n",
      "\n",
      "Engineers often write the imaginary unit as $j$, and that is the convention that\n",
      "Python uses:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "1j * 1j"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "8"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An imaginary number is a real number multiplied by the imaginary unit.  For\n",
      "example $3i$ is an imaginary number.\n",
      "\n",
      "A *complex* number is a number that has a *real* part and an *imaginary* part.\n",
      "The real part is a real number.  The imaginary part is an imaginary number.\n",
      "\n",
      "For example, consider the complex number $a = (4 + 3i)$.  $a$ has two parts.\n",
      "The first is a real number, often written as $\\R{a}$, and the second is an\n",
      "imaginary number $\\I{a}$:\n",
      "\n",
      "$$\n",
      "\\R{a} = 4 \\\\\n",
      "\\I{a} = 3i\n",
      "$$\n",
      "\n",
      "Now consider the complex number $c = (p + qi)$.  $\\R{c} = p, \\I{c} = qi$.\n",
      "\n",
      "To multiply a complex number $c$ by a real number $r$, we multiply both real and\n",
      "imaginary parts by the real number:\n",
      "\n",
      "$$\n",
      "r (p + qi) = (rp + rqi)\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = (4 + 3j)\n",
      "2 * a"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "9"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Multiplying a complex number by an imaginary number follows the same logic, but\n",
      "remembering that $i^2 = -1$. Let us say $s i$ is our imaginary number:\n",
      "\n",
      "$$\n",
      "si (p + qi) = (s i p + s i q i) \\\\\n",
      "= (s p i + s q i^2) \\\\\n",
      "= (s p i - s q) \\\\\n",
      "= (-s q + s p i)\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "3j * a"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "10"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Back to scrubbing the $e^i$\n",
      "\n",
      "Our first tool in this enterprise is Euler's formula:\n",
      "\n",
      "$$\n",
      "e^{ix} = \\cos x + i\\sin x\n",
      "$$\n",
      "\n",
      "This is the basis for thinking of $e^{ix}$ as being rotation around a circle, of\n",
      "which you will hear no more in this page.  In our case, it allows us to rewrite\n",
      "the forward and inverse Fourier transforms:\n",
      "\n",
      "First let's define a new value $f$, that depends on $N$ - the number of\n",
      "observations in our vector $\\vec{x}$.\n",
      "\n",
      "$$\n",
      "D \\triangleq \\frac{2 \\pi}{N}\n",
      "$$\n",
      "\n",
      "With that value:\n",
      "\n",
      "$$\n",
      "X_k = \\sum_{n=0}^{N-1} x_n \\cdot \\cos(-k n D) + i \\sum_{n=0}^{N-1} x_n \\cdot\n",
      "\\sin(-k n D) \\\\\n",
      "x_n = \\frac{1}{N} \\sum_{k=0}^{N-1} X_k \\cdot \\cos(k n D) + i \\frac{1}{N}\n",
      "\\sum_{k=0}^{N-1} X_k \\cdot \\sin(k n D)\n",
      "$$\n",
      "\n",
      "We can simplify this a bit further, because, for any angle $\\alpha$:\n",
      "\n",
      "$$\n",
      "\\cos(-\\alpha) = cos(\\alpha) \\\\\n",
      "\\sin(-\\alpha) = -sin(\\alpha)\n",
      "$$\n",
      "\n",
      "$$\n",
      "X_k = \\sum_{n=0}^{N-1} x_n \\cdot \\cos(k n D) - i \\sum_{n=0}^{N-1} x_n \\cdot\n",
      "\\sin(k n D)\n",
      "\\\\\n",
      "x_n = \\frac{1}{N} \\sum_{k=0}^{N-1} X_k \\cdot \\cos(k n D) + i \\frac{1}{N}\n",
      "\\sum_{k=0}^{N-1} X_k \\cdot \\sin(k n D)\n",
      "$$\n",
      "\n",
      "### Rewriting as dot products\n",
      "\n",
      "We can simplify the notation, and maybe make the process clearer, by rewriting\n",
      "these sums in terms of dot products.\n",
      "\n",
      "As y'all remember, the dot product of two length $N$ vectors $\\vec{v}, \\vec{w}$\n",
      "is given by:\n",
      "\n",
      "$$\n",
      "\\vec{v} \\cdot \\vec{w} \\triangleq \\sum_{i=0}^{N-1} v_i w_i\n",
      "$$\n",
      "\n",
      "Clearly, because $v_i w_i = w_i v_i$:\n",
      "\n",
      "$$\n",
      "\\vec{v} \\cdot \\vec{w} = \\vec{w} \\cdot \\vec{v}\n",
      "$$\n",
      "\n",
      "For the moment, let us concentrate on the forward transform.\n",
      "\n",
      "$$\n",
      "\\vec{t_k} \\triangleq \\left[ k 2 \\pi \\frac{n}{N} \\;\\mathrm{for}\\; n \\in\n",
      "0,1,\\ldots,N-1 \\right] \\\\\n",
      "\\vec{c_k} \\triangleq \\left[ \\cos(t_{k,n}) \\;\\mathrm{for}\\; n \\in 0,1,\\ldots,N-1\n",
      "\\right] \\\\\n",
      "\\vec{s_k} \\triangleq \\left[ \\sin(t_{k,n}) \\;\\mathrm{for}\\; n \\in 0,1,\\ldots,N-1\n",
      "\\right]\n",
      "$$\n",
      "\n",
      "Now we can rewrite the sums in the forward transform as the sum of two dot\n",
      "products:\n",
      "\n",
      "$$\n",
      "X_k = \\vec{x} \\cdot \\vec{c_k} - i \\vec{x} \\cdot \\vec{s_k}\n",
      "$$\n",
      "\n",
      "The vector $\\vec{t_k}$ is key to understanding what is going on.  $t_k$ sets up\n",
      "the horizontal axis values to sample a $\\sin$ or $\\cos$ function so the function\n",
      "gives us $k$ cycles over the indices $0 .. N-1$.\n",
      "\n",
      "In the formulae above, $n / N$ is the proportion of the whole signal width $N$,\n",
      "so it varies between 0 and $(N-1) / N$.  The $2 \\pi$ corresponds to one cycle of\n",
      "the cosine or sine function.\n",
      "\n",
      "So, $\\vec{t_0}$ gives a vector of zeros corresponding to $k=0$ cycles across $0\n",
      "... N-1$.  $\\vec{t_1}$ gives us $0$ up to (not including) $2 \\pi$ - one cycle\n",
      "across the indices $0 .. N-1$.  $\\vec{t_2}$ gives us $0$ up to (not including)\n",
      "$4 \\pi$ - two cycles.\n",
      "\n",
      "Here are some plots of $\\vec{c_k}$, $\\vec{s_k}$ for $k \\in 0, 1, 2, 3, 4, 5$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(6, 1, figsize=(8, 5))\n",
      "ns = np.arange(N)\n",
      "one_cycle = 2 * np.pi * ns / N\n",
      "for k in range(6):\n",
      "    t_k = k * one_cycle\n",
      "    axes[k].plot(ns, np.cos(t_k), label='cos')\n",
      "    axes[k].plot(ns, np.sin(t_k), label='sin')\n",
      "    axes[k].set_xlim(0, N-1)\n",
      "    axes[k].set_ylim(-1.1, 1.1)\n",
      "axes[0].legend()\n",
      "plt.tight_layout()"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "11"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To rephrase: $\\vec{c_k}, \\vec{s_k}$ are cosine / sine waves with $k$ cycles over\n",
      "the $N$ indices.\n",
      "\n",
      "So, the $X_k$ value is the dot product of the $\\vec{x}$ with a cosine wave of\n",
      "$k$ cycles minus $i$ times the dot product of $\\vec{x}$ with the sine wave of\n",
      "$k$ cycles.\n",
      "\n",
      "While this is all fresh in our minds, let us fill out the equivalent notation\n",
      "for the inverse transform.\n",
      "\n",
      "$$\n",
      "\\vec{t_n} \\triangleq \\left[ n 2 \\pi \\frac{k}{N} \\;\\mathrm{for}\\; k \\in\n",
      "=0,1,\\ldots,N-1 \\right] \\\\\n",
      "\\vec{c_n} \\triangleq \\left[ \\cos(t_{n,k}) \\;\\mathrm{for}\\; k \\in 0,1,\\ldots,N-1\n",
      "\\right] \\\\\n",
      "\\vec{s_n} \\triangleq \\left[ \\sin(t_{n,k}) \\;\\mathrm{for}\\; k \\in 0,1,\\ldots,N-1\n",
      "\\right]\n",
      "$$\n",
      "\n",
      "Because both $n$ and $k$ have indices from $0 .. N-1$:\n",
      "$$\n",
      "\\vec{t_k} = \\vec{t_n} \\;\\mathrm{where}\\; k=n \\\\\n",
      "\\vec{c_k} = \\vec{c_n} \\;\\mathrm{where}\\; k=n \\\\\n",
      "\\vec{s_k} = \\vec{s_n} \\;\\mathrm{where}\\; k=n\n",
      "$$\n",
      "\n",
      "We will return to this point fairly soon.\n",
      "\n",
      "The inverse transform is now:\n",
      "\n",
      "$$\n",
      "x_n = \\frac{1}{N} \\vec{X} \\cdot \\vec{c_n} + i \\frac{1}{N} \\vec{X} \\cdot\n",
      "\\vec{s_n}\n",
      "$$\n",
      "\n",
      "### Rewriting the DFT with cosine and sine basis matrices\n",
      "\n",
      "Instead of writing the formulae for the individual elements $X_k$ and $x_n$, we\n",
      "can use matrices to express our formulae in terms of the vectors $\\vec{X},\n",
      "\\vec{x}$.\n",
      "\n",
      "$\\newcommand{C}{\\mathbf{C}} \\newcommand{S}{\\mathbf{S}}$\n",
      "Define a matrix $\\C$ that has rows $[\\vec{c_0}, \\vec{c_1}, ..., \\vec{c_{N-1}}]$:\n",
      "\n",
      "$$\n",
      "\\C \\triangleq\n",
      " \\begin{bmatrix}\n",
      "    c_{0,0}, c_{0, 1}, ..., c_{0, N-1} \\\\\n",
      "    c_{1,0}, c_{1, 1}, ..., c_{1, N-1} \\\\\n",
      "    ... \\\\\n",
      "    c_{N-1,0}, c_{N-1, 1}, ..., c_{N-1, N-1} \\\\\n",
      " \\end{bmatrix}\n",
      "$$\n",
      "\n",
      "Call $\\C$ the *cosine basis matrix*.\n",
      "\n",
      "Define a matrix $\\S$ that has rows $[\\vec{s_0}, \\vec{s_1}, ..., \\vec{s_{N-1}}]$:\n",
      "\n",
      "$$\n",
      "\\S \\triangleq\n",
      " \\begin{bmatrix}\n",
      "    s_{0,0}, s_{0, 1}, ..., s_{0, N-1} \\\\\n",
      "    s_{1,0}, s_{1, 1}, ..., s_{1, N-1} \\\\\n",
      "    ... \\\\\n",
      "    s_{N-1,0}, s_{N-1, 1}, ..., s_{N-1, N-1} \\\\\n",
      " \\end{bmatrix}\n",
      "$$\n",
      "\n",
      "Call $\\S$ the *sine basis matrix*.\n",
      "\n",
      "Now we can rewrite the forward and inverse DFT as matrix products:\n",
      "\n",
      "$$\n",
      "\\vec{X} = \\C \\cdot \\vec{x} - i \\S \\cdot \\vec{x} \\\\\n",
      "\\vec{x} = \\frac{1}{N} \\C \\cdot \\vec{X} + i \\frac{1}{N} \\S \\cdot \\vec{X}\n",
      "$$\n",
      "\n",
      "This gives us the same calculation for $X_k$ and $x_n$ as we have above using\n",
      "the vector dot products.  Write row $k$ of $\\C$ as $C_{k,:}$. Row $k$ of $\\S$ is\n",
      "$S_{k,:}$.  Thus, from the rules of matrix multiplication:\n",
      "\n",
      "$$\n",
      "X_k = C_{k,:} \\cdot \\vec{x} - i S_{k,:} \\cdot \\vec{x} \\\\\n",
      "= \\vec{x} \\cdot \\vec{c_k} - i \\vec{x} \\cdot \\vec{s_k}\n",
      "$$\n",
      "\n",
      "and (inverse transform):\n",
      "\n",
      "$$\n",
      "x_n = \\frac{1}{N} C_{n,:} \\cdot \\vec{X} + i \\frac{1}{N} S_{n,:} \\cdot \\vec{X} \\\\\n",
      "= \\frac{1}{N} \\vec{X} \\cdot \\vec{c_n} + i \\frac{1}{N} \\vec{X} \\cdot \\vec{s_n}\n",
      "$$\n",
      "\n",
      "We can build $\\C$ and $\\S$ for our case with $N=32$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.zeros((N, N))\n",
      "S = np.zeros((N, N))\n",
      "ns = np.arange(N)\n",
      "one_cycle = 2 * np.pi * ns / N\n",
      "for k in range(N):\n",
      "    t_k = k * one_cycle\n",
      "    C[k, :] = np.cos(t_k)\n",
      "    S[k, :] = np.sin(t_k)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "12"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We get the same result using this matrix formula, as we do using the canned DFT:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Recalculate the forward transform with C and S\n",
      "X_again = C.dot(x) - 1j * S.dot(x)\n",
      "assert np.allclose(X, X_again)  # same result as for np.fft.fft\n",
      "# Recalculate the inverse transform\n",
      "x_again = 1. / N * C.dot(X) + 1j / N * S.dot(X)\n",
      "assert np.allclose(x, x_again)  # as for np.fft.ifft, we get x back"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "13"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Displaying the DFT transform\n",
      "\n",
      "We can show the matrix calculation of the DFT as images.   To do this we will\n",
      "use some specialized code.  If you are running this notebook yourself, download\n",
      "[dft_plots.py](dft_plots.py) to the directory containing the notebook."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Import the custom DFT plotting code for this notebook\n",
      "import dft_plots as dftp"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "14"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here we show the forward DFT given by the formula:\n",
      "\n",
      "$$\n",
      "\\vec{X} = \\C \\cdot \\vec{x} - i \\S \\cdot \\vec{x}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Show image sketch for forward DFT\n",
      "sketch = dftp.DFTSketch(x)\n",
      "sketch.sketch(figsize=(12, 5))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "15"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The plot shows each matrix and vector as grayscale, where mid gray corresponds\n",
      "to 0, black corresponds to the most negative value and white to the most\n",
      "positive value.  For example the first four values in the vector $\\vec{x}$ are:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[:4]"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "16"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can see $\\vec{x}$ shown at the right of the graphic as a column vector.  The\n",
      "grayscale of the top four values in the graphic are light gray, mid gray, light\n",
      "gray, and near white, corresponding to the values above.\n",
      "\n",
      "$\\vec{X}$ is a vector of complex numbers.\n",
      "\n",
      "On the left of the equals sign you see the complex vector $\\vec{X}$ displayed as\n",
      "two columns.\n",
      "\n",
      "Define $\\R{\\vec{X}}$ to be the vector containing the real parts of the complex\n",
      "values in $\\vec{X}$.  Define $\\I{\\vec{X}}$ to be the vector containing the\n",
      "imaginary parts of $\\vec{X}$:\n",
      "\n",
      "$$\n",
      "\\R{\\vec{X}} \\triangleq [\\R{X_0}, \\R{X_1}, ..., \\R{X_{N-1}}] \\\\\n",
      "\\I{\\vec{X}} \\triangleq [\\I{X_0}, \\I{X_1}, ..., \\I{X_{N-1}}]\n",
      "$$\n",
      "\n",
      "The left hand column in the graphic shows $\\R{\\vec{X}}$, and the column to the\n",
      "right of that shows $\\I{\\vec{X}}$.\n",
      "\n",
      "To the right of the equals sign we see the representation of $\\C \\cdot \\vec{x}$\n",
      "and $\\S \\cdot \\vec{x}$, with $\\vec{x}$ displayed as a column vector.\n",
      "\n",
      "$\\C$ and $\\S$ have some interesting patterns which we will explore in the next\n",
      "section.\n",
      "\n",
      "We can show the inverse DFT in the same way:\n",
      "\n",
      "$$\n",
      "\\vec{x} = \\frac{1}{N} \\C \\cdot \\vec{X} + i \\frac{1}{N} \\S \\cdot \\vec{X}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sketch.sketch(inverse=True, figsize=(12, 5))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "17"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The output from the inverse transform is a complex vector, but in our case,\n",
      "where the input to the DFT was a vector of real numbers, the imaginary parts are\n",
      "all zero, and the real part is equal to our input to the forward DFT :\n",
      "$\\vec{x}$.  We will see why the imaginary parts are all zero in the following\n",
      "sections.\n",
      "\n",
      "### Real and complex input to the DFT\n",
      "\n",
      "This page is mostly concerned with the common case where the input to the\n",
      "forward DFT is a vector of real numbers.  The mathematics also works for the\n",
      "case where the input to the forward DFT is a vector of complex numbers:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "complex_x = np.array(  # A Random array of complex numbers\n",
      "      [ 0.61-0.83j, -0.82-0.12j, -0.50+1.14j,  2.37+1.67j,  1.62+0.69j,\n",
      "        1.61-0.06j,  0.54-0.73j,  0.89-1.j  ,  0.17-0.71j,  0.75-0.01j,\n",
      "       -1.06-0.14j, -2.53-0.33j,  1.74+0.83j,  1.34-0.64j,  1.47+0.71j,\n",
      "        0.82+0.4j , -1.59-0.58j,  0.13-1.02j,  0.47-0.73j,  1.45+1.31j,\n",
      "        1.32-0.28j,  1.58-2.13j,  0.75-0.43j,  1.24+0.4j ,  0.02+1.08j,\n",
      "        0.07-0.57j, -1.21+1.08j,  1.38+0.54j, -1.35+0.3j , -0.61+1.08j,\n",
      "       -0.96+1.81j, -1.95+1.64j])\n",
      "complex_X = np.fft.fft(complex_x)  # Canned DFT\n",
      "complex_X_again = C.dot(complex_x) - 1j * S.dot(complex_x)  # Our DFT\n",
      "# We get the same result as the canned DFT\n",
      "assert np.allclose(complex_X, complex_X_again)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "18"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The sketch of the complex forward DFT looks like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sketch = dftp.DFTSketch(complex_x)\n",
      "sketch.sketch(figsize=(12, 5))\n",
      "sketch.title('Forward DFT for complex input vector')"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "19"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The input $\\vec{x}$ vectors following $\\C$ and $\\S$ are now complex, with a real\n",
      "and a complex column for the real and complex vectors in $\\vec{x}$.\n",
      "\n",
      "For what follows, unless we say otherwise, we will always be talking about real\n",
      "number input to the DFT.\n",
      "\n",
      "## Some properties of the cosine and sine basis matrices\n",
      "\n",
      "First we note that $\\C$ and $\\S$ are always real matrices, regardless of the\n",
      "input $\\vec{x}$ or $\\vec{X}$.\n",
      "\n",
      "Let's show $\\C$ and $\\S$ as grayscale images again:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
      "dftp.show_array(axes[0], dftp.scale_array(C))\n",
      "axes[0].set_title(\"$\\mathbf{C}$\")\n",
      "dftp.show_array(axes[1], dftp.scale_array(S))\n",
      "axes[1].set_title(\"$\\mathbf{S}$\")"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "20"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Mirror symmetry\n",
      "\n",
      "From the images we see that the bottom half of $\\C$ looks like a mirror image\n",
      "of the top half of $\\C$.  The bottom half of $\\S$ looks like a sign flipped\n",
      "(black $\\Leftrightarrow$ white) mirror image of the top half of $\\S$.  In fact\n",
      "this is correct:\n",
      "\n",
      "$$\n",
      "C_{p,:} = C_{N-p,:} \\; \\mathrm{for} \\; p > 0 \\\\\n",
      "S_{p,:} = -S_{N-p,:} \\; \\mathrm{for} \\; p > 0\n",
      "$$\n",
      "\n",
      "Why is this?  Let's look at lines from the center of $\\C$.  Here we are plotting\n",
      "the continuous cosine function with dotted lines, with filled circles to\n",
      "represent the discrete samples we took to fill the row of $\\C$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "center_rows = [N / 2. - 1, N / 2., N / 2. + 1]\n",
      "fig = dftp.plot_cs_rows('C', N, center_rows)\n",
      "fig.suptitle('Rows $N / 2 - 1$ through $N / 2 + 1$ of $\\mathbf{C}$',\n",
      "             fontsize=20)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "21"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first plot in this grid is for row $k = N / 2 - 1$.  This row starts\n",
      "sampling just before the peak and trough of the cosine.\n",
      "In the center is row $k = N / 2$ of $\\C$.  This is sampling the cosine wave\n",
      "exactly at the peak and trough.  When we get to next row, at $k = N / 2 + 1$, we\n",
      "start sampling after the peak and trough of the cosine, and these samples are\n",
      "identical to the samples just before the peak and trough, at row $k = N / 2 -\n",
      "1$.  Row $k = N / 2$ is sampling at the Nyquist sampling frequency, and row $k =\n",
      "N / 2 + 1$ is sampling at a frequency lower than Nyquist and therefore it is\n",
      "being *aliased* to the same apparent frequency as row $k = N / 2 - 1$.\n",
      "\n",
      "This might be more obvious plotting rows 1 and N-1 of $\\C$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = dftp.plot_cs_rows('C', N, [1, N-1])\n",
      "fig.suptitle('Rows $1$ and $N - 1$ of $\\mathbf{C}$',\n",
      "             fontsize=20)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "22"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Of course we get the same kind of effect for $\\S$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = dftp.plot_cs_rows('S', N, center_rows)\n",
      "fig.suptitle('Rows $N / 2 - 1$ through $N / 2 + 1$ of $\\mathbf{S}$',\n",
      "             fontsize=20)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "23"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = dftp.plot_cs_rows('S', N, [1, N-1])\n",
      "fig.suptitle('Rows $1$ and $N - 1$ of $\\mathbf{S}$',\n",
      "             fontsize=20)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "24"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that for $\\S$, the sine waves after $k = N / 2$ are sign-flipped relative\n",
      "to their matching rows before $k = N / 2$.  Thus row $k = N / 2 + 1$ will be\n",
      "aliased to the same frequency as for row $k = N / 2 - 1$, but with a negative\n",
      "sign.\n",
      "\n",
      "It is this sign-flip that leads to the concept of *negative frequency* in the\n",
      "DFT, and to the property of *conjugate symmetry* from the DFT on a vector of\n",
      "real numbers.  We will hear more about these later.\n",
      "\n",
      "### Matrix symmetry\n",
      "\n",
      "The next thing we notice about $\\C$ and $\\S$ is that they are transpose\n",
      "*symmetric* matrices:\n",
      "\n",
      "$$\n",
      "\\C = \\C^T \\\\\n",
      "\\S = \\S^T \\\\\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "assert np.allclose(C, C.T)\n",
      "assert np.allclose(S, S.T)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "25"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Why is this?  Consider the first *column* of $\\C$.  This is given by $\\cos(k 2\n",
      "\\pi 0 / N) = \\cos(0)$, and thus, like the first *row* of $\\C$, is always = 1.\n",
      "\n",
      "Now consider the second row of $\\C$.  This is a cosine sampled at horizontal\n",
      "axis values:\n",
      "\n",
      "$$\n",
      "\\vec{t_1} \\triangleq \\left[ 2 \\pi \\frac{n}{N} \\;\\mathrm{for}\\; n \\in\n",
      "0,1,\\ldots,N-1 \\right]\n",
      "$$\n",
      "\n",
      "Call $t_{k, n}$ the value of $\\vec{t_k}$ at index $n$.  Now consider the second\n",
      "*column* of $\\C$.  This is a cosine sampled at horizontal axis values for $n\n",
      "= 1$:\n",
      "\n",
      "$$\n",
      "t_{0,1} = (0) 2 \\pi \\frac{1}{N} \\\\\n",
      "t_{1,1} = (1) 2 \\pi \\frac{1}{N} \\\\\n",
      "... \\\\\n",
      "t_{N-1,1} = (N-1) 2 \\pi \\frac{1}{N} \\\\\n",
      "$$\n",
      "\n",
      "In general, because the sequence $k \\in 0,1,\\ldots,N-1 $ is equal to the\n",
      "sequence $n \\in 0,1,\\ldots,N-1$, this means that the column sampling positions\n",
      "for row $n \\in t_{0, n}, t_{1, n}, ... , t_{N-1, n}$ are equal to the row\n",
      "sampling positions for corresponding ($k = n$) row $k \\in t_{k, 0}, t_{k, 1},\n",
      "... , t_{k, N-1}$.  Write column $z$ of $\\C$ as $C_{:,z}$; column $z$ of $\\S$ is\n",
      "$S_{:, z}$.  Therefore $C_{z, :} = C_{:, z}, S_{z, :} = S_{:, z}$.\n",
      "\n",
      "### Row dot products and lengths\n",
      "\n",
      "It is useful to look at the dot products of the rows of $\\C$ and $\\S$.  The dot\n",
      "product of each row with itself gives the squared *length* of the vector in that\n",
      "row.\n",
      "\n",
      "The vector length of a vector $\\vec{v}$ with $N$ elements is written as $\\|\n",
      "\\vec{v} \\|$, and defined as:\n",
      "\n",
      "$$\n",
      "\\| \\vec{v} \\| \\triangleq \\sqrt{\\sum_{n=0}^{N-1} v_n^2}\n",
      "= \\sqrt{ \\vec{v} \\cdot \\vec{v} }\n",
      "$$\n",
      "\n",
      "The dot products of different rows of $\\C$ and $\\S$ give an index of the\n",
      "strength of the relationship between the rows.  We can look at the dot products\n",
      "of all the rows of $\\C$ with all other rows with the matrix multiplication\n",
      "$\\C^T \\C$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dftp.show_array(plt.gca(), dftp.scale_array(C.T.dot(C)))\n",
      "plt.title(\"$\\mathbf{C^TC}$\")"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "26"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The image shows us that the dot product between the rows of $\\C$ is\n",
      "0 everywhere except:\n",
      "\n",
      "* the dot products of the rows with themselves (the squared vector lengths);\n",
      "* the dot products of the mirror image vectors such as $\\vec{c_1}$ and\n",
      "  $\\vec{c_{N-1}}$.  Because $\\vec{c_n} = \\vec{c_{N-n}}$, these dot products are\n",
      "  the same as the $\\| \\vec{c_n} \\|^2$.\n",
      "\n",
      "The squared row lengths are:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.diag(C.T.dot(C))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "27"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that the rows $\\vec{c_0}$ and $\\vec{c_{N / 2}}$ have squared length $N$,\n",
      "and the other rows have squared length $N / 2$.\n",
      "\n",
      "We can do the same for $\\S$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dftp.show_array(plt.gca(), dftp.scale_array(S.T.dot(S)))\n",
      "plt.title(\"$\\mathbf{S^TS}$\")"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "28"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remember that $\\vec{s_0}$ and $\\vec{s_{n/2}}$ are all 0 vectors.  The dot\n",
      "product of these rows with any other row, including themselves, is 0.  All other\n",
      "entries in this $\\S^T \\S$ matrix are zero except:\n",
      "\n",
      "* the dot products of rows with themselves (other than $\\vec{s_0}$ and\n",
      "$\\vec{s_{n/2}}$);\n",
      "* the dot products of the flipped mirror image vectors such as $\\vec{s_1}$ and\n",
      "  $\\vec{s_{N-1}}$.  Because $\\vec{s_n} = -\\vec{s_{N-n}}$, these dot products\n",
      "  are the same as $-\\| \\vec{s_n} \\|^2$.\n",
      "\n",
      "The squared row lengths are:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.diag(S.T.dot(S))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "29"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The rows $\\vec{s_0}$ and $\\vec{s_{N / 2}}$ have squared length $0$,\n",
      "and the other rows have squared length $N / 2$.\n",
      "\n",
      "Finally, let's look at the relationship between the rows of $\\C$ and the rows of\n",
      "$\\S$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.allclose(C.T.dot(S), 0)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "30"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The rows of $\\C$ and $\\S$ are completely orthogonal.\n",
      "\n",
      "In fact these relationships hold for $\\C$ and $\\S$ for any $N$.\n",
      "\n",
      "#### Proof for $\\C, \\S$ dot products\n",
      "\n",
      "We can show these relationships with some more or less basic trigonometry.\n",
      "\n",
      "Let's start by looking at the dot product of two rows from $\\C$.  We will take\n",
      "rows $\\vec{c_p} =\\C_{p,:}$ and $\\vec{c_q} = \\C_{q,:}$. As we remember, these\n",
      "vectors are:\n",
      "\n",
      "$$\n",
      "\\vec{c_p} = \\left[ \\cos(p n \\frac{2 \\pi}{N}) \\;\\mathrm{for}\\;\n",
      "n \\in 0,1,\\ldots,N-1 \\right] \\\\\n",
      "\\vec{c_q} = \\left[ \\cos(q n \\frac{2 \\pi}{N}) \\;\\mathrm{for}\\;\n",
      "n \\in 0,1,\\ldots,N-1 \\right]\n",
      "$$\n",
      "\n",
      "So:\n",
      "\n",
      "$$\n",
      "\\vec{c_p} \\cdot \\vec{c_q} = \\sum_{n=0}^{N-1} \\cos(p n \\frac{2 \\pi}{N}) \\cos(q\n",
      "n \\frac{2 \\pi}{N})\n",
      "$$\n",
      "\n",
      "Our trigonometry tells us that:\n",
      "\n",
      "$$\n",
      "\\cos \\alpha \\cos \\beta = \\frac{1}{2} [ \\cos(\\alpha + \\beta) - \\cos(\\alpha -\n",
      "\\beta) ]\n",
      "$$\n",
      "\n",
      "We can rewrite the dot product as the addition of two sums of cosines:\n",
      "\n",
      "$$\n",
      "\\vec{c_p} \\cdot \\vec{c_q} =\n",
      "\\frac{1}{2} \\sum_{n=0}^{N-1} \\cos((p + q) n \\frac{2 \\pi}{N}) +\n",
      "\\frac{1}{2} \\sum_{n=0}^{N-1} \\cos((p - q) n \\frac{2 \\pi}{N})\n",
      "$$\n",
      "\n",
      "Now we can use the formulae for sums of [arithmetic progressions of cosines and\n",
      "sines](sum_of_cosines.html) to solve these equations.  Here are the formulae:\n",
      "\n",
      "$$\n",
      "R \\triangleq \\frac{\\sin(N \\frac{1}{2}d)}{\\sin(\\frac{1}{2} d)} \\\\\n",
      "\\sum_{n=0}^{N-1} \\cos(a + nd) =\n",
      "\\begin{cases}\n",
      "N \\cos a & \\text{if } \\sin(\\frac{1}{2}d) = 0 \\\\\n",
      "R \\cos ( a + (N - 1) \\frac{1}{2} d) & \\text{otherwise}\n",
      "\\end{cases}\n",
      "\\\\\n",
      "\\sum_{n=0}^{N-1} \\sin(a + nd) =\n",
      "\\begin{cases}\n",
      "N \\sin a & \\text{if } \\sin(\\frac{1}{2}d) = 0 \\\\\n",
      "R \\sin ( a + (N - 1) \\frac{1}{2} d) & \\text{otherwise}\n",
      "\\end{cases}\n",
      "$$\n",
      "\n",
      "For our $\\C, \\S$ row dot product sums, starting angle $a$ is always 0, and the\n",
      "$d$ value in the formulae are always integer multiples of $\\frac{2 \\pi}{N}$.\n",
      "For example, $d = (p \\pm q) \\frac{2 \\pi}{N}$ in the equations above.  For our\n",
      "case, we can write $d = g \\frac{2 \\pi}{N}$ where $g$ is an integer.\n",
      "\n",
      "$$\n",
      "R = \\frac{\n",
      "\\sin( g N \\frac{1}{2} \\frac{2 \\pi}{N} )\n",
      "}\n",
      "{\n",
      "\\sin( g \\frac{1}{2} \\frac{2 \\pi}{N} )\n",
      "} \\\\\n",
      "= \\frac{ \\sin( g \\pi ) } { \\sin( \\frac{g}{N} \\pi ) }\n",
      "$$\n",
      "\n",
      "Because $g$ is an integer, the numerator of $R$ will always be 0, so the\n",
      "resulting sum is zero unless the denominator of $R$ is zero.  The denominator\n",
      "is zero only if $g$ is a multiple of N, including 0.  When the denominator is\n",
      "zero, the sum will be equal to $N \\cos(a) = N \\cos(0) = N$ for a cosine series\n",
      "or $N \\sin(a) = N \\sin(0) = 0$ for a sine series.\n",
      "\n",
      "Now we can calculate our dot product:\n",
      "\n",
      "$$\n",
      "\\vec{c_p} \\cdot \\vec{c_q} =\n",
      "\\begin{cases}\n",
      "\\frac{1}{2} N + \\frac{1}{2} N = N & \\text{if } p = q, p \\in 0, N/2 \\\\\n",
      "\\frac{1}{2} N & \\text{if } p = q, p \\notin 0, N/2 \\\\\n",
      "\\frac{1}{2} N & \\text{if } p + q = N, p \\ne N/2 \\\\\n",
      "0 & \\text{otherwise}\n",
      "\\end{cases}\n",
      "$$\n",
      "\n",
      "We can apply the same kind of logic to the rows of $\\S$:\n",
      "\n",
      "$$\n",
      "\\sin \\alpha \\sin \\beta = \\frac{1}{2} [ \\cos(\\alpha - \\beta) - \\cos(\\alpha +\n",
      "\\beta) ]\n",
      "$$\n",
      "\n",
      "So:\n",
      "\n",
      "$$\n",
      "\\vec{s_p} \\cdot \\vec{s_q} =\n",
      "\\frac{1}{2} \\sum_{n=0}^{N-1} \\cos((p - q) n \\frac{2 \\pi}{N}) -\n",
      "\\frac{1}{2} \\sum_{n=0}^{N-1} \\cos((p + q) n \\frac{2 \\pi}{N})\n",
      "$$\n",
      "\n",
      "This gives:\n",
      "\n",
      "$$\n",
      "\\vec{s_p} \\cdot \\vec{s_q} =\n",
      "\\begin{cases}\n",
      "0 & \\text{if } p = q, p \\in 0, N/2 \\\\\n",
      "\\frac{1}{2} N & \\text{if } p = q, p \\notin 0, N/2 \\\\\n",
      "-\\frac{1}{2} N & \\text{if } p + q = N, p \\ne N/2 \\\\\n",
      "0 & \\text{otherwise}\n",
      "\\end{cases}\n",
      "$$\n",
      "\n",
      "\n",
      "### Introducing vector projection\n",
      "\n",
      "If you are not familiar with projection, I highly recommend the tutorials over\n",
      "at [Khan academy](https://www.khanacademy.org/math/linear-\n",
      "algebra/matrix_transformations/lin_trans_examples/v/introduction-to-\n",
      "projections).\n",
      "\n",
      "If you know projection, you may think of a dot product like $\\vec{x} \\cdot\n",
      "\\vec{c_k}$ as part of the projection of our input signal $\\vec{x}$ onto the\n",
      "cosine vector $\\vec{c_k}$.\n",
      "\n",
      "Projection involves calculating the amount of a particular signal vector (such\n",
      "as a cosine wave) in another signal vector (such as our input data $\\vec{x}$).\n",
      "\n",
      "The Pearson product-moment correlation coefficient uses the dot product to test\n",
      "for relationship between two variables.  In our case, except for the first\n",
      "cosine vector $\\vec{c_0} = \\vec{1}$, the dot products $\\vec{x} \\cdot \\vec{c_k}$\n",
      "and $\\vec{x} \\cdot \\vec{s_k}$ are proportional to the Pearson product-moment\n",
      "correlation coefficient between $\\vec{c_k}$ and $\\vec{x}$ or $\\vec{s_k}$ and\n",
      "$\\vec{x}$, respectively.\n",
      "\n",
      "The projection of a vector $\\vec{a}$ onto a vector $\\vec{b}$ is given by:\n",
      "\n",
      "$$\n",
      "proj_{\\vec{b}}\\vec{a} \\triangleq g \\vec{b}\n",
      "$$\n",
      "\n",
      "where $g$ is a scalar that we will call the *projection coefficient*:\n",
      "\n",
      "$$\n",
      "g = \\frac{\\vec{a} \\cdot \\vec{b}}{\\vec{b} \\cdot \\vec{b}}\n",
      "$$\n",
      "\n",
      "Note that $\\vec{b} \\cdot \\vec{b}$ is also $\\| \\vec{b} \\|^2$, so we can also\n",
      "write:\n",
      "\n",
      "$$\n",
      "g = \\frac{\\vec{a} \\cdot \\vec{b}}{\\| \\vec{b} \\|^2}\n",
      "$$\n",
      "\n",
      "The result of the projection $proj_{\\vec{b}}\\vec{a}$ is a copy of $\\vec{b}$\n",
      "scaled by $g$ - the scalar amount of $\\vec{a}$ present in $\\vec{b}$.\n",
      "\n",
      "### Forward and inverse DFT as vector projection\n",
      "\n",
      "### Projection and the DFT\n",
      "\n",
      "The principle of the DFT on real input is the following.\n",
      "\n",
      "In the forward transform:\n",
      "\n",
      "* We calculate the data we need to form the projection coefficients for\n",
      "  projecting the input data onto the cosines and sine waves in the rows of\n",
      "  $\\C$ and $\\S$.\n",
      "* The projection data for the cosines goes into the real part of $\\vec{X}$ :\n",
      "  $\\R{\\vec{X}}$.  The projection data for the sines goes into the imaginary\n",
      "  part $\\I{\\vec{X}}$;\n",
      "\n",
      "In the inverse transform:\n",
      "\n",
      "* We complete the calculation of the projection coefficients $g$ for each\n",
      "  cosine and sine wave in $\\C, \\S$;\n",
      "* We use the projection coefficients to project the original data $\\vec{X}$ onto\n",
      "the set of cosines and sines in $\\C$, $\\S$. Each projection forms a new output\n",
      "vector, to give projection vectors $[proj_{\\vec{c_0}} \\vec{x}, proj_{\\vec{c_1}}\n",
      "\\vec{x}, ..., proj_{\\vec{s_0}} \\vec{x}, proj_{\\vec{s_1}} \\vec{x}, ...]$;\n",
      "* We sum up all the projection vectors to reconstruct the original data\n",
      "  $\\vec{X}$.\n",
      "\n",
      "This is how it works in principle.  There are some complications to the way it\n",
      "works in practice, due to the generality of the DFT in accepting real and\n",
      "complex input.  In the next sections we will go through some examples to show\n",
      "how\n",
      "the forward and inverse transform work in detail."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Does it actually work?\n",
      "unique_Cs = C[:N/2+1, :]\n",
      "unique_Ss = S[1:N/2, :]\n",
      "small_n = len(unique_Ss)\n",
      "cos_dots = unique_Cs.dot(x)\n",
      "sin_dots = unique_Ss.dot(x)\n",
      "cos_gs = cos_dots / ([N] + [N/2] * small_n + [N])\n",
      "sin_gs = sin_dots / ([N/2] * small_n)\n",
      "cos_projections = cos_gs[:, None] * unique_Cs\n",
      "sin_projections = sin_gs[:, None] * unique_Ss\n",
      "x_back = np.sum(np.vstack((cos_projections, sin_projections)), axis=0)\n",
      "x_back - x"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "31"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### The first element in $\\vec{X}$ for real input\n",
      "\n",
      "From our matrix multiplication, we know the first element of $\\vec{X}$ comes\n",
      "from:\n",
      "\n",
      "$$\n",
      "X_0 = C_{0, :} \\cdot \\vec{x} - i S_{0, :} \\cdot \\vec{x}\n",
      "    = \\vec{c_0} \\cdot \\vec{x} - i \\vec{s_0} \\cdot \\vec{x}\n",
      "$$\n",
      "\n",
      "We can represent this by highlighting the relevant parts of the matrix\n",
      "multiplication:\n",
      "\n",
      "We can simplify further because we know what $\\vec{c_0}$ and $\\vec{s_0}$ are:\n",
      "\n",
      "$$\n",
      "X_0 = \\vec{1} \\cdot \\vec{x} - i \\vec{0} \\cdot \\vec{x}\n",
      "    = \\vec{1} \\cdot \\vec{x}\n",
      "$$\n",
      "\n",
      "This final dot product can also be written as:\n",
      "\n",
      "$$\n",
      "X_0 = \\sum_{n=0}^{N-1}{x_n}\n",
      "$$\n",
      "\n",
      "That is, $X_0$ is a complex number with imaginary part = 0, where the real\n",
      "part contains the sum of the elements in $\\vec{x}$.\n",
      "\n",
      "Is this true of our original input vector $\\vec{x}$?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Sum of x', np.sum(x))\n",
      "print('First DFT coefficient X[0]', X[0])"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "32"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can show how $X_0$ comes about in the matrix multiplication by highlighting\n",
      "\n",
      "* $X_0$;\n",
      "* the relevant row of $\\C$ : $C_{0,:}$;\n",
      "* the vector $\\vec{x}$.\n",
      "\n",
      "We can leave out the relevant row of $\\S$ : $S_{0,:}$ because it is all 0."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sketch = dftp.DFTSketch(x)\n",
      "sketch.sketch(figsize=(12, 5))\n",
      "sketch.highlight('X_real', [0])\n",
      "sketch.highlight('C', [[0, ':']])\n",
      "sketch.highlight('x_c', [':'])"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "33"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### DFT of a constant input vector\n",
      "\n",
      "Next we will consider the forward and inverse DFT of an input vector that is\n",
      "constant.\n",
      "\n",
      "Our input is vector with $N$ elements, where every element = 2:\n",
      "\n",
      "$$\n",
      "\\vec{w} = [2, 2, ... 2]\n",
      "$$\n",
      "\n",
      "We could also write $\\vec{w}$ as $\\vec{2}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w = np.ones(N) * 2\n",
      "w"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "34"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What DFT output $\\vec{W}$ will this generate?\n",
      "\n",
      "We already know that $W_0$ must be the sum of $\\vec{w}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "W = np.fft.fft(w)\n",
      "print('Sum of w', np.sum(w))\n",
      "print('First DFT coefficient W[0]', W[0])"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "35"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How about the rest of $\\vec{W}$?  All the remaining cosine and sine waves in\n",
      "$\\C, \\S$ sum to zero over the rows (and columns):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Sums over rows of C after first', np.sum(C[1:], axis=1))\n",
      "print('Sums over rows of S', np.sum(S, axis=1))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "36"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For any vector $\\vec{v}$ that sums to zero, the dot product $\\vec{2} \\cdot\n",
      "\\vec{v}$ will be $\\sum_{n=0}^{N-1} 2 v_n = 2 \\sum_{n=0}^{N-1} v_n = 0$.\n",
      "\n",
      "So, we predict that all the rest of $W_0$, real and imaginary, will be 0:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "W"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "37"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us show this in the matrix form:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sketch = dftp.DFTSketch(w)\n",
      "sketch.sketch(figsize=(12, 5))\n",
      "sketch.highlight('X_real', [0])\n",
      "sketch.highlight('C', [[0, ':']])\n",
      "sketch.highlight('x_c', [':'])"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "38"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Cosines in the real part, sines in the imaginary part\n",
      "\n",
      "The following only applies to real input to the DFT.\n",
      "\n",
      "From the forward DFT formula on a vector of real numbers, we see that the\n",
      "$\\R{X}$ will contain the dot product of $\\vec{x}$ with the cosine basis, and\n",
      "$\\I{X}$ will contain the dot product of $\\vec{x}$ with the sine basis.\n",
      "\n",
      "Imagine, for simplicity, that $\\vec{s_k} \\cdot \\vec{x} = 0$ for every $k$, or\n",
      "(saying the same thing in a different way) $\\S \\cdot \\vec{x} = \\vec{0}$.\n",
      "\n",
      "In that case our forward DFT would be:\n",
      "\n",
      "$$\n",
      "X = \\C \\cdot \\vec{x}\n",
      "$$\n",
      "\n",
      "and the inverse DFT would be:\n",
      "\n",
      "$$\n",
      "X = \\frac{1}{N} \\C \\cdot \\vec{X}\n",
      "$$\n",
      "\n",
      "In that case, $\\vec{X}$ would be a vector of real numbers, each expressing the\n",
      "amount of the corresponding cosine vector is present in the data.\n",
      "\n",
      "We could then perfectly reconstruct our original data by summing up the result\n",
      "of projecting onto each cosine vector.\n",
      "\n",
      "In the case of our constant input vector $\\vec{w}$, this is the case - there are\n",
      "no sine components in $\\vec{w}$ and $\\S \\cdot \\vec{x} = \\vec{0}$.\n",
      "\n",
      "So, $\\R{\\vec{X}}$ contains all the information in $\\vec{w}$.  In fact, as we\n",
      "have seen, $\\R{X_0}$ contains all the information in $\\vec{w}$.\n",
      "\n",
      "Rephrasing in terms of projection, $W_0$ comes from $\\vec{1} \\cdot \\vec{w}$.\n",
      "This the top half of the $g$ value for projecting the signal $\\vec{w}$ onto a\n",
      "vector of ones $\\vec{c_0}$ : $g = \\frac{\\vec{w} \\cdot \\vec{1}}{\\vec{1} \\cdot\n",
      "\\vec{1}}$. We know $\\vec{1} \\cdot \\vec{1} = N$ so the projection of $\\vec{w}$\n",
      "onto $\\vec{1}$ is $\\frac{\\vec{w} \\cdot \\vec{1}}{\\vec{1} \\cdot \\vec{1}} \\vec{1} =\n",
      "\\frac{1}{N} \\vec{w} \\cdot \\vec{1}$, and this is precisely what the inverse DFT\n",
      "will do:\n",
      "\n",
      "$$\n",
      "w_0 = \\frac{1}{N} \\vec{1} \\cdot \\vec{W} - i \\frac{1}{N} \\vec{0} \\cdot \\vec{W} =\n",
      "\\frac{1}{N} \\vec{1} \\cdot \\vec{W} \\\\\n",
      "w_1 = \\frac{1}{N} \\vec{1} \\cdot \\vec{W} \\\\\n",
      "... \\\\\n",
      "w_{N-1} = \\frac{1}{N} \\vec{1} \\cdot \\vec{W}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w_again = np.zeros(w.shape, dtype=np.complex)\n",
      "c_0 = np.ones(N)\n",
      "for n in np.arange(N):\n",
      "    w_again[n] = 1. / N * c_0.dot(W)\n",
      "w_again"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "39"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In matrix form:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "1. / N * C.dot(W)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "40"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sketch = dftp.DFTSketch(w)\n",
      "sketch.sketch(inverse=True, figsize=(12, 5))\n",
      "sketch.highlight('x_real', [':'])\n",
      "sketch.highlight('C', [[':', 0]])\n",
      "sketch.highlight('X_c_real', [0])"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "41"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### DFT on a signal with a single cosine\n",
      "\n",
      "Now let us look at the second coefficient, $X_1$.\n",
      "\n",
      "This was formed by dot products of the signal with cosine and sine waves having\n",
      "a single cycle across the whole signal:\n",
      "\n",
      "$$\n",
      "\\vec{t_1} \\triangleq \\left[ 2 \\pi \\frac{n}{N} \\;\\mathrm{for}\\; n \\in\n",
      "0,1,\\ldots,N-1 \\right] \\\\\n",
      "\\vec{c_1} \\triangleq \\left[ \\cos(t_{1,n}) \\;\\mathrm{for}\\; n \\in 0,1,\\ldots,N-1\n",
      "\\right] \\\\\n",
      "\\vec{s_1} \\triangleq \\left[ \\sin(t_{1,n}) \\;\\mathrm{for}\\; n \\in 0,1,\\ldots,N-1\n",
      "\\right]\n",
      "$$\n",
      "\n",
      "Here are plots of $\\vec{c_1}, \\vec{s_1}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ns = np.arange(N)\n",
      "t_1 = 2 * np.pi * ns / N\n",
      "plt.plot(ns, np.cos(t_1), 'o:')\n",
      "plt.plot(ns, np.sin(t_1), 'o:')\n",
      "plt.xlim(0, N-1)\n",
      "plt.xlabel('n')"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "42"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If the input signal is a single cosine wave of amplitude 3, with one cycle over\n",
      "the signal, then we can predict $X_1$.  It will be the dot product of the input\n",
      "signal with $c_1$, which is the same as $3 c_1 \\cdot c_1$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "t_1 = 2 * np.pi * ns / N\n",
      "cos_x = 3 * np.cos(t_1)\n",
      "c_1 = np.cos(t_1)\n",
      "X = np.fft.fft(cos_x)\n",
      "print('First DFT coefficient for single cosine', X[1])\n",
      "print('Dot product of single cosine with c_1', cos_x.dot(c_1))\n",
      "print('3 * dot product of c_1 with itself', 3 * c_1.T.dot(c_1))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "43"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Fitting all cosine phases with an added sine\n",
      "\n",
      "Now it is time to bring the $i \\vec{x} \\cdot \\vec{s_k}$ part of the DFT into\n",
      "play.\n",
      "\n",
      "By calculating the dot product of our input vector with a cosine wave of a given\n",
      "frequency, we detect any signal that matches that cosine with the given phase\n",
      "and the given frequency. In our example above, we used the DFT $\\vec{c_1}$ dot\n",
      "product to detect a cosine with phase offset 0 - the cosine starts at $n = 0$.\n",
      "\n",
      "What happens if the cosine in the signal has a different phase?  For example,\n",
      "what happens to the dot product if the cosine wave in our data is shifted by\n",
      "0.8,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cos_x_shifted = 3 * np.cos(t_1 + 0.8)\n",
      "plt.plot(t_1, cos_x_shifted)\n",
      "print('Dot product of shifted cosine with c_1', cos_x_shifted.dot(c_1))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "44"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "When the cosine wave is shifted in our data, relative to the $\\vec{c_1}$, then\n",
      "the dot product of the signal against $\\vec{c_1}$ drops in value, and is\n",
      "therefore less successful at detecting this cosine wave.\n",
      "\n",
      "This is the role of the $\\vec{s_k}$ vectors in the DFT.  By calculating dot\n",
      "products with the $\\vec{s_k}$ vectors, we can detect cosine waves of any phase.\n",
      "\n",
      "Let us see that in action first, and then explain why this is so.\n",
      "\n",
      "First, here is what happens to the dot products for the shifted and unshifted\n",
      "cosine waves:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s_1 = np.sin(t_1)\n",
      "plt.plot(t_1, cos_x, label='3 * cos wave')\n",
      "plt.plot(t_1, cos_x_shifted, label='3 * cos wave, shifted')\n",
      "plt.legend()\n",
      "print('Dot product of unshifted cosine with c_1', cos_x.dot(c_1))\n",
      "print('Dot product of unshifted cosine with s_1', cos_x.dot(s_1))\n",
      "print('Dot product of shifted cosine with c_1', cos_x_shifted.dot(c_1))\n",
      "print('Dot product of shifted cosine with s_1', cos_x_shifted.dot(s_1))"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "45"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that the dot product with $\\vec{s_1}$ is effectively zero in the\n",
      "unshifted case, and goes up to around 34 in the shifted case.\n",
      "\n",
      "Now let us use the projections from these dot products to reconstruct the\n",
      "original vector (as we will soon do using the inverse DFT).\n",
      "\n",
      "First we use the dot product with $\\vec{c_1}$ to reconstruct the unshifted\n",
      "cosine (the dot product with $\\vec{s_1}$ is zero, so we do not need it)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reconstruct unshifted cos from dot product projection\n",
      "c_unshifted = cos_x.dot(c_1) / c_1.dot(c_1)\n",
      "proj_onto_c1 = c_unshifted * c_1\n",
      "plt.plot(ns, proj_onto_c1)\n",
      "plt.title('Reconstructed unshifted cosine')"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "46"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can use the cosine and sine dot product to reconstruct the shifted cosine\n",
      "vector:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reconstruct shifted cos from dot product projection\n",
      "c_cos_shifted = cos_x_shifted.dot(c_1) / c_1.dot(c_1)\n",
      "c_sin_shifted = cos_x_shifted.dot(s_1) / s_1.dot(s_1)\n",
      "proj_onto_c1 = c_cos_shifted * c_1\n",
      "proj_onto_s1 = c_sin_shifted * s_1\n",
      "reconstructed = proj_onto_c1 + proj_onto_s1\n",
      "plt.plot(ns, reconstructed)\n",
      "plt.title('Reconstructed shifted cosine')\n",
      "assert np.allclose(reconstructed, cos_x_shifted)"
     ],
     "language": "python",
     "metadata": {
      "attributes": {
       "classes": [],
       "id": "",
       "n": "47"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The reason that this works for any phase shift is the angle sum rule.\n",
      "\n",
      "The angle sum rule is:\n",
      "\n",
      "$$\n",
      "\\cos(\\alpha \\pm \\beta) = \\cos \\alpha \\cos \\beta \\mp \\sin \\alpha \\sin \\beta\n",
      "$$\n",
      "\n",
      "To unpack the $\\pm, \\mp$:\n",
      "\n",
      "$$\n",
      "\\cos(\\alpha + \\beta) = \\cos \\alpha \\cos \\beta - \\sin \\alpha \\sin \\beta \\\\\n",
      "\\cos(\\alpha - \\beta) = \\cos \\alpha \\cos \\beta + \\sin \\alpha \\sin \\beta\\\n",
      "$$\n",
      "\n",
      "See [angle sum proof](https://perrin.dynevor.org/angle_sum.html) for a visual\n",
      "proof in the case of real angles $\\alpha, \\beta$."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}