{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 6. Time series and data smoothing\n",
"\n",
"**[Dataset download](https://s3.amazonaws.com/bebi103.caltech.edu/data/retina_trace.zip)**\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"\n",
"(function(root) {\n",
" function now() {\n",
" return new Date();\n",
" }\n",
"\n",
" var force = true;\n",
"\n",
" if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
" root._bokeh_onload_callbacks = [];\n",
" root._bokeh_is_loading = undefined;\n",
" }\n",
"\n",
" var JS_MIME_TYPE = 'application/javascript';\n",
" var HTML_MIME_TYPE = 'text/html';\n",
" var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
" var CLASS_NAME = 'output_bokeh rendered_html';\n",
"\n",
" /**\n",
" * Render data to the DOM node\n",
" */\n",
" function render(props, node) {\n",
" var script = document.createElement(\"script\");\n",
" node.appendChild(script);\n",
" }\n",
"\n",
" /**\n",
" * Handle when an output is cleared or removed\n",
" */\n",
" function handleClearOutput(event, handle) {\n",
" var cell = handle.cell;\n",
"\n",
" var id = cell.output_area._bokeh_element_id;\n",
" var server_id = cell.output_area._bokeh_server_id;\n",
" // Clean up Bokeh references\n",
" if (id != null && id in Bokeh.index) {\n",
" Bokeh.index[id].model.document.clear();\n",
" delete Bokeh.index[id];\n",
" }\n",
"\n",
" if (server_id !== undefined) {\n",
" // Clean up Bokeh references\n",
" var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
" cell.notebook.kernel.execute(cmd, {\n",
" iopub: {\n",
" output: function(msg) {\n",
" var id = msg.content.text.trim();\n",
" if (id in Bokeh.index) {\n",
" Bokeh.index[id].model.document.clear();\n",
" delete Bokeh.index[id];\n",
" }\n",
" }\n",
" }\n",
" });\n",
" // Destroy server and session\n",
" var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
" cell.notebook.kernel.execute(cmd);\n",
" }\n",
" }\n",
"\n",
" /**\n",
" * Handle when a new output is added\n",
" */\n",
" function handleAddOutput(event, handle) {\n",
" var output_area = handle.output_area;\n",
" var output = handle.output;\n",
"\n",
" // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
" if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
" return\n",
" }\n",
"\n",
" var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
"\n",
" if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
" toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
" // store reference to embed id on output_area\n",
" output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
" }\n",
" if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
" var bk_div = document.createElement(\"div\");\n",
" bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
" var script_attrs = bk_div.children[0].attributes;\n",
" for (var i = 0; i < script_attrs.length; i++) {\n",
" toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
" }\n",
" // store reference to server id on output_area\n",
" output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
" }\n",
" }\n",
"\n",
" function register_renderer(events, OutputArea) {\n",
"\n",
" function append_mime(data, metadata, element) {\n",
" // create a DOM node to render to\n",
" var toinsert = this.create_output_subarea(\n",
" metadata,\n",
" CLASS_NAME,\n",
" EXEC_MIME_TYPE\n",
" );\n",
" this.keyboard_manager.register_events(toinsert);\n",
" // Render to node\n",
" var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
" render(props, toinsert[toinsert.length - 1]);\n",
" element.append(toinsert);\n",
" return toinsert\n",
" }\n",
"\n",
" /* Handle when an output is cleared or removed */\n",
" events.on('clear_output.CodeCell', handleClearOutput);\n",
" events.on('delete.Cell', handleClearOutput);\n",
"\n",
" /* Handle when a new output is added */\n",
" events.on('output_added.OutputArea', handleAddOutput);\n",
"\n",
" /**\n",
" * Register the mime type and append_mime function with output_area\n",
" */\n",
" OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
" /* Is output safe? */\n",
" safe: true,\n",
" /* Index of renderer in `output_area.display_order` */\n",
" index: 0\n",
" });\n",
" }\n",
"\n",
" // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
" if (root.Jupyter !== undefined) {\n",
" var events = require('base/js/events');\n",
" var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
"\n",
" if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
" register_renderer(events, OutputArea);\n",
" }\n",
" }\n",
"\n",
" \n",
" if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
" root._bokeh_timeout = Date.now() + 5000;\n",
" root._bokeh_failed_load = false;\n",
" }\n",
"\n",
" var NB_LOAD_WARNING = {'data': {'text/html':\n",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
" \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"
re-rerun `output_notebook()` to attempt to load from CDN again, or
\"}};\n",
"\n",
" function display_loaded() {\n",
" var el = document.getElementById(null);\n",
" if (el != null) {\n",
" el.textContent = \"BokehJS is loading...\";\n",
" }\n",
" if (root.Bokeh !== undefined) {\n",
" if (el != null) {\n",
" el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
" }\n",
" } else if (Date.now() < root._bokeh_timeout) {\n",
" setTimeout(display_loaded, 100)\n",
" }\n",
" }\n",
"\n",
"\n",
" function run_callbacks() {\n",
" try {\n",
" root._bokeh_onload_callbacks.forEach(function(callback) {\n",
" if (callback != null)\n",
" callback();\n",
" });\n",
" } finally {\n",
" delete root._bokeh_onload_callbacks\n",
" }\n",
" console.debug(\"Bokeh: all callbacks have finished\");\n",
" }\n",
"\n",
" function load_libs(css_urls, js_urls, callback) {\n",
" if (css_urls == null) css_urls = [];\n",
" if (js_urls == null) js_urls = [];\n",
"\n",
" root._bokeh_onload_callbacks.push(callback);\n",
" if (root._bokeh_is_loading > 0) {\n",
" console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
" return null;\n",
" }\n",
" if (js_urls == null || js_urls.length === 0) {\n",
" run_callbacks();\n",
" return null;\n",
" }\n",
" console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
" root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
"\n",
" function on_load() {\n",
" root._bokeh_is_loading--;\n",
" if (root._bokeh_is_loading === 0) {\n",
" console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
" run_callbacks()\n",
" }\n",
" }\n",
"\n",
" function on_error() {\n",
" console.error(\"failed to load \" + url);\n",
" }\n",
"\n",
" for (var i = 0; i < css_urls.length; i++) {\n",
" var url = css_urls[i];\n",
" const element = document.createElement(\"link\");\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\n",
" element.rel = \"stylesheet\";\n",
" element.type = \"text/css\";\n",
" element.href = url;\n",
" console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
" document.body.appendChild(element);\n",
" }\n",
"\n",
" for (var i = 0; i < js_urls.length; i++) {\n",
" var url = js_urls[i];\n",
" var element = document.createElement('script');\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\n",
" element.async = false;\n",
" element.src = url;\n",
" console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
" document.head.appendChild(element);\n",
" }\n",
" };\n",
"\n",
" function inject_raw_css(css) {\n",
" const element = document.createElement(\"style\");\n",
" element.appendChild(document.createTextNode(css));\n",
" document.body.appendChild(element);\n",
" }\n",
"\n",
" \n",
" var js_urls = [];\n",
" var css_urls = [];\n",
" \n",
"\n",
" var inline_js = [\n",
" function(Bokeh) {\n",
" /* BEGIN bokeh.min.js */\n",
" /*!\n",
" * Copyright (c) 2012 - 2019, Anaconda, Inc., and Bokeh Contributors\n",
" * All rights reserved.\n",
" * \n",
" * Redistribution and use in source and binary forms, with or without modification,\n",
" * are permitted provided that the following conditions are met:\n",
" * \n",
" * Redistributions of source code must retain the above copyright notice,\n",
" * this list of conditions and the following disclaimer.\n",
" * \n",
" * Redistributions in binary form must reproduce the above copyright notice,\n",
" * this list of conditions and the following disclaimer in the documentation\n",
" * and/or other materials provided with the distribution.\n",
" * \n",
" * Neither the name of Anaconda nor the names of any contributors\n",
" * may be used to endorse or promote products derived from this software\n",
" * without specific prior written permission.\n",
" * \n",
" * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n",
" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n",
" * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\n",
" * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE\n",
" * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR\n",
" * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF\n",
" * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS\n",
" * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN\n",
" * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)\n",
" * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF\n",
" * THE POSSIBILITY OF SUCH DAMAGE.\n",
" */\n",
" (function(root, factory) {\n",
" root[\"Bokeh\"] = factory();\n",
" })(this, function() {\n",
" var define;\n",
" var parent_require = typeof require === \"function\" && require\n",
" return (function(modules, entry, aliases, externals) {\n",
" if (aliases === undefined) aliases = {};\n",
" if (externals === undefined) externals = {};\n",
"\n",
" var cache = {};\n",
"\n",
" var normalize = function(name) {\n",
" if (typeof name === \"number\")\n",
" return name;\n",
"\n",
" if (name === \"bokehjs\")\n",
" return entry;\n",
"\n",
" var prefix = \"@bokehjs/\"\n",
" if (name.slice(0, prefix.length) === prefix)\n",
" name = name.slice(prefix.length)\n",
"\n",
" var alias = aliases[name]\n",
" if (alias != null)\n",
" return alias;\n",
"\n",
" var trailing = name.length > 0 && name[name.lenght-1] === \"/\";\n",
" var index = aliases[name + (trailing ? \"\" : \"/\") + \"index\"];\n",
" if (index != null)\n",
" return index;\n",
"\n",
" return name;\n",
" }\n",
"\n",
" var require = function(name) {\n",
" var mod = cache[name];\n",
" if (!mod) {\n",
" var id = normalize(name);\n",
"\n",
" mod = cache[id];\n",
" if (!mod) {\n",
" if (!modules[id]) {\n",
" if (parent_require && externals[id]) {\n",
" try {\n",
" mod = {exports: parent_require(id)};\n",
" cache[id] = cache[name] = mod;\n",
" return mod.exports;\n",
" } catch (e) {}\n",
" }\n",
"\n",
" var err = new Error(\"Cannot find module '\" + name + \"'\");\n",
" err.code = 'MODULE_NOT_FOUND';\n",
" throw err;\n",
" }\n",
"\n",
" mod = {exports: {}};\n",
" cache[id] = cache[name] = mod;\n",
" modules[id].call(mod.exports, require, mod, mod.exports);\n",
" } else\n",
" cache[name] = mod;\n",
" }\n",
"\n",
" return mod.exports;\n",
" }\n",
"\n",
" var main = require(entry);\n",
" main.require = require;\n",
"\n",
" main.register_plugin = function(plugin_modules, plugin_entry, plugin_aliases, plugin_externals) {\n",
" if (plugin_aliases === undefined) plugin_aliases = {};\n",
" if (plugin_externals === undefined) plugin_externals = {};\n",
"\n",
" for (var name in plugin_modules) {\n",
" modules[name] = plugin_modules[name];\n",
" }\n",
"\n",
" for (var name in plugin_aliases) {\n",
" aliases[name] = plugin_aliases[name];\n",
" }\n",
"\n",
" for (var name in plugin_externals) {\n",
" externals[name] = plugin_externals[name];\n",
" }\n",
"\n",
" var plugin = require(plugin_entry);\n",
"\n",
" for (var name in plugin) {\n",
" main[name] = plugin[name];\n",
" }\n",
"\n",
" return plugin;\n",
" }\n",
"\n",
" return main;\n",
" })\n",
" ([\n",
" function _(n,o,r){n(1),function(n){for(var o in n)r.hasOwnProperty(o)||(r[o]=n[o])}(n(102))},\n",
" function _(n,c,f){n(2),n(11),n(14),n(21),n(49),n(52),n(87),n(94),n(100)},\n",
" function _(e,n,a){e(3)()||Object.defineProperty(Object,\"assign\",{value:e(4),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(r,t,o){t.exports=function(){var r,t=Object.assign;return\"function\"==typeof t&&(t(r={foo:\"raz\"},{bar:\"dwa\"},{trzy:\"trzy\"}),r.foo+r.bar+r.trzy===\"razdwatrzy\")}},\n",
" function _(t,r,n){var o=t(5),c=t(10),a=Math.max;r.exports=function(t,r){var n,f,h,i=a(arguments.length,2);for(t=Object(c(t)),h=function(o){try{t[o]=r[o]}catch(t){n||(n=t)}},f=1;f= 0\");if(!isFinite(r))throw new RangeError(\"Count must be < ∞\");for(n=\"\";r;)r%2&&(n+=t),r>1&&(t+=t),r>>=1;return n}},\n",
" function _(t,i,n){var r=t(18),a=Math.abs,o=Math.floor;i.exports=function(t){return isNaN(t)?0:0!==(t=Number(t))&&isFinite(t)?r(t)*o(a(t)):t}},\n",
" function _(n,t,i){t.exports=n(19)()?Math.sign:n(20)},\n",
" function _(n,t,o){t.exports=function(){var n=Math.sign;return\"function\"==typeof n&&(1===n(10)&&-1===n(-20))}},\n",
" function _(n,r,t){r.exports=function(n){return n=Number(n),isNaN(n)||0===n?n:n>0?1:-1}},\n",
" function _(e,r,a){e(22)()||Object.defineProperty(Array,\"from\",{value:e(23),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(n,o,r){o.exports=function(){var n,o,r=Array.from;return\"function\"==typeof r&&(o=r(n=[\"raz\",\"dwa\"]),Boolean(o&&o!==n&&\"dwa\"===o[1]))}},\n",
" function _(e,l,r){var n=e(24).iterator,t=e(44),a=e(45),i=e(46),u=e(47),o=e(10),f=e(8),c=e(48),v=Array.isArray,h=Function.prototype.call,y={configurable:!0,enumerable:!0,writable:!0,value:null},s=Object.defineProperty;l.exports=function(e){var l,r,A,g,p,w,b,d,x,j,O=arguments[1],m=arguments[2];if(e=Object(o(e)),f(O)&&u(O),this&&this!==Array&&a(this))l=this;else{if(!O){if(t(e))return 1!==(p=e.length)?Array.apply(null,e):((g=new Array(1))[0]=e[0],g);if(v(e)){for(g=new Array(p=e.length),r=0;r
\n",
""
],
"text/plain": [
":Curve [t (s)] (V (µV))"
]
},
"execution_count": 9,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1116"
}
},
"output_type": "execute_result"
}
],
"source": [
"hv.Curve(\n",
" data=df.iloc[10000:30000, :],\n",
" kdims=['t (s)'],\n",
" vdims=['V (µV)']\n",
").opts(\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" frame_height=300,\n",
" frame_width=625,\n",
" show_grid=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the zoomed out view, we can see spiking events as obvious drops in the voltage. When we zoom in, we can see what an individual spike looks like. We'll reproduce it here."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Curve [t (s)] (V (µV))"
]
},
"execution_count": 10,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1226"
}
},
"output_type": "execute_result"
}
],
"source": [
"inds = (0.675 < df['t (s)']) & (df['t (s)'] < 0.69)\n",
"\n",
"hv.Curve(\n",
" data=df.loc[inds, :],\n",
" kdims=['t (s)'],\n",
" vdims=['V (µV)']\n",
").opts(\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" frame_height=300,\n",
" frame_width=625,\n",
" show_grid=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The signal in the frequency domain\n",
"\n",
"We are used to looking at time series data in the time domain, but this is only one way to look at it. Another way is in the **frequency domain**. We will not go into the mathematical details here, but talk very loosely to get a feel for what it means to represent data in the frequency domain. You can imagine that a signal can be represented by the sum of sine waves of different frequency. High frequency aspects of the signal oscillate rapidly, while low frequency aspects of the signal vary slowly in time. We can represent a signal by the magnitudes that each frequency contributes to the total signal. So, instead of plotting a signal versus $t$, we can plot the magnitudes, or **Fourier coefficients**, of each frequency versus the frequency, $f$. A plot like this is called a **periodogram**, and is useful to see what is present in the signal.\n",
"\n",
"We can do much more with these magnitudes corresponding to frequencies. Say we wanted to filter out high frequency parts of the signal because they constitute noise. We could just go to the frequency domain, set the magnitudes for high frequencies to zero, and then go back to the time domain.\n",
"\n",
"The **Fourier transform** is a useful mathematical tool to take us to the frequency domain and back, via the inverse Fourier transform. The **fast Fourier transform**, or FFT, is a very fast algorithm to do this, and it is implemented in NumPy and SciPy. Furthermore, the `scipy.signal` module has many convenient functions for filtering in the frequency domain. We will make use of a **Butterworth filter**, which is one of the simplest frequency domain filters available."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The periodogram\n",
"\n",
"Let's see what frequencies are present in our signal by computing the periodogram. We do this by getting the frequencies associated with the Fourier coefficients using `np.fft.fftfreq()` and then computing the **power spectral density**, which is the square of the magnitude of the Fourier coefficients. We can then plot the result. For reasons we will not get in to here, we will only consider positive frequencies, since the Fourier coefficients for the negative frequencies are symmetric. Note that we could generate the periodogram automatically using\n",
"\n",
" f, psd = scipy.signal.periodogram(df['V (µV)'], fs=sampling_freq)\n",
" \n",
"but this is much slower to compute, due to all the bells and whistles included in `scipy.signal`.\n",
"\n",
"Computing an FFT for a long signal that does not have many prime factors is very inefficient. This data set has 2682401 time points, which has a prime factorization of 141179×19, so the performance of the FFT algorithm may be very slow. For the purposes of this periodogram, we will only consider the first 2097152 points, since this is a power of two. (I timed this a couple years ago, and computing the periodogram for the array of length 2097152 takes about 160 ms, and computing it for the array of length 2682401 takes 8.5 minutes! That's a 3000-fold difference!)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
""
],
"text/plain": [
":DynamicMap []\n",
" :RGB [f (Hz),PSD] (R,G,B,A)"
]
},
"execution_count": 11,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1336"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Truncated signal for faster calculation\n",
"max_ind = 2**(int(np.floor(np.log2(len(df)))))\n",
"V_trunc = df['V (µV)'].values[:max_ind]\n",
"\n",
"# Determine frequencies\n",
"f = np.fft.fftfreq(len(V_trunc)) * sampling_freq\n",
"\n",
"# Compute power spectral density\n",
"psd = np.abs(np.fft.fft(V_trunc))**2 / len(V_trunc)\n",
"\n",
"# Make into DataFrame to enable DataShaded viewing\n",
"df_psd = pd.DataFrame(data={'f (Hz)': f[f>=0], 'PSD': psd[f>=0]})\n",
"\n",
"# Make plot\n",
"hv.operation.datashader.datashade(\n",
" hv.Curve(\n",
" data=df_psd,\n",
" kdims=['f (Hz)'],\n",
" vdims=['PSD']\n",
" ),\n",
" aggregator='any',\n",
").opts(\n",
" frame_height=300,\n",
" frame_width=625,\n",
" show_grid=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we have frequencies up to 12,500 Hz. This is the **Nyquist frequency**, half the sampling frequency, which is the highest frequency that can be coded at a given sampling rate (in our case, 25,000 Hz).\n",
"\n",
"There is a big peak at 8000 Hz, which corresponds to the dominant frequency in the signal. This is high-frequency noise, but it is important to keep these frequencies because the signal in spikes undergoes a rapid change in just a fraction of a millisecond. So, we probably want to keep the high frequencies.\n",
"\n",
"The very low frequency peak at $f = 0$ corresponds to the constant, or DC, signal. That obviously also does not include spikes. There is also some low frequency oscillation, which are clearly not related to spiking. If we want to detect spikes by setting a threshold voltage, though, it is important to filter out these low frequency variations.\n",
"\n",
"As we saw from looking at the plots, a spike occurs over about 10 or 20 ms, which would be captured by sampling at about 200 Hz (2/0.01 s, according to the Nyquist sampling theorem). So, to make sure we get good resolution on the spikes, we should be sure not to filter out frequencies that are too close to 100 Hz."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Butterworth filter: a frequency filter\n",
"\n",
"There are lots and lots of filters available, but we will use a simple [Butterworth filter](https://en.wikipedia.org/wiki/Butterworth_filter). We won't go into the theory of filters, transfer functions, etc., here, but will get a graphical idea of how a Butterworth filter works. A Butterworth filter is a **band pass** filter, meaning that it lets certain bands of frequency \"pass\" and be included in the filtered signal. We want to construct a highpass filter because we just want to filter out the low frequencies.\n",
"\n",
"We construct a Butterworth filter using `scipy.signal.butter()`. This function returns coefficients on polynomials for the numerator and denominator of a Butterworth filter, which you can read about on the [Wikipedia page](https://en.wikipedia.org/wiki/Butterworth_filter). We have to specify the order of the polynomials (we'll use third order) and the cutoff frequency, which defines the frequency below which modes are cut off. The cutoff frequency has to be in units of the Nyquist frequency. We will put our cutoff at 50 Hz."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Compute Nyquist frequency\n",
"nyquist_freq = sampling_freq / 2\n",
"\n",
"# Design a butterworth filter\n",
"b, a = scipy.signal.butter(3, 100 / nyquist_freq, btype='high')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's look at what frequencies get through this filter (the so-called **frequency response**). We can get the frequencies (in units of $\\pi/f_\\text{Nyquist}$, or radians/sample) and the the response using `scipy.signal.freqz()`."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Curve [freq. (Hz)] (response)"
]
},
"execution_count": 13,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1448"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Get frequency response curve\n",
"w, h = scipy.signal.freqz(b, a, worN=2000)\n",
"\n",
"# Make plot\n",
"hv.Curve(\n",
" data=((nyquist_freq / np.pi) * w, abs(h)),\n",
" kdims=['freq. (Hz)'],\n",
" vdims=['response'],\n",
").opts(\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" frame_height=250,\n",
" frame_width=400,\n",
" logx=True,\n",
" padding=0.05,\n",
" show_grid=True,\n",
" xlim=(10, 1e4),\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, we will block out frequencies below about 50 Hz and let high frequencies pass through. Let's filter the voltages and look at the result. It's easiest to see how the filter worked if we again just look at a subset of the whole signal."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.Unfiltered :Curve [t (s)] (V (µV))\n",
" .Curve.Filtered :Curve [t (s)] (V filtered (µV))"
]
},
"execution_count": 14,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1558"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Filter data\n",
"df['V filtered (µV)'] = scipy.signal.lfilter(b, a, df['V (µV)'])\n",
"\n",
"# Plot the result, overlaying filtered and unfiltered\n",
"curves = [\n",
" hv.Curve(\n",
" data=df.iloc[10000:30000, :],\n",
" kdims=['t (s)'],\n",
" vdims=V,\n",
" label=filtered,\n",
" ).opts(\n",
" color=hv.Cycle(bebi103.hv.default_categorical_cmap),\n",
" frame_height=300,\n",
" frame_width=575,\n",
" show_grid=True,\n",
" )\n",
" for filtered, V in zip(['unfiltered', 'filtered'], ['V (µV)', 'V filtered (µV)'])\n",
"]\n",
"\n",
"hv.Overlay(curves)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The filtering did a good job of limited longer time scale variations while keeping the spikes intact."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding the spikes\n",
"\n",
"We'll now use our filtered trace to locate the spikes. Our goal is to find the spikes, store a neighborhood around the spike, and then find the position of the spike to sub-sampling accuracy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Locating spikes\n",
"To locate the spikes, we will search for places where the signal dips below a threshold and then rises up above it again. This will mark boundaries of the spike, and we will then keep half a millisecond of samples on the left side of the threshold crossing and 2 milliseconds on the right."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# Threshold below which spike has to drop (in µV)\n",
"thresh = -40\n",
"\n",
"# Amount of time in seconds to keep on either side of spike\n",
"time_window_left = 0.0005\n",
"time_window_right = 0.002\n",
"\n",
"# Number of samples to keep on either side of spike\n",
"n_window_left = int(time_window_left * sampling_freq)\n",
"n_window_right = int(time_window_right * sampling_freq)\n",
"\n",
"# DataFrame to store spikes\n",
"df_spike = pd.DataFrame(columns=['spike', 't (s)', 'V (µV)'])\n",
"\n",
"# Use a NumPy array for speed in looping\n",
"V = df['V filtered (µV)'].values\n",
"\n",
"# Initialize while loop\n",
"i = 1\n",
"spike = 0\n",
"while i < len(df):\n",
" if V[i] < thresh:\n",
" # Found a spike, get crossings\n",
" cross_1 = i\n",
" while V[i] < thresh:\n",
" i += 1\n",
" cross_2 = i\n",
" \n",
" # Store perintent quantities in DataFrame\n",
" t_in = df['t (s)'][cross_1-n_window_left:cross_2+n_window_right]\n",
" V_in = df['V (µV)'][cross_1-n_window_left:cross_2+n_window_right]\n",
" data={'t (s)': t_in,\n",
" 'V (µV)': V_in,\n",
" 'spike': spike * np.ones_like(t_in, dtype=int)}\n",
" df_add = pd.DataFrame(data)\n",
" df_spike = pd.concat((df_spike, df_add), sort=True)\n",
"\n",
" spike += 1\n",
"\n",
" i += 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have found the spikes, let's plot them all on top of each other to ascertain their shape. First, we should make a set of time points to plot them all on top of each other."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"df_spike[\"relative time (s)\"] = df_spike.groupby(\"spike\")[\"t (s)\"].transform(\n",
" lambda x: np.arange(len(x)) / sampling_freq\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can make the plots."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":NdOverlay [spike]\n",
" :Curve [relative time (s)] (V (µV),spike)"
]
},
"execution_count": 17,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1779"
}
},
"output_type": "execute_result"
}
],
"source": [
"hv.Curve(\n",
" data=df_spike,\n",
" kdims=['relative time (s)'],\n",
" vdims=['V (µV)', 'spike']\n",
").groupby(\n",
" 'spike'\n",
").opts(\n",
" hv.opts.Curve(\n",
" alpha=0.05,\n",
" frame_height=300,\n",
" frame_width=575,\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" line_width=1,\n",
" show_grid=True,\n",
" )\n",
").overlay()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A smooth interpolant for the spike potential\n",
"\n",
"We would like to describe the spike with a smooth function. In doing this, we can estimate the location of the bottom of the spike to greater accuracy than dictated by the sampling rate. We call a smooth curve that passes through the noisy data and **interpolant**.\n",
"\n",
"Let's look at a single spike as we develop the interpolant."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.I :Curve [t (s)] (V (µV))\n",
" .Scatter.I :Scatter [t (s)] (V (µV))"
]
},
"execution_count": 18,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1945"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Pull out 0th spike\n",
"spike = 0\n",
"df_0 = df_spike.loc[df_spike['spike']==spike, :]\n",
"\n",
"# Plot the spike\n",
"def base_spike_plot():\n",
" line_plot = hv.Curve(\n",
" data=df_0,\n",
" kdims=['t (s)'],\n",
" vdims=['V (µV)']\n",
" ).opts(\n",
" frame_height=300,\n",
" frame_width=575,\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" line_width=1,\n",
" padding=0.05,\n",
" show_grid=True, \n",
" )\n",
"\n",
" scatter_plot = hv.Scatter(\n",
" data=df_0,\n",
" kdims=['t (s)'],\n",
" vdims=['V (µV)']\n",
" ).opts(\n",
" color=bebi103.hv.default_categorical_cmap[0],\n",
" size=5,\n",
" )\n",
" \n",
" return line_plot * scatter_plot\n",
"\n",
"base_spike_plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will consider several ways to get a smooth interpolant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Smoothing with basis functions (splines)\n",
"\n",
"A function $f(x)$ (or data) can generally be written as an expansion of $M$ orthogonal **basis functions** $h_m$, such that\n",
"\n",
"\\begin{align}\n",
"f(x) = \\sum_{m=1}^M \\beta_m\\,h_m(x),\n",
"\\end{align}\n",
"\n",
"where $\\beta_m$ are coefficients. Examples of basis functions you may be familiar with are sines/cosines (used to construct a Fourier series), Legendre polynomials, spherical harmonics, etc. The idea behind smoothing with basis functions is to control the complexity of $f(x)$ by restricting the number of basis functions or selecting certain basis functions to use. This restricted $f(x)$, being less complex than the data, is smoother.\n",
"\n",
"A commonly used set of basis function are cubic polynomials. An approximation of this procedure is called a **spline**. If we have $k$ joining points of cubic functions (called *knots*), the data are represented with $k$ basis functions.\n",
"\n",
"The `scipy.interpolate` module has spline tools. Importantly, it selects the number of knots to use until a smoothing condition is met. This condition is defined by\n",
"\n",
"\\begin{align}\n",
"\\sum_{i\\in D} \\left[w_i(y_i - \\hat{y}(x_i))\\right]^2 \\le s,\n",
"\\end{align}\n",
"\n",
"where $w_i$ is the weight we wish to assign to data point $i$, and $\\hat{y}(x_i)$ is the smoothed function at $x_i$. We will take all $w_i$'s to be equal in our analysis. A good rule of thumb is to choose\n",
"\n",
"\\begin{align}\n",
"s = f\\sum_{i\\in D} y_i^2,\n",
"\\end{align}\n",
"\n",
"where $f$ is roughly the fractional difference between the rough and smooth curves."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.I :Curve [t (s)] (V (µV))\n",
" .Scatter.I :Scatter [t (s)] (V (µV))\n",
" .Curve.II :Curve [x] (y)"
]
},
"execution_count": 19,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "2149"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Determine smoothing factor from rule of thumb (use f = 0.05)\n",
"smooth_factor = 0.05 * (df_0['V (µV)']**2).sum()\n",
"\n",
"# Set up an scipy.interpolate.UnivariateSpline instance\n",
"ts = df_0['t (s)'].values\n",
"Vs = df_0['V (µV)'].values\n",
"spl = scipy.interpolate.UnivariateSpline(ts, Vs, s=smooth_factor)\n",
"\n",
"# spl is now a callable function\n",
"ts_spline = np.linspace(ts[0], ts[-1], 400)\n",
"Vs_spline = spl(ts_spline)\n",
"\n",
"# Plot results\n",
"base_spike_plot() * hv.Curve(\n",
" data=(ts_spline, Vs_spline)\n",
").opts(\n",
" color='#ff7f0e', \n",
" line_width=2\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Perhaps we chose our smoothing factor to be too small. But a better explanation for the windy path of the spline curve is that the spike is highly localized and then imposes a large variation on the data. It might be better to only fit a spline in the close neighborhood of the spike."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.I :Curve [t (s)] (V (µV))\n",
" .Scatter.I :Scatter [t (s)] (V (µV))\n",
" .Curve.II :Curve [x] (y)"
]
},
"execution_count": 20,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "2391"
}
},
"output_type": "execute_result"
}
],
"source": [
"minind = np.argmin(Vs)\n",
"ts_short = ts[minind-7:minind+7]\n",
"Vs_short = Vs[minind-7:minind+7]\n",
"\n",
"smooth_factor = 0.05 * (Vs_short**2).sum()\n",
"spl = scipy.interpolate.UnivariateSpline(ts_short, Vs_short, s=smooth_factor)\n",
"ts_spline = np.linspace(ts_short[0], ts_short[-1], 400)\n",
"Vs_spline = spl(ts_spline)\n",
"\n",
"# Plot results\n",
"base_spike_plot() * hv.Curve(\n",
" data=(ts_spline, Vs_spline)\n",
").opts(\n",
" color='#ff7f0e', \n",
" line_width=2\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Kernel-based smoothing\n",
"\n",
"Another option for smoothing is a **kernel-based method**. We will use a **Nadaraya-Watson kernel estimator**, which is essentially a weighted average of the data points. \n",
"\n",
"When we have a high sampling rate, the extra fine resolution often just gives the noise. We could do a sliding average over the data. We already performed this operation in Tutorial 2 when we were looking at the fish sleep data. I will write this process more concretely. For ease of notation, we arrange our data points in time, i.e., index $i+1$ follows $i$ in time, with a total of $n$ data points. For data point $(x_i, y_i)$, we get a smoothed value $\\hat{x}_i$ by\n",
"\n",
"\\begin{align}\n",
"\\hat{y}(x_0) = \\frac{\\sum_{i=1}^n K_\\lambda(x_0, x_i)\\, y_i}{\\sum_{i=1}^n K_\\lambda(x_0,x_i)}\n",
"\\end{align}\n",
"\n",
"where $K_\\lambda(x_0, x_i)$ is the **kernel**, or weight. If we took each data point to be the mean of itself and its 10 nearest-neighbors, we would have (neglecting end effects)\n",
"\n",
"\\begin{align}\n",
"K_\\lambda(x_j, x_i) = \\left\\{\\begin{array}{cl}\n",
"\\frac{1}{11} & \\left|i - j\\right| \\le 5 \\\\[1mm]\n",
"0 & \\text{otherwise}.\n",
"\\end{array}\n",
"\\right.\n",
"\\end{align}\n",
"\n",
"Commonly used Kernels are the Epanechnikov quadratic kernel, the tri-cube kernel, and the Gaussian kernel. We can write them generically as\n",
"\n",
"\\begin{align}\n",
"K_\\lambda(x_0, x_i) &= D\\left(\\frac{\\left|x_i-x_0\\right|}{\\lambda}\\right)\n",
"\\equiv D(t), \\\\[1mm]\n",
"\\text{where } t &= \\frac{\\left|x_i-x_0\\right|}{\\lambda}.\n",
"\\end{align}\n",
"\n",
"These kernels are\n",
"\n",
"\\begin{align}\n",
"\\text{Epanechnikov:}\\;\\;\\;\\;\\; D &= \\left\\{\\begin{array}{cl}\n",
"\\frac{3}{4}\\left(1 - t^2\\right) & |t| \\le 1 \\\\[1mm]\n",
"0 & \\text{otherwise}\n",
"\\end{array}\\right. \\\\[5mm]\n",
"\\text{tri-cube:}\\;\\;\\;\\;\\;D&=\\left\\{\\begin{array}{cl}\n",
"\\left(1 - \\left|t\\right|^3\\right)^3 & |t| \\le 1 \\\\[1mm]\n",
"0 & \\text{otherwise}\n",
"\\end{array}\\right. \\\\[5mm]\n",
"\\text{Gaussian:}\\;\\;\\;\\;\\;D &= \\mathrm{e}^{-t^2}.\n",
"\\end{align}\n",
"\n",
"Let's code up functions for these kernels."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# Define our kernels\n",
"def epan_kernel(t):\n",
" \"\"\"\n",
" Epanechnikov kernel.\n",
" \"\"\"\n",
" return np.logical_and(t > -1.0, t < 1.0) * 3.0 * (1.0 - t**2) / 4.0\n",
"\n",
"\n",
"def tri_cube_kernel(t):\n",
" \"\"\"\n",
" Tri-cube kernel.\n",
" \"\"\"\n",
" return np.logical_and(t > -1.0, t < 1.0) * (1.0 - abs(t**3))**3\n",
"\n",
"\n",
"def gauss_kernel(t):\n",
" \"\"\"\n",
" Gaussian kernel.\n",
" \"\"\"\n",
" return np.exp(-t**2 / 2.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can compare the respective kernels to get a feel for how the averaging goes."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"(function(root) {\n",
" function embed_document(root) {\n",
" \n",
" var docs_json = {\"7826c7a6-4eb7-4e8b-8cd8-f26d537e7e31\":{\"roots\":{\"references\":[{\"attributes\":{\"below\":[{\"id\":\"2640\",\"type\":\"LinearAxis\"}],\"center\":[{\"id\":\"2644\",\"type\":\"Grid\"},{\"id\":\"2649\",\"type\":\"Grid\"},{\"id\":\"2675\",\"type\":\"Legend\"}],\"left\":[{\"id\":\"2645\",\"type\":\"LinearAxis\"}],\"plot_height\":300,\"plot_width\":450,\"renderers\":[{\"id\":\"2666\",\"type\":\"GlyphRenderer\"},{\"id\":\"2680\",\"type\":\"GlyphRenderer\"},{\"id\":\"2695\",\"type\":\"GlyphRenderer\"}],\"title\":{\"id\":\"2668\",\"type\":\"Title\"},\"toolbar\":{\"id\":\"2656\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"2632\",\"type\":\"DataRange1d\"},\"x_scale\":{\"id\":\"2636\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"2634\",\"type\":\"DataRange1d\"},\"y_scale\":{\"id\":\"2638\",\"type\":\"LinearScale\"}},\"id\":\"2631\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"AAAAAAAAEMC+NHEYq60PwH1p4jBWWw/AO55TSQEJD8D60sRhrLYOwLgHNnpXZA7AdjynkgISDsA1cRirrb8NwPOlicNYbQ3Astr62wMbDcBwD2z0rsgMwC5E3QxadgzA7XhOJQUkDMCrrb89sNELwGriMFZbfwvAKBeibgYtC8DmSxOHsdoKwKWAhJ9ciArAY7X1twc2CsAi6mbQsuMJwOAe2OhdkQnAnlNJAQk/CcBciLoZtOwIwBu9KzJfmgjA2vGcSgpICMCYJg5jtfUHwFZbf3tgowfAFJDwkwtRB8DTxGGstv4GwJL50sRhrAbAUC5E3QxaBsAOY7X1twcGwMyXJg5jtQXAi8yXJg5jBcBKAQk/uRAFwAg2eldkvgTAxmrrbw9sBMCEn1yIuhkEwEPUzaBlxwPAAQk/uRB1A8DAPbDRuyIDwH5yIepm0ALAPKeSAhJ+AsD72wMbvSsCwLkQdTNo2QHAeEXmSxOHAcA2eldkvjQBwPSuyHxp4gDAs+M5lRSQAMBxGKutvz0AwGCaOIzV1v+/3AMbvSsy/79Ybf3tgY3+v9bW3x7Y6P2/UkDCTy5E/b/OqaSAhJ/8v0wTh7Ha+vu/yHxp4jBW+79G5ksTh7H6v8JPLkTdDPq/PrkQdTNo+b+8IvOlicP4vziM1dbfHvi/tvW3BzZ6978yX5o4jNX2v67IfGniMPa/LDJfmjiM9b+om0HLjuf0vyYFJPzkQvS/om4GLTue878e2Ohdkfnyv5xBy47nVPK/GKutvz2w8b+WFJDwkwvxvxJ+ciHqZvC/HM+ppICE778Yom4GLTvuvxB1M2jZ8ey/CEj4yYWo678EG70rMl/qv/ztgY3eFem/+MBG74rM57/wkwtRN4Pmv+hm0LLjOeW/5DmVFJDw47/cDFp2PKfiv9jfHtjoXeG/0LLjOZUU4L+QC1E3g5bdv4ix2vrbA9u/eFdkvjRx2L9w/e2Bjd7Vv2Cjd0XmS9O/UEkBCT+50L+Q3hWZL03Mv3AqKSDhJ8e/YHY8p5ICwr+AhJ9ciLq5v4A4jNXW366/gNCy4zmVlL8A0LLjOZWUP4A4jNXW364/gISfXIi6uT9gdjynkgLCP4AqKSDhJ8c/gN4VmS9NzD9QSQEJP7nQP2Cjd0XmS9M/cP3tgY3e1T+AV2S+NHHYP5Cx2vrbA9s/kAtRN4OW3T/QsuM5lRTgP9jfHtjoXeE/4Axadjyn4j/oOZUUkPDjP+hm0LLjOeU/8JMLUTeD5j/4wEbvisznPwDugY3eFek/CBu9KzJf6j8ISPjJhajrPxB1M2jZ8ew/GKJuBi077j8gz6mkgITvPxR+ciHqZvA/lBSQ8JML8T8Yq62/PbDxP5xBy47nVPI/INjoXZH58j+kbgYtO57zPyQFJPzkQvQ/qJtBy47n9D8sMl+aOIz1P7DIfGniMPY/NF+aOIzV9j+09bcHNnr3PziM1dbfHvg/vCLzpYnD+D9AuRB1M2j5P8RPLkTdDPo/ROZLE4ex+j/IfGniMFb7P0wTh7Ha+vs/0KmkgISf/D9UQMJPLkT9P9TW3x7Y6P0/WG397YGN/j/cAxu9KzL/P2CaOIzV1v8/chirrb89AEC04zmVFJAAQPSuyHxp4gBANnpXZL40AUB4ReZLE4cBQLoQdTNo2QFA/NsDG70rAkA8p5ICEn4CQH5yIepm0AJAwD2w0bsiA0ACCT+5EHUDQETUzaBlxwNAhJ9ciLoZBEDGautvD2wEQAg2eldkvgRASgEJP7kQBUCMzJcmDmMFQMyXJg5jtQVADmO19bcHBkBQLkTdDFoGQJL50sRhrAZA1MRhrLb+BkAUkPCTC1EHQFZbf3tgowdAmCYOY7X1B0Da8ZxKCkgIQBy9KzJfmghAXIi6GbTsCECeU0kBCT8JQOAe2OhdkQlAIupm0LLjCUBktfW3BzYKQKSAhJ9ciApA5ksTh7HaCkAoF6JuBi0LQGriMFZbfwtArK2/PbDRC0DseE4lBSQMQC5E3QxadgxAcA9s9K7IDECy2vrbAxsNQPSlicNYbQ1ANHEYq62/DUB2PKeSAhIOQLgHNnpXZA5A+tLEYay2DkA8nlNJAQkPQH5p4jBWWw9AvjRxGKutD0AAAAAAAAAQQA==\",\"dtype\":\"float64\",\"shape\":[200]},\"y\":{\"__ndarray__\":\"AAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAzKjrzTn6lj/UIo6d1aK0P/jxmj0idME/qhY5hm1HyD9tf6GoTMvOPzcWatLff9I/wo5oPWNy1T9oKUyVMD3YPyHmFNpH4No/5cTCC6lb3T/GxVUqVK/fP1r05pqk7eA/3ZYVF8Tv4T9pyraJCN7iP/2OyvJxuOM/neRQUgB/5D9Ey0moszHlP/dCtfSL0OU/tEuTN4lb5j965eNwq9LmP0kQp6DyNec/Iszcxl6F5z8FGYXj78DnP/H2n/al6Oc/6GUtAIH85z/oZS0AgfznP/H2n/al6Oc/BRmF4+/A5z8izNzGXoXnP0gQp6DyNec/e+XjcKvS5j+0S5M3iVvmP/dCtfSL0OU/RMtJqLMx5T+b5FBSAH/kP/qOyvJxuOM/acq2iQje4j/dlhUXxO/hP1r05pqk7eA/v8VVKlSv3z/fxMILqVvdPyHmFNpH4No/aClMlTA92D/Cjmg9Y3LVPy0WatLff9I/W3+hqEzLzj+qFjmGbUfIP/jxmj0idME/1CKOndWitD8MqOvNOfqWPwAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgA==\",\"dtype\":\"float64\",\"shape\":[200]}},\"selected\":{\"id\":\"2689\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"2688\",\"type\":\"UnionRenderers\"}},\"id\":\"2663\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"2672\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"2654\",\"type\":\"ResetTool\"},{\"attributes\":{},\"id\":\"2651\",\"type\":\"WheelZoomTool\"},{\"attributes\":{},\"id\":\"2705\",\"type\":\"UnionRenderers\"},{\"attributes\":{\"data_source\":{\"id\":\"2692\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"2693\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"2694\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"2696\",\"type\":\"CDSView\"}},\"id\":\"2695\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"callback\":null},\"id\":\"2632\",\"type\":\"DataRange1d\"},{\"attributes\":{},\"id\":\"2653\",\"type\":\"SaveTool\"},{\"attributes\":{\"source\":{\"id\":\"2692\",\"type\":\"ColumnDataSource\"}},\"id\":\"2696\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"2688\",\"type\":\"UnionRenderers\"},{\"attributes\":{},\"id\":\"2650\",\"type\":\"PanTool\"},{\"attributes\":{},\"id\":\"2670\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{\"label\":{\"value\":\"Tri-cube\"},\"renderers\":[{\"id\":\"2680\",\"type\":\"GlyphRenderer\"}]},\"id\":\"2691\",\"type\":\"LegendItem\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"AAAAAAAAEMC+NHEYq60PwH1p4jBWWw/AO55TSQEJD8D60sRhrLYOwLgHNnpXZA7AdjynkgISDsA1cRirrb8NwPOlicNYbQ3Astr62wMbDcBwD2z0rsgMwC5E3QxadgzA7XhOJQUkDMCrrb89sNELwGriMFZbfwvAKBeibgYtC8DmSxOHsdoKwKWAhJ9ciArAY7X1twc2CsAi6mbQsuMJwOAe2OhdkQnAnlNJAQk/CcBciLoZtOwIwBu9KzJfmgjA2vGcSgpICMCYJg5jtfUHwFZbf3tgowfAFJDwkwtRB8DTxGGstv4GwJL50sRhrAbAUC5E3QxaBsAOY7X1twcGwMyXJg5jtQXAi8yXJg5jBcBKAQk/uRAFwAg2eldkvgTAxmrrbw9sBMCEn1yIuhkEwEPUzaBlxwPAAQk/uRB1A8DAPbDRuyIDwH5yIepm0ALAPKeSAhJ+AsD72wMbvSsCwLkQdTNo2QHAeEXmSxOHAcA2eldkvjQBwPSuyHxp4gDAs+M5lRSQAMBxGKutvz0AwGCaOIzV1v+/3AMbvSsy/79Ybf3tgY3+v9bW3x7Y6P2/UkDCTy5E/b/OqaSAhJ/8v0wTh7Ha+vu/yHxp4jBW+79G5ksTh7H6v8JPLkTdDPq/PrkQdTNo+b+8IvOlicP4vziM1dbfHvi/tvW3BzZ6978yX5o4jNX2v67IfGniMPa/LDJfmjiM9b+om0HLjuf0vyYFJPzkQvS/om4GLTue878e2Ohdkfnyv5xBy47nVPK/GKutvz2w8b+WFJDwkwvxvxJ+ciHqZvC/HM+ppICE778Yom4GLTvuvxB1M2jZ8ey/CEj4yYWo678EG70rMl/qv/ztgY3eFem/+MBG74rM57/wkwtRN4Pmv+hm0LLjOeW/5DmVFJDw47/cDFp2PKfiv9jfHtjoXeG/0LLjOZUU4L+QC1E3g5bdv4ix2vrbA9u/eFdkvjRx2L9w/e2Bjd7Vv2Cjd0XmS9O/UEkBCT+50L+Q3hWZL03Mv3AqKSDhJ8e/YHY8p5ICwr+AhJ9ciLq5v4A4jNXW366/gNCy4zmVlL8A0LLjOZWUP4A4jNXW364/gISfXIi6uT9gdjynkgLCP4AqKSDhJ8c/gN4VmS9NzD9QSQEJP7nQP2Cjd0XmS9M/cP3tgY3e1T+AV2S+NHHYP5Cx2vrbA9s/kAtRN4OW3T/QsuM5lRTgP9jfHtjoXeE/4Axadjyn4j/oOZUUkPDjP+hm0LLjOeU/8JMLUTeD5j/4wEbvisznPwDugY3eFek/CBu9KzJf6j8ISPjJhajrPxB1M2jZ8ew/GKJuBi077j8gz6mkgITvPxR+ciHqZvA/lBSQ8JML8T8Yq62/PbDxP5xBy47nVPI/INjoXZH58j+kbgYtO57zPyQFJPzkQvQ/qJtBy47n9D8sMl+aOIz1P7DIfGniMPY/NF+aOIzV9j+09bcHNnr3PziM1dbfHvg/vCLzpYnD+D9AuRB1M2j5P8RPLkTdDPo/ROZLE4ex+j/IfGniMFb7P0wTh7Ha+vs/0KmkgISf/D9UQMJPLkT9P9TW3x7Y6P0/WG397YGN/j/cAxu9KzL/P2CaOIzV1v8/chirrb89AEC04zmVFJAAQPSuyHxp4gBANnpXZL40AUB4ReZLE4cBQLoQdTNo2QFA/NsDG70rAkA8p5ICEn4CQH5yIepm0AJAwD2w0bsiA0ACCT+5EHUDQETUzaBlxwNAhJ9ciLoZBEDGautvD2wEQAg2eldkvgRASgEJP7kQBUCMzJcmDmMFQMyXJg5jtQVADmO19bcHBkBQLkTdDFoGQJL50sRhrAZA1MRhrLb+BkAUkPCTC1EHQFZbf3tgowdAmCYOY7X1B0Da8ZxKCkgIQBy9KzJfmghAXIi6GbTsCECeU0kBCT8JQOAe2OhdkQlAIupm0LLjCUBktfW3BzYKQKSAhJ9ciApA5ksTh7HaCkAoF6JuBi0LQGriMFZbfwtArK2/PbDRC0DseE4lBSQMQC5E3QxadgxAcA9s9K7IDECy2vrbAxsNQPSlicNYbQ1ANHEYq62/DUB2PKeSAhIOQLgHNnpXZA5A+tLEYay2DkA8nlNJAQkPQH5p4jBWWw9AvjRxGKutD0AAAAAAAAAQQA==\",\"dtype\":\"float64\",\"shape\":[200]},\"y\":{\"__ndarray__\":\"AAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAp2NJMxYHFD8GioU3E09rP+IcVuFTF48/xltrEQyuoz8eC04kpOGyP8Y7se23y74/Ur/YkrKPxj/6hxLRvaHOP+Z8lzIepNM/PDlO/gQY2D9qltf7D4bcPwimB0hOZuA/VXCyOv1n4j/64gDpcj3kP56G65tG3+U/495IEiFJ5z9UH5QKd3noP1ZdXw0xcek/dPGBikkz6j8tgXGkaMTqPw32vXyDKus/TXleeoNs6z+nY1/W+JHrP1Pvjr3aous/vl4Ahlan6z++XgCGVqfrP1Pvjr3aous/p2Nf1viR6z9NeV56g2zrPw32vXyDKus/L4FxpGjE6j908YGKSTPqP1ZdXw0xcek/VB+UCnd56D/f3kgSIUnnP5qG65tG3+U/+uIA6XI95D9VcLI6/WfiPwimB0hOZuA/XJbX+w+G3D8uOU7+BBjYP+Z8lzIepNM/+ocS0b2hzj9Sv9iSso/GP547se23y74//gpOJKThsj/GW2sRDK6jP+IcVuFTF48/BoqFNxNPaz+tYUkzFgcUPwAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAgA==\",\"dtype\":\"float64\",\"shape\":[200]}},\"selected\":{\"id\":\"2706\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"2705\",\"type\":\"UnionRenderers\"}},\"id\":\"2677\",\"type\":\"ColumnDataSource\"},{\"attributes\":{\"overlay\":{\"id\":\"2674\",\"type\":\"BoxAnnotation\"}},\"id\":\"2652\",\"type\":\"BoxZoomTool\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2679\",\"type\":\"Line\"},{\"attributes\":{\"data_source\":{\"id\":\"2663\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"2664\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"2665\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"2667\",\"type\":\"CDSView\"}},\"id\":\"2666\",\"type\":\"GlyphRenderer\"},{\"attributes\":{\"callback\":null},\"id\":\"2634\",\"type\":\"DataRange1d\"},{\"attributes\":{},\"id\":\"2641\",\"type\":\"BasicTicker\"},{\"attributes\":{\"source\":{\"id\":\"2677\",\"type\":\"ColumnDataSource\"}},\"id\":\"2681\",\"type\":\"CDSView\"},{\"attributes\":{\"line_color\":\"#9467bd\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2693\",\"type\":\"Line\"},{\"attributes\":{\"label\":{\"value\":\"Epanechnikov\"},\"renderers\":[{\"id\":\"2666\",\"type\":\"GlyphRenderer\"}]},\"id\":\"2676\",\"type\":\"LegendItem\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2665\",\"type\":\"Line\"},{\"attributes\":{\"label\":{\"value\":\"Gaussian\"},\"renderers\":[{\"id\":\"2695\",\"type\":\"GlyphRenderer\"}]},\"id\":\"2708\",\"type\":\"LegendItem\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"2674\",\"type\":\"BoxAnnotation\"},{\"attributes\":{\"ticker\":{\"id\":\"2641\",\"type\":\"BasicTicker\"}},\"id\":\"2644\",\"type\":\"Grid\"},{\"attributes\":{\"text\":\"\",\"text_color\":{\"value\":\"black\"},\"text_font_size\":{\"value\":\"12pt\"}},\"id\":\"2668\",\"type\":\"Title\"},{\"attributes\":{\"items\":[{\"id\":\"2676\",\"type\":\"LegendItem\"},{\"id\":\"2691\",\"type\":\"LegendItem\"},{\"id\":\"2708\",\"type\":\"LegendItem\"}]},\"id\":\"2675\",\"type\":\"Legend\"},{\"attributes\":{},\"id\":\"2636\",\"type\":\"LinearScale\"},{\"attributes\":{\"callback\":null,\"data\":{\"x\":{\"__ndarray__\":\"AAAAAAAAEMC+NHEYq60PwH1p4jBWWw/AO55TSQEJD8D60sRhrLYOwLgHNnpXZA7AdjynkgISDsA1cRirrb8NwPOlicNYbQ3Astr62wMbDcBwD2z0rsgMwC5E3QxadgzA7XhOJQUkDMCrrb89sNELwGriMFZbfwvAKBeibgYtC8DmSxOHsdoKwKWAhJ9ciArAY7X1twc2CsAi6mbQsuMJwOAe2OhdkQnAnlNJAQk/CcBciLoZtOwIwBu9KzJfmgjA2vGcSgpICMCYJg5jtfUHwFZbf3tgowfAFJDwkwtRB8DTxGGstv4GwJL50sRhrAbAUC5E3QxaBsAOY7X1twcGwMyXJg5jtQXAi8yXJg5jBcBKAQk/uRAFwAg2eldkvgTAxmrrbw9sBMCEn1yIuhkEwEPUzaBlxwPAAQk/uRB1A8DAPbDRuyIDwH5yIepm0ALAPKeSAhJ+AsD72wMbvSsCwLkQdTNo2QHAeEXmSxOHAcA2eldkvjQBwPSuyHxp4gDAs+M5lRSQAMBxGKutvz0AwGCaOIzV1v+/3AMbvSsy/79Ybf3tgY3+v9bW3x7Y6P2/UkDCTy5E/b/OqaSAhJ/8v0wTh7Ha+vu/yHxp4jBW+79G5ksTh7H6v8JPLkTdDPq/PrkQdTNo+b+8IvOlicP4vziM1dbfHvi/tvW3BzZ6978yX5o4jNX2v67IfGniMPa/LDJfmjiM9b+om0HLjuf0vyYFJPzkQvS/om4GLTue878e2Ohdkfnyv5xBy47nVPK/GKutvz2w8b+WFJDwkwvxvxJ+ciHqZvC/HM+ppICE778Yom4GLTvuvxB1M2jZ8ey/CEj4yYWo678EG70rMl/qv/ztgY3eFem/+MBG74rM57/wkwtRN4Pmv+hm0LLjOeW/5DmVFJDw47/cDFp2PKfiv9jfHtjoXeG/0LLjOZUU4L+QC1E3g5bdv4ix2vrbA9u/eFdkvjRx2L9w/e2Bjd7Vv2Cjd0XmS9O/UEkBCT+50L+Q3hWZL03Mv3AqKSDhJ8e/YHY8p5ICwr+AhJ9ciLq5v4A4jNXW366/gNCy4zmVlL8A0LLjOZWUP4A4jNXW364/gISfXIi6uT9gdjynkgLCP4AqKSDhJ8c/gN4VmS9NzD9QSQEJP7nQP2Cjd0XmS9M/cP3tgY3e1T+AV2S+NHHYP5Cx2vrbA9s/kAtRN4OW3T/QsuM5lRTgP9jfHtjoXeE/4Axadjyn4j/oOZUUkPDjP+hm0LLjOeU/8JMLUTeD5j/4wEbvisznPwDugY3eFek/CBu9KzJf6j8ISPjJhajrPxB1M2jZ8ew/GKJuBi077j8gz6mkgITvPxR+ciHqZvA/lBSQ8JML8T8Yq62/PbDxP5xBy47nVPI/INjoXZH58j+kbgYtO57zPyQFJPzkQvQ/qJtBy47n9D8sMl+aOIz1P7DIfGniMPY/NF+aOIzV9j+09bcHNnr3PziM1dbfHvg/vCLzpYnD+D9AuRB1M2j5P8RPLkTdDPo/ROZLE4ex+j/IfGniMFb7P0wTh7Ha+vs/0KmkgISf/D9UQMJPLkT9P9TW3x7Y6P0/WG397YGN/j/cAxu9KzL/P2CaOIzV1v8/chirrb89AEC04zmVFJAAQPSuyHxp4gBANnpXZL40AUB4ReZLE4cBQLoQdTNo2QFA/NsDG70rAkA8p5ICEn4CQH5yIepm0AJAwD2w0bsiA0ACCT+5EHUDQETUzaBlxwNAhJ9ciLoZBEDGautvD2wEQAg2eldkvgRASgEJP7kQBUCMzJcmDmMFQMyXJg5jtQVADmO19bcHBkBQLkTdDFoGQJL50sRhrAZA1MRhrLb+BkAUkPCTC1EHQFZbf3tgowdAmCYOY7X1B0Da8ZxKCkgIQBy9KzJfmghAXIi6GbTsCECeU0kBCT8JQOAe2OhdkQlAIupm0LLjCUBktfW3BzYKQKSAhJ9ciApA5ksTh7HaCkAoF6JuBi0LQGriMFZbfwtArK2/PbDRC0DseE4lBSQMQC5E3QxadgxAcA9s9K7IDECy2vrbAxsNQPSlicNYbQ1ANHEYq62/DUB2PKeSAhIOQLgHNnpXZA5A+tLEYay2DkA8nlNJAQkPQH5p4jBWWw9AvjRxGKutD0AAAAAAAAAQQA==\",\"dtype\":\"float64\",\"shape\":[200]},\"y\":{\"__ndarray__\":\"fpA/sNuKIT9QmgKxDJYkP62R0eFyHig//90cCl82LD/5xnmGQHkwP51GYO8UNTM/rx0wXclbNj/7XPnbFPw5P6jaXO5xJj4/+FOEK6V2QT+OTIHQkjJEP8udgC1tUkc/9i3m9+DiSj+AegJ99/FOP8SRMHmZx1E/5EMJ59VlVD9VjVeAF11XP0KLLXGut1o/i999TeSAXj9WXwrxhWJhPzU/UX7IyGM/RsPh/4N6Zj8PmUGJnn9pP3jfqFqh4Gw/2eU24F9TcD+A122z7m1yP3MKpP9JxXQ/WzFZCBxfdz9ZFfJlbkF6P3zfeFyrcn0/WvSMRM98gD9bTNtnum6CP+cgt7jdkoQ/9LhiGi/thj/52hNE0oGJPzlw88MWVYw/E9/OaHVrjz9LlACExmSRP9f78EoPOpM/hGwERAQ4lT9ce/a8IGGXP1q47wbpt5k/Ly0IA+c+nD87HJ5JpviePy8KMvbX86A/48gc6kKHoj/J1ZXezjekP52cLuqtBqY/ZCmGfgb1pz8qJblh8AOqP6Bu54JxNKw/jC/FrnqHrj9+3P0U8n6wP1Fc9Rs1zLE/VpyQytQrsz/I79H2DZ60PyW84ysKI7Y/Rh+E7t26tz+eCGoNh2W5Pz/x2QLrIrs/XfvNbNXyvD/yuTOh9tS+PyvU5zFxZMA/hwqpYQdnwT+T5YiR6XHCPzeXWruzhMM/cxRp2vGexD8f7ZvCH8DFP0M7lRep58Y/RUCtZekUyD+UoGZeLEfJP/oQozmufco/zmqEPJy3yz8HaYJmFfTMP53d0kQrMs4/NjzS6+Jwzz8lJ9MKm1fQPy0mbzIK9tA/UAuuZDKT0T+MzK2Cgy7SP6ep/Kdqx9I/tasUCVNd0z8+VlPdpu/TPzLOgFLQfdQ/NebNhzoH1T9wVAqOUovVP5mQtGqICdY/9R9sG1CB1j9e9T2XIvLWP254Oct+W9c/ixG8j+q81z9b9OiT8xXYP2Vd0zwwZtg/0n36dUCt2D/yqddwzurYPx6vaFGPHtk/3v7QxUNI2T+p/2OHuGfZP+qSqMPGfNk/8OIobFSH2T/w4ihsVIfZP+qSqMPGfNk/qf9jh7hn2T/e/tDFQ0jZPx6vaFGPHtk/8qnXcM7q2D/Sffp1QK3YP2Vd0zwwZtg/W/Tok/MV2D+JEbyP6rzXP2x4Oct+W9c/XvU9lyLy1j/1H2wbUIHWP5mQtGqICdY/blQKjlKL1T8z5s2HOgfVPzLOgFLQfdQ/PlZT3abv0z+1qxQJU13TP6Wp/Kdqx9I/isytgoMu0j9QC65kMpPRPy0mbzIK9tA/JSfTCptX0D8xPNLr4nDPP5jd0kQrMs4/DGmCZhX0zD/OaoQ8nLfLP/oQozmufco/kaBmXixHyT9BQK1l6RTIP0c7lRep58Y/H+2bwh/AxT9zFGna8Z7EPzOXWruzhMM/j+WIkelxwj+KCqlhB2fBPyvU5zFxZMA/8rkzofbUvj9Y+81s1fK8Pzrx2QLrIrs/owhqDYdluT9GH4Tu3bq3PyW84ysKI7Y/xO/R9g2etD9RnJDK1CuzP1Zc9Rs1zLE/ftz9FPJ+sD+ML8WueoeuP6Bu54JxNKw/JCW5YfADqj9dKYZ+BvWnP52cLuqtBqY/ydWV3s43pD/jyBzqQoeiPysKMvbX86A/MxyeSab4nj8vLQgD5z6cP1q47wbpt5k/XHv2vCBhlz9/bAREBDiVP9L78EoPOpM/S5QAhMZkkT8T385odWuPPzlw88MWVYw/+doTRNKBiT/ruGIaL+2GP+cgt7jdkoQ/W0zbZ7pugj9a9IxEz3yAP3zfeFyrcn0/TRXyZW5Bej9bMVkIHF93P3MKpP9JxXQ/gNdts+5tcj/Z5TbgX1NwP2nfqFqh4Gw/D5lBiZ5/aT9Gw+H/g3pmPzU/UX7IyGM/Vl8K8YViYT97331N5IBeP0+LLXGut1o/VY1XgBddVz/kQwnn1WVUP8SRMHmZx1E/eHoCfffxTj8FLub34OJKP8udgC1tUkc/jkyB0JIyRD/4U4QrpXZBP5jaXO5xJj4/CF352xT8OT+vHTBdyVs2P51GYO8UNTM/+cZ5hkB5MD/y3RwKXzYsP6GR0eFyHig/UJoCsQyWJD9+kD+w24ohPw==\",\"dtype\":\"float64\",\"shape\":[200]}},\"selected\":{\"id\":\"2720\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"2719\",\"type\":\"UnionRenderers\"}},\"id\":\"2692\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"2720\",\"type\":\"Selection\"},{\"attributes\":{\"line_color\":\"#ff7f0e\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2664\",\"type\":\"Line\"},{\"attributes\":{\"source\":{\"id\":\"2663\",\"type\":\"ColumnDataSource\"}},\"id\":\"2667\",\"type\":\"CDSView\"},{\"attributes\":{\"line_alpha\":0.1,\"line_color\":\"#1f77b4\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2694\",\"type\":\"Line\"},{\"attributes\":{},\"id\":\"2646\",\"type\":\"BasicTicker\"},{\"attributes\":{\"axis_label\":\"K(x0, x)\",\"formatter\":{\"id\":\"2670\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"2646\",\"type\":\"BasicTicker\"}},\"id\":\"2645\",\"type\":\"LinearAxis\"},{\"attributes\":{\"axis_label\":\"t = |x-x0|/\\u03bb\",\"formatter\":{\"id\":\"2672\",\"type\":\"BasicTickFormatter\"},\"ticker\":{\"id\":\"2641\",\"type\":\"BasicTicker\"}},\"id\":\"2640\",\"type\":\"LinearAxis\"},{\"attributes\":{\"line_color\":\"#2ca02c\",\"line_width\":2,\"x\":{\"field\":\"x\"},\"y\":{\"field\":\"y\"}},\"id\":\"2678\",\"type\":\"Line\"},{\"attributes\":{},\"id\":\"2719\",\"type\":\"UnionRenderers\"},{\"attributes\":{},\"id\":\"2655\",\"type\":\"HelpTool\"},{\"attributes\":{\"data_source\":{\"id\":\"2677\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"2678\",\"type\":\"Line\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"2679\",\"type\":\"Line\"},\"selection_glyph\":null,\"view\":{\"id\":\"2681\",\"type\":\"CDSView\"}},\"id\":\"2680\",\"type\":\"GlyphRenderer\"},{\"attributes\":{},\"id\":\"2638\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"2689\",\"type\":\"Selection\"},{\"attributes\":{},\"id\":\"2706\",\"type\":\"Selection\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"2650\",\"type\":\"PanTool\"},{\"id\":\"2651\",\"type\":\"WheelZoomTool\"},{\"id\":\"2652\",\"type\":\"BoxZoomTool\"},{\"id\":\"2653\",\"type\":\"SaveTool\"},{\"id\":\"2654\",\"type\":\"ResetTool\"},{\"id\":\"2655\",\"type\":\"HelpTool\"}]},\"id\":\"2656\",\"type\":\"Toolbar\"},{\"attributes\":{\"dimension\":1,\"ticker\":{\"id\":\"2646\",\"type\":\"BasicTicker\"}},\"id\":\"2649\",\"type\":\"Grid\"}],\"root_ids\":[\"2631\"]},\"title\":\"Bokeh Application\",\"version\":\"1.4.0\"}};\n",
" var render_items = [{\"docid\":\"7826c7a6-4eb7-4e8b-8cd8-f26d537e7e31\",\"roots\":{\"2631\":\"c336d6e8-23a6-46b4-8b2b-16c72c46e1ed\"}}];\n",
" root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
"\n",
" }\n",
" if (root.Bokeh !== undefined) {\n",
" embed_document(root);\n",
" } else {\n",
" var attempts = 0;\n",
" var timer = setInterval(function(root) {\n",
" if (root.Bokeh !== undefined) {\n",
" clearInterval(timer);\n",
" embed_document(root);\n",
" } else {\n",
" attempts++;\n",
" if (attempts > 100) {\n",
" clearInterval(timer);\n",
" console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
" }\n",
" }\n",
" }, 10, root)\n",
" }\n",
"})(window);"
],
"application/vnd.bokehjs_exec.v0+json": ""
},
"metadata": {
"application/vnd.bokehjs_exec.v0+json": {
"id": "2631"
}
},
"output_type": "display_data"
}
],
"source": [
"# Set up time points\n",
"t = np.linspace(-4.0, 4.0, 200)\n",
"dt = t[1] - t[0]\n",
"\n",
"# Compute unnormalized kernels\n",
"K_epan = epan_kernel(t)\n",
"K_tri_cube = tri_cube_kernel(t)\n",
"K_gauss = gauss_kernel(t)\n",
"\n",
"# Adjust to approximately integrate to unity for ease of comparison\n",
"K_epan /= K_epan.sum() * dt\n",
"K_tri_cube /= K_tri_cube.sum() * dt\n",
"K_gauss /= K_gauss.sum() * dt\n",
"\n",
"# Plot kernels; easier in base BOkeh\n",
"p = bokeh.plotting.figure(plot_height=300,\n",
" plot_width=450,\n",
" x_axis_label='t = |x-x0|/λ',\n",
" y_axis_label='K(x0, x)')\n",
"p.line(t, K_epan, color='#ff7f0e', line_width=2, legend_label='Epanechnikov')\n",
"p.line(t, K_tri_cube, color='#2ca02c', line_width=2, legend_label='Tri-cube')\n",
"p.line(t, K_gauss, color='#9467bd', line_width=2, legend_label='Gaussian')\n",
"bokeh.io.show(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the scale of the $x$-axis in this plot, we see that the parameter $\\lambda$ sets the width of the smoothing kernel. The bigger $\\lambda$ is, the stronger the smoothing and the more fine detail will be averaged over. For this reason, $\\lambda$ is sometimes called the **bandwidth** of the kernel. The Gaussian kernel is broad and has no cutoff like the Epanechnikov or tri-cube kernels. If therefore uses all of the data points.\n",
"\n",
"We will now write a function to give the smoothed value of a function for a given kernel using Nadaraya-Watson."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def nw_kernel_smooth(x_0, x, y, kernel_fun, lam):\n",
" \"\"\"\n",
" Gives smoothed data at points x_0 using a Nadaraya-Watson kernel \n",
" estimator. The data points are given by NumPy arrays x, y.\n",
" \n",
" kernel_fun must be of the form\n",
" kernel_fun(t), \n",
" where t = |x - x_0| / lam\n",
" \n",
" This is not a fast way to do it, but it simply implemented!\n",
" \"\"\"\n",
" \n",
" # Function to give estimate of smoothed curve at single point.\n",
" def single_point_estimate(x_0_single):\n",
" \"\"\"\n",
" Estimate at a single point x_0_single.\n",
" \"\"\"\n",
" t = np.abs(x_0_single - x) / lam\n",
" return np.dot(kernel_fun(t), y) / kernel_fun(t).sum()\n",
" \n",
" # If we only want an estimate at a single data point\n",
" if np.isscalar(x_0):\n",
" return single_point_estimate(x_0)\n",
" else: # Get estimate at all points\n",
" return np.array([single_point_estimate(x) for x in x_0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All right, we're all set. Let's try applying some kernels to our spike"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.I :Curve [t (s)] (V (µV))\n",
" .Scatter.I :Scatter [t (s)] (V (µV))\n",
" .Curve.Epanechnikov :Curve [x] (y)\n",
" .Curve.Tri_hyphen_minus_cube :Curve [x] (y)\n",
" .Curve.Gaussian :Curve [x] (y)"
]
},
"execution_count": 24,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "2802"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Specify time points\n",
"ts_kernel = np.linspace(ts[0], ts[-1], 400)\n",
"\n",
"# Averaging window (0.2 ms)\n",
"lam = 0.0002\n",
"\n",
"# Compute smooth curves with kernel\n",
"Vs_epan = nw_kernel_smooth(ts_kernel, ts, Vs, epan_kernel, lam)\n",
"Vs_tri_cube = nw_kernel_smooth(ts_kernel, ts, Vs, tri_cube_kernel, lam)\n",
"Vs_gauss = nw_kernel_smooth(ts_kernel, ts, Vs, gauss_kernel, lam)\n",
"\n",
"# Plot results\n",
"plot = (\n",
" base_spike_plot()\n",
" * hv.Curve((ts_kernel, Vs_epan), label=\"Epanechnikov\").opts(\n",
" color=\"#ff7f0e\", line_width=2\n",
" )\n",
" * hv.Curve((ts_kernel, Vs_tri_cube), label=\"Tri-cube\").opts(\n",
" color=\"#2ca02c\", line_width=2\n",
" )\n",
" * hv.Curve((ts_kernel, Vs_gauss), label=\"Gaussian\").opts(\n",
" color=\"#9467bd\", line_width=2\n",
" )\n",
")\n",
"\n",
"plot.opts(legend_position='bottom_right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All of these dull the spike. In fact, this is generally what you will find with smoothing methods."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The simple option\n",
"\n",
"Why don't we just fit a parabola to the bottom three points and then compute the minimum that way? If we approximate $V(t) = at^2 + bt + c$ at the bottom of the spike, the minimum is at $t = -b/2a$. We can do this easily using `np.polyfit()`."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"