{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "def css_styling():\n", " styles = open(\"../styles/custom.css\", \"r\").read()\n", " return HTML(styles)\n", "css_styling()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "# [Physics 411](http://jklymak.github.io/Phy411/) Time Series Analysis\n", "*Jody Klymak*\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Week 4: More spectral properties and discrete spectra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\"Turbulence\n", "
\n", "\n", "As we saw in the [Introduction](Lecture-00-Intro-Python), it possible to have two time series that have completely different character, but the same univariate statistics. i.e. their underlying probability distribution densities can be exactly the same, but the ordering of the data matters. \n", "\n", "Here we define the meaning of a **Stationary Time Series** and discuss two of the most straight forward ways of characterizing the relation of that time series in time, the **Lag Covariance** and the **Spectral Density**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bandwidth considerations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An important property of spectra is that the broader the signal is in time, the narrower it is in spectral space, or vice-versa. This is somewhat non-intuitive, so lets look at a couple of examples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examples:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Ex 1** First, consider the power spectrum of a delta function in time $x(t)=\\delta(t-t_0)$. This is as localized a signal in time as possible. The Fourier Transform is:\n", "\n", "\\begin{align}\n", " X(f,T) &= \\int_{0}^{T} \\mathrm{e}^{-j2\\pi f t} \\delta(t-t_0)\\ \\mathrm{d}t\\\\\n", " &=\\mathrm{e}^{-j2\\pi f t_0}\n", "\\end{align}\n", "\n", "So the power spectrum is approximated by:\n", "\n", "\\begin{align}\n", " S_{xx}(f,T) &= \\frac{1}{T}X^*(f,T)X(f,T)\\\\\n", " &=\\frac{1}{T} \\mathrm{e}^{+j2\\pi f t_0}\\mathrm{e}^{-j2\\pi f t_0}\\\\\n", " &=\\frac{1}{T}\n", "\\end{align}\n", "\n", "This is just a constant for all $f$. So, a very narrow signal in $t$ leads to a very broad signal in $f$. \n", "\n", "**Ex 2** No consider a Gaussian in $t$ of standard deviation $b$:\n", "\n", "\\begin{equation}\n", " x(t)=\\exp\\left({\\frac{-t^2}{2b^2}}\\right)\n", "\\end{equation}\n", "\n", "Then the Fourier Transform is just:\n", "\n", "\\begin{align}\n", " X(f,T) &= \\int_{-T}^{T} \\mathrm{e}^{\\frac{-t^2}{2b^2}} \\mathrm{e}^{-j2\\pi f t}\\ \\mathrm{d}t\\\\\n", " &=\\sqrt{\\frac{\\pi}{2}}b\\ \\mathrm{e}^{-\\frac{4\\pi^2f^2b^2}{2}}\n", "\\end{align}\n", "\n", "or in terms of the spectral density:\n", "\n", "\\begin{align}\n", " S_{xx}(f) &= \\frac{\\pi b^2}{2}\\ \\mathrm{e}^{-4\\pi^2f^2b^2}\n", "\\end{align}\n", "\n", "which is a Gaussian with width $\\frac{1}{\\sqrt{8\\pi a}}$. Therefore, a narrow Gaussian signal in time will yield a wide Gaussian signal in frequency, and vice-versa." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Formal Proof:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, normalize $x(t)$ so that \n", "\n", "\\begin{equation}\n", " \\int x^2(t)\\ \\mathrm{d} t = \\int X^2(f)\\ \\mathrm{d} f = 1\n", "\\end{equation}\n", "\n", "And define the bandwidth of the signals in time and frequency as: \n", "\n", "\\begin{align}\n", " T_o^2 & \\equiv \\int_{-\\infty}^{\\infty} t^2 x^2(t)\\ \\mathrm{d}t\\\\\n", " B_o^2 & \\equiv \\int_{-\\infty}^{\\infty} f^2 X^2(f)\\ \\mathrm{d}f\\\\ \n", "\\end{align}\n", "\n", "Then $T_o^2B_o^2\\geq\\frac{1}{\\left(4\\pi\\right)^2}$, which says that as the badnwidth in time goes down, it has to go up in frequency, or vice-versa. \n", "\n", "The choice of the form of the bandwidth functions, $T_o$ and $B_o$ is somewhat arbitrary, but clearly both $x(t)$ and $X(f)$ have to fall off faster than $t^{-3}$ and $f^{-3}$ respectively for the bandwidths to be limited. \n", "\n", "We next note that $B_o$ is related to the integral of $\\dot{x}(t)=dx/dt$ squared:\n", "\n", "\\begin{align}\n", " \\left(2\\pi\\right)^2B_o^2 & = \\int_{-\\infty}^{\\infty} \\left(2\\pi f\\right)^2\\left|X(f)\\right|^2\\ \\mathrm{d}f\\\\\n", " & = \\int_{-\\infty}^{\\infty} \\left|FT (\\dot{x}) \\right|^2\\ \\mathrm{d}f\\\\\n", " & = \\int_{-\\infty}^{\\infty} \\dot{x}^2(t) \\mathrm{d}t\\\\ \n", "\\end{align}\n", "\n", "So, then we have $B_o$ in terms of a function of $x$ instead of $X$. We use the [Cauchy-Schwarz inequality](http://en.wikipedia.org/wiki/Cauchy–Schwarz_inequality) which says that \n", "\n", "\\begin{equation}\n", " \\int \\left| f(t) \\right|^2\\ \\mathrm{d}t \\int \\left| g(t) \\right|^2\\ \\mathrm{d}t \\geq \\left(\\int f\\,g\\ \\mathrm{d}t \\right)^2\n", "\\end{equation}\n", "\n", "to argue that:\n", "\n", "\\begin{align}\n", " T_o B_o & = \\frac{1}{2\\pi} \\left(\\int_{-\\infty}^{\\infty} t^2x^2\\ \\mathrm{d}t \\right)^{1/2} \\left( \\int_{-\\infty}^{\\infty} \\dot{x}^2\\ \\mathrm{d}t\\right)^{1/2} \\\\\n", " & \\geq \\frac{1}{2\\pi} \\int_{-\\infty}^{\\infty} t x \\dot{x}\\ \\mathrm{d}t\\\\\n", " & \\geq \\frac{1}{4\\pi} \\int x^2\\ \\mathrm{d}t\\\\\n", " & \\geq \\frac{1}{4\\pi}\n", "\\end{align}\n", "\n", "where the third step follows from integration by parts, and the fourth from the fact that we normalized $x(t)$ so that the integral here was $1$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cross spectrum and Coherence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far we have shown a few results with the cross spectrum $G_{xy}(f)$ derived from the cross lag-correlation:\n", "\n", "\\begin{equation}\n", " G_{xy}(f)\\equiv 2 \\int_0^{\\infty}R_{xy}(\\tau)\\ \\mathrm{e}^{-j2\\pi f\\tau}\\ \\mathrm{d}\\tau\n", "\\end{equation}\n", "\n", "This defines the frequency dependence of the relations ship between $x(t)$ and $y(t)$. \n", "\n", "This is more powerful than a straight linear regression between the two time series, in that even quite burried signals can be seen to be significantly correlated at certain frequencies. Consider the following examples where there is a sine wave with frequency of $f=0.033\\ \\mathrm{Hz}$ burried in each signal that is completely swamped by noise:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Close figure', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import mlab\n", "Gxy,f = mlab.csd(x,y,NFFT=400,Fs=1.) # Don't use mlab.csd until you prove it works!\n", "fig,ax=plt.subplots(1,1)\n", "ax.loglog(f,np.abs(Gxy))\n", "ax.set_xlabel('f [Hz]');ax.set_ylabel('$|G_{xy}(f)| [V^2 Hz^{-1}]$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We find a strong peak at the frequency of the sine waves, which indicates the two waves are correlated. \n", "\n", "The functions we input had a $0.6$ radian phase shift, and the cross-sepctrum gives information about that as well. Since $G_{xy}$ is complex, the phase is given by $Im(G_{xy})/Re(G_{xy}) = \\tan \\phi$. Looking at our estimate, we see that" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.71326976]\n" ] } ], "source": [ "ind = np.where(abs(Gxy)>1.e2)[0]\n", "\n", "pha = np.arctan2(np.imag(Gxy[ind]),np.real(Gxy[ind]))\n", "print pha" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we see that the phase shift is relatively well captured. \n", "\n", "More formally. If $G_{xy}(f)$ is complex, then we can split into the co-incident and quadrature parts of the cross-spectrum:\n", "\n", "\\begin{equation}\n", " G_{xy}(f)=C_{xy}(f)-jQ_{xy}(f)\n", "\\end{equation}\n", "\n", "or in terms of amplitude and phase:\n", "\n", "\\begin{equation}\n", " G_{xy}(f)=|G_{xy}(f)|\\mathrm{e}^{-j\\theta_{xy}(f)}\n", "\\end{equation}\n", "\n", "where $\\tan\\left(\\theta_{xy}(f)\\right)=\\frac{Q_{xy}(f)}{C_{xy}(f)}$.\n", "\n", "Like for the Correlation Coefficient, we define a **Coherence Squared Function** by normalizing by the spectral desnity of our two signals:\n", "\n", "\\begin{equation}\n", " \\gamma_{xy}^2 \\equiv \\frac{\\left|G_{xy}(f)\\right|}{G_{xx}(f)G_{yy}(f)} \\leq 1\n", "\\end{equation}\n", "\n", "We can prove the inequality by the same quadratic trick we used for the Correlation Co-efficient. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Time delay" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As another example, consider a function $y(t)$ that is a noisy and delayed version of $x(t)$, i.e.\n", "\n", "\\begin{equation}\n", " y(t)=a \\, x(t-\\tau_0) + n(t)\n", "\\end{equation}\n", "\n", "Where $n(t)$ is noise that is not correlated with $x(t)$. Then we can calculate $R_{xy}(\\tau)$:\n", "\n", "\\begin{align}\n", " R_{xy}(\\tau)&=E\\left[ x(t)y(t+\\tau) \\right]\\\\\n", " &=E\\left[ a x(t)x(t-\\tau_0+\\tau) +x(t)n(t)\\right]\\\\\n", " &=aE\\left[x(t)x(t-\\tau_0+\\tau)\\right]\\\\\n", " &=a\\,R_{xx}(\\tau-\\tau_0)\n", "\\end{align}\n", "\n", "So by a simple change of variables:\n", "\n", "\\begin{align}\n", " G_{xy}(f)=a\\, G_{xx}(f)\\ \\mathrm{e}^{-j2\\pi f\\tau}\n", "\\end{align}\n", "\n", "and we see that \n", "\n", "\\begin{align}\n", " \\left|G_{xy}(f)\\right|=a\\, G_{xx}(f)\n", "\\end{align}\n", "\n", "and\n", "\n", "\\begin{align}\n", " \\theta_{xy}(f)&=\\arctan \\left(\\frac{-\\sin(2\\pi f\\tau_0)}{\\cos(2\\pi f\\tau_0}\\right)\\\\\n", " & = 2\\pi f\\tau_0\n", "\\end{align}\n", "\n", "We also note that the coherence squared is\n", "\n", "\\begin{align}\n", " \\gamma_{xy}^2 &= \\frac{\\left| G_{xy}\\right|^2}{G_{xx}G_{yy}}\\\\\n", " &=\\frac{\\left|G_{xx}\\right|^2}{G_{xx}^2+G_{xx}G_{nn}}\n", "\\end{align}\n", "\n", "So as the noise gets smaller, this approaches 1, whereas a stronger noise signal reduces the coherence.\n", "\n", "Here we test this numerically, with $a=1$, and a delay of $\\tau_0=10\\ \\mathrm{s}$:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Close figure', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/jklymak/anaconda/lib/python2.7/site-packages/matplotlib/mathtext.py:860: MathTextWarning: Substituting with a symbol from Computer Modern.\n", " MathTextWarning)\n" ] } ], "source": [ "Gxy,f=mlab.csd(x,y,Fs=1,NFFT=200)\n", "Gxx,f=mlab.psd(x,Fs=1,NFFT=200)\n", "Gyy,f=mlab.psd(y,Fs=1,NFFT=200)\n", "fig,ax=plt.subplots(3,1,sharex=True,figsize=(5,9))\n", "ax[0].loglog(f,Gxx,label='x')\n", "ax[0].loglog(f,Gyy,label='y: delayed')\n", "ax[0].legend(loc=3)\n", "ax[0].set_ylabel(r'$G_{xx} [V^2 Hz^{-1}]$')\n", "ax[1].semilogx(f,np.abs(Gxy)**2/Gxx/Gyy)\n", "ax[1].set_ylabel(r'$\\gamma_{xy} $')\n", "ax[2].semilogx(f,np.arctan2(-np.imag(Gxy),np.real(Gxy)))\n", "ax[2].semilogx(f,np.mod(2*np.pi*f*10.+np.pi,2*np.pi)-np.pi,'k--')\n", "ax[2].set_ylabel(r'$\\theta_{xy} [rad]$')\n", "ax[2].set_xlabel('f [Hz]')\n", "ax[2].set_yticks(np.arange(-np.pi,np.pi,np.pi/4))\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, here we see a fair bit of information. The spectra of the two time series look the same until the noise spectrum starts to dominate at high frequencies. This leads to a corresponding drop in the coherence squared at high frequencies (second panel). The time offset is readily apparent in the third panel (note the wrapping at $\\pm\\pi$). The dashed line is what we would expect from a 10-s phase lag between $x$ and $y$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discrete Power Spectra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course in the plots I've made above, I have been using discrete power spectra, but there are a number of caveats to be kept in mind when dealing with discrete data. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Discrete time series" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we have a real time series $x(t)$, then in real life we discretely sample it, lets say every $\\Delta t$ seconds to get a series of data at \n", "\n", "\\begin{align}\n", " t_n &\\equiv n\\,\\Delta t\\\\\n", " x_n &\\equiv x(t_n)\n", "\\end{align}\n", "\n", "Sampling at even intervals is completely a convenience, and not all time series have this luxury! However, you can usually bin or interpolate your data to get it onto a regular grid without losing too much fidelity.\n", "\n", "We also can only sample for a finite amount of time, say $T$ seconds, and therefore our time series will be of finite length:\n", "\n", "\\begin{align}\n", " T=N\\, \\Delta t\n", "\\end{align}\n", "\n", "We define $f_S =\\frac{1}{\\Delta t}$ as the **sampling frequency**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Digital Fourier Transform" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Formally the Fourier Transform of $x(t)$ is \n", "\n", "\\begin{equation}\n", " X(f)=\\lim_{T\\to\\infty} \\int_0^{T} x(t)\\ \\mathrm{e}^{-j2\\pi f t}\\ \\mathrm{d}t\n", "\\end{equation}\n", "\n", "The digital approximation is simply\n", "\n", "\\begin{align}\n", " X(f_k)&=\\sum_{n=1}^{N} x_n\\ \\mathrm{e}^{-j2\\pi f_k t_n}\\ \\Delta t\\\\\n", " &=\\frac{T}{N}\\sum_{n=1}^{N} x_n\\ \\mathrm{e}^{\\frac{-j2\\pi f_k n T}{N}}\n", "\\end{align}\n", "\n", "Now, if we are smart about how we choose our values of $f_k$ we can efficiently decompose $x_n$. To do so, we choose: \n", "\\begin{align}\n", " f_k &= \\frac{k}{T} && \\text{for } k=0,1,...,N/2\n", "\\end{align}\n", "\n", "then\n", "\n", "\\begin{align}\n", " X(f_k)&=\\frac{T}{N}\\sum_{n=1}^{N} x_n\\ \\mathrm{e}^{\\frac{-j2\\pi k n }{N}}\\\\\n", " &= \\frac{T}{N}\\left(\\sum_{n=1}^{N} x_n\\ \\cos\\left(\\frac{2\\pi k n }{N}\\right) -j \\sum_{n=1}^{N} x_n\\ \\sin\\left(\\frac{2\\pi k n }{N}\\right)\\right)\\\\\n", " &= \\frac{T}{N}\\sum_{n=1}^{N} x_n\\ C_{k}(n) - j\\frac{T}{N}\\sum_{n=1}^{N} x_n\\ S_{k}(n)\n", "\\end{align}\n", "\n", "There are a couple of edge cases to consider in this. For $k=0$ and $N/2$, $S_k(n)=0$, so we can ignore those. Also note that $C_0(n)=1$ and $f_0=0$, so \n", "\n", "\\begin{align}\n", " X(f_0)&=\\frac{T}{N}\\sum_{n=1}^{N} x_n\\\\\n", " &= T\\ \\overline{x}\n", "\\end{align}\n", "\n", "represents the *mean* of $x_m$ (times $T$). \n", "\n", "Why do we only use $N/2$ frequencies? There are two reasons. First, there are $N/2-1$ values for $S_k(n)$, and there are $N/2+1$ values for $C_k(n)$, so in total this choice of $f_k$ yields us $N$ Fourier co-efficients. From an information point of view, if each of those $N$ co-efficients is *independent* then there is as much information in the Fourier decomposition as there is in the original time series, i.e. our tranform is lossless.\n", "\n", "Second, the highest frequency $f_{N/2}=\\frac{1}{2\\Delta t}$ is a sine wave that will only have two samples in it. That is the minimum needed to define the sine wave, and therefore is the highest frequency we can resolve. $f_{N/2}=\\frac{1}{2}f_S$ is fundamental to time-series analysis and is called the **Nyquist Frequency**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sin and Cos form a basis that spans the time series" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is possible to choose any $f$ to calculate a Fourier Component, but we use the discrete $f_k$ as defined because the resulting $C_k(n)$ and $S_k(n)$ are vectors that as a set form an orthogonal basis that spans the space $x_n \\in \\mathbb{R}^N$. It is relatively straightforward to show that the basis vectors are orthogonal:\n", "\n", "\\begin{align}\n", " \\mathbf{C}_k\\cdot \\mathbf{C}_l & = \\sum_{n=1}^N \\cos\\left(\\frac{2\\pi k n}{N} \\right) \\cos\\left(\\frac{2\\pi l n}{N} \\right)\\\\\n", " &= \\sum_{n=1}^N \\frac{1}{2} \\left(\n", " \\cos\\left(\\frac{2\\pi n}{N} (k+l)\\right) + \n", " \\cos\\left(\\frac{2\\pi n}{N} (k-l)\\right) \n", " \\right)\\\\\n", " &= \\frac{N}{2}\\, \\delta(k-l) \\\\\n", " \\mathbf{S}_k\\cdot \\mathbf{S}_l &= \\frac{N}{2}\\,\\delta(k-l)\\\\\n", " \\mathbf{C}_k\\cdot \\mathbf{S}_l &= 0\\\\\n", "\\end{align}\n", "\n", "This then lets us write $x_n$ in a matrix equation:\n", "\n", "\\begin{align}\n", " \\mathbf{x}& = \\frac{2}{N}\\mathrm{F} \\mathbf{X^0}\n", "\\end{align}\n", "\n", "where the columns of the matrix $\\mathrm{F}$ are \n", "\n", "\\begin{align}\n", " \\mathrm{F}&=[ \\mathbf{C}_{0}, \\mathbf{C}_{1}, ... \\mathbf{C}_{N/2}, \\mathbf{S}_1,\\mathbf{S}_2,...\\mathbf{S}_{N/2-1}]\n", "\\end{align}\n", "\n", "Note that $\\mathrm{F}$ is $NxN$. Since the vectors in the matrix $\\mathrm{F}$ are independent the matrix is invertable, and we can solve the matrix equation to get the Fourier co-efficients:\n", "\n", "\\begin{align}\n", " \\mathbf{X^0}& = \\frac{N}{2}\\mathrm{F}^{-1} \\mathbf{x}\n", "\\end{align}\n", "\n", "The advantage of thinking about these as matrices are three-fold. First, there are a number of computer algorithms for inverting matrices that are quite efficient compared to perfoming the sums by hand. Second, once you invert $\\mathrm{F}$ once, it can be applied to other vectors $\\mathbf{x}$ that you may need to calculate the Fourier transform of. Third, sines and cosines are not the only possible orthogonal basis of $\\mathbb{R}^n$, so this technique may be more applicable. For instance there are cases where the waves that propagate in time (or space) may have a different (known) form than a sine wave. For instance if you have a wave propagating through an inhomogenous media, the sine wave will be distorted, but there may exist a basis function to decompose your signal into.\n", "\n", "To get the Fourier transform from the product of the inversion, we form complex numbers out of the appropriate entries in the matrix, and we apply the proper normalization:\n", "\n", "\\begin{align}\n", " X_k& = \\frac{T}{N}X^0_k-j\\frac{T}{N} X^0_{k+N/2} && \\text{for } k=1,2,...,N/2-1\\\\\n", " X_0&=\\frac{T}{N}X^0_0 && \\text{for } k=0\\\\\n", " X_{N/2}&=\\frac{T}{N}X^0_{N/2} && \\text{for } k=N/2\\\\ \n", "\\end{align}\n", "\n", "giving $N/2+1$ co-efficients, $N/2-1$ of them imaginary, the other two real. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using computing packages to compute the Fourier Transform" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Many computing packages provide routines to calculate the Fourier transform, often using a very quick algorithm called the *Fast Fourier Transform*. In general, you can use this algorithm. In python it is called as:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(400,)\n", "[ 6939.81683831 +0.j 337.56386658+2293.61434898j\n", " 38.43233748+2091.71264929j 287.73663901 +775.72214725j\n", " 50.60860807+1047.02881806j 182.44177424 +843.38411605j]\n", "[ 182.44177424 -843.38411605j 50.60860807-1047.02881806j\n", " 287.73663901 -775.72214725j 38.43233748-2091.71264929j\n", " 337.56386658-2293.61434898j]\n", "X[N/2]= -3.704e+01 +4.086e-14i \n", "sum x: 6939.817\n" ] } ], "source": [ "x = np.cumsum(np.random.randn(400))\n", "X = np.fft.fft(x)\n", "print np.shape(X)\n", "print X[:6]\n", "print X[-5:]\n", "print 'X[N/2]= %1.3e %+1.3ei '%(X[200].real,X[200].imag)\n", "print 'sum x: %1.3f'%np.sum(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, what happened here? Its a good idea to check the [manual](http://docs.scipy.org/doc/numpy/reference/routines.fft.html), but their implimentation of the FFT returns $N$ entries from $k=0...N-1$. However, above, I said that we only go from $k=0...N/2$. However, look at the first few entries for $X$ and the last few. They are mirror images of one another, and therefore this implimentation of FFT has returned extraneous information.\n", "\n", "You can also note that $X_{N/2}$ is (to machine precision) real, as we expected from above. \n", "\n", "Finally note that $X_0=\\sum x_n$, also as expected. \n", "\n", "There is usually an inverse FFT that we call as below. Note that it returns the correct answer, again within machine precision." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-0.55145917 +1.73638881e-15j -1.66816803 -1.34559031e-15j\n", " -0.99621976 -5.19154609e-15j -1.37562181 +4.10782519e-15j\n", " -2.00659662 +1.93045580e-14j -3.38316970 -1.07425180e-14j\n", " -2.60799393 +3.33525319e-15j -3.44388475 -2.11874962e-14j\n", " -4.33911080 +2.20223839e-14j -5.61199721 +9.45465928e-15j]\n", "[-0.55145917 -1.66816803 -0.99621976 -1.37562181 -2.00659662 -3.3831697\n", " -2.60799393 -3.44388475 -4.3391108 -5.61199721]\n" ] } ], "source": [ "xinv=np.fft.ifft(X)\n", "print xinv[:10]\n", "print x[:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note also that `fft` does *not* include the proper units in front of $\\Delta t=\\frac{T}{N}$, so if your time series is sampled at a different rate than $\\Delta t=1$, and you use `fft`, you need to add the appropriate normalization. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discrete Power Spectra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We rarely consider just the Fourier transform, but rather the power spectra or cross-spectra, $G_{xx}(f)$ and $G_{xy}(f)$ Recall that \n", "\n", "\\begin{equation}\n", " G_{xx}(f) = 2 \\lim_{T\\to\\infty} \\frac{1}{T} E\\left[ \\left| X(f,T) \\right|^2 \\right]\n", "\\end{equation}\n", "\n", "so it is natural to approximate the discrete spectrum using the discrete Fourier Transform:\n", "\n", "\\begin{equation}\n", " \\hat{G}_{xx}(f_k) = \\frac{2}{T} X_k^*X_k\n", "\\end{equation}\n", "\n", "Note again the units: $X_k$ has units of $\\mathrm{V^2\\,s^2}$, so $\\hat{G}_{xx}(f_k)$ has units of \n", "$\\mathrm{V^2\\,s}=\\mathrm{V^2\\,Hz^{-1}}$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Statistics of the spectral estimate?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimate of $\\hat{G}_{xx}(f_k)$ is a random variable, and of course has a probability distribution that we can estimate. Note that \n", "\n", "\\begin{equation}\n", " X_k^*X_k = \\mathbb{R}\\left(X_k\\right)^2+\\mathbb{I}\\left(X_k\\right)^2\n", "\\end{equation}\n", "\n", "If we assume that $\\mathbb{R}(X_k)$ and $\\mathbb{I}(X_k)$ are random independent variables, each with a normal distribution and a varaince estimate of $\\sigma^2=\\hat{G}_{xx}(f_k)/2$ (this not always true, but as we see below we will end up averaging a number of estimates together, which due to the Central Limit Theorem will make this closer to true in all cases), then $(X_k^*X_k)/\\sigma^2$ follows a [$\\chi^2_n$ distribution](http://en.wikipedia.org/wiki/Chi-squared_distribution) with $n=2$ degrees of freedom.\n", "\n", "If we want the 95% confidence interval that denotes where we expect the true value of $G_{xx}(f)$ to be then:\n", "\n", "\\begin{equation}\n", " \\frac{\\left. \\chi^2_n \\right|_{0.025}}{n} \\leq \\frac{\\hat{G}_{xx}/2}{G_{xx}/2} \\leq \\frac{\\left. \\chi^2_n \\right|_{0.975}}{n}\n", "\\end{equation}\n", "\n", "or in a manner that is what we might plot as error bars on a plot:\n", "\n", "\\begin{equation}\n", " \\frac{n\\hat{G}_{xx}/2}{\\left. \\chi^2_n \\right|_{0.975}} \\leq {G_{xx}/2} \\leq \\frac{n\\hat{G}_{xx}/2}{\\left. \\chi^2_n \\right|_{0.025}}\n", "\\end{equation}\n", "\n", "where again, $\\hat{G}_{kk}$ is our best guess at the variance of $\\mathbb{R}(X_k)^2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example error estimate:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('