fourier_transform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison between Discrete and Fast Fourier Transform.\n",
    "#### Following JVP's online tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def DFT_slow(x):\n",
    "    \"\"\"\n",
    "    Compute the discrete Fourier Transform of the 1D array x.\n",
    "    This is a slow-running code on purpose.\n",
    "    \"\"\"\n",
    "    \n",
    "    x = np.asarray(x, dtype=float)\n",
    "    N = x.shape[0]\n",
    "    n = np.arange(N)\n",
    "    k = n.reshape((N, 1))\n",
    "    M = np.exp(-2j * np.pi * k * n / N) #Standard Fourier Transform equation\n",
    "    \n",
    "    return np.dot(M, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#I can check that function works against numpy's built in fft function\n",
    "x = np.random.random(1024)\n",
    "np.allclose(DFT_slow(x), np.fft.fft(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can speed up the DFT process by exploiting symmetry in what's known as FFT.\n",
    "\n",
    "Essentially, $X_{N+k}=X_{k}$. \n",
    "\n",
    "We can further the argument to complex equivalencies: $X_{iN+k}=X_{k}$\n",
    "\n",
    "This allows us to split Fourier Transform integrals into two. One for even integers and one for odd integers. We can see from symmetry that we would only need one of these, and we can therefore cut our computation in half."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def FFT(x):\n",
    "    \"\"\"\n",
    "    A recursive implementation of the 1D Cooley-Tukey FFT by JVP.\n",
    "    \"\"\"\n",
    "    x = np.asarray(x, dtype=float)\n",
    "    N = x.shape[0]\n",
    "    \n",
    "    if N % 2 > 0:\n",
    "        raise ValueError(\"size of x must be a power of 2\")\n",
    "    elif N <= 30:  # This cutoff has ben optimpized for the random array, but can be switched for our data\n",
    "        return DFT_slow(x)\n",
    "    else:\n",
    "        X_even = FFT(x[::2]) #This is a very matlab logic argument\n",
    "        X_odd = FFT(x[1::2])\n",
    "        factor = np.exp(-2j * np.pi * np.arange(N) / N) #The fourier transform equation\n",
    "        return np.concatenate([X_even + factor[:N // 2] * X_odd,\n",
    "                               X_even + factor[N // 2:] * X_odd]) #The equation solves N/2 for the odd numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Checking that this functions works\n",
    "np.allclose(FFT(x), np.fft.fft(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 84.4 ms per loop\n"
     ]
    }
   ],
   "source": [
    "#Now I can check if the timing changes!\n",
    "\n",
    "%timeit DFT_slow(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 4.56 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit FFT(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 23.66 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "10000 loops, best of 3: 43.2 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit np.fft.fft(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy's fft is clearly the fastest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fft(a, n=None, axis=-1, norm=None):\n",
    "    \"\"\"\n",
    "    Compute the one-dimensional discrete Fourier Transform.\n",
    "    This function computes the one-dimensional *n*-point discrete Fourier\n",
    "    Transform (DFT) with the efficient Fast Fourier Transform (FFT)\n",
    "    algorithm [CT].\n",
    "    Parameters\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    FFT (Fast Fourier Transform) refers to a way the discrete Fourier\n",
    "    Transform (DFT) can be calculated efficiently, by using symmetries in the\n",
    "    calculated terms.  The symmetry is highest when `n` is a power of 2, and\n",
    "    the transform is therefore most efficient for these sizes.\n",
    "    The DFT is defined, with the conventions used in this implementation, in\n",
    "    the documentation for the `numpy.fft` module.\n",
    "    References\n",
    "    ----------\n",
    "    .. [CT] Cooley, James W., and John W. Tukey, 1965, \"An algorithm for the\n",
    "            machine calculation of complex Fourier series,\" *Math. Comput.*\n",
    "            19: 297-301.\n",
    "            \n",
    "    Examples\n",
    "    --------\n",
    "    >>> np.fft.fft(np.exp(2j * np.pi * np.arange(8) / 8))\n",
    "    array([ -3.44505240e-16 +1.14383329e-17j,\n",
    "             8.00000000e+00 -5.71092652e-15j,\n",
    "             2.33482938e-16 +1.22460635e-16j,\n",
    "             1.64863782e-15 +1.77635684e-15j,\n",
    "             9.95839695e-17 +2.33482938e-16j,\n",
    "             0.00000000e+00 +1.66837030e-15j,\n",
    "             1.14383329e-17 +1.22460635e-16j,\n",
    "             -1.64863782e-15 +1.77635684e-15j])\n",
    "    In this example, real input has an FFT which is Hermitian, i.e., symmetric\n",
    "    in the real part and anti-symmetric in the imaginary part, as described in\n",
    "    the `numpy.fft` documentation:\n",
    "    >>> import matplotlib.pyplot as plt\n",
    "    >>> t = np.arange(256)\n",
    "    >>> sp = np.fft.fft(np.sin(t))\n",
    "    >>> freq = np.fft.fftfreq(t.shape[-1])\n",
    "    >>> plt.plot(freq, sp.real, freq, sp.imag)\n",
    "    [<matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>]\n",
    "    >>> plt.show()\n",
    "    \"\"\"\n",
    "\n",
    "    a = np.asarray(a).astype(complex, copy=False)\n",
    "    if n is None:\n",
    "        n = a.shape[axis]\n",
    "    output = _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache) #must look closer at these packs\n",
    "    if _unitary(norm):\n",
    "        output *= 1 / sqrt(n)\n",
    "    return output"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}