{
 "metadata": {
  "name": "",
  "signature": "sha256:dbe8957aa0bb05f517cde33eaf63b12966e21a229ac09508d5f0b187237c3037"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Colored Noise\n",
      "\n",
      "In this problem we generate some white noise in the time domain, color it in the frequency domain, and then inverse Fourier transform it back into the frequency domain.  We also investigate the properties of the power spectral density."
     ]
    },
    {
     "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",
      "rcParams['figure.figsize'] = (10.0, 8.0)\n",
      "rcParams['font.size'] = 14\n",
      "from __future__ import division"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Data generation\n",
      "\n",
      "We generate $T=16\\,\\textrm{sec}$ of data sampled at $\\frac{1}{\\delta t}=1024\\,\\textrm{Hz}$, so that $N=16\\,384$.  We do this by generating an $N$-point series $\\{x_j\\}$ of random samples drawn from a distribution with zero mean and unit variance, i.e., $\\langle x_j\\rangle=0$ and $\\langle x_jx_l\\rangle=\\delta_{jl}$.  Before calling the random number generator, we seed it so that the results will be reproducible."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N=16384; dt_sec=1./1024.;\n",
      "T_sec=N*dt_sec;\n",
      "seed(20140204);\n",
      "x=randn(N);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## White noise in time domain\n",
      "\n",
      "Now we plot this in the time domain:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "t_sec = arange(N) * dt_sec\n",
      "plot(t_sec,x,'k.')\n",
      "title(\"White noise in the time domain\")\n",
      "xlabel(r'$t_j$ (sec)')\n",
      "ylabel(r'$x_j$')\n",
      "ymax = max(abs(array(ylim())))\n",
      "ylim([-ymax,ymax])\n",
      "yticks(arange(-ceil(ymax),ceil(ymax)+1))\n",
      "grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## White noise in the frequency domain\n",
      "\n",
      "Now we take the discrete Fourier transform $\\widehat{x}_k$ and plot $|\\widehat{x}_k|$ versus $f_k=k\\delta f=k/T$ for $f_k$ from $-512\\,\\textrm{Hz}$ to $512\\,\\textrm{Hz}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df_Hz = 1/T_sec; f_Hz = arange(-N/2,N/2) * df_Hz; fNy_Hz = 0.5/dt_sec\n",
      "xhat = fftshift(fft.fft(x))\n",
      "plot(f_Hz,abs(xhat),'k.')\n",
      "title('White noise in the frequency domain')\n",
      "xlabel(r'$f_k$ (Hz)')\n",
      "ylabel(r'$|\\hat{x}_k|$ (Hz$^{-1}$)')\n",
      "xlim([-fNy_Hz,fNy_Hz])\n",
      "xticks(128*arange(-4,5))\n",
      "grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Mean and standard deviation of the periodogram\n",
      "\n",
      "We calculate the mean and standard deviation of the periodogram values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pgram = abs(xhat)**2;\n",
      "print \"Observed mean of |xhat|^2 is %g\" % mean(pgram);\n",
      "print \"Observed standard deviation of |xhat|^2 is %g\" % std(pgram);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**QUESTION**: How does this compare to your theoretical expectations for $\\langle|\\widehat{x}_k|^2\\rangle$?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Colored noise in the frequency domain\n",
      "\n",
      "Next we \"color\" the noise in the frequency domain, by constructing $\\widehat{y}_k=\\widehat{x}_k\\,e^{-f_k^2/2\\sigma_f^2}$, where $\\sigma_f=16\\,\\textrm{Hz}$, for $k\\in[-N/2,N/2-1]$ and plot $|\\widehat{y}_k|$ versus $f_k$ for $f_k$ from $-512\\,\\textrm{Hz}$ to $512\\,\\textrm{Hz}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigmaf_Hz = 16.\n",
      "yhat = xhat * exp(-f_Hz**2/(2*sigmaf_Hz**2))\n",
      "plot(f_Hz,abs(yhat),'k.')\n",
      "title('Red noise in the frequency domain')\n",
      "xlabel(r'$f_k$ (Hz)')\n",
      "ylabel(r'$|\\hat{y}_k|$ (Hz$^{-1}$)')\n",
      "xlim([-fNy_Hz,fNy_Hz])\n",
      "xticks(128*arange(-4,5))\n",
      "grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is called \"red noise\" because it has more power at low frequencies\n",
      "\n",
      "## Colored noise in time domain\n",
      "\n",
      "Now we take the inverse discrete Fourier Transform to get $y_j$.  Note that our construction should give us data which is real to machine precision:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = ifft(fftshift(yhat)); print max(abs(imag(y)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We plot $y_j$ versus $t_j-t_0$ for $t_j-t_0$ from $0$ to $16\\,\\textrm{sec}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(t_sec,real(y),'k.'); title('Red noise in the time domain');\n",
      "xlabel(r'$t_j-t_0$ (sec)'); ylabel(r'$y_j$'); grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Power spectrum estimation\n",
      "\n",
      "Our plots of $|\\widehat{x}_k|^2$ and $|\\widehat{y}_k|^2$ are, up to constants, periodograms, which are estimates of the power spectral density.  We can use numpy's `psd()` function to get a more accurate estimate of the power spectral density.  First, check out the documentation:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "?psd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll use the defaults, except that we'll specify the sampling frequency `Fs`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psd(x,Fs=1/dt_sec); xlim([0,fNy_Hz]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And for the colored noise:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "psd(y,Fs=1/dt_sec); xlim([0,fNy_Hz]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We get a semilog plot by default, but we can instead save the values and compare them to our periodogram.  For simplicity we'll ask for two-sided psds.  We'll also call `mlab.psd()` which is the function under the hood of `pylab.psd()`, so we don't get the auto-generated plot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xPgram = (dt_sec/N)*abs(xhat)**2; (xPsd,fPsd_Hz) = mlab.psd(x,Fs=1/dt_sec,sides='twosided'); figure();\n",
      "plot(f_Hz,xPgram,'k.',color=(0.3,0.3,0.3),label='Periodogram');\n",
      "plot(fPsd_Hz,xPsd,'b-',label='Welch PSD estimate',linewidth=2);\n",
      "xlabel(r'$f_k$ (Hz)')\n",
      "ylabel(r'$S_x(f_k)$ (Hz$^{-1}$)')\n",
      "xlim([-fNy_Hz,fNy_Hz]); legend();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice the standard deviation of the PSD estimate using Welch's method is reduced relative to the periodogram:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'White noise periodogram has mean %e and standard deviation %e' % (mean(xPgram), std(xPgram));\n",
      "print \"White noise PSD estimated with Welch's method has mean %e and standard deviation %e\" % (mean(xPsd), std(xPsd))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we do for the same thing for the colored noise, and also include the theoretical PSD from our model:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "yPgram = (dt_sec/N)*abs(yhat)**2; (yPsd,fPsd_Hz) = mlab.psd(y,Fs=1/dt_sec,sides='twosided'); figure();\n",
      "yTheor = dt_sec * exp(-(f_Hz/sigmaf_Hz)**2)\n",
      "plot(f_Hz,yPgram,'k.',color=(0.3,0.3,0.3),label='Periodogram');\n",
      "plot(fPsd_Hz,yPsd,'b-',label='Welch PSD estimate',linewidth=2);\n",
      "plot(f_Hz,yTheor,'r-',label='Theory')\n",
      "title('PSD estimates for red noise');\n",
      "xlabel(r'$f_k$ (Hz)')\n",
      "ylabel(r'$S_y(f_k)$ (Hz$^{-1}$)');\n",
      "xlim([-5*sigmaf_Hz,5*sigmaf_Hz]); legend();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}