Fourier Intro -- 4 pt FT

gistfile1.txt

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<script type=\"text/javascript\">\n",
      "    var div_cell = \"div.cell.text_cell{ width:1000px; margin-left:auto; margin-right:auto;}\"\n",
      "    $('head').append( \"<style type='text/css'>\" + div_cell + \"</style>\");\n",
      "</script>\n",
      "\n",
      "# Fourier Transform Introduction\n",
      "\n",
      "For a continuous signal, the [Fourier transform](http://en.wikipedia.org/wiki/Fourier_transform) (FT, Equation 1) is simply the integral of a time-domain (complex) function, $f(t)$, multiplied by an imaginary exponential function.\n",
      "\n",
      "$$\\mathrm{Equation\\ 1:}\\ \\ \\ F( \\nu ) = \\int_{-\\infty}^{+\\infty} f(t) e^{- i 2 \\pi \\nu t} \\mathrm{d}t$$\n",
      "\n",
      "There are two things to keep in mind here. First of all, the frequency term in the exponential is equivalent to radians (Equation 2). Also, an exponential of an imaginary number can be expanded by [Euler's Formula](http://en.wikipedia.org/wiki/Euler%27s_formula) (Equation 3).\n",
      "\n",
      "$$\\mathrm{Equation\\ 2:}\\ \\ \\ \\omega=2 \\pi \\nu$$\n",
      "\n",
      "$$\\mathrm{Equation\\ 3:}\\ \\ \\ e^{- i 2 \\pi \\nu t}=\\cos(2 \\pi \\nu t) + i \\sin(2 \\pi \\nu t)$$\n",
      "\n",
      "In this case, Equation 1 can be expanded to a form in terms of sines and cosines (Equation 4).\n",
      "\n",
      "$$\\mathrm{Equation\\ 4:}\\ \\ \\ F( \\nu ) = \\int_{-\\infty}^{+\\infty} f(t) \\cos(2 \\pi \\nu t) \\mathrm{d}t + i \\int_{-\\infty}^{+\\infty} f(t) \\sin(2 \\pi \\nu t) \\mathrm{d}t$$\n",
      "\n",
      "For functions that are easily integrable and infinite, this equation can be solved exactly. However, in most real situations (such as a spectroscopic signal) this time domain function is a very complex, finite signal that can not be easily integrated. In this case, we need to use a [Discrete Fourier Transform](http://en.wikipedia.org/wiki/Discrete_Fourier_transform) (Equation 5).\n",
      "\n",
      "$$\\mathrm{Equation\\ 5:}\\ \\ \\ F( \\nu_{m} ) = \\sum_{n=0}^{N-1} x_{n} e^{- i 2 \\pi n m / N}$$\n",
      "\n",
      "Where the *x* values are the *N* total discrete data points, and *m/N* are the frequencies for the complex sinusoidal functions with the following *m* values: \n",
      "\n",
      "$$m( \\mathrm{even} ) = \\frac{-N}{2}, \\frac{-N}{2}+1, \\cdots , \\frac{-N}{2} + (N-1)$$\n",
      "\n",
      "$$m( \\mathrm{odd} ) = \\frac{-(N-1)}{2}, \\frac{-(N-1)}{2}+1, \\cdots , \\frac{-(N-1)}{2} + N$$\n",
      "\n",
      "Again this can be expanded using Euler's Formula. \n",
      "\n",
      "$$\\mathrm{Equation\\ 6:}\\ \\ \\ F( \\nu_{m} ) = \\sum_{n=0}^{N-1} x_{n} [\\cos(2 \\pi n m / N) + i \\sin(2 \\pi n m / N)]$$\n",
      "\n",
      "Fortunately, efficient computational alogrithms were determined for equation many years ago and are collectively known as [Fast Fourier Transformations](http://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT). It is not an overstatement to say that these algorithms revolutionized spectrocopic methods. \n",
      "\n",
      "Below, we will first look at a very simple, contrived example FT just to make sure that we understand how it works in general. Next, we will explore some of the properties of this transformation by looking at some simulated data that more accurately reflects real world situations.\n",
      "\n",
      "First of all, be sure to import Numpy and Matplotlib, which have the code for FT and plotting, respecively. The 'inline' statement will put all of the pictures into the web browswer window."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "# This statement below sets some printing options so we \n",
      "# don't get so many digits after the decimal (precision) and lots of very small floats (suppress)\n",
      "np.set_printoptions(precision=5, suppress=True) "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## \"Simple\" FT Example\n",
      "\n",
      "First of all, let's try to perform a 4-point FT on an arbitrary set of data. In this case, we'll focus on Equation 6 from above, which for 4 data points will have the following values \n",
      "\n",
      "* *n* values will range from 0 to 3\n",
      "* *m* values will range from (-4/2)+0,...,(-4/2)+3 = -2, -1, 0, 1"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = np.arange(4)\n",
      "m = np.array([0., -1., -2., 1.]) # This goofy ordering has important consequences as seen below"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we are going to need to calculate our angular 2*&#960;*n*m/N terms. We need to calculate the product for each n value at every freqency, which means we will need to ultimately get a 2D array of numbers as shown below. \n",
      "\n",
      "* theta1 = (m0,n0), (m0, n1), (m0, n2), (m0, n3)\n",
      "* theta2 = (m1,n0), (m1, n1), (m1, n2), (m1, n3)\n",
      "* theta3 = (m2,n0), (m2, n1), (m2, n2), (m2, n3)\n",
      "* theta4 = (m3,n0), (m3, n1), (m3, n2), (m3, n3)\n",
      "\n",
      "To efficiently accomplish this, rotate the *m* array so that it will be a column rather than a row vector. Then we can calculate the product, which generates a 2D array. The rows of our 2D array will be individual frequencies (m/N values), and the columns will be each *n* value.\n",
      "\n",
      "Print out some of the values so that we can check some intermediate steps of the calculation. Also, we will print the cosines and sines of these values to see what our 2D matrices look like."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = m.reshape(-1, 1)\n",
      "theta = 2*np.pi*n*m/4.\n",
      "\n",
      "print '2*pi*n Values'\n",
      "print 2*np.pi*n\n",
      "print 'm/4 Values'\n",
      "print m/4.\n",
      "print '2D angular array.'\n",
      "print theta\n",
      "print 'Cosines of the angle array.'\n",
      "print np.cos(theta)\n",
      "print 'Sines of the angle array.'\n",
      "print np.sin(theta)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2*pi*n Values\n",
        "[  0.        6.28319  12.56637  18.84956]\n",
        "m/4 Values\n",
        "[[ 0.  ]\n",
        " [-0.25]\n",
        " [-0.5 ]\n",
        " [ 0.25]]\n",
        "2D angular array.\n",
        "[[ 0.       0.       0.       0.     ]\n",
        " [-0.      -1.5708  -3.14159 -4.71239]\n",
        " [-0.      -3.14159 -6.28319 -9.42478]\n",
        " [ 0.       1.5708   3.14159  4.71239]]\n",
        "Cosines of the angle array.\n",
        "[[ 1.  1.  1.  1.]\n",
        " [ 1.  0. -1. -0.]\n",
        " [ 1. -1.  1. -1.]\n",
        " [ 1.  0. -1. -0.]]\n",
        "Sines of the angle array.\n",
        "[[ 0.  0.  0.  0.]\n",
        " [-0. -1. -0.  1.]\n",
        " [-0. -0.  0. -0.]\n",
        " [ 0.  1.  0. -1.]]\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The second term in Equation 6 is the complex sum of the cosine and sine functions, so let's take a look at that array. \n",
      "\n",
      "Notice anything interesting about this array? It is symmetric about the diagonal. It turns out that this will always be the case if the number of data points is a power of two (e.g. 2**n). This symmetry allows for considerable reductions in the number of computations that need to be performed, and it is the basis of the original Fast Fourier Transform algorithms."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "euler = np.cos(theta) + np.sin(theta)*1j\n",
      "\n",
      "print euler"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]\n",
        " [ 1.+0.j  0.-1.j -1.-0.j -0.+1.j]\n",
        " [ 1.+0.j -1.-0.j  1.+0.j -1.-0.j]\n",
        " [ 1.+0.j  0.+1.j -1.+0.j -0.-1.j]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To finish things out, this complex array needs to be multiplied by every discrete data point in order to finish the FT calculation, so let's make up some imaginary data to see what this looks like. The frequency spectrum will be the sum of the rows after this multiplication. We can also compare our data with a FFT algorithm, specifically the one implemented in Numpy, [*numpy.fft.fft*](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fake = np.array([5., 7., 1., 3.])\n",
      "ft_fake = fake*euler\n",
      "\n",
      "print ft_fake\n",
      "print\n",
      "print 'Our FT analysis:', ft_fake.sum(axis=1) # Sum the 2D array across the rows to get the final FT.\n",
      "print 'FFT algorithm:', np.fft.fft( fake )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 5.+0.j  7.+0.j  1.+0.j  3.+0.j]\n",
        " [ 5.+0.j  0.-7.j -1.-0.j -0.+3.j]\n",
        " [ 5.+0.j -7.-0.j  1.+0.j -3.-0.j]\n",
        " [ 5.+0.j  0.+7.j -1.+0.j -0.-3.j]]\n",
        "\n",
        "Our FT analysis: [ 16.+0.j   4.-4.j  -4.-0.j   4.+4.j]\n",
        "FFT algorithm: [ 16.+0.j   4.-4.j  -4.+0.j   4.+4.j]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}