{
 "metadata": {
  "name": "ps01"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Discrete Fourier Transform of Cosine-Gaussian signal\n",
      "\n",
      "On this problem set you showed that the (continuous) Fourier transform of\n",
      "$$\n",
      "h(t) = h_0\\\\,e^{-\\alpha t^2} \\cos 2\\pi f_1 t\n",
      "$$\n",
      "is\n",
      "$$\n",
      "\\widetilde{h}(f) = \\frac{h_0}{2}\\sqrt{\\frac{\\pi}{\\alpha}}\n",
      "\\left(\n",
      "e^{-\\pi^2 (f+f_1)^2/\\alpha} + e^{-\\pi^2 (f-f_1)^2/\\alpha}\n",
      "\\right)\n",
      "$$\n",
      "You will now investigate the discrete Fourier transform as an approximation to this result."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Preliminaries\n",
      "\n",
      "We turn on matlab-like syntax for numpy and also allow inline matplotlib plots.  And we import the Python3 behavior of doing all division in floating point by default."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from __future__ import division"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Setting constants\n",
      "\n",
      "We consider what happens if we sample this signal with some discrete spacing $\\delta t$, and limit ourselves to some time duration $T$ centered on $t_0$.\n",
      "Note that every number handled by a numerical code is a dimensionless quantity, so we don't actually need to define $\\delta t$, and can instead define variables which are equal to $\\frac{T}{\\delta t}$, $f_1\\\\,\\delta t$ etc.  But the dimensinless parameter $N=\\frac{T}{\\delta t}$ does mean something, and in particular, $\\delta f=\\frac{N}{\\delta t}$, so we might get confused about whether we're keeping track of $\\frac{f_1}{\\delta f}$ or $f_1\\\\,\\delta t$. so we'll define both $\\delta t$ and $\\delta f$ to be safe, even though this is a bit redundant.  We also define the Nyquist frequency $f_{\\textrm{Ny}}=\\frac{1}{2\\delta t}$ for convenience."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deltat=1.; N=16384; T=N*deltat; deltaf=1./T; fNy=0.5/deltat;\n",
      "f1=1./(128.*deltat)\n",
      "alpha=1./(256.*deltat)**2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Defining the time variable\n",
      "\n",
      "We want to have N time points from $-\\frac{T}{2}$ to $\\frac{T}{2}-\\delta t$, at a spacing of $\\delta t$.  The variable `t` is actually $\\frac{t}{\\delta t}$,since we're measuring everything in units of $\\delta t$, but we'll pretend we don't know that, in case we later decide to measure everything in units of $T$ or something. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "t = arange(-N//2,N//2) * deltat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The time series\n",
      "\n",
      "We define a variable `ht` which is $h(t)/h_0$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ht = exp(-alpha*t**2)*cos(2*pi*f1*t)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can plot this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(t/deltat,ht); xlim([-0.5*T/deltat,0.5*T/deltat]); xlabel(r'$t/\\delta t$'); ylabel(r'$h(t)/h_0$'); xticks(linspace(-0.5*T/deltat,0.5*T/deltat,9)); grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can zoom in and see the cosine-Gaussian structure:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(t*sqrt(alpha),ht); xlim([-3,3]); xlabel(r'$t\\sqrt{\\alpha}$'); ylabel(r'$h(t)/h_0$'); grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Calculate the discrete Fourier transform\n",
      "\n",
      "We know that $\\widehat{h}_k\\\\,\\delta t$ should be an approximation to $\\widetilde{h}(f_k)$, so we calculate this using numpy's `fft` function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hf=fftshift(fft.fft(fftshift(ht)))*deltat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A couple of notes here:\n",
      "\n",
      "  * Numpy's `fft` routine already has the same conventions as we use in class for the dft, which we can see by either looking up\n",
      "    [the `numpy.fft` documentation online](http://docs.scipy.org/doc/numpy-1.7.0/reference/routines.fft.html#implementation-details)\n",
      "    or checking the documentation via ipython:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "?fft"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* We have used [the `fftshift` function](docs.scipy.org/doc/numpy-1.7.0/reference/generated/numpy.fft.fftshift.html) to exchange the two halves of the data array before and after taking the discrete Fourier transform.  That's because numpy's `fft` routine naturally operates on indices from `0` to `N-1`, but it makes more sense to define both our time and frequency arrays to range from indices `-N/2` to `N/2-1`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "?fftshift"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Define the frequency vector"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f=np.arange(-N/2,N/2)*deltaf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Plot the discrete approximation to the Fourier transform\n",
      "\n",
      "Note that in general the Fourier transform is complex, but the symmetry of the cosine-Gaussian means that the imaginary part is zero up to machine precision:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "max(abs(imag(hf)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "so we can just plot the real part:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(f/f1,real(hf)/T); xlim([-2*f1/f1,2*f1/f1]); xlabel(r'$f/f_1$'); ylabel(r'$\\widetilde{h}(f)/(h_0 T)$'); grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that we've plotted $\\frac{\\widetilde{h}(f)}{h_0T}$ vs $\\frac{f}{f_1}$ so that both quantities are dimensionless.\n",
      "\n",
      "Define and plot the continuous Fourier transform from the analytical calculation:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hftheor=0.5*sqrt(pi/alpha)*(exp(-pi**2/alpha*(f+f1)**2)+exp(-(pi)**2/alpha*(f-f1)**2))\n",
      "plot(f/f1,hftheor/T); xlim([-2*f1/f1,2*f1/f1]); xlabel(r'$f/f_1$'); ylabel(r'$\\widetilde{h}(f)/(h_0 T)$'); grid(True)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The plots look pretty similar, but we can plot one on top of the other to verify they look identical:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(f/f1,hftheor/T,'r--',linewidth='2',label=\"theor\");plot(f/f1,real(hf)/T,'b-',label=\"numer\");\n",
      "xlim([-2*f1/f1,2*f1/f1]); xlabel(r'$f/f_1$'); ylabel(r'$\\widetilde{h}(f)/(h_0 T)$'); legend(loc=\"upper center\"); grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}