{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Numpy DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy has lots of functions for executing the DFT.  Let's take a look at what they all do.\n",
    "\n",
    "Functions explored in this tutorial: **`fft.fft`**, **`fftshift`**, **`fft.ifft`**, **`fft.rfft`**, **`fft.irfft`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import packages\n",
    "import csv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import *\n",
    "\n",
    "# Display plots in the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll construct and plot a simple function `x(t)`.  Note that the time series consists of 201 data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t = np.linspace(0,200,201)\n",
    "x = 1+.5*np.cos(t*2*np.pi/20)\n",
    "\n",
    "plt.plot(t,x)\n",
    "xlabel('t')\n",
    "ylabel('x(t)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll take the Discrete Fourier Transform (DFT) of the time-series data.  All numpy DFT commands live in the `fft` library (an abbreviation for \"Fast Fourier Transform,\" a set of algorithms to speed up the DFT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_fft = np.fft.fft(x)\n",
    "k = np.linspace(0,200,201)\n",
    "\n",
    "plt.plot(k,x_fft)\n",
    "xlabel('k')\n",
    "ylabel('X_k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hm, something's not quite right... the **`fft.fft`** command returns 201 coefficients (as you would expect from the DFT algorithm), but just as if you evaluated the sums by hand, those coefficients are often complex!  IPython throws a warning because we are asking it to plot complex numbers on a real graph.  (It deals with this by plotting only the real component of each value, but that's not really what we want.)  Let's fix this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(k,abs(x_fft))\n",
    "xlabel('k')\n",
    "ylabel('X_k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, that's better.  There are three peaks in the spectrum (why?), but we can hardly see the DC term at `k=0` because it is right on the edge of the graph.  Fortunately, python provides a built-in function **`fftshift`** that centers the frequency-domain plot on `k=0`.  (This is totally OK because the DFT is periodic in N.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_fftshift=fftshift(abs(x_fft))\n",
    "k_fftshift=np.linspace(-100,100,201)\n",
    "\n",
    "plt.plot(k_fftshift,x_fftshift)\n",
    "xlabel('k')\n",
    "ylabel('X_k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use **`fft.ifft`** to take the inverse transform of this frequency spectrum, which should return the time-series data that we started with.  (Note that we invert the original fft coefficients, **not** the versions resequenced with fftshift.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_ifft = np.fft.ifft((x_fft))\n",
    "plt.plot(k,x_ifft)\n",
    "xlabel('t')\n",
    "ylabel('x(t)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like it worked, but we got another warning about plotting complex numbers!  The inverse transform should just return the real-valued time-series data we started with... let's figure out what's going on here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"The first time-series value is:\", x_ifft[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like because of rounding errors, our \"real\" data points have accumulated tiny imaginary components during the inverse transform.  This won't affect the result much (considering the real portion of each data point is an extremely good approximation in this case), but to make the warning disappear we can again find the magnitude of the points before plotting them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(t,abs(x_ifft))\n",
    "xlabel('t')\n",
    "ylabel('x(t)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's examine another FFT function, **`fft.rfft`**.  Short for \"real FFT\", this function assumes that the input is real-valued, which means that the resulting spectrum has a well-known relationship between the first half of its coefficients and the second half.  (How are the coefficients in the two halves related?)  `fft.rfft` therefore only returns the first half of the coefficients - it assumes you are clever enough to figure out the second half if you need to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_rfft=np.fft.rfft(x)\n",
    "print(\"For an input of\",len(x),\"data points, rfft returned only\",len(x_rfft),\"coefficients.\")\n",
    "plt.plot(k[0:len(x_rfft)],abs(x_rfft))\n",
    "xlabel('k')\n",
    "ylabel('X_k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy also supplies an **`fft.irfft`** function, which takes as input half of a frequency spectrum of a real-valued signal (such as we just found using `fft.rfft`).  `fft.irfft` infers the other half of the spectrum and then calculates the inverse DFT to return time-series data.\n",
    "\n",
    "`fft.rfft` has to make certain assumptions about how the final frequency coefficient value relates to the original time-series data (see the documentation of this function for more information).  The result of this is that it sometimes returns a time series that is shorter by one point than you might expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_irfft=np.fft.irfft(x_rfft)\n",
    "print(\"Length of original time series:\",len(x))\n",
    "print(\"Length of the irfft series:\",len(x_irfft))\n",
    "plt.plot(t[:-1],x_irfft)\n",
    "xlabel('t')\n",
    "ylabel('x(t)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, note that because `fft.rfft` and `fft.irfft` are intended to deal only with real-valued time series, the results of `fft.rfft` are real-valued, and we do not need an extra magnitude calculation to clean up rounding errors."
   ]
  }
 ],
 "metadata": {
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}