Tutorial_063_Exercises_Answers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In class exercises for 11.21.2018 (Filtering)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First import data from eeg_data00.npz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first import the modules that you'll need. \n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# import our sinc filter implementation function...\n",
    "import my_sinc_filt as msf\n",
    "\n",
    "# also define the default font we'll use for figures. \n",
    "fig_font = {'fontname':'Arial', 'size':'20'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Just like in the first FFT tutorial load eeg_data00.npz so that we have an EEG data set to play with\n",
    "* you should find that it has three arrays: data, tx, and sr (with 'data' a #trials x # timepoints array containing the EEG data, the time of each data acquisition in tx, and the sample rate stored as sr).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of the big eeg data set:  (960, 876)\n",
      "Sample rate:  250\n"
     ]
    }
   ],
   "source": [
    "# load the data...\n",
    "eeg = np.load('eeg_data00.npz')\n",
    "\n",
    "# get the different arrays like this...kind of like a structure in matlab. \n",
    "eeg['data']\n",
    "eeg['sr']\n",
    "eeg['tx']\n",
    "\n",
    "# and can query the attributes of the data like this...which will tell us that there 960 trials and 876 timepoints per trial sampled\n",
    "# at 250Hz\n",
    "print('Shape of the big eeg data set: ', eeg['data'].shape)\n",
    "print('Sample rate: ', eeg['sr'])\n",
    "\n",
    "# and if you want to save some typing, especially because we only have a few variables, you reassign the different arrays like this\n",
    "data = eeg['data']\n",
    "sr = eeg['sr']\n",
    "tx = eeg['tx']\n",
    "# now you can just write data instead of eeg['data'] (and can for tx and sr as well...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract the data from time 0 until time 2000\n",
    "* hint: you can use boolean indexing by selecting the columns from the 'data' matrix that meet: (tx>=0)  & (tx<2000)\n",
    "* while you're at it, go ahead and grab the timepoints corresponding to that interval as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "seg_data = data[:,(tx>=0)  & (tx<2000)]\n",
    "seg_time = tx[(tx>=0) & (tx<2000)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then FFT each trial and take a look at the average power across frequencies\n",
    "* FFT each trial \n",
    "* Then compute the amplitude on each trial\n",
    "* Then compute the average amplitude across all trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAETCAYAAADH1SqlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcVNWd9/HPr7uBbgFlDaKIoMElijTY44YG1GhkUeNEXGISTMxDoskzZHklE8dJBOOTx6iJieM2uIHLaNQoOq4xjoi7gCIgIIKAogiIIKDQ0PRv/ji3oGmr6nZ119bV3/frdV9VdZdzT3V117fPufeea+6OiIhIOmWFroCIiBQ/hYWIiMRSWIiISCyFhYiIxFJYiIhILIWFiIjEUliIiEgshYWIiMRSWIiISKyKQlcgW3r06OH9+vUrdDVERFqVWbNmfezuPePWK5mw6NevHzNnzix0NUREWhUzW96U9dQNJSIisRQWIiISS2EhIiKxSuaYhYgUzrZt21ixYgVbtmwpdFUkhcrKSvr06UO7du2atb3CQkRabMWKFXTu3Jl+/fphZoWujjTi7qxdu5YVK1bQv3//ZpVR0G4oM6s0s9fM7E0ze8vMJkbz+5vZq2b2jpn91czaF7KeIpLeli1b6N69u4KiSJkZ3bt3b1HLr9DHLGqBE9x9EFANnGJmRwF/AK5x9wHAOuCCAtZRRJpAQVHcWvr5FDQsPNgUvWwXTQ6cADwQzZ8CfCNXdVi1CubOzVXpIiKlodAtC8ys3MxmA6uBp4ElwHp3r4tWWQHsnWLbcWY208xmrlmzpln7nz4dLrusWZuKiOwwfPjwHRcGjxw5kvXr17eovGnTpjF69Oiky84991wOO+wwrrnmGhYuXEh1dTWDBw9m1qxZ3HDDDS3abyoFDwt33+7u1UAf4Ajg4GSrpdh2krvXuHtNz56xV6snVVYG9fXN2lRESkRdXV38Shl4/PHH6dKlS1bLTPjoo4946aWXmDNnDj/72c+YOnUqp59+Om+88Qbdu3fPWVgUzdlQ7r7ezKYBRwFdzKwial30AT7M1X4VFiKt37JlyxgxYgTHHnssL730EnvvvTcPP/wwVVVVzJ49mx/96Ed8/vnn7L///tx222107dqV4cOHc8wxx/Diiy9y2mmnMXfuXKqqqli4cCHLly/n9ttvZ8qUKbz88ssceeSRTJ48GYALL7yQGTNmsHnzZs4880wmTpz4hfokhh964IEHuOmmmwD49NNP6devH88++yx///vfufTSS6mtrWX//ffn9ttvp1OnTjz55JP89Kc/pUePHgwZMiTpez355JNZvXo11dXVnHHGGdx4442Ul5czffp0evXqxZIlS6iuruakk07iqquuyt4P2d0LNgE9gS7R8yrgeWA0cD9wTjT/JuCiuLIOP/xwb46pU91PPbVZm4pIZP78+QXd/9KlS728vNzfeOMNd3cfM2aM33nnne7uPnDgQJ82bZq7u//mN7/x8ePHu7v7sGHD/MILL9xRxtixY/3ss8/2+vp6nzp1qnfu3NnnzJnj27dv9yFDhuwoe+3ate7uXldX58OGDfM333xzR3kzZsxwd/d9993X16xZs6PsrVu3+rHHHuuPPPKIr1mzxo877jjftGmTu7tfccUVPnHiRN+8ebP36dPHFy1a5PX19T5mzBgfNWpU0vd6yCGH7Hh96aWX+lVXXZV0WWPJPidgpjfh+7rQ3VC9gWfNbA4wA3ja3R8F/hX4uZktBroDt+aqAmpZiGSfWfanOP3796e6uhqAww8/nGXLlvHpp5+yfv16hg0bBsDYsWOZPn36jm3OPvvsXco49dRTMTMGDhxIr169GDhwIGVlZRxyyCEsW7YMgPvuu48hQ4YwePBg3nrrLebPnx9bt/Hjx3PCCSdw6qmn8sorrzB//nyGDh1KdXU1U6ZMYfny5SxcuJD+/fszYMAAzIxvf/vbTfxp50dBu6HcfQ4wOMn8dwnHL3KuvFxhIZJtnvQoY2516NBhx/Py8nI2b94cu03Hjh2TllFWVrZLeWVlZdTV1bF06VKuvvpqZsyYQdeuXTn//PNjr12YPHkyy5cv57rrrgNCb85JJ53EPffcs8t6s2fPLurTjwvdsig4tSxEStcee+xB165def755wG48847d7QymmPDhg107NiRPfbYg1WrVvHEE0+kXX/WrFlcffXV3HXXXZSVha/bo446ihdffJHFixcD8Pnnn7No0SIOOuggli5dypIlSwC+ECZN0blzZzZu3Jjxdk1RNAe4C0VhIVLapkyZsuMA93777cftt9/e7LIGDRrE4MGDOeSQQ9hvv/0YOnRo2vWvu+46PvnkE44//ngAampquOWWW5g8eTLnnnsutbW1AFx++eUccMABTJo0iVGjRtGjRw+OPfZY5s2bl1H9unfvztChQzn00EMZMWJEVg9wmxeivZgDNTU13pybH/3jH3DFFeFRRJpnwYIFHHxwsrPepZgk+5zMbJa718Rtq24otSxERGIpLBQWIiKxFBYKC5GsKJUu7VLV0s9HYVEG27cXuhYirVtlZSVr165VYBQpj+5nUVlZ2ewydDaUWhYiLdanTx9WrFhBcwf0lNxL3Cmvudp8WOiiPJGWa9euXbPvwCatg7qh1LIQEYmlsFBYiIjEUlgoLEREYiksFBYiIrEUFgoLEZFYCguFhYhILIWFLsoTEYmlsFDLQkQkVpsPC12UJyISr82HhVoWIiLxFBYKCxGRWAoLhYWISCyFhcJCRCSWwkJhISISS2GhsBARiVXQsDCzfczsWTNbYGZvmdn4aP4EM/vAzGZH08hc1UEX5YmIxCv0zY/qgF+4++tm1hmYZWZPR8uucferc10BtSxEROIVNCzcfSWwMnq+0cwWAHvnsw66KE9EJF7RHLMws37AYODVaNZPzGyOmd1mZl1TbDPOzGaa2czm3vtXLQsRkXhFERZm1gn4G/BTd98A3AjsD1QTWh5/TLadu09y9xp3r+nZs2ez9q2wEBGJV/CwMLN2hKC4290fBHD3Ve6+3d3rgZuBI3K1f4WFiEi8Qp8NZcCtwAJ3/1OD+b0brHYGMC9XdVBYiIjEK/TZUEOB7wBzzWx2NO/fgHPNrBpwYBnww1xVQGEhIhKv0GdDvQBYkkWP56sOCgsRkXgFP2ZRaLooT0QknsKiDNzDJCIiybX5sDALk8JCRCS1Nh8WoOMWIiJxFBYoLERE4igsUFiIiMRRWKCwEBGJo7BAYSEiEkdhgcJCRCSOwgJdmCciEkdhgW6AJCISR2GBuqFEROIoLFBYiIjEUVigsBARiaOwQGEhIhJHYYHCQkQkjsIChYWISJxmhYWZdTKzo81sVPS6a3arlV8KCxGR9DIKCzPrZWb3Ap8ALwAPR4suMrPFZnZctiuYD7ooT0QkvSaHhZn1BF4GzgJeBV5n5/2zPwP2BZ4ws4HZrmSu6aI8EZH0MmlZTAT6Aqe5+3HAo4kF7v5n4CSgAvj3rNYwD9QNJSKSXiZhcRrwoLs/mmyhu08DHgSOzkK98kphISKSXiZh0QN4N2adFUDP5lenMBQWIiLpZRIWK4AhMescGa3XqigsRETSyyQsHgBONLMfJltoZj8HjgUeamqBZraPmT1rZgvM7C0zGx/N72ZmT5vZO9FjTk/NVViIiKSXSVj8P2A+cIOZzSGcFYWZTTazucBVwGLg9xmUWQf8wt0PBo4CfmxmXwF+DTzj7gOAZ6LXOaOwEBFJr8lh4e4bgaHATUA/4GDCqbPfBb4M3AkMdff1GZS50t1fb1D+AmBv4HRgSrTaFOAbTS2zORQWIiLpVWSysrtvIPz3/y/AgUAXYBPwtrvXtqQiZtYPGEy4hqOXu6+M9rnSzL7UkrLj6KI8EZH0MgqLBHffTuiSygoz6wT8Dfipu28ws7hNEtuNA8YB9O3bt9n710V5IiLppQwLM/ttM8t0d/9dU1c2s3aEoLjb3R+MZq8ys95Rq6I3sDrFjiYBkwBqamq8mfVVN5SISIx0LYsJSeYlvpCT/evv0XwHmhQWFpoQtwIL3P1PDRY9AowFrogeH06yedYoLERE0ksXFmc0el0F/IkQCNcALxEGFOwE/BPwc2A78P0M9j8U+A4w18xmR/P+jRAS95nZBcB7wJgMysyYwkJEJL2UYeHuu/w3b2ZXA+2AI9x9aaPVXzOzBwmDC36HMCJtLHd/geStFIATm1JGNigsRETSy+Q6i28TxoZqHBRAOGuJMDbUmdmoWD4pLERE0sskLCqJP3tqd1K3FIqWwkJEJL1MwmIW8E0zOzjZQjM7mtCqmJ6NiuWTwkJEJL1MrrOYAPwDeMXMJhPCYyOwB+FA9XnAVlrp/Sx0UZ6ISGpNDgt3f97MTgduAP4vO0+jhdD1NB/4vrvPy24Vc08X5YmIpJfpcB9PmtmXCUORDwK6AuuAWe4+Iwf1ywt1Q4mIpJfxcB/uXk+4F/fL2a9OYSgsRETSa3JYRIMHNoW7+380sz4FobAQEUkvk5bFn9k5pEdjDYcBcUBhISJSQjIJi++lmL8b4X4W3yEc5P5ZSyuVbwoLEZH0Mjkbakq65WZ2LfAGcHz02GooLERE0svkory03H05cB/wo2yVmS8KCxGR9LIWFpEtwD5ZLjPndFGeiEh6WQuLaBiQbwHLslVmvuiiPBGR9DI5dfb1VIsI97ToTwifJt8lr1ioG0pEJL1MzoaqTrNsKzAPuNndr29ZlfJPYSEikl4mZ0Nl+/hG0VBYiIik1+QAMLPvmtlhMesMNbPftLxa+aWwEBFJL5PWwmTg9Jh1/hm4uNm1KRCFhYhIeim7oczsfOC0RrPPNbPBKTZpT7ggb212qpY/CgsRkfTSHbN4EriWcKYThDGfDoqmVGqB32anavmjsBARSS9lWLj7R2a2P2HsJwPeJQwm+JdkqwPbgI/dfVsuKppLuihPRCS9tGdDufuaxHMz+x4wOxrWo6ToojwRkfSyNpBga6ZuKBGR9FKeDWVmn5jZLxu9bsrU5APcZnabma02s3kN5k0wsw/MbHY0jWz+22sahYWISHrpWhYbCAesG772FOs212TgOuCORvOvcfers7yvlBQWIiLppTvA3S/d62xw9+lmlvVyM6WwEBFJr1iH8PiJmc2Juqm65npnCgsRkfTSXZTX+IK8JnP3R5q7LXAjYeRajx7/CHw/2YpmNg4YB9C3b99m71BhISKSXrpjFlPJ/BiFRduUN7dC7r5qR2FmNwOPpll3EjAJoKamptnHUxQWIiLppQuLy8j+Ae1YZtbb3VdGL88gDH2eU2VlsK3VXUooIpI/6Q5wT8j1zs3sHmA40MPMVgCXAsPNrJoQVMuAH+a6HuXlUFsbv56ISFuVyc2PADAzA44CBgF7AB8Dr7p7xi0Adz83yexbMy2npdQNJSKSXkZhYWZfIxyA3i8xK3p0M5sJjHP3N7NYv7xQWIiIpJfJPbiPIRxsLgPuBl4DPgK6AMcA3wL+x8yOcfe3c1DXnFFYiIikl0nLYiJQBxzv7jMaLbvFzG4BniEcGD87S/XLC4WFiEh6mVyU90/A/UmCAgB3fwl4ADghGxXLJ4WFiEh6mYTFVmBjzDobMiyzKCgsRETSy+SL/R7gvFRjOZnZnoR7cN/b8mrll8JCRCS9TI5Z3AwMBd4ws2uB6cAHQBWhi+qXhCu3n2s8VEgLh//IOd0pT0QkvUzCYg7hQjkDfsMXr+5OnEZ7T6N5LRr+Ix90pzwRkfQyCYuCDP+RD+qGEhFJL5Pbqk7IYT0KSmEhIpJeqztzKRcUFiIi6WU63MephCu1+wEdUqzm7n54C+uVVwoLEZH0Mhnu4yLgP9h5IDuVVndcQ2EhIpJeJt1Q4wkjzA4Hqty9LMVU1Gc+JaOwEBFJL5Ow6APc5e7T3b2k7v6gsBARSS+TsFgAfClXFSkkXZQnIpJeJmFxGXCWmZ2Sq8oUii7KExFJL5PrLB6Jhvl4zMzeBpYCybqj3N2/ma0K5oO6oURE0svkbKjzgJ8TzoY6KJqS0dlQjTz2GIwcCRZ3HpmISJHK5DqLfwc+JwTGi8BnOalRAeQyLOrrYfRo+Phj6N49N/sQEcm1TMKiH3CLu9+co7oUTC7D4rMoUt9/X2EhIq1XJge4lwLtc1WRQsplWGzaFB7ffz835YuI5EMmYXENcI6ZVeeqMoWSr5aFiEhrlUk31GfAEuA1M3seWEzy4xbu7r/IRuXyRS0LEZH0MgmLuxo8Pz6aknGgSWFhZrcBo4HV7n5oNK8b8FfCMZJlwFnuvi6DemYslxflJcLivfdyU76ISD5kEhapwqGxQzMoczJwHXBHg3m/Bp5x9yvM7NfR63/NoMyM5fKivM8+gw4d1LIQkdYtk4vynku1zMw6AucAPyDcj/v6JpY53cz6NZp9OmGwQoApwDRyHBa57oY68ECFhYi0bi26+ZGZHWlmNwMrgUnAkcCGFtapl7uvBIgecz4eVT7C4sMPdZW4iLReGYeFmXUzs/FmNhd4CbgA2A14BjgP6J3dKqatyzgzm2lmM9esWdPscnJ9NlSPHrDHHrBqVW72ISKSa5kM93EioZvpG4TrLRKDVzwLnO/u2epoWWVmvd19pZn1BlanWtHdJxFaNNTU1DR7mJFctyw6dYJ99gldUb3zFqUiItmTtmVhZnuZ2SVmtgT4O3A28ClwAzA0Wu3tLAYFwCPA2Oj5WODhLJadVK7DomNH2H33nWdGiYi0NilbFmb238DXo3U+JZyxdC/wtLvXR+u0aOdmdg/hYHYPM1sBXApcAdxnZhcA7wFjWrSTJsh1N9Ree0FVFWzenJt9iIjkWrpuqFGEi+5+B1yZi7vjufu5KRadmO19pZOPbqjKStiyJTf7EBHJtXTdUE8BlcAE4AMzu8fMvmFmJTc+VK4vyuvYMbQsFBYi0lqlDAt3H0G47/avgA8Ixyv+Bqw2s9vM7OT8VDH3cn1RXqJloW4oEWmt0h7gdvdV7v5Hdx8EDAL+TLinxfnAE4ShPYaY2dG5rmgu5eMAt7qhRKQ1a/J1Fu4+NxogsA8wkjB+0xbgCOAFM1tiZpeZ2YG5qWru5OOYhbqhRKQ1y/iiPHevd/cn3f1bQC/CtRfPEwb++3fgrazWMA9yfTaUuqFEpLVr0XAf7r7J3W9z9+FAf+C3wDvZqFg+qRtKRCS9FoVFQ+7+nrtf7u4HZ6vMfFE3lIhIelkLi9Ys191QiZaFuqFEpLVSWJC76yy2bw+tiaoqdUOJSOumsCB3YZFoVZSVqRtKRFo3hQWw22656SJKhAWoG0pEWjeFBdClC6xfD97sQc6TSxzcBnVDiUjrprAA2rUL3UQbN2a33MRps6BuKBFp3RQWka5dYd267Ja5fn0oF9SyKHXLl8OECYWuhUjuKCwiuQiLTz7ZNSx0zKJ0vfMOPPFEoWshkjsKi0guwmLdOujWLTxXN1Rpq60Nk0ipUlhEctWySISFuqFKm8JCSp3CIpKPsFA3VOlSWEipU1hEunTJbVioG6q0KSyk1CksIuqGkpZQWEipU1hEcn2Au3172LYtd/f6lsLaskVhIaVNYRHJdcvCLLQu9IVSmtSykFKnsIjk+joLUFdUKautDS3HXA11L1JoCotIrlsWoDOiSlmiVaHWhZSqikJXIBUzWwZsBLYDde5ek8v9ZTsstm0LrYjOnXfO0xlRpathWFRVFbYuIrlQtGEROd7dP87HjrIdFuvWhTLNds5TN1TpUstCSp26oSJdu2Z3mPLGXVCgbqhSprCQUlfMYeHA381slpmNy/XO2rcP06ZN2SkvWVioG6p0JT5XhYWUqmLuhhrq7h+a2ZeAp81sobtPb7hCFCLjAPr27dviHfbqBStX7nqcoblStSwUFqVJLQspdUXbsnD3D6PH1cBDwBFJ1pnk7jXuXtOzZ88W73PgQJgzp8XFAKlbFuqGKk0KCyl1RRkWZtbRzDonngMnA/Nyvd/qapg9OztlrVoFjfNLLYvSpbCQUleUYQH0Al4wszeB14DH3P3JXO900CB4883slLVkCey//67zFBalKxES+nylVBXlMQt3fxcYlO/9ZrNlsXgx/PM/7zpPB7hLV21tuN+6WhZSqoq1ZVEQ/frBhg2wdm3Ly1q8GL785V3n6dTZ0lVbC7vvrrCQ0qWwaKCsDA47DN54o2Xl1NaGs6oan6DVrRusWdOysqU4bdmisJDSprBo5LjjYNq0lpWxbFkIiopGnXwDBsA777SsbClOallIqVNYNHLKKfBkCw+lJ+uCAjjgAFi0qGVlS3FSWEipU1g0cvTR4Uym1aubX0aqsEi0LLI1pIgUj9pa2GMPhYWULoVFI+3awQknwFNPNb+MxYu/eNoshPt8V1XBRx81v2wpTmpZSKlTWCQxZgzceGPzWwAvvQRDhiRfpq6o0pQIC50aLaVKYZHEmDHw2WcwdWrm265eHbqxjj46+fIDDtBB7lLjrm4oKX0KiyTKy+EPf4CLL4a6usy2ffppGD48dGclM2CAWhalJvE7sttuCgspXQqLFL7+ddh7b7j11sy2e+qpsG0qhxwCr7/esrpJcamthQ4dwqSwkFKlsEjBDK68EiZObPoV3R9+CI8/DiNHpl7nxBNhxgz4OC/3/5N8qK0NV+crLKSUKSzSOPxw+Na34Pzzdx7s3rgxBEJj7nDhhXDRRbDvvqnL7NgxXMvx4IM5qbIUgFoW0hYoLGL8/vehFXDZZeFMlzPOgNNOg7ff3nW9xx4L8y65JL7Ms86Ce+/NTX0l/xQW0hYoLGK0bw8PPQS33w69e8Nee8FvfxvCI6GuDn75S7j66vCFEWfkSFi6NBwMl9ZPYSFtQVEOUV5s9twTZs0KXwR77RW6og46CB54AEaNgh/8IIwFNWpU08qrqoLrrw/dVrNnQ6dOua2/5JbCQtoCtSyaqHv3EBQQ7tH92GPh+ESvXmG02oceCgfFm2rkyHCl+Nln73p67mefhfA580z4xz+y+x4kN7Zs2RkWuihPSpXCopmqq2HuXFi+HO68M5xjn6nrrw8j01ZXwx13wHPPhQCpqwutlPPOg/vuy37dJbt0NpS0BeqGaoFevVq2fbt24Srxxx+Hm24KxzEuuih0T5mFEDnzTLj/fvjqV+ErXwlXhjcnmCR3Et1QlZUKCyldCosCMwutiGTHOwYPDsc07rgD5s2Dv/41PA4eHLrCtm2DrVvD40knwc9/Hk7NlfzSMQtpCxQWRa5zZ/jxj3e+Xr06BMbGjeFMrcSwIv/5n3DEEeFYSr9+Balqm6WwkLZAYdHKfOlL4bhGY1/7GvzlL6Gb6rbbYMSI/NetrSrWsLjyShg9OnRfirSUwqKEjB8fjnN873thBNSxY8MpvhUVMHBg6FMvKwsDJSaeS8sVY1i4wx//GLooFRaSDQqLEjNsWLj50rRp4VjHE0/A5s3w1lvhi6O+PkzbtoVTgQ8/PEx9+4YRcQ86KHR9lZcX+p20Hg1PnS2WsFi0KHRZTp/etFEFROIoLEpQWVnoqkrWXZWwbRssWxZGwH399XDf8WuuCUHz2WfhjKuuXaFHj9D1VVkZBlTca68w9tW++8I++0C3buFU306dwp0AKyrCQfe6unA8pX37L07t2pVWq6YYT52dPj1cy/PCC+GzqNBfurRQ0f4KmdkpwF+AcuAWd7+iwFUqKe3ahZbEgAHhwsCG6uth0yb45JMQEKtWhf+eu3WDlSvDtSXz5oVTftetC2Vt2hSeb98eAqGiInxJbd2afKqoSB4kyYJl69bwuPvuYd7GjSFsKip2Tu3apX+9fXvomuvZM2zbcNq+PbS+tm0LX/hVVbtO7duH97d+/c7la9aEei1aFFpw114b5kO4D8qwYTv3VV6+62OqeYmpoiIEr1noTnLPPFynTw/jmC1fHs6oq6nJzu+NtF3mzb13aA6ZWTmwCDgJWAHMAM519/mptqmpqfGZM2fmqYbSEu7pg6ThtG1b+OLftg02bNh5+9JEGXV1YVnieePXiefl5eHL/uOPw7aJ7rj6+vBFXFUV9lNbG4Kj4bR1a2g5Je6Et3lzCIL27UMLa8yYcAdECINJ3nprGB5m/foQRPX1X3xMPN++/Yv12bo1hLNZ2J9Z2H8iLBuGYIcOoeW3aVMoq3PnMO/ZZ8NQ+H/+M8ycCQceGFqLVVVf7GI0Cy2jysqw3Cy0Ljt3DqFVXr6zvu47u9wqK8O6DbvhKivDz9Es+VRWtvN54h+Gioqwj4ZTsnmJYE3l00/D6eVDh4b7xkjTmNksd4/9d6JYw+JoYIK7fz16fTGAu///VNsoLKSUJIYNad8+PG7aFMIyEX6Jxy1b4KOPwhd7RUVYb8uW0IU4fHhY9uKL8PnnO6f6+l33VV8fttmyJQShe7heZ+PGna3FxJe12c4wq60N21ZW7py3eXOoV6JFlGyqr98Z9tu2hSkRnA2nurovzoOdgdPwsaIilHniifDqq6FuiRZmotszbtq+PYRkeXn4uSdCNTGMT7LHdMvysU6HDtCnT/inoLIy9tcqqaaGRbF2Q+0NvN/g9QrgyALVRSTvGv/h7757mDK1557wzW9mp07FIBE0DR/r60OwlJWFFlhdXQjWhmHUuPWWbCorCyFZX78zCBP/Syd7TLcsX+ts3gwrVjRttOuWKtawSDYk3xeaQGY2DhgH0Ldv31zXSUQKLNENle5svYqKcHxNsqtYz0lZAezT4HUf4MPGK7n7JHevcfeanj175q1yIiJtTbGGxQxggJn1N7P2wDnAIwWuk4hIm1WU3VDuXmdmPwGeIpw6e5u7v1XgaomItFlFGRYA7v448Hih6yEiIsXbDSUiIkVEYSEiIrEUFiIiEkthISIisYpyuI/mMLM1wPJmbt4D+DiL1SkEvYfCa+31B72HYpDv+u/r7rEXqpVMWLSEmc1sytgoxUzvofBae/1B76EYFGv91Q0lIiKxFBYiIhJLYRFMKnQFskDvofBae/1B76EYFGX9dcxCRERiqWUhIiKx2nxYmNkpZva2mS02s18Xuj5xzGwfM3vWzBaY2VtmNj6a383Mnjazd6LHroWuaxwzKzezN8zs0eh1fzN7NXoPf41GHC5aZtbFzB4ws4XR53F0a/qomLRZAAAJHElEQVQczOxn0e/QPDO7x8wqi/0zMLPbzGy1mc1rMC/pz9yCa6O/7TlmNqRwNd8pxXu4Kvo9mmNmD5lZlwbLLo7ew9tm9vXC1LqNh0V0r+/rgRHAV4Bzzewrha1VrDrgF+5+MHAU8OOozr8GnnH3AcAz0etiNx5Y0OD1H4BrovewDrigILVqur8AT7r7QcAgwntpFZ+Dme0N/AtQ4+6HEkZ3Pofi/wwmA6c0mpfqZz4CGBBN44Ab81THOJP54nt4GjjU3Q8DFgEXA0R/2+cAh0Tb3BB9b+Vdmw4L4Ahgsbu/6+5bgXuB0wtcp7TcfaW7vx4930j4gtqbUO8p0WpTgG8UpoZNY2Z9gFHALdFrA04AHohWKer3YGa7A18FbgVw963uvp7W9TlUAFVmVgHsBqykyD8Dd58OfNJodqqf+enAHR68AnQxs975qWlqyd6Du//d3euil68QbvgG4T3c6+617r4UWEz43sq7th4Wye71vXeB6pIxM+sHDAZeBXq5+0oIgQJ8qXA1a5I/A78C6qPX3YH1Df5giv2z2A9YA9wedaXdYmYdaSWfg7t/AFwNvEcIiU+BWbSuzyAh1c+8tf59fx94InpeNO+hrYdFk+71XYzMrBPwN+Cn7r6h0PXJhJmNBla7+6yGs5OsWsyfRQUwBLjR3QcDn1GkXU7JRP36pwP9gb2AjoRum8aK+TOI09p+pzCzSwhdzXcnZiVZrSDvoa2HRZPu9V1szKwdISjudvcHo9mrEk3s6HF1oerXBEOB08xsGaHr7wRCS6NL1CUCxf9ZrABWuPur0esHCOHRWj6HrwFL3X2Nu28DHgSOoXV9Bgmpfuat6u/bzMYCo4HzfOc1DUXzHtp6WLS6e31Hffu3Agvc/U8NFj0CjI2ejwUeznfdmsrdL3b3Pu7ej/Az/x93Pw94FjgzWq3Y38NHwPtmdmA060RgPq3nc3gPOMrMdot+pxL1bzWfQQOpfuaPAN+Nzoo6Cvg00V1VbMzsFOBfgdPc/fMGix4BzjGzDmbWn3Cw/rVC1BF3b9MTMJJw9sES4JJC16cJ9T2W0AydA8yOppGEPv9ngHeix26FrmsT389w4NHo+X6EP4TFwP1Ah0LXL6bu1cDM6LOYCnRtTZ8DMBFYCMwD7gQ6FPtnANxDOMayjfBf9wWpfuaELpzro7/tuYQzv4r1PSwmHJtI/E3f1GD9S6L38DYwolD11hXcIiISq613Q4mISBMoLEREJJbCQkREYiksREQklsJCRERiKSykVTGz883MmzAV1ZhGIq1dRfwqIkXpOWBamuUL81QPkTZBYSGt1TR3n1DoSoi0FeqGEhGRWAoLKVlmNiE6fnFidPe32uhuY52i5Xua2Q1mtiJattTM/mBmnZOU1dfM7ojucLbRzKZG85aZ2bQG602O9lmdpAw3s9mN5pmZ/cjMXjezzWa2zsweMbPBjdYbHm1/vpl938zmmtmWqO5Xm9luSfZ3uoW7Kq43s4/N7B9m9tVo2W5mtsHM3ovGhmq87W3R/gZk8COXEqawkLbgbmAz8B+E7qtNZtaXMJDkjwj3cbiGMPbOr4DnontTACEoCPcM+Q7hxjQ3Ee6s+CxhaO+WmEK4g1v7qNz7CTdVesnMTkiy/k+i9eYB1wJbgF9Ez3cws4sJ41V9hTAi7r2EUXGfMbOTPQxW9zfCiKbHNdq2Evgm8Iq7v9PC9yelotCDamnSlMkEnE8YSHEaMCHF1C9ad0K07mtAWaNyHiPceGl0o/n/Em1zZYN590TzLmgwrwp4OVGXBvMnR/Oqk9TdgdkNXo+J5t0NVDSY3x9YSxhkrn00b3i0bh1wdIN19yAMyb0Z6BjNO4AwSN0CYM8G634Z2ATMiV4fH5V5Y6N6Jup1UaE/b03FMxW8Apo0ZTI1CIt00/Bo3URYXNyojN5RUDyapPwywvDda6LXVcDnwPwk6x7TwrB4MprXM8m6l0fLRkWvE2HxTJJ1p0bLDopeXxK9/laSdX8I/AxoRxiVdRnwMdCuwToPA1uB7oX+vDUVz6SzoaS1muhNPxtqWaPXQwhflN3NLFkZW4F9zGxvwrDjVYQuq8ZeIfxH31yHE7qRfpzksMFB0WM1oRWUsChJOZ9Gjx2ix0HR48uNV3T3/2z42szuIoTLycBjZtYNOAV43N3XNu1tSFugsJC2oPEXepfo8ahoSqUbsHv0fFPjhe5eb2brWlCvLoS/wUtj6tBQbZJ1EvcZSCRO1+ixKbfbvYMQFucSQukswvGTO5uwrbQhOsAtbVHii/937m5pprlAIgx6Ny4kOouoe6PZjb+4E+t+4WylqB7vx9ThFy14f8nO6qoysx1/9+6+iNBCOt3MOhDCYh3waDP2KyVMYSFt0ZzosSbZQjObaGa/jm61uwjYCNRE9z5v6FB2dv0kbI0eOzWav3+KevQxsz2T1GG0mV1uZoOSbBdnbvR4RJJl1wKfR7foTLgjqu/phDsx3ufuyVow0oYpLKTNcfelwHRghJmd2XCZmX0H+C1wirtvdfc64HbCKaa/arBee+APSYpPDDMyusG6ZYSunsYmE1og10XlJdbvTTid9mKSdH81wX8RDuBfYmY7Wj5mtj+h5fBu9DNIuJcQclcSDnyrC0q+QMcspK0aBzwP3G9mTxCuWziQ8CX/CXBRg3V/A5wAXG5mxxP+cz8R6JOk3P8Cfgf8Mvpyfpdw8Lgr4SyrhiYDpxGuaZhrZk8R/ibPInRv/drdl2T6xtx9YXTg/jLgTTP7b0IonQNUAt9rtP66aJ1vEoLkxUz3KaVPLQtpk9z9bcLZSDcDhwHjCWcR3Qn8k7vPb7DuBsKFa9cRup5+CHxAOGuocbmrCKe5PgOMAP4P4XqHoew8/pFY14Ezo31/DvwAOBuYD5zh7slaLk19f7+Lynof+C7wbcL1JsPc/dUkm9wXPd7V3H1KabPw+yoimTKzLoQAeM7dhxe4Oi1iZr8ndHsNcPfFha6PFB+1LETaODPrA1wAPKugkFQUFiJtlJmdZ2ZvEMbE6km44l0kKYWFSNu1AuhLuAL8h+4+vcD1kSKmYxYiIhJLLQsREYmlsBARkVgKCxERiaWwEBGRWAoLERGJpbAQEZFY/wvucZfb44DVXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29096037e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r,c = seg_data.shape\n",
    "\n",
    "sp = np.fft.rfft(seg_data, axis=1)                    # compute FFT\n",
    "\n",
    "# compute amplitude\n",
    "amp = (np.abs(sp) / c) *2\n",
    "\n",
    "mean_amp = np.mean(amp, axis=0)\n",
    "\n",
    "# get freq axis\n",
    "freq = np.fft.rfftfreq(c, 1/sr)   # Sample spacing is 2nd arg\n",
    "\n",
    "plt.plot(freq, mean_amp, 'b', linewidth=1)          # plot amplitude (power = np.absolute(sp)**2)\n",
    "\n",
    "plt.legend(['normalized fft'], loc='best', fancybox=1)  \n",
    "plt.xlabel('Frequency', **fig_font)\n",
    "plt.ylabel('Amplitude', **fig_font)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notice that there are two prominent features in the freq domain representation. Make a lp filter that gets rid of power > 20Hz (you'll know if it worked if you get rid of the spike at 60Hz)\n",
    "* Because our apply_filter function isn't optimized to work on matrices (we could do that, but haven't) you will need to loop and apply filter to each trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc=0.1\n",
    "b=0.1\n",
    "win = 'blackman'\n",
    "type_filt = 'lp'\n",
    "\n",
    "# make our function call to generate a lp filter\n",
    "n, sinc_func = msf.define_filt(fc=fc, b=b, window=win, type_filt=type_filt)\n",
    "\n",
    "# filter the signal\n",
    "filt_sig = np.zeros((r,c))\n",
    "for i in np.arange(r):\n",
    "    filt_sig[i,:] = msf.apply_filt(input_sig=seg_data[i,:], input_filter=sinc_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then FFT it again to see if you got rid of the hf components of the signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAETCAYAAADH1SqlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYVNWd//H3lwYBgQgCQQIqS3AJKg22SgQjanTEjZi4RhNNzBBNMsHEZ4zGMWrib8YtMWaM+kMj4DIYowaJEaNxIBi3AIqggIoC2gahRfZ9+c4f5zaUbdW9Vd21dn9ez1NPVd26de6pLro/nHPuOdfcHRERkTitSl0BEREpfwoLERFJpLAQEZFECgsREUmksBARkUQKCxERSaSwEBGRRAoLERFJpLAQEZFErUtdgXzp1q2b9+nTp9TVEBGpKLNmzfrI3bsn7ddswqJPnz7MnDmz1NUQEakoZrYkm/3UDSUiIokUFiIikkhhISIiiZrNmIWIlM7WrVupra1l06ZNpa6KZNCuXTt69+5NmzZtGvV+hYWINFltbS2dOnWiT58+mFmpqyMNuDsrVqygtraWvn37NqqMknZDmVk7M/uHmb1mZm+Y2XXR9r5m9rKZvW1mvzez3UpZTxGJt2nTJrp27aqgKFNmRteuXZvU8iv1mMVm4Fh3HwRUAyea2VDgRuBWdx8ArAQuKmEdRSQLCory1tTvp6Rh4cG66Gmb6ObAscAj0fYJwFcKVYdly2Du3EKVLiLSPJS6ZYGZVZnZbGA58AzwDrDK3bdFu9QCvTK8d7SZzTSzmXV1dY06/vTp8POfN+qtIiI7jRgxYufE4JNOOolVq1Y1qbxp06ZxyimnpH3t3HPP5ZBDDuHWW29lwYIFVFdXM3jwYGbNmsUdd9zRpONmUvKwcPft7l4N9AYOBw5Mt1uG94519xp3r+nePXG2elqtWsGOHY16q4g0E9u2bUveKQdPPvkknTt3zmuZ9T788ENeeOEF5syZw49+9CMmTZrEqFGjePXVV+natWvBwqJszoZy91VmNg0YCnQ2s9ZR66I38M9CHVdhIVL5Fi9ezMiRIxk+fDgvvPACvXr14vHHH6d9+/bMnj2biy++mA0bNtC/f3/uvfdeunTpwogRIzjyyCN5/vnnOe2005g7dy7t27dnwYIFLFmyhHHjxjFhwgRefPFFjjjiCMaPHw/AJZdcwowZM9i4cSNnnHEG11133afqU7/80COPPMJdd90FwOrVq+nTpw9Tp07l6aef5pprrmHz5s3079+fcePG0bFjR5566ikuvfRSunXrxpAhQ9J+1hNOOIHly5dTXV3N6aefzp133klVVRXTp0+nR48evPPOO1RXV3P88cdz88035++H7O4luwHdgc7R4/bAc8ApwB+Ac6LtdwHfSyrr0EMP9caYNMn91FMb9VYRicybN6+kx1+0aJFXVVX5q6++6u7uZ555pt9///3u7n7wwQf7tGnT3N396quv9jFjxri7+9FHH+2XXHLJzjIuuOACP/vss33Hjh0+adIk79Spk8+ZM8e3b9/uQ4YM2Vn2ihUr3N1927ZtfvTRR/trr722s7wZM2a4u/u+++7rdXV1O8vesmWLDx8+3CdPnux1dXV+1FFH+bp169zd/YYbbvDrrrvON27c6L179/a33nrLd+zY4WeeeaaffPLJaT/rwIEDdz6/5ppr/Oabb077WkPpvidgpmfx97rU3VA9galmNgeYATzj7k8APwF+bGYLga7A7wpVgaoqtSxE8s0s/7ckffv2pbq6GoBDDz2UxYsXs3r1alatWsXRRx8NwAUXXMD06dN3vufss8/+RBmnnnoqZsbBBx9Mjx49OPjgg2nVqhUDBw5k8eLFADz88MMMGTKEwYMH88YbbzBv3rzEuo0ZM4Zjjz2WU089lZdeeol58+YxbNgwqqurmTBhAkuWLGHBggX07duXAQMGYGacf/75Wf60i6Ok3VDuPgcYnGb7u4Txi4JTN5RI/nnaUcbCatu27c7HVVVVbNy4MfE9HTp0SFtGq1atPlFeq1at2LZtG4sWLeKWW25hxowZdOnShQsvvDBx7sL48eNZsmQJt99+OxB6c44//ngmTpz4if1mz55d1qcfl7plUXIKC5Hma4899qBLly4899xzANx///07WxmNsWbNGjp06MAee+zBsmXLmDJlSuz+s2bN4pZbbuGBBx6gVavw53bo0KE8//zzLFy4EIANGzbw1ltvccABB7Bo0SLeeecdgE+FSTY6derE2rVrc35fNspmgLtUFBYizduECRN2DnD369ePcePGNbqsQYMGMXjwYAYOHEi/fv0YNmxY7P633347H3/8MccccwwANTU13HPPPYwfP55zzz2XzZs3A3D99dez3377MXbsWE4++WS6devG8OHDef3113OqX9euXRk2bBgHHXQQI0eOzOsAt3kp2osFUFNT4425+NFf/wo33BDuRaRx5s+fz4EHpjvrXcpJuu/JzGa5e03Se9UNpZaFiEgihYXCQkQkkcKiFWzfXupaiFS+5tKl3Vw19ftRWKhlIdJk7dq1Y8WKFQqMMuXR9SzatWvX6DJa/NlQmpQn0nS9e/emtraWxi7oKYVXf6W8xmrxYaGWhUjTtWnTptFXYJPKoG4ohYWISCKFhcJCRCSRwkJhISKSSGGhsBARSaSwUFiIiCRSWGhSnohIohYfFppnISKSrMWHhbqhRESSKSwUFiIiiRQWCgsRkUQKC4WFiEgihYXCQkQkkcJCYSEikkhhoXkWIiKJWnxYaJ6FiEiykoaFme1tZlPNbL6ZvWFmY6Lt15rZB2Y2O7qdVKg6qBtKRCRZqS9+tA24zN1fMbNOwCwzeyZ67VZ3v6XQFVBYiIgkK2lYuPtSYGn0eK2ZzQd6FbMOCgsRkWRlM2ZhZn2AwcDL0aYfmNkcM7vXzLpkeM9oM5tpZjMbe+1fhYWISLKyCAsz6wg8Clzq7muAO4H+QDWh5fHLdO9z97HuXuPuNd27d2/UsRUWIiLJSh4WZtaGEBQPuvtjAO6+zN23u/sO4G7g8EIdX2EhIpKs1GdDGfA7YL67/yple8+U3U4HXi9UHTTPQkQkWanPhhoGfAOYa2azo20/Bc41s2rAgcXAdwtVAc2zEBFJVuqzof4OWJqXnixWHdQNJSKSrORjFqWmsBARSdbiw8IM3MNNRETSU1jYrsAQEZH0WnxYgLqiRESSKCxQWIiIJFFYoLkWIiJJFBaoZSEikkRhgSbmiYgkUVigloWISBKFBQoLEZEkCgsUFiIiSRQWKCxERJIoLFBYiIgkUVigsBARSaKwQJPyRESSKCzQPAsRkSQKC9QNJSKSpFFhYWYdzeyLZnZy9LxLfqtVXAoLEZF4OYWFmfUws4eAj4G/A49HL33PzBaa2VH5rmAxKCxEROJlHRZm1h14ETgLeBl4hV3Xz14P7AtMMbOD813JQlNYiIjEy6VlcR2wD3Caux8FPFH/grv/GjgeaA38R15rWAQKCxGReLmExWnAY+7+RLoX3X0a8BjwxTzUq6gUFiIi8XIJi27Auwn71ALdG1+d0tA8CxGReLmERS0wJGGfI6L9KormWYiIxMslLB4BjjOz76Z70cx+DAwH/phtgWa2t5lNNbP5ZvaGmY2Jtu9pZs+Y2dvRfUFPzVU3lIhIvFzC4v8B84A7zGwO4awozGy8mc0FbgYWAv+ZQ5nbgMvc/UBgKPB9M/sCcAXwrLsPAJ6NnheMwkJEJF7WYeHua4FhwF1AH+BAwqmz3wQ+D9wPDHP3VTmUudTdX0kpfz7QCxgFTIh2mwB8JdsyG0NhISISr3UuO7v7GsL//n8I7A90BtYBb7r75qZUxMz6AIMJczh6uPvS6JhLzeyzGd4zGhgNsM8++zT62AoLEZF4OYVFPXffTuiSygsz6wg8Clzq7mvMLOkt9fUYC4wFqKmp8cYeX2EhIhIvY1iY2c8aWaa7+y+y3dnM2hCC4kF3fyzavMzMekatip7A8kbWJSsKCxGReHEti2vTbKv/33u6//p7tN2BrMLCQhPid8B8d/9VykuTgQuAG6L7x9O8PW80z0JEJF5cWJze4Hl74FeEQLgVeIGwoGBH4DDgx8B24Ns5HH8Y8A1grpnNjrb9lBASD5vZRcB7wJk5lJkztSxEROJlDAt3/8T/5s3sFqANcLi7L2qw+z/M7DHC4oLfIKxIm8jd/076VgrAcdmUkQ+alCciEi+XeRbnE9aGahgUQDhribA21Bn5qFgxqWUhIhIvl7BoR/LZU58hc0uhbCksRETi5RIWs4CvmdmB6V40sy8SWhXT81GxYlJYiIjEy2WexbXAX4GXzGw8ITzWAnsQBqrPA7ag61mIiDQ7WYeFuz9nZqOAO4B/Y9dptBC6nuYB33b31/NbxcJTWIiIxMt1uY+nzOzzhKXIBwFdgJXALHefUYD6FYXmWYiIxMt5uQ9330G4FveL+a9OaahlISISL+uwiBYPzIa7+383sj4loXkWIiLxcmlZ/JpdS3o0lLoMiAMVFRZqWYiIxMslLL6VYfvuhOtZfIMwyP2jplaq2BQWIiLxcjkbakLc62b2G+BV4JjovmIoLERE4uUyKS+Wuy8BHgYuzleZxaKwEBGJl7ewiGwC9s5zmQWnsBARiZe3sIiWAfk6sDhfZRaLwkJEJF4up86+kuklwjUt+hLCJ+ur5JULTcoTEYmXy9lQ1TGvbQFeB+529982rUrFp3kWIiLxcjkbKt/jG2VD3VAiIvGyDgAz+6aZHZKwzzAzu7rp1SouhYWISLxcWgvjgVEJ+3wVuLLRtSkRhYWISLyM3VBmdiFwWoPN55rZ4Axv2Y0wIW9FfqpWPAoLEZF4cWMWTwG/IZzpBGHNpwOiWyabgZ/lp2rFo7AQEYmXMSzc/UMz609Y+8mAdwmLCd6WbndgK/CRu28tREULSWEhIhIv9mwod6+rf2xm3wJmR8t6NCuaZyEiEi/rAW53n+Dur+Xz4GZ2r5ktN7PXU7Zda2YfmNns6HZSPo+ZjuZZiIjEyxgWZvaxmf17g+fZ3HIZ4B4PnJhm+63uXh3dnsyhvEZRN5SISLy4bqg1hAHr1OeeYd9GcffpZtYnn2U2hsJCRCRe3AB3n7jnBfYDM/smMBO4zN1XFvJgCgsRkXjluITHnUB/wlpUS4FfZtrRzEab2Uwzm1lXV5dpt0QKCxGReHGT8hpOyMuau09uwnuXpdThbuCJmH3HAmMBampqGt1FprAQEYkXN2YxidzHKCx6T1VjK2RmPd19afT0dMJqtgWlsBARiRcXFj8nzwPaDZnZRGAE0M3MaoFrgBFmVh0dezHw3ULWAUJYbK24qYQiIsUTN8B9baEP7u7nptn8u0IftyG1LERE4uVy8SMAzMyAocAgYA/gI+Bldy94d1GhaFKeiEi8nMLCzL5MOFupX/2m6N7NbCYwOt+zvItBLQsRkXi5XIP7SMKZSa2AB4F/AB8CnYEjga8D/2tmR7r7mwWoa8EoLERE4uXSsrgO2AYc4+4zGrx2j5ndAzxLGBg/O0/1KwqFhYhIvFwm5R0G/CFNUADg7i8AjwDH5qNixaSwEBGJl0tYbAHWJuyzJscyy4LCQkQkXi5/2CcC52Va+M/M9iJcg/uhpleruHQ9CxGReLmMWdwNDANeNbPfANOBD4D2hC6qfyfM3P5bw6VCmrL8RzGoZSEiEi+XsJhDmFVtwNV8enZ3/Wm0Extsa9LyH8WgeRYiIvFyCYuCL/9RKmpZiIjEyzosirH8R6koLERE4lXcmUuFoLAQEYmX63IfpxJmavcB2mbYzd390CbWq6gUFiIi8XJZ7uN7wH+zayA7k4ob11BYiIjEy6UbagxhhdkRQHt3b5XhVtZnPqWjsBARiZdLWPQGHnD36e6+uVAVKgVNyhMRiZdLWMwHPluoipSS5lmIiMTLJSx+DpxlZicWqjKlom4oEZF4ucyzmBwt8/FnM3sTWASk645yd/9avipYDAoLEZF4uZwNdR7wY8LZUAdEt3R0NlQDX/0qPPwwtM75IrYiIuUhlz9f/wFsIATG88D6gtSoBAoZFlu2wB//CP/8J+yzT2GOISJSaLmERR/gHne/u0B1KZlChsW6deH+vfcUFiJSuXIZ4F4E7FaoipRSIcNifdT+ev/9wpQvIlIMuYTFrcA5ZlZdqMqUSiHnWdS3LBQWIlLJcumGWg+8A/zDzJ4DFpJ+3MLd/bJsCjSze4FTgOXuflC0bU/g94Rur8XAWe6+Mod65qyQ8ywUFiLSHOQSFg+kPD4muqXjQFZhAYwHbgfuS9l2BfCsu99gZldEz3+SQz1zVowxC4WFiFSyXMIiUzg0dFC2Bbr79DTX9B5FWH8KYAIwjQoOi/Xr4bOfVViISGXLZVLe3zK9ZmYdgHOA7xCux/3bJtSph7svjY651MwKvsRIoVsWBx4Ib7xRmPJFRIqhSRc/MrMjzOxuYCkwFjgCWJOPimV5/NFmNtPMZtbV1TW6nEKHRb9+sGYNbNxYmGOIiBRazmFhZnua2Rgzmwu8AFwE7A48C5wH9GxinZaZWc/oWD2B5Zl2dPex7l7j7jXdu3dv9AELHRadOkGvXlBbW5hjiIgUWtZhYWbHmdlE4APgV8BAwtIfU4G+7n6Cu0/Mw/Llk4ELoscXAI83sbxEhR6z6NgRevaEDz8szDFERAotNizM7HNmdpWZvQM8DZwNrAbuAIZFu73p7o0avo3C50VgfzOrNbOLgBuA483sbeD46HlBFXqeRceOsPvu6oYSkcqVcYDbzP4E/Eu0z2rC6a0PAc+4+45onyYd3N3PzfDScU0qOEeF7obq2RPatYNNmwpzDBGRQos7G+pkwqS7XwA3Nber46Uq9KS8jh2hfXu1LESkcsV1Q/0FaAdcC3xgZhPN7Ctm1uzWhyrGmIVaFiJSyTKGhbuPJFx3+3LCoPbZwKPAcjO718xOKE4VC6/Q3VD1LQuFhYhUqtgBbndf5u6/dPdBwCDg14RrWlwITCEs7THEzL5Y6IoWUqHDokOH0LJQN5SIVKqsT51197nRAoG9gZMIi/1tAg4H/m5m75jZz81s/8JUtXCK0bJQN5SIVLKcJ+W5+w53f8rdvw70ICzx8Rxhldj/ACpuYYtijFmoG0pEKlmTlvtw93Xufq+7jwD6Aj8D3s5HxYqpGPMs1A0lIpWsSWGRyt3fc/fr3f3AfJVZLMUYs1DLQkQqWd7CopIVap6Fuwa4RaR5UFhQuJbFpk3Qpk24aYBbRCqZwoLwx3xzAeanr18fWhWgbigRqWwKC6BzZ1i1Kv/l1g9ug7qhRKSyKSwI/+t3z/8f89SwUMtCRCqZwgIwgy5dYOXK/Ja7Zk248BFozEJEKpvCIlKIsFi5EvbcMzxWN5SIVDKFRaQQYfHxx7vCQt1QIlLJFBaRQoeFWhYiUskUFhG1LEREMlNYRIrRslBYiEilUlhE1A0lIpKZwiJSrJaFe36PISJSDAqLSKHDoqoKWreGrVvzewwRkWJQWEQKHRYQBrnVFSUilUhhESlGWGiQW0QqVetSVyATM1sMrAW2A9vcvaaQx8t3WOzYAatXh0UK6yksRKRSlW1YRI5x94+KcaB8h8Xq1WERwdYpP2F1Q4lIpVI3VCTfYdGwCwrUshCRylXOYeHA02Y2y8xGF/pgu+8euo7y9cc8dRHBempZiEilKuewGObuQ4CRwPfN7EsNdzCz0WY208xm1tXVNelgZtC1KzSxmJ3UshCR5qRsw8Ld/xndLwf+CByeZp+x7l7j7jXdu3dv8jEPPBDmzWtyMYDCQkSal7IMCzPrYGad6h8DJwCvF/q4gwbB7Nn5KauuLrRUUqkbSkQqVbmeDdUD+KOZQajj/7j7U4U+aHU1TJmSn7IWLYK+fT+5TS0LEalUZRkW7v4uMKjYx62uhv/6r/yUtXAhfKnBKIuWKReRSlWW3VClcsABsGQJrF/f9LIWLoTPf/6T27TyrIhUKoVFit12C4Exd27Tytm+PXRD9ev3ye0dOsCaNU0rW0SkFBQWDRxxBDz/fNPK+OCDMLi9++6f3N6vH7z7btPKFhEpBYVFAyee2PRB7nRdUAD77QdvvdW0skVESkFh0cCxx8LLL8O6dY0vIy4s3n678eWKiJSKwqKBTp3gsMNg6tTGl5EpLD73ubDA4Nq1jS9bRKQUFBZpnHYa3H9/49//yiswcOCnt7dqBQMGqHUhIpVHYZHGv/5rGOSeMSP3965fH7qxRoxI//qAARq3EJHKo7BIo0MHuPZauPxycM/tvdOmwaGHhu6sdDTILSKVSGGRwbe+BR9+mPuZUX/5SzijKpP998/fYoUiIsWisMigdWu48cbQutiyJbv3bNgAkybByJGZ9znuOHj66ezLFBEpBwqLGKeeCv37w5VX7tq2YwcsXpx+/6uugqOOCqvXZtKrVxj8fuaZvFZVRKSgFBYxzGDcOHj0UXj44TB+8cMfhj/2y5d/ct9XX4WJE+G225LLPeusUJ6ISKUwz3UEt0zV1NT4zJkzC1L2a6/Bl78c5k5s2QIHHxyW8/jlL8Pr7qF76ayz4OKLk8tbuhQOOiicYrvvvgWpsohIVsxslrvXJO1XlkuUl5tBg8LA9bJlYYb3ihUhMM47DwYPDmMby5bBd76TXXk9e8Kll8K//Rs8/nhowYiIlDOFRZaGDNn1+HOfg3vugZNOCo8hnDXVOoef5uWXh5nit90WgiPVhg3Qti1UVTW93iIi+aCwaKTTT4cePWDbNhg6NCxvnou2beFPf4IjjwzLll92WZjfMWUKnH9+uEjST34CV1+tloeIlJ7CogmOPLJp7993X/j73+GKK2DvvcNlWJcvDyHSrx8cf3xYS+rGG3NrtYiI5Jv+BJVY377w+99DbW2Y2f2lL+0KhqlTw7jIEUfA174WBtEPOyysMSUiUkwKizLRu3e4perWLXRLTZkCf/1rmFW+fTucc04YQxkwILRA2rUrTZ1FpOXQ/1HLXKtWcPLJcOut8MYbYWB961YYOxa++lXo3Dl0Z/30p2FgXESkENSyqCBmYYb4UUft2rZ9e7h+xs9+Frqo6sc7RETySS2LCldVFRYnfOghuOSScGbWxIlhWRIRkXxRWDQTZvCDH8Cf/ww33RQC5KabwjIkq1eXunYiUunKthvKzE4EbgOqgHvc/YYSV6kiHHZYWEbk5ZfDuMaDD8KiRXDAAfCFL4QQ2W+/8Hi//TTxT0SyU5ZhYWZVwG+B44FaYIaZTXZ3XQkiC2ahO2ro0PB80yaYORPefDPc7rsvDJYvXRoWRayuDkuaDBoEhxwCn/lMaesvIuWnLMMCOBxY6O7vApjZQ8AoQGHRCO3awfDh4ZZq7VqYOxdmzw6LJdaHSNeuYSn1vfYK61jttdenH3fpAu3ba86HSEtRrmHRC3g/5XktcESJ6tJsdeoUZqGnzkTfvh2WLAmtjg8/3HX/0kuf3LZ6dWixtG8PHTuGpUrqb+3bh4mFDW9VVZ98Hhc0cUucNOY1LZkizdltt0GbNoU9RrmGRbpf7U+tpW5mo4HRAPvss0+h69QiVFWFU2+zOf12+/Ywt2P9+l23detg48bw2vbtYe2s1Fv9tq1bM1/fPG7V/Hy+R6S5KMZ/hso1LGqBvVOe9wb+2XAndx8LjIVwPYviVE3qVVWF1kmnTqWuiYgUWrn2OM8ABphZXzPbDTgHmFziOomItFhl2bJw921m9gPgL4RTZ+919zdKXC0RkRarLMMCwN2fBJ4sdT1ERKR8u6FERKSMKCxERCSRwkJERBIpLEREJJHCQkREEpk3k+mtZlYHLGnk27sBH+WxOqWgz1B6lV5/0GcoB8Wu/77u3j1pp2YTFk1hZjPdvabU9WgKfYbSq/T6gz5DOSjX+qsbSkREEiksREQkkcIiGFvqCuSBPkPpVXr9QZ+hHJRl/TVmISIiidSyEBGRRC0+LMzsRDN708wWmtkVpa5PEjPb28ymmtl8M3vDzMZE2/c0s2fM7O3ovkup65rEzKrM7FUzeyJ63tfMXo4+w++j5enLlpl1NrNHzGxB9H18sZK+BzP7UfRv6HUzm2hm7cr9OzCze81suZm9nrIt7c/cgt9Ev9tzzGxI6Wq+S4bPcHP072iOmf3RzDqnvHZl9BneNLN/KU2tW3hYmFkV8FtgJPAF4Fwz+0Jpa5VoG3CZux8IDAW+H9X5CuBZdx8APBs9L3djgPkpz28Ebo0+w0rgopLUKnu3AU+5+wHAIMJnqYjvwcx6AT8Eatz9IMKlAM6h/L+D8cCJDbZl+pmPBAZEt9HAnUWqY5LxfPozPAMc5O6HAG8BVwJEv9vnAAOj99wR/d0quhYdFsDhwEJ3f9fdtwAPAaNKXKdY7r7U3V+JHq8l/IHqRaj3hGi3CcBXSlPD7JhZb+Bk4J7ouQHHAo9Eu5T1ZzCzzwBfAn4H4O5b3H0VlfU9tAbam1lrYHdgKWX+Hbj7dODjBpsz/cxHAfd58BLQ2cx6FqemmaX7DO7+tLtvi56+RLg6KITP8JC7b3b3RcBCwt+tomvpYdELeD/leW20rSKYWR9gMPAy0MPdl0IIFOCzpatZVn4NXA7siJ53BVal/MKU+3fRD6gDxkVdafeYWQcq5Htw9w+AW4D3CCGxGphFZX0H9TL9zCv19/vbwJTocdl8hpYeFukuc14Rp4eZWUfgUeBSd19T6vrkwsxOAZa7+6zUzWl2LefvojUwBLjT3QcD6ynTLqd0on79UUBf4HNAB0K3TUPl/B0kqbR/U5jZVYSu5gfrN6XZrSSfoaWHRS2wd8rz3sA/S1SXrJlZG0JQPOjuj0Wbl9U3saP75aWqXxaGAaeZ2WJC19+xhJZG56hLBMr/u6gFat395ej5I4TwqJTv4cvAInevc/etwGPAkVTWd1Av08+8on6/zewC4BTgPN81p6FsPkNLD4sZwIDoDJDdCANJk0tcp1hR3/7vgPnu/quUlyYDF0SPLwAeL3bdsuXuV7p7b3fvQ/iZ/6+7nwdMBc6Idiv3z/Ah8L6Z7R9tOg6YR+V8D+8BQ81s9+jfVH39K+Y7SJHpZz4Z+GZ0VtRQYHV9d1W5MbMTgZ8Ap7n7hpSXJgPnmFmvMhHwAAAGJUlEQVRbM+tLGKz/RynqiLu36BtwEuHsg3eAq0pdnyzqO5zQDJ0DzI5uJxH6/J8F3o7u9yx1XbP8PCOAJ6LH/Qi/CAuBPwBtS12/hLpXAzOj72IS0KWSvgfgOmAB8DpwP9C23L8DYCJhjGUr4X/dF2X6mRO6cH4b/W7PJZz5Va6fYSFhbKL+d/qulP2vij7Dm8DIUtVbM7hFRCRRS++GEhGRLCgsREQkkcJCREQSKSxERCSRwkJERBIpLKSimNmFZuZZ3MpqTSORStc6eReRsvQ3YFrM6wuKVA+RFkFhIZVqmrtfW+pKiLQU6oYSEZFECgtptszs2mj84rjo6m+bo6uNdYxe38vM7jCz2ui1RWZ2o5l1SlPWPmZ2X3SFs7VmNinattjMpqXsNz46ZnWaMtzMZjfYZmZ2sZm9YmYbzWylmU02s8EN9hsRvf9CM/u2mc01s01R3W8xs93THG+UhasqrjKzj8zsr2b2pei13c1sjZm9F60N1fC990bHG5DDj1yaMYWFtAQPAhuB/yZ0X60zs30IC0leTLiOw62EtXcuB/4WXZsCCEFBuGbINwgXprmLcGXFqYSlvZtiAuEKbrtF5f6BcFGlF8zs2DT7/yDa73XgN8Am4LLo8U5mdiVhvaovEFbEfYiwKu6zZnaCh8XqHiWsaHpUg/e2A74GvOTubzfx80lzUepFtXTTLZcbcCFhIcVpwLUZbn2ifa+N9v0H0KpBOX8mXHjplAbbfxi956aUbROjbRelbGsPvFhfl5Tt46Nt1Wnq7sDslOdnRtseBFqnbO8LrCAsMrdbtG1EtO824Isp++5BWJJ7I9Ah2rYfYZG6+cBeKft+HlgHzImeHxOVeWeDetbX63ul/r51K59bySugm2653FLCIu42Itq3PiyubFBGzygonkhTfivC8t110fP2wAZgXpp9j2xiWDwVbeueZt/ro9dOjp7Xh8WzafadFL12QPT8quj519Ps+13gR0Abwqqsi4GPgDYp+zwObAG6lvr71q18bjobSirVdZ792VCLGzwfQvhD2dXM0pWxBdjbzHoRlh1vT+iyauglwv/oG+tQQjfS99MMGxwQ3VcTWkH13kpTzurovm10Pyi6f7Hhju7+/1Ofm9kDhHA5Afizme0JnAg86e4rsvsY0hIoLKQlaPgHvXN0PzS6ZbIn8Jno8bqGL7r7DjNb2YR6dSb8Dl6TUIdUm9PsU3+dgfrE6RLdZ3O53fsIYXEuIZTOIoyf3J/Fe6UF0QC3tET1f/h/4e4Wc5sL1IdBz4aFRGcRdW2wueEf7vp9P3W2UlSP9xPqcFkTPl+6s7ram9nO33t3f4vQQhplZm0JYbESeKIRx5VmTGEhLdGc6L4m3Ytmdp2ZXRFdavctYC1QE137PNVB7Or6qbcluu/YYHv/DPXobWZ7panDKWZ2vZkNSvO+JHOj+8PTvPYbYEN0ic5690X1HUW4EuPD7p6uBSMtmMJCWhx3XwRMB0aa2Rmpr5nZN4CfASe6+xZ33waMI5xiennKfrsBN6Ypvn6ZkVNS9m1F6OppaDyhBXJ7VF79/j0Jp9NeSZruryz8D2EA/yoz29nyMbP+hJbDu9HPoN5DhJC7iTDwrS4o+RSNWUhLNRp4DviDmU0hzFvYn/BH/mPgeyn7Xg0cC1xvZscQ/ud+HNA7Tbn/A/wC+Pfoj/O7hMHjLoSzrFKNB04jzGmYa2Z/IfxOnkXo3rrC3d/J9YO5+4Jo4P7nwGtm9idCKJ0DtAO+1WD/ldE+XyMEyfO5HlOaP7UspEVy9zcJZyPdDRwCjCGcRXQ/cJi7z0vZdw1h4trthK6n7wIfEM4aaljuMsJprs8CI4F/Jcx3GMau8Y/6fR04Izr2BuA7wNnAPOB0d0/Xcsn28/0iKut94JvA+YT5Jke7+8tp3vJwdP9AY48pzZuFf68ikisz60wIgL+5+4gSV6dJzOw/Cd1eA9x9YanrI+VHLQuRFs7MegMXAVMVFJKJwkKkhTKz88zsVcKaWN0JM95F0lJYiLRctcA+hBng33X36SWuj5QxjVmIiEgitSxERCSRwkJERBIpLEREJJHCQkREEiksREQkkcJCREQS/R+e+mI/N/L3OwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29096037cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r,c = filt_sig.shape\n",
    "\n",
    "# compute FFT - all trials at once\n",
    "sp = np.fft.rfft(filt_sig, axis=1)                    \n",
    "\n",
    "# compute amplitude\n",
    "amp = (np.abs(sp) / c) * 2\n",
    "\n",
    "mean_amp = np.mean(amp, axis=0)\n",
    "\n",
    "# get freq axis\n",
    "freq = np.fft.rfftfreq(c, 1/sr)   # Sample spacing is 2nd arg\n",
    "\n",
    "plt.plot(freq, mean_amp, 'b', linewidth=1)          # plot amplitude (power = np.absolute(sp)**2)\n",
    "\n",
    "plt.legend(['normalized fft'], loc='best', fancybox=1)  \n",
    "plt.xlabel('Frequency', **fig_font)\n",
    "plt.ylabel('Amplitude', **fig_font)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now design a band pass filter that gets rid of everything outside of the 6-14Hz range (where that big bump is...)\n",
    "* To design a band pass filter, make a lp filter (with an appropriate fc), then make a hp filter (with an appropriate fc)\n",
    "* For fun use np.convolve to combine the two filters (instead of multiply in frequ domain) \n",
    "* Then call the apply_filt function\n",
    "* Then you should be left with only those frequencies that pass the lp and the hp filters (i.e. a band pass filter)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.024 0.056]\n",
      "500\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8nPV94PHPd2Z0X9ZlS9ZhyRf4whyOzZlkQ4hNCDjbTVJSktIufbGvlibpkTaw7abbpOwmm21Is03S0hxNUxogNE0cQgADIYANBgM22JYP+ZRs65at+5iZ7/4xz9iK0DHSHM8z0vf9eumlmWd+zzNfaUb6zu8WVcUYY4yZjs/tAIwxxqQHSxjGGGNiYgnDGGNMTCxhGGOMiYklDGOMMTGxhGGMMSYmljCMMcbExBKGMcaYmFjCMMYYE5OA2wEkUllZmdbV1bkdhjHGpJXXX3+9Q1XLpys3pxJGXV0du3fvdjsMY4xJKyJyMpZy1iRljDEmJpYwjDHGxMQShjHGmJhYwjDGGBMTSxjGGGNiYgnDGGNMTCxhGGOMiYklDGOSrKNvmB/tbiIctu2QTXqbUxP3jPGa0VCY//aD13n9ZDct54f41I0r3A7JmFmzGoYxSfR/njzI6ye7WVtVyAPPHGZnY4fbIRkza5YwjEmSp/e38E8vHue3r1nCI3dfw9LyfD798Ju09Qy5HZoxs2IJw5gkONU5wJ/+aC/rqor4i1tWkZcV4Ft3XEn/cIhP/fBNgqGw2yEaM2OWMIxJsJFgmHv+7Q0AvnnHlWQF/ACsWFTA/f95LbuOd/G1Z464GaIxs2IJw5gEe/FIO2+fPs/ffHgtNSW5v/bYb1xZzW3rF/NPLx5jaDTkUoTGzI4lDGMSbEdjJ1kBH5vXVEz4+NbLFzMcDPPGqe4UR2ZMfCxhGJNgO492sKGumOwM/4SPb6wvwe8TdjZ2pjgyY+JjCcOYBGrvHeZgSy/XLiubtExBdgbrq4vYcdSG2Jr0YgnDmATa6SSB65ZPnjCij+9tOkfP0GgqwjImISxhGJNAOxs7KcgOsK6qaMpy1y0vI6yw61hXiiIzJn6WMIxJoB1HO7hmaSl+n0xZ7oraBWRn+NhhM79NGrGEYUyCnOocoLl7cNrmKICsgJ931ZVcaMIyJh1YwjAmQV5qjPZflMZU/rrlZRxu7bOlQkzasIRhTILsONrBosIslpXnx1T+eqcmsvOoDa816cEShjEJEA4rLx/t5LplZYhM3X8RtbqykAW5GdaPYdJGQhKGiGwRkUMi0igi907weJaIPOI8vktE6sY8dp9z/JCIbB5zfIGIPCYiB0WkQUSuSUSsxiTDwZZeuvpHuDaG/oson0+4ZmkpO492omqbKxnvizthiIgf+AZwM7Aa+LiIrB5X7C6gW1WXAw8AX3bOXQ3cDqwBtgDfdK4H8HfAk6p6KbAeaIg3VmOS5eL8i9j6L6KuXV7G6XODnOwcSEZYxiRUImoYG4FGVT2mqiPAw8DWcWW2At93bj8G3CiRevtW4GFVHVbV40AjsFFECoF3A98BUNURVT2XgFiNSYqXGjtYWp5HZVHOjM6L9mO8ZM1SJg0kImFUAU1j7jc7xyYso6pB4DxQOsW5S4F24Hsi8qaIfFtE8hIQqzEJFworrx7v4roplgOZTF1pLouLsnn5mHV8G+9LRMKYqIdvfIPsZGUmOx4ArgS+papXAP3AO/pGAETkbhHZLSK729vbY4/amAQ52dnPwEiIy6qnnt09ERHhsuoFNJztSUJkxiRWIhJGM1Az5n41cGayMiISAIqArinObQaaVXWXc/wxIgnkHVT1QVXdoKobysvL4/xRjJm5w629AFxSUTCr81dWFHCio9/2xzCel4iE8RqwQkTqRSSTSCf2tnFltgF3Orc/AjynkWEh24DbnVFU9cAK4FVVbQGaROQS55wbgQMJiNWYhDvU0gfA8oWxzb8Yb+WifMIKR9v7EhmWMQkXiPcCqhoUkT8EngL8wHdVdb+IfAHYrarbiHRe/0BEGonULG53zt0vIo8SSQZB4B5VjX7M+hTwkJOEjgG/G2+sxiTD4bZeaktyyc2c3Z/TJYsiNZMjrX2sWTzzZi1jUiXuhAGgqk8AT4w79vkxt4eAj05y7v3A/RMc3wNsSER8xiTT4ZZeVi6aXXMUQF1ZHhl+4ZDTtGWMV9lMb2PiMBIMc7yjn5WLZtccBZDh97G0LJ8jljCMx1nCMCYOxzv6CYZ11h3eUSsrCqyGYTzPEoYxcYj+k1+xMM6EsTCfpq5B+oeDiQjLmKSwhGFMHI609uL3CUvL45tXutKpoTS22Ugp412WMIyJw6GWXupKc8nO8E9feArRkVLWLGW8zBKGMXE40tYX1wipqJqSXLICPuv4Np5mCcOYWRoaDXGisz8hCcPvE1YsyudQqzVJGe+yhGHMLDW29aFKQhIGwMqFBRxusRqG8S5LGMbM0sU1pGY/B2OslRUFtPQMcX5wNCHXMybRLGEYM0uHWnvJ9PtYUpqYlfcvLhFitQzjTZYwjJmlI619LC3PI8OfmD+jFc5s8cPWj2E8yhKGMbN0KM41pMarWpBDXqb/QlOXMV5jCcOYWegbDnL63GBca0iNJyKsWFTAIev4Nh5lCcOYWYj2MySyhgGRfowjbZYwjDdZwjBmFuLdZW8yKysK6OgbobNvOKHXNSYRLGEYMwuHW/vIzvBRU5yb0OuutI5v42GWMIyZhcOtvaxYWIDPJwm97oU1pVp6EnpdYxLBEoYxs3CsvZ9lca5QO5HygiwKsgMc7+hP+LWNiZclDGNmaDgY4sz5wYRN2BtLRFhSmsuJzoGEX9uYeFnCMGaGmroGUYUlpYntv4haUpLHqS5LGMZ7LGEYM0OnuiLNRclKGLWluTR3DxAKa1Kub8xsWcIwZoZOOs1FtSWJb5ICWFKSy2hIOXNuMCnXN2a2LGEYM0MnOwfIzfRTlp+ZlOvXOjUXa5YyXmMJw5gZOtU1QG1JLiKJHVIbFe1MP2kd38ZjLGEYM0MnO/upS8IIqaiKwmwy/T5OdtnQWuMtljCMmYFQWGnqGkxahzdEtmutLsnhZIfVMIy3WMIwZgZaeoYYCYUv9DMky5KSXE5aH4bxmIQkDBHZIiKHRKRRRO6d4PEsEXnEeXyXiNSNeew+5/ghEdk87jy/iLwpIo8nIk5j4nWy0xlSm6QRUlFLSvM41dmPqg2tNd4Rd8IQET/wDeBmYDXwcRFZPa7YXUC3qi4HHgC+7Jy7GrgdWANsAb7pXC/qM0BDvDEakyinnI7oZDZJRa/fPxKis38kqc9jzEwkooaxEWhU1WOqOgI8DGwdV2Yr8H3n9mPAjRIZYrIVeFhVh1X1ONDoXA8RqQZuAb6dgBiNSYiTXQMEfEJlUXZSnyeakGyklPGSRCSMKqBpzP1m59iEZVQ1CJwHSqc592vAnwPhBMRoTEKc6hygujiHQIL28Z5MdFLgKRspZTwkEe/6iQajj294nazMhMdF5ENAm6q+Pu2Ti9wtIrtFZHd7e/v00RoTh5Nd/UlZdHC8mpIcRKyGYbwlEQmjGagZc78aODNZGREJAEVA1xTnXgfcJiIniDRxvU9E/nWiJ1fVB1V1g6puKC8vj/+nMWYSqsrJzoGk918AZAX8VBZmX+gzMcYLEpEwXgNWiEi9iGQS6cTeNq7MNuBO5/ZHgOc0MvxjG3C7M4qqHlgBvKqq96lqtarWOdd7TlU/kYBYjZm17oFReoeC1JYkP2FAZImQE53WJGW8IxDvBVQ1KCJ/CDwF+IHvqup+EfkCsFtVtwHfAX4gIo1Eaha3O+fuF5FHgQNAELhHVUPxxmRMMlwYUpuCJimIDN199mBrSp7LmFjEnTAAVPUJ4Ilxxz4/5vYQ8NFJzr0fuH+Kaz8PPJ+IOI2JR3QxwFQ0SUGkhtHRN0LfcJD8rIT8qRoTF5vpbUyMLi5rnpqEEU1M1o9hvMIShjExOtk5wKLCLLIz/NMXToDoAoc2tNZ4hSUMY2J0qqs/6UuCjFVrk/eMx1jCMCZGqRpSG1WYnUFxboYtQmg8wxKGMTEYGAnS1juc0oQBUFuad2F0ljFus4RhTAyiI6RqUzSkNmpJSa41SRnPsIRhTAyi/7SXpGiEVNSS0lzOnBtkJGhLqhn3WcIwJgapWtZ8vNqSXMIKp88NpvR5jZmIJQxjYnCyq5/C7AALcjNT+rzRWeXWj2G8wBKGMTFo6hpM+rasE4lOEmzqthqGcZ8lDGNi0Nw9QPWC1CeMhQVZZPp9NHdbx7dxnyUMY6ahqjR3D1JTkpPy5/b5hKriHJqthmE8wBKGMdPo6BthOBimujj1NQyAaksYxiMsYRgzjSanOai6OPU1jOjzNttsb+MBljCMmUb00717NYxcOvtHGBgJuvL8xkRZwjBmGs0eqGEAnLZmKeMySxjGTKO5e5CSvEzyXNrEKFqzsX4M4zZLGMZMo6lrwLXaBUCN89xNNrTWuMwShjHTON096GrCKMvPIjPgsxqGcZ0lDGOmEA4rzecGqXGpwxsiczEiQ2uthmHcZQnDmCl09A0zEgy7WsOASD+G1TCM2yxhGDOFi3Mw3KthRJ4/hyabi2FcZgnDmClcnIPhdg0jh+6BUfqGbS6GcY8lDGOm4PakvahoH4rNxTBusoRhzBSauwcoy88kJ9PvahzRGo51fBs3WcIwZgpNXYNUuVy7gIs1HOvHMG6yhGHMFJq73Z20F1WWn0l2hs3FMO6yhGHMJMJh5bTLczCiRMSG1hrXJSRhiMgWETkkIo0icu8Ej2eJyCPO47tEpG7MY/c5xw+JyGbnWI2I/FJEGkRkv4h8JhFxGjMTbb3DjIbUEzUMcJY5P2dNUsY9cScMEfED3wBuBlYDHxeR1eOK3QV0q+py4AHgy865q4HbgTXAFuCbzvWCwJ+q6irgauCeCa5pTFK5vQ/GeJG5GFbDMO5JRA1jI9CoqsdUdQR4GNg6rsxW4PvO7ceAG0VEnOMPq+qwqh4HGoGNqnpWVd8AUNVeoAGoSkCsxsSs2SOT9qKqi3M5PzhKz9Co26GYeSoRCaMKaBpzv5l3/nO/UEZVg8B5oDSWc53mqyuAXQmI1ZiYNXd5Y9JelM3FMG5LRMKQCY5pjGWmPFdE8oF/B/5IVXsmfHKRu0Vkt4jsbm9vjzFkY6bX1D1AeUEW2RnuzsGIiiYuG1pr3JKIhNEM1Iy5Xw2cmayMiASAIqBrqnNFJINIsnhIVX882ZOr6oOqukFVN5SXl8f5oxhzUbPLy5qPd3HyntUwjDsSkTBeA1aISL2IZBLpxN42rsw24E7n9keA51RVneO3O6Oo6oEVwKtO/8Z3gAZV/WoCYjRmxpq7vTGkNqokL5PcTL8lDOOauBOG0yfxh8BTRDqnH1XV/SLyBRG5zSn2HaBURBqBPwHudc7dDzwKHACeBO5R1RBwHfBJ4H0issf5+mC8sRoTq1BYOXPOWzWMyFwM2xfDuCchmxSr6hPAE+OOfX7M7SHgo5Ocez9w/7hjLzFx/4YxKdHSM0QwrJ4ZIRVVXZxLk9UwjEtsprcxE2ju8tYcjCirYRg3WcIwZgLRfoKaEm/VMGqKc+kdCnJ+wOZimNSzhGHMBJq6BxCBxQuy3Q7l11wYWmu1DOMCSxjGTKCpa5BFBdlkBbwxByMqWuOxZinjBksYxkygqXuAmhJv9V/AxdnetqaUcYMlDGMm0Nw14LkRUgCFOQEKsgLWJGVcYQnDmHFGgmFaeoao8dgIKXDmYpTk2vIgxhWWMIwZ5+z5QcIK1R4bIRVVU5xjs72NKyxhGDNOtH/AS8uCjFVTEtl5L7K6jjGpYwnDmHGi/QNe7PSGSA1jcDRER9+I26GYecYShjHjNHUNEPAJlUUeTRhOU5l1fJtUs4RhzDhN3YMsXpCD3+fN5cwuJAzr+DYpZgnDmHGaurw5ByPK9sUwbrGEYcw4zd2DVC/wZoc3QG5mgNK8TKthmJSzhGHMGIMjITr6hj1dw4DIkF/rwzCpZgnDmDGaL4yQ8m4NAyIjpWx5EJNqljCMGSP6qd2Ly4KMVVOSy5lzg4TCNhfDpI4lDGPGuDBpz+NNUjXFuQTDSkvPkNuhmHnEEoYxYzR1DZCd4aM8P8vtUKYUTWjW8W1SyRKGMWM0dUdWqRXx5hyMqIvLnFvCMKljCcOYMZq7Bz25Su14ixfkIBKZZGhMqljCMGaMyKQ9b3d4A2QGfFQWZtNsNQyTQpYwjHGcHxylZyh4YSa111UX21wMk1qWMIxxRPsDvLqs+XjVJTYXw6SWJQxjHOkyaS+qpjiX1t4hhoMht0Mx80TA7QCMd6gqJzsH2Nt8jv1neujuH6F/JEj/cIhQWLlySTHvWVnO+uoiAv6591nD6xsnjVdTkosqnO4eZGl5vtvhJFzfcJCdjR28cKSdo2395Gb6yc0KkJ/lp6Iwh8tqilhfvYCSvEy3Q503LGHMY+GwcuBsDzuPdrDzaCd7ms5xbmAUiHSqluZlkpvpJz8rQDCs/L/njvD1Z49QmB3gptUV/OUtqyieQ3+sTd0DFGQHKMrNcDuUmNSMWbV2LiWMvU3n+N+/aGD3iW6CYSU308+lFQX0DI0yMBKibzhIR98w0Q0Ha0pyeNeSEq5ZVsp1y8tYvCA9+qDSUUIShohsAf4O8APfVtUvjXs8C/gX4CqgE/hNVT3hPHYfcBcQAj6tqk/Fck0zc+Gw0tDSw65jXew63smu410XEsSy8jy2rKlgfc0C1lcvYOWi/HfUIs4NjPBSYwe/OtTOT/ec4ZVjnXzzjitZX7PAjR8n4SJDatOjdgFzbyMlVeVfXznJFx4/QGleFr93w1LevbKMDUtKyAz8+nuxbzjI283neav5HHuazvH84XZ+/OZpAOpKc7l6aSmblpawqb7UEkgCxZ0wRMQPfAO4CWgGXhORbap6YEyxu4BuVV0uIrcDXwZ+U0RWA7cDa4DFwDMistI5Z7prmimoKmfOD/F283n2nT7PW6fPs+dUNz1DQSCyp8L7Vy3iuuWlXLO0jIqi7GmvuSA3kw9dtpgPXbaYT1y9hD946A0++g8v8z9uXc0nNtV6frLbdJq6Blhanud2GDFbVJhNhl/mRMf3wEiQ//7jt/nJnjO895JyHvjY5VPWXvOzAlyzrJRrlpUCkQ9Dh1p72Xm0k52NHfz87bM8/FoTEHmvr69ZwGVVRayrLmLN4iKKctKjFuk1iahhbAQaVfUYgIg8DGwFxv5z3wr8T+f2Y8DfS+S/y1bgYVUdBo6LSKNzPWK45rylqgyMhDg/OEpH3zCtPcO09gzR2jPE8Y7+C18DI5HOUL9PWLmogA+uq2RjfQmblpZSFeenrvU1C3j8U9fzx4/u4X/8ZB8NZ3u4/8Nr0zZpqCrN3YO8Z2W526HEzO8TqhbkpH0No7t/hNsffIUjbb189gMr+YP3Lsc3w90OfT5hVWUhqyoLuev6ekJhpeFsD7uOd7H7RBd7Tp3j52+dvVC+LD+LpWV51JflUVuay8KCLBYVZrOoMJuSvEwKcwJkBfyJ/lHTXiISRhXQNOZ+M7BpsjKqGhSR80Cpc/yVcedWObenu2bC/Mebzexs7EzY9cauH6o6/jG9UEiJ/KMKK4RVCasSCivBkDIaVoKhMKOhMEOjYYZGQwwFQ/QPh+gZHCU4wSqlPomMza8vy2NjfQlLy/NZuzjyR5Sdkfg3f3FeJt+98118+cmD/OMLx1hVWcgnr16S8OdJhY6+EQZHQ2kzByOqujg3rSfvhcPKHz+6h+Md/fzz727k3QlK2H6fsLaqiLVVRdx1fT0AXf0jvH36PA1nezjeHvlQ9ezBVjr6Ria8RlbAR2FOBrmZfnIy/GRl+MkO+Mjw+wj4hYDPR4Zf8Ing8wk+AQFEBAEQcG5dMP7zVCI/Xt11Qz2XVhQm8IrvlIiEMdHPPP6/2WRlJjs+0RCcCddxFpG7gbsBamtrJ49yCkfb+tnR2DGrcycz1Sft6EPivKF8gvOGi9yOvhED/sj3svxMcjL9ZAf85Gb5KcrJoDA7g8KcDMrys1hUGPl0VJqXmfLRSz6f8Lktl3KotZcv/uwA66uLuKw6/fo0mtJsSG1UTUkuT+47O31Bj/rWr47y/KF2vrh1TcKSxWRK8jJ5z8ryd9Qih0ZDtPUM09o7RMv5Ic4NjNAzFKRncPRCR/vQaIhB54PbwEiQYFgZDUU+1IVV0Qsf+iIfClXf+WFxPJ2uwAz9xpXVCb3eRBKRMJqBmjH3q4Ezk5RpFpEAUAR0TXPudNcEQFUfBB4E2LBhw6xegc9uvoTPbr5kNqcaIknjgY9dzi1ff5E/eOgNfv6pG9JmpFHUiY5+AJaUpk8fBsCS0ly6B0Y5Pziadu3yO4928LdPH+LW9ZE+MbdkZ/ipLc2ltjS9Piy4IREfR18DVohIvYhkEunE3jauzDbgTuf2R4DnNJJetwG3i0iWiNQDK4BXY7ym8ZDivEz+/o4rae0Z4k9/tDfhn56S7URHPz6B2jSrYdQ5CS6a8NJFW88Qn/7hHurK8vjfv7Eubfu+5pu4E4aqBoE/BJ4CGoBHVXW/iHxBRG5zin0HKHU6tf8EuNc5dz/wKJHO7CeBe1Q1NNk1443VJNeVtcXcd/Mqnmlo5Xs7Trgdzowc74wsaz5++KbXRUd1nehMn4ShGum36Bse5Vt3XEV+lk0HSxcJeaVU9QngiXHHPj/m9hDw0UnOvR+4P5ZrGu/73evq+NXhdr72zGH+y5XVadM0daKjn7qy9GqOgkiNSASOp1EN49mGNnY0dvLFrWu4pKLA7XDMDKTXxynjeSLCvTdfSu9wkH944ajb4cREVTnR0U99GrZhZ2f4WVyUkzZNUqGw8pWnDlFflsfHN85ukIpxjyUMk3CrKgu5bf1ivrfjOG1psOd0Z/8IvcPBtKxhANSV5XK8Mz2G1m7be5pDrb38yU0r5+R6ZHOdvWImKf7kppUEQ8rXnzvidijTin46T9uEUZqXFjWMkWCYr24/zJrFhdyyrtLtcMwsWMIwSbGkNI/bN9bw8KtNnPR4h+wx559tfZoNqY2qL8vj/OAo3f0TT0DziodfO0VT1yB/tvmSGc/kNt5gCcMkzafft4KAX/jq9sNuhzKlEx39BHySdrO8o6JDa497ODEPjAT5+rONbKwvSavlV8yvs4RhkmZhYTa/c2092/aeoeFsj9vhTOpEZz81Jblp26YebUrzcrPU93acoKNvmM9tucTmXKSx9PwLMWnj99+zjNwMP9956bjboUzqeMcAdWk4QiqqtiQXn3g3YYwEw3xvxwnes7Kcq5aUuB2OiYMlDJNURbkZbL2iip/tPcN5Z+8NL4nsMpieczCiMgM+qopzPDtSavuBVjr6hvmda+vcDsXEyRKGSbo7NtUyHAzz2BvNbofyDm29wwyMhKhP44QB3h4p9a+vnKRqQU7SFxc0yWcJwyTdmsVFXFG7gId2nfTcGlPRGdJ1aTpCKqq+LI/jHf2e+/02tvXx8rFOfmtTLX4bGZX2LGGYlPjEpiUca+/n5WOJ23ckEaKfyudCDSOy17W3htb+265TZPiFj22omb6w8TxLGCYlbrmskqKcDB565ZTbofya4539ZPp9ab/vczTheWkRwsGREI+93sTmNRWUF2S5HY5JAEsYJiWyM/x89Kpqntrf4qnlQk509FNbmpv2zSXRhOGlRQh/9tYZeoaCru51YRLLEoZJmd/aVEswrDy6u2n6wilyomMg7fsvAKqLcwj4xFMd3w/tOsXyhflsqrehtHOFJQyTMkvL87lueSk/fLWJ0AR7kqdaOKyc6Oynvix952BEBfw+akpyPdMkte/0efY2neOOTbU2UW8OsYRhUuqOTUs4fW6QF4+0ux0KZ3uGGA6G03oOxlh1pbkc7/DGXIxHXmsiK+BLyT7TJnUsYZiUunHVQgqyA/xs71m3Q7k4QmoONElBZImQk53uD60NhsL8/O2zvH/1orTbZ9xMzRKGSamsgJ8tayp4en8LQ6MhV2M5nubLmo9XX5bHwEiItt5hV+PYcbSTrv4Rblu/2NU4TOJZwjApd+v6xfQOB3n+UJurcZzo6Ccr4KOiMNvVOBLlwqq1Lnd8b9tzhoKsgK1KOwdZwjApd+2yUsryM11vljrR2U9dad6c2Zuh3gOr1g6Nhnh6fwub11aQneF3LQ6THJYwTMoF/D4+uK6SZxpa6RsOuhbH8Y5+6ubACKmoxQtyyPT7XN0X4/lD7fQOB605ao6yhGFccev6xQwHw2w/0OLK84fCSlPX4JzpvwDw+4SakhyOt7uXMH629wyleZlcu6zUtRhM8ljCMK64qraYxUXZrjVLHe/oZyQUZsXCAleeP1lWLCzgSFufK8/dNxzk2YOtfHBdZdpuRmWmZq+qcYXPJ9y6fjEvHG53ZS/q6A6AqyrnVsJYVVnIic5+BkZS39T3zIFWhkbD3Ha5NUfNVZYwjGtuXb+YYFh5cn/qm6UazvYQ8AnLF+an/LmT6dLKAlThYEtvyp/7Z3vPUFmUzVW1xSl/bpMaljCMa9YsLmRpWR7b9pxJ+XMfbOllWXk+WYG5NZJndWUhAAfPpjZhnBsY4YUj7dy6fvGcGXVm3skShnGNSKRZ6pXjnbT1pnYF24azPXOuOQoiixAWZAUuNLmlylP7WxgNKbdeZs1Rc5klDOOqD66rRBWe3t+asuc8NzDC2fNDrHI+jc8lIsKllQUpTxi/2NdCTUkOa6vm3u/UXBRXwhCREhHZLiJHnO8TNl6KyJ1OmSMicueY41eJyNsi0igiXxdnWUsR+YqIHBSRt0TkP0RkQTxxGu9auSif+rI8nkphP8aBCx3ec/Of26rKQg629BJO0YrAPUOj7GjsYMuaCluZdo6Lt4ZxL/Csqq4AnnXu/xoRKQH+CtgEbAT+akxi+RZwN7DC+driHN8OrFXVy4DDwH1xxmk8SkTYsraCl492cm4gNaOlGpz2/UvnYJMUwKUVhfQNB2nuHkzJ8z3X0MZoSNmytjIlz2fcE2/C2Ap837n9feDDE5TZDGxX1S5V7SaSDLaISCVQqKova2R5zX+Jnq+qT6tqdFzgK4CtkTyH3by2gmBY2X4gNc1SB8/2UJZfdNskAAAPIUlEQVSfycKCubGG1HjRvpmGltQ0S/1i31kWFWZxRY01BMx18SaMRap6FsD5vnCCMlXA2C3Wmp1jVc7t8cfH+6/AL+KM03jYuqoiqhbkpKxZqqGlZ842RwFcUlGACCnpxxgYCfKrw+1sXlNho6PmgWkThog8IyL7JvjaGuNzTPQu0imOj33uvwCCwENTxHe3iOwWkd3t7e5vymNmTkTYvKaCF450JH1tqWAozOHWvjmdMHIzA9SX5qUkYfzqUDtDo2G2rK1I+nMZ902bMFT1/aq6doKvnwKtTtMSzveJ1qtuBmrG3K8GzjjHqyc4jnO9O4EPAXfoFDvCqOqDqrpBVTeUl9tyyunq5nUVjATD/PJgcpc8P9bRz0gwPCeH1I61qrLwQl9NMj25v4WSvEw21tm+3fNBvE1S24DoqKc7gZ9OUOYp4AMiUux0dn8AeMppwuoVkaud0VG/HT1fRLYAnwNuU1Vv7DlpkurK2mLK8rN4cl9ym6Win7ovrZi7NQyASysKONU1QO/QaNKeYzgY4rmGNm5atcjWjpon4n2VvwTcJCJHgJuc+4jIBhH5NoCqdgFfBF5zvr7gHAP4feDbQCNwlIt9FX8PFADbRWSPiPxDnHEaj/P7hM1rFvHLQ21J3Ymv4WwvGX5hWfncWhJkvGiT2+HW5NUydjR20DscZMs6a46aLwLxnKyqncCNExzfDfzemPvfBb47Sbm1ExxfHk9cJj1tWVvBQ7tO8cLhdj6wJjn/hBrO9rB8YQGZgbn9iXjV4kjCOHC2l6uWJKe56Ml9LRRkBWwp83lkbv/VmLRy9dJSinIyktosNVeXBBlvcVE2hdnJWyIkGAqz/UArN65aOOfW4zKTs4RhPCPD7+Om1YvY3tDKcDDxzVKdfcO09Q6zao73X0B0iZDCpCWMV4510T0wapP15hlLGMZTbrmskt6hIC8e7kj4taOjhubykNqxVlcWcihJS4Q8/tYZ8jL9vPcSG5k4n1jCMJ5y/fIyinIy+Pnbid+J72DL3Nw0aTKrKgsYGAlxqiuxAw1HQ2Ge3N/CTasXkZ1hzVHziSUM4ykZfh9b1lSw/UBrwkdLHTjbw8KCLErzsxJ6Xa+K1qQS3Sy182gn5wZGucWWMp93LGEYz7nlskr6hiNLTiTSgTNze0mQ8VYuKsAnsP9MYhPG43vPUJAV4N0ryxJ6XeN9ljCM51y7rJTi3Ax+/lbimqU6+oY52NLLu+rmz/ah2Rl+1lUvYOfRxPUHjQTDPLW/hZvWLLLRUfOQJQzjOQG/jy1rK3mmoZXBkcQ0S+1ojPzTvH7F/OqkvWF5GXubz3N+MDEzvl9qbKdnKMiHLrPRUfORJQzjSR+6rJKBkRDPH0rM2lIvHemgKCeDdVVFCbleurhhRRmhsPLy0c6EXO/xt85SmB3g+uXzK/GaCEsYxpM21ZdQlp/J4wkYLaWqvHikg+uWl+KfZ0twX1FbTG6mn5ca4+8PGg6G2L6/lc1rKub8THkzMXvVjSdFmqUqeK6hjYGR+JY8P9reR0vP0Lz8VJwZ8HH10lJeOhJ/P8YLhyNrR91izVHzliUM41m3rFvM4GiI5+Jc8vxF55/lDSvm56ieG1aUcaJzgKY452P8/K0zLMjN4Lrl8/P3aCxhGA/bWF9CeUEW//HG6biu89KRDupKc6kpyU1QZOklmihfjKOW0Ts0yvYDrWxZU0GGLWU+b9krbzzL7xM+tqGa5w61zfrT8UgwzMvHOrl+ntYuAJaV51NRmB1XP8ZP3jxN/0iI33xXzfSFzZxlCcN42m9tWoIA//bqqVmd/+apbgZGQtwwz4bTjiUi3LCijB2NnYRmsa6UqvKDV06yrqqIy2sWJCFCky4sYRhPq1qQw/tXLeKR15pmtYLtS40d+H3CNfN8z4brV5RxfnCUt0+fn/G5u453cbi1j09es4TI5phmvrKEYTzvk9csoat/hCdmMcT2hSMdrK8uojA7IwmRpY9oR/VLR2beLPWDV05SlJPBrbZ21LxnCcN43nXLylhalse/vHxyRuedHxjl7eZz87o5KqosP4s1iwt5YYYd3209Qzy1r4WPbagmJ9OWApnvLGEYz/P5hDuuXsKbp86xbwZNKjuPdhDW+TucdrzrV5Tx5qlu+odjn9fyw1ebCIaVOzYtSWJkJl1YwjBp4SNXVZOd4eMHM6hlbG9oJT8rwHrrqAXghuXljIaU5w/F1iw1Ggrzb6+e5D0ry6kry0tydCYdWMIwaaEoJ4MPX17FT/ee5vzA9AvpHWnt5SdvnuYjV1XbvAHHpqUl1Jfl8cAzhwmGwtOWf+ZAK609w3zyaqtdmAj7SzJp45PXLGFoNMyDLx6dtuyXfnGQvMwAn75xRQoiSw8Zfh+f23IJjW19PLq7ecqyo6Ew33i+kaoFOfynSxemKELjdZYwTNpYs7iI/3JlNd96/ii7jk2++urOox08e7CNe963nJK8zBRG6H2b11SwYUkxX91+eMq+jL99+jD7Tvfwl7esmncLNprJWcIwaeWvt66htiSXP3pkD+cGRt7xeDis/K8nGqhakMPvXFuX+gA9TkT477esoqNvmAdfODZhmZeOdPCPLxzl4xtruHmdLTRoLrKEYdJKflaA//fxK+noG+bef38b1V+fubxt7xn2ne7hs5tXkp1hw0AncmVtMbesq+TBF47R1jP0a4919g3zx4/uYVl5Pp//0BqXIjReZQnDpJ111UX82eZLeHJ/Cz98tenC8f7hIF956hBrqwrZur7KxQi978+3XEIwHOZvnz58YbkQVeXPHnuL8wOjfP32K2zehXmHgNsBGDMbv3f9Ul480sFfbdvHA88cpmdwlOFgZOTPVz5yGT5rd5/SktI8PnH1Er634wSP7G4iPytAbqaftt5h/uetq1m9uNDtEI0HxZUwRKQEeASoA04AH1PV7gnK3Qn8pXP3b1T1+87xq4B/BnKAJ4DP6Jg2BhH5LPAVoFxVE7eTvUl7Pp/wwG9ezteeiXxCLszOoDAng0srCrjW9muIyee2XMryhfm09w7TMxikZ2iUutJc7rS+HzMJGd8GPKOTRf4P0KWqXxKRe4FiVf3cuDIlwG5gA6DA68BVqtotIq8CnwFeIZIwvq6qv3DOqwG+DVzqlJ82YWzYsEF3794965/HGGPmIxF5XVU3TFcu3j6MrcD3ndvfBz48QZnNwHZV7XJqH9uBLSJSCRSq6stOreJfxp3/APDnRJKMMcYYl8WbMBap6lkA5/tEM3yqgKYx95udY1XO7fHHEZHbgNOqujfO+IwxxiTItH0YIvIMUDHBQ38R43NM1Puokx0XkVzn2h+I6eIidwN3A9TW1sYYkjHGmJmaNmGo6vsne0xEWkWkUlXPOk1MbRMUawbeO+Z+NfC8c7x63PEzwDKgHtjrbNZSDbwhIhtVtWWC+B4EHoRIH8Z0P48xxpjZibdJahtwp3P7TuCnE5R5CviAiBSLSDGRmsNTThNWr4hcLZHM8NvAT1X1bVVdqKp1qlpHJLFcOVGyMMYYkzrxJowvATeJyBHgJuc+IrJBRL4NoKpdwBeB15yvLzjHAH6fyEioRuAo8Is44zHGGJMkcQ2r9RobVmuMMTOXqmG1xhhj5ok5VcMQkXZgZhs/X1QGeHE2ucU1MxbXzFhcM+PVuCC+2Jaoavl0heZUwoiHiOyOpUqWahbXzFhcM2NxzYxX44LUxGZNUsYYY2JiCcMYY0xMLGFc9KDbAUzC4poZi2tmLK6Z8WpckILYrA/DGGNMTKyGYYwxJiaWMAAR2SIih0Sk0dnXw604visibSKyb8yxEhHZLiJHnO/FLsRVIyK/FJEGEdkvIp/xQmwiki0ir4rIXieuv3aO14vILieuR0QkM5VxjYnPLyJvisjjXolLRE6IyNsiskdEdjvHvPAeWyAij4nIQed9do3bcYnIJc7vKfrVIyJ/5HZcTmx/7Lzn94nID52/haS/v+Z9whARP/AN4GZgNfBxEVntUjj/DGwZd+xe4FlVXQE869xPtSDwp6q6CrgauMf5Hbkd2zDwPlVdD1xOZJ+Vq4EvAw84cXUDd6U4rqjPAA1j7nslrv+kqpePGYLp9usI8HfAk6p6KbCeyO/N1bhU9ZDze7ocuAoYAP7D7bhEpAr4NLBBVdcCfuB2UvH+UtV5/QVcQ2QxxOj9+4D7XIynDtg35v4hoNK5XQkc8sDv7KdE1g7zTGxALvAGsInI5KXARK9vCuOpJvLP5H3A40SW8/dCXCeAsnHHXH0dgULgOE6fqlfiGhfLB4AdXoiLi3sMlRBZcfxxIhvVJf39Ne9rGEy+wZNXxLJJVcqISB1wBbALD8TmNPvsIbK0/nYii1ieU9WgU8St1/NrRHaMDDv3Sz0SlwJPi8jrzl4y4P7ruBRoB77nNOF9W0TyPBDXWLcDP3RuuxqXqp4G/i9wCjgLnCey9XXS31+WMCbf4MmMIyL5wL8Df6SqPW7HA6CqIY00GVQDG4FVExVLZUwi8iGgTVVfH3t4gqJuvM+uU9UriTTB3iMi73YhhvECwJXAt1T1CqAfd5rFJuT0BdwG/MjtWACcPpOtRPYNWgzkEXk9x0v4+8sSRiQT14y5H93IyStanc2pmGKTqqQTkQwiyeIhVf2xl2IDUNVzRDbmuhpYICLRzcHceD2vA24TkRPAw0Sapb7mgbhQ1TPO9zYi7fEbcf91bAaaVXWXc/8xIgnE7biibgbeUNVW577bcb0fOK6q7ao6CvwYuJYUvL8sYUT26FjhjDDIJFL13OZyTGPFsklVUomIAN8BGlT1q16JTUTKRWSBczuHyB9SA/BL4CNuxaWq96lqtUY2ALsdeE5V73A7LhHJE5GC6G0i7fL7cPl11MjmaE0icolz6EbggNtxjfFxLjZHgftxnQKuFpFc528z+vtK/vvLrU4kL30BHwQOE2n//gsX4/ghkTbJUSKfuu4i0vb9LHDE+V7iQlzXE6nevgXscb4+6HZswGXAm05c+4DPO8eXAq8S2ZjrR0CWi6/pe4HHvRCX8/x7na/90fe626+jE8PlwG7ntfwJUOyRuHKBTqBozDEvxPXXwEHnff8DICsV7y+b6W2MMSYm1iRljDEmJpYwjDHGxMQShjHGmJhYwjDGGBMTSxjGGGNiYgnDGGNMTCxhGGOMiYklDGOMMTH5/9Yuao3q3z55AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29096ba87f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# figure out our target fc given a desired band pass range\n",
    "bp_range = [6, 14]\n",
    "# convert from Hz to fc\n",
    "fcs = bp_range/sr  \n",
    "print(fcs)\n",
    "\n",
    "# build the low pass filter\n",
    "fc=fcs[0]\n",
    "b=0.1\n",
    "win = 'blackman'\n",
    "type_filt = 'lp'\n",
    "\n",
    "# make our function call to generate a lp filter\n",
    "n, sinc_func_lp = msf.define_filt(fc=fc, b=b, window=win, type_filt=type_filt)\n",
    "\n",
    "# now make the hp filter\n",
    "fc=fcs[1]\n",
    "b=0.1\n",
    "win = 'blackman'\n",
    "type_filt = 'hp'\n",
    "\n",
    "# make our function call to generate a hp filter\n",
    "n, sinc_func_hp = msf.define_filt(fc=fc, b=b, window=win, type_filt=type_filt)\n",
    "\n",
    "# now convolve the two filters\n",
    "bp_filt = np.convolve(sinc_func_lp,sinc_func_hp)\n",
    "plt.plot(bp_filt)\n",
    "\n",
    "# filter the signal\n",
    "filt_sig = np.zeros((r,c))\n",
    "for i in np.arange(r):\n",
    "    filt_sig[i,:] = msf.apply_filt(input_sig=seg_data[i,:], input_filter=bp_filt)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then FFT and look at filtered signal in the frequency domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAETCAYAAAAh/OHhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYVOWZ9/Hv3c0qIiAQTUACGhJFUcCOGtEB14gLjBON4obRhGiSicZM8poxUVHfTBJ9R+O4hRgB0eCoSZRxcBsj6rgCYkARDApqu4EICLK23O8fz2koi+rldFfVqXP697muuqrqnKdO3dVF982zm7sjIiLSXFVJByAiIumixCEiIrEocYiISCxKHCIiEosSh4iIxKLEISIisShxiIhILEocIiISixKHiIjE0i7pAEqhV69e3r9//6TDEBFJlTlz5nzo7r2bKpfJxNG/f39mz56ddBgiIqliZm82p1yiTVVmdpuZLTOzlxs4b2Z2vZktNrN5Zjas3DGKiMhnJd3HMRk4ppHzo4CB0W08cHMZYhIRkUYkmjjc/Ungo0aKjAFu9+A5oLuZfb480YmISCGV3sfRB3g753ltdOy9ZMIRkaZs3ryZ2tpaNmzYkHQo0oBOnTrRt29f2rdv36LXV3risALHCm4gYmbjCc1Z9OvXr5QxiUgjamtr6dq1K/3798es0K+wJMndWbFiBbW1tQwYMKBF10i6j6MptcBuOc/7Au8WKujuE929xt1revducjSZiJTIhg0b6Nmzp5JGhTIzevbs2aoaYaUnjunAWdHoqoOA1e6uZiqRCqekUdla+/0k2lRlZtOAkUAvM6sFLgPaA7j7LcAM4FhgMbAO+FYykRbPW2+BWtJEJM2SHlU11t0/7+7t3b2vu//B3W+JkgbRaKrvu/se7j7Y3VM/q+/II2Hx4qSjEJFSGjly5NZJyMceeyyrVq1q1fVmzpzJ8ccfX/Dc2LFj2Xfffbn22mtZuHAhQ4YMYejQocyZM4ebbrqpVe/bkEpvqsqcDRvg44+TjkJEGlJXV1fU682YMYPu3bsX9Zr13n//fZ555hnmzZvHj370I+677z7GjBnD3Llz6dmzZ8kSR6WPqsqcujr45JOkoxDJrqVLlzJq1CgOOeQQnnnmGfr06cP9999P586deemllzjvvPNYt24de+yxB7fddhs9evRg5MiRHHzwwTz99NOMHj2a+fPn07lzZxYuXMibb77JpEmTmDJlCs8++ywHHnggkydPBuD8889n1qxZrF+/npNOOokJEyZsF0/9Ekj33nsvt9xyCwCrV6+mf//+PP744zzyyCNcdtllbNy4kT322INJkyax44478tBDD3HhhRfSq1cvhg0rvGjG0UcfzbJlyxgyZAgnnngiN998M9XV1Tz55JPssssuvP766wwZMoSjjjqKq6++ung/ZHfP3G3//ff3StWrl/uDDyYdhUjpLFiwINH3X7JkiVdXV/vcuXPd3f3kk0/2qVOnurv74MGDfebMme7u/otf/MIvuOACd3cfMWKEn3/++VuvMW7cOD/llFN8y5Ytft9993nXrl193rx5/umnn/qwYcO2XnvFihXu7l5XV+cjRozwv/3tb1uvN2vWLHd3/+IXv+jLly/feu1Nmzb5IYcc4tOnT/fly5f7oYce6mvXrnV391/96lc+YcIEX79+vfft29dfe+0137Jli5988sl+3HHHFfyse++999bnl112mV999dUFz+Ur9D0Bs70Zf2PVVFVmqnFIW2NW/FtTBgwYwJAhQwDYf//9Wbp0KatXr2bVqlWMGDECgHHjxvHkk09ufc0pp5zymWuccMIJmBmDBw9ml112YfDgwVRVVbH33nuzdOlSAO6++26GDRvG0KFDeeWVV1iwYEGTsV1wwQUcfvjhnHDCCTz33HMsWLCA4cOHM2TIEKZMmcKbb77JwoULGTBgAAMHDsTMOOOMM5r50y4PNVWVmRKHtDVecMpuaXXs2HHr4+rqatavX9/ka7p06VLwGlVVVZ+5XlVVFXV1dSxZsoRrrrmGWbNm0aNHD84+++wm50ZMnjyZN998kxtuuAEILT5HHXUU06ZN+0y5l156qaKHNKvGUWabN8O6dUlHIdL2dOvWjR49evDUU08BMHXq1K21j5b4+OOP6dKlC926deODDz7gwQcfbLT8nDlzuOaaa7jjjjuoqgp/eg866CCefvppFkdDLdetW8drr73GnnvuyZIlS3j99dcBtksszdG1a1fWrFkT+3XNoRpHmanGIZKcKVOmbO0c33333Zk0aVKLr7XffvsxdOhQ9t57b3bffXeGDx/eaPkbbriBjz76iMMOOwyAmpoabr31ViZPnszYsWPZuHEjAFdddRVf/vKXmThxIscddxy9evXikEMO4eWXC+4+0aCePXsyfPhw9tlnH0aNGlXUznHzJOqRJVZTU+OVuJGTO1RVwYQJcOmlSUcjUhqvvvoqe+21V9JhSBMKfU9mNsfda5p6rZqqyqh+eLhqHCKSZkocZaTEISJZoMRRRkoc0lZksQk8S1r7/ShxlNHmzeFeiUOyrFOnTqxYsULJo0J5tB9Hp06dWnwNjaoqI9U4pC3o27cvtbW1LF++POlQpAH1OwC2lBJHGSlxSFvQvn37Fu8sJ+mgpqoyUuIQkSxQ4igj9XGISBYocZRRXR20b6/EISLppsRRRnV10L271qoSkXRT4iijzZthp51U4xCRdFPiKKO6OujaNWwfu2VL0tGIiLSMEkcZ1dVBhw7QubOaq0QkvZQ4ymjz5tA53qWLmqtEJL2UOMqorg7atVPiEJF0U+IoIyUOEckCJY4yUlOViGSBEkcZqcYhIlmgxFFG9Yljhx2UOEQkvbQ6bhnVLznSoYMSh4iklxJHGW3eHGoc1dXbFjwUEUkbJY4yqm+qqqqCTz9NOhoRkZZR4iij+qYqUOIQkfRKvHPczI4xs0VmttjMLi5wvp+ZPW5mc81snpkdm0ScxVDfVNWunRKHiKRXoonDzKqBG4FRwCBgrJkNyiv2c+Budx8KnArcVN4oi6e+qaq6ettugCIiaZN0jeMAYLG7v+Hum4C7gDF5ZRzYKXrcDXi3jPEVVX1TVXW1ahwikl5J93H0Ad7OeV4LHJhX5nLgETP7Z6ALcGR5Qiu++qYqUOIQkfRKusZhBY553vOxwGR37wscC0w1s+3iNrPxZjbbzGYvX768BKG2Xm5TlRKHiKRV0omjFtgt53lftm+KOhe4G8DdnwU6Ab3yL+TuE929xt1revfuXaJwW0dNVSKSBUknjlnAQDMbYGYdCJ3f0/PKvAUcAWBmexESR2VWKZqQOwFQiUNE0irRxOHudcAPgIeBVwmjp14xsyvMbHRU7MfAd8zsb8A04Gx3z2/OSgU1VYlIFiTdOY67zwBm5B27NOfxAmB4ueMqhfrEoXkcIpJmSTdVtSm5fRyaxyEiaaXEUUbq4xCRLFDiKCP1cYhIFihxlJGG44pIFihxlJGaqkQkC5Q4ykhNVSKSBUocZaSmKhHJAiWOMtJ+HCKSBUocZaT9OEQkC5Q4ykhNVSKSBUocZaRRVSKSBUocZaRRVSKSBUocZaSmKhHJAiWOMlKNQ0SyQImjjNTHISJZ0KLEYWY7mtnXzOy46HmP4oaVTfVNVe3aaTiuiKRXrMRhZruY2V3AR8D/AvdHp75nZovN7NBiB5glaqoSkSxoduIws97As8A3geeBFwGLTn8CfBF40MwGFzvIrFBTlYhkQZwaxwSgHzDa3Q8FHqg/4e7XAUcRtqL9eVEjzBCNqhKRLIiTOEYDf3b3BwqddPeZwJ+BrxUhrkxSU5WIZEGcxNELeKOJMrVA75aHk21qqhKRLIiTOGqBYU2UOTAqJwWoxiEiWRAncdwLHGFm3y100swuAg4B/lKMwLJIfRwikgXtYpT9v8BxwE1m9n2gGsDMJgP7A4OAxcAvixxjZuTux6F5HCKSVs2ucbj7GmA4cAvQH9iLMBz3LOBLwFRguLuvKn6Y2aCmKhHJgjg1Dtz9Y+D7ZvZD4CtAd2AtsMjdN5YgvkxRU5WIZEGsxFHP3T8FFhQ5lsxTjUNEsqDBxGFml7bwmu7uV7bwtZlVnyiqqpQ4RCTdGqtxXF7gmEf31sA5i+6VOPLUN1OBEoeIpFtjiePEvOedgX8nJIdrgWcIix3uCHwVuAj4FDin+GGmX30zFShxiEi6NZg43P3+3Odmdg3QHjjA3ZfkFX/BzP5MWPjwTMLKuZKjfiguKHGISLrFmQB4BmGtqvykAYC7v0dYq+qkYgSWNbk1Ds3jEJE0i5M4OtH0KKydKNz/0SAzO8bMFkX7eVzcQJlvmtkCM3vFzP4Y5/qVYuNG6NgxPFaNQ0TSLE7imAN8w8z2KnTSzL5GqG082dwLmlk1cCMwijDzfKyZDcorMxD4GWFy4d7AhTFirhibNilxiEg2xJnHcTnwP8Bz0TIjc4A1QDfCjPLTgU3E24/jAGCxu78BEO0uOIbPzhH5DnCju68EcPdlMa5fMTZuhA4dwmMlDhFJs2YnDnd/yszGADcB/8y2obkQmqcWAOe4+8sx3r8P8HbO81rCCru5vgxgZk8T1se63N0fyr+QmY0HxgP069cvRgjloaYqEcmKuEuOPGRmXyL8cd8P6AGsBOa4+6wWvH9D80HyYxwIjAT6Ak+Z2T75a2K5+0RgIkBNTU3+NRK3aZNqHCKSDbGXHHH3LYS9x58twvvXArvlPO8LvFugzHPuvhlYYmaLCImkJYkqMapxiEhWNDtxRAsbNoe7+380s+wsYKCZDQDeAU4FTssrcx8wFphsZr0ITVdN7URYcVTjEJGsiFPjuI5ty4rky12KxIFmJQ53rzOzHwAPE/ovbnP3V8zsCmC2u0+Pzh1tZgsIM9N/4u4rYsRdEXJrHJrHISJpFidxfKuB4zsQ9uM4k9BB/qM4Abj7DGBG3rFLcx47YTmTi+Jct9LkDsetqgL3cLNYs15ERJIXZ1TVlMbOm9n1wFzgsOhecuQOxzULyePTT7fNJhcRSYs4EwAb5e5vAncD5xXrmlmSW+MA9XOISHoVLXFENvDZUVISya1xgBKHiKRX0RJHtBTJacDSYl0zS1TjEJGsiDMc98WGThH25BhASETaxKkA1ThEJCvidM0OaeTcJuBl4PfufmPrQsqm3OG4oMQhIukVZ1RVsftD2pTcCYCguRwikl7NTgZmdpaZ7dtEmeFm9ovWh5U9qnGISFbEqUVMJix53ph/IuydIXnyaxxKHCKSVg02VZnZ2cDovMNjzWxoAy/pQJj8l7rlQMpBNQ4RyYrG+jgeAq4njJiCsAbVntGtIRuBSxs532ZpOK6IZEWDicPd3zezPQhrURlhRdrrgN8WKg5sBj6Mlj+XPBqOKyJZ0eioKndfXv/YzL4FvBQtLSIxqcYhIllRtEUOpXGqcYhIVjTWOf4R8G/ufnXO8+Zwd+9ZjOCyJL/GoXkcIpJWjdU4PiZ0duc+r7i9vNNCNQ4RyYrGOsf7N/Zc4lEfh4hkhZYRKRPVOEQkKxrr48if/Nds0V7hkkMTAEUkKxrr47iP+H0aFr2musURZZSWHBGRrGgscVyBOsOLRjUOEcmKxjrHLy9jHJmnGoeIZEWcjZwAMDMDDgL2A7oBHwLPu/vLRY4tU/JrHJrHISJpFStxmNmRwM3A7vWHons3s9nAeHf/WxHjywwNxxWRrIiz5/jBwAOEIbx3Ai8A7wPdgYOB04C/mtnB7r6oBLGmmobjikhWxKlxTADqgMPcfVbeuVvN7FbgMUKn+ilFii8zVOMQkayIMwHwq8A9BZIGAO7+DHAvcHgxAssS95A42rffdkyJQ0TSKk7i2ASsaaLMxzGv2SZs3hySRlXOT0aJQ0TSKs4f+WnA6WbWv9BJM9uVsOf4Xa0PK1vy+zdAiUNE0itOH8fvgeHAXDO7HngSeAfoTGjG+glhxvgT+cuVtPUlSPKH4oISh4ikV5zEMY8wk9yAX7D9rPL6obnT8o61+SVI8if/geZxiEh6xUkcJVmCxMyOIexjXg3c6u6/aqDcScA9wFfdfXax4ygl1ThEJEvibB17ebHf3MyqgRuBo4BaYJaZTXf3BXnlugI/BJ4vdgzlUKjGocQhImmV9AioA4DF7v6Gu28idKyPKVDuSuA3wIZyBlcsqnGISJbEXXLkBMIM8f5AxwaKubvv38xL9gHeznleCxyY955Dgd3c/QEz+5c48VYK1ThEJEviLDnyPeA/2NYJ3pA4/SCFrrX19WZWBVwLnN2M+MYD4wH69esXI4TSU41DRLIkTlPVBYSVcEcCnd29qoFbnBFUtcBuOc/7Au/mPO8K7APMNLOlhFV5p5tZTf6F3H2iu9e4e03v3r1jhFB6+cuNgBKHiKRXnMTRF7jD3Z90941Fev9ZwEAzG2BmHYBTga1zPtx9tbv3cvf+7t4feA4YncZRVYWG4ypxiEgaxUkcrwKfK+abu3sd8APg4ej6d7v7K2Z2RWv2PK8069dD586fPVZdrXkcIpJOcedx3GNmd7r7Q8UKwN1nADPyjl3aQNmRxXrfclqzBrp2/ewxNVWJSFrFmccxPVpq5L/NbBGwBCjUZOXu/o1iBZgFa9bAjjt+9pgSh4ikVZxRVacDFxFGQu0Z3Qop+uzytFu7tnCNY2OxeopERMooTlPVz4F1hOTxNPBJSSLKIDVViUiWxEkc/QlrSf2+RLFk1po10KvXZ48pcYhIWsUZVbUE6NBkKdmOahwikiVxEse1wKlmNqRUwWRVocShZdVFJK3iNFV9ArwOvGBmTwGLKdzP4e7+42IElxVr124/qkqJQ0TSKk7iuCPn8WHRrRAHlDhyFKpxdOkCn2h4gYikUJzE0VCiyLdPSwLJskKJY6ed4OOPk4lHRKQ14kwAfKKhc2bWhbDO1LcJ+4/f2PrQsqNQ4ujaNRwXEUmbWPtx5DOzAwnJ4hSgC2Fy4KoixJUpqnGISJbEThxmtjNwJiFhDCIkiy3AY8Ak4M/FDDALCnWO77STahwikk5xlhw5gpAs/pEwn6N+E6bHgbPd/e2GXtuW1dXB5s3br47btatqHCKSTo0mDjP7AvAt4BzCzHEDlgH3AH8kLD2ySEmjYfULHFreXodqqhKRtGowcZjZfwFfj8qsBm4H7gIedfctUZlyxJhqhfo3ICSTTz6BLVugKs40TBGRhDVW4ziOMMHvSuA3Rdz1r00ptDIuhCVHOncOyaPQeRGRStXY/3UfBjoBlwPvmNk0M/vHaItXaaZCe3HUU3OViKRRg4nD3UcR9hn/KfAOYcjtn4BlZnabmR1dnhDTraGmKtBcDhFJp0Zb1939A3f/f+6+H7AfcB1hT46zgQcJy4sMM7OvlTrQtGoscajGISJp1OxuWXefHy1e2Bc4FvhPYANwAPC/Zva6mV1hZl8pTajp1FTiUI1DRNIm9nged9/i7g+5+2nALoS5HU8Rhuv+HHilqBGmXEOd46C5HCKSTq0aCOrua939NncfCQwALgX+XozAskKd4yKSNUWbQeDub7n7Ve6+V7GumQWrV4cEUYiaqkQkjTT1rMQ++AB23bXwOTVViUgaKXGU2PvvN5w4VOMQkTRS4iixxhKHahwikkZKHCXWVI1DiUNE0kaJo4S2bIHly+Fznyt8Xk1VIpJGShwltGJFSA4dGljdS01VIpJGShwl1FgzFUD37rByZfniEREpBiWOEmoqcXzhC/DOO+WLR0SkGJQ4SqipxLHLLmGC4IYN5YtJRKS1Ek8cZnaMmS0ys8VmdnGB8xeZ2QIzm2dmj5nZF5OIsyXefz8kh4ZUVUGfPlBbW76YRERaK9HEYWbVwI3AKGAQMNbMBuUVmwvUuPu+wL3Ab8obZcs1VeMA2G03eFs7totIiiRd4zgAWOzub7j7JsKe5mNyC7j74+6+Lnr6HGFZ91RobuJ4663yxCMiUgxJJ44+QO7/t2ujYw05l7CB1HbMbLyZzTaz2cuXLy9iiC1XWxs6wBujGoeIpE3SicMKHPOCBc3OAGqAqwudd/eJ7l7j7jW9e/cuYogt4w4vvwx77914OSUOEUmbpBNHLbBbzvO+wLv5hczsSOASYLS7byxTbK3y3nuh87uxznFQ4hCR9Ek6ccwCBprZADPrAJwKTM8tYGZDgd8RksayBGJskXnzYN99wQrVqXIocYhI2iSaONy9DvgB8DDwKnC3u78S7V0+Oip2NbAjcI+ZvWRm0xu4XEWZPx8GD266nDrHRSRt2iUdgLvPAGbkHbs05/GRZQ+qCObPhxEjmi63886hVrJsWcOLIYqIVJKkm6oyq7k1DjMYMgTmzi19TCIixaDEUQJ1dbBoUdMjquoNHarEISLpocRRAq+9FpYS6dKleeWHDoWXXiptTCIixaLEUQLNbaaqN2yYahwikh5KHCUQN3HsuWeYZb52beliEhEpFiWOEpg/P8zhaK527WCffeDFF0sXk4hIsShxlMC8efFqHACHHAJPPVWaeEREikmJo8jWrAlzMvbYI97rRoyAJ54oTUwiIsWkxFFkL7wQmqmqq+O97tBD4dlnYfPm0sQlIlIsShxF9vDDcPTR8V/Xo0eopcyeXfyYRESKSYmjyB56CI45pmWvPfzw8HoRkUqmxFFE774bhtV+9aste/0ZZ8Dtt8OWLcWNS0SkmJQ4iuiBB+DII8Pw2pYYOhS6dYOZM4salohIUSlxFNHEiXD22S1/vRmccw78/vdFC0lEpOiUOIpk1ixYsQK+/vXWXWfcOHjkEW3uJCKVS4mjlR59FP7yFxg/Hs4/P/4w3HzduoXk8dvfFic+EZFiS3wjpzT7+9/hpJPggAPgwgvhzDOLc90LL4SaGth/fxg7tjjXFBEpFtU4WuHRR+HEE8P9uHFQVaSfZr9+8NhjcMklISktWFCc64qIFIMSRys88kjLJvs1x+DBYV+P0aPhX/6lNO8hItISShwttHlzGDZ7ZAl3RG/XDn7yE3jlFXj66dK9j4hIHEocLfT00zBgAHzuc6V9n44d4eKL4brrSvs+IiLNpcTRQpMmFa8zvCmnnhqaxVauLM/7iYg0RomjBVauhPvvh7POKs/79egBRx0F99xTnvcTEWmMEkcL3HwzjBoFvXqV7z3POgv+8AdwL997iogUosQR06xZob/hl78s7/sedxysWqV1rEQkeUocMaxaBaecArfcEjrGy6m6OnSSX3WVah0ikiwljmZyh29/O/zP/5/+KZkYzjgDli8PiUtEJClacqSZbr0VXn8d7rwzuRjat4c//QmGDw9LsB90UHKxiEjbpRpHMyxcCP/6rzBtWphXkaSBA8Oy69/8JixblmwsItI2KXE0YePGsNDgVVfBnnsmHU0wZkzY9+OEE+CTT5KORkTaGiWORrjD974Hu+8elk2vJBMmwKBBoc/lww+TjkZE2pLEE4eZHWNmi8xssZldXOB8RzP7z+j882bWv1yx/du/wdy5MGVK2J2vkpiFfpcDD4Rhw8LkQI22EpFySDRxmFk1cCMwChgEjDWzQXnFzgVWuvuXgGuBX5cjtmnT4He/C/uI77hjOd4xvupq+PWvYepUuPLKsODiU08pgYhIaSVd4zgAWOzub7j7JuAuYExemTHAlOjxvcARZqX7/7873HADXHRRSBpf+EKp3ql4RoyAF1+Ek0+G73wH9tsPrr8+jALbsiXp6EQka5IejtsHyN1duxY4sKEy7l5nZquBnkDRW/aXLIHTT4cNG8Lqt7vvXux3KJ127eC88+C734W//hVuvz3URlatgi99CXr3hi5dYIcdwi1306n6NJybjht6LFKI/o1UjksvhV13Le17JJ04Cv1zy29oaU4ZzGw8MB6gX79+LQqmV6+wb/hpp7V+7/CkmMERR4QbwJo1YYvbFStg3bowCmvdum3NWfn3jT0WKUT/RipLOaYMJJ04aoHdcp73Bd5toEytmbUDugEf5V/I3ScCEwFqampa9E+5a9fyLZVeLl27hs5zEZFiSbqPYxYw0MwGmFkH4FRgel6Z6cC46PFJwF/d9X8cEZGkJFrjiPosfgA8DFQDt7n7K2Z2BTDb3acDfwCmmtliQk3j1OQiFhGRpJuqcPcZwIy8Y5fmPN4AnFzuuEREpLCkm6pERCRllDhERCQWJQ4REYlFiUNERGJR4hARkVgsi1MizGw58GYLX96LEixnUmb6DMlLe/ygz1AJyh3/F929d1OFMpk4WsPMZrt7TdJxtIY+Q/LSHj/oM1SCSo1fTVUiIhKLEoeIiMSixLG9iUkHUAT6DMlLe/ygz1AJKjJ+9XGIiEgsqnGIiEgsShw5zOwYM1tkZovN7OKk42mKme1mZo+b2atm9oqZXRAd39nMHjWzv0f3PZKOtSlmVm1mc83sgej5ADN7PvoM/xktu1+xzKy7md1rZguj7+NrafoezOxH0b+hl81smpl1qvTvwMxuM7NlZvZyzrGCP3MLro9+t+eZWUXsUtPAZ7g6+nc0z8z+Ymbdc879LPoMi8zs68lErcSxlZlVAzcCo4BBwFgzG5RsVE2qA37s7nsBBwHfj2K+GHjM3QcCj0XPK90FwKs5z38NXBt9hpXAuYlE1Xy/BR5y9z2B/QifJRXfg5n1AX4I1Lj7PoQtDk6l8r+DycAxecca+pmPAgZGt/HAzWWKsSmT2f4zPArs4+77Aq8BPwOIfrdPBfaOXnNT9Her7JQ4tjkAWOzub7j7JuAuYEzCMTXK3d9z9xejx2sIf6z6EOKeEhWbAvxjMhE2j5n1BY4Dbo2eG3A4cG9UpKI/g5ntBPwDYe8Y3H2Tu68iXd9DO6BztMvmDsB7VPh34O5Psv1uoA39zMcAt3vwHNDdzD5fnkgbVugzuPsj7l4XPX2OsDMqhM9wl7tvdPclwGLC362yU+LYpg/wds7z2uhYKphZf2Ao8Dywi7u/ByG5AJ9LLrJmuQ74KbAlet4TWJXzy1Pp38XuwHJgUtTcdquZdSEl34O7vwNcA7xFSBirgTmk6zuo19DPPK2/3+cAD0aPK+YzKHFsYwWOpWLImZntCPwJuNDdP046njjM7HhgmbvPyT1coGglfxftgGHAze4+FPiECm2WKiTqBxgDDAC+AHQv3U3JAAAHJklEQVQhNO3kq+TvoClp+zeFmV1CaI6+s/5QgWKJfAYljm1qgd1ynvcF3k0olmYzs/aEpHGnu/85OvxBfTU8ul+WVHzNMBwYbWZLCc2DhxNqIN2jZhOo/O+iFqh19+ej5/cSEklavocjgSXuvtzdNwN/Bg4mXd9BvYZ+5qn6/TazccDxwOm+bc5ExXwGJY5tZgEDo5EkHQidUNMTjqlRUV/AH4BX3f3fc05NB8ZFj8cB95c7tuZy95+5e19370/4mf/V3U8HHgdOiopV+md4H3jbzL4SHToCWEB6voe3gIPMbIfo31R9/Kn5DnI09DOfDpwVja46CFhd36RVaczsGOD/AKPdfV3OqenAqWbW0cwGEDr6X0giRtxdt+gGHEsYxfA6cEnS8TQj3kMIVdV5wEvR7VhCH8FjwN+j+52TjrWZn2ck8ED0eHfCL8Vi4B6gY9LxNRH7EGB29F3cB/RI0/cATAAWAi8DU4GOlf4dANMIfTKbCf8bP7ehnzmhmefG6Hd7PmEEWaV+hsWEvoz63+lbcspfEn2GRcCopOLWzHEREYlFTVUiIhKLEoeIiMSixCEiIrEocYiISCxKHCIiEosSh6SWmZ1tZt6MW0WtsSSSdu2aLiJS8Z4AZjZyfmGZ4hBpE5Q4JAtmuvvlSQch0laoqUpERGJR4pA2wcwuj/o7joh2tdsY7aK2Y3R+VzO7ycxqo3NLzOzXZta1wLX6mdnt0c5ta8zsvujYUjObmVNucvSeQwpcw83spbxjZmbnmdmLZrbezFaa2XQzG5pXbmT0+rPN7Bwzm29mG6LYrzGzHQq83xgLu0WuMrMPzex/zOwfonM7mNnHZvZWtFZV/mtvi95vYIwfuWSYEoe0NXcC64H/IDRxrTWzfoRFLs8j7ENxLWEtoJ8CT0R7awAhaRD2PDmTsMnOLYQdIx8nLEfeGlMIO9N1iK57D2GDqGfM7PAC5X8QlXsZuB7YAPw4eryVmf2MsH7WIMLKvXcRVu99zMyO9rCQ3p8IK68emvfaTsA3gOfc/e+t/HySFUkv8qWbbi29AWcTFnmcCVzewK1/VPbyqOwLQFXedf6bsInU8XnHfxi95jc5x6ZFx87NOdYZeLY+lpzjk6NjQwrE7sBLOc9Pjo7dCbTLOT4AWEFYAK9DdGxkVLYO+FpO2W6EZcTXA12iY18mLKD3KrBrTtkvAWuBedHzw6Jr3pwXZ31c30v6+9atcm6JB6Cbbi295SSOxm4jo7L1ieNnedf4fJQ0Hihw/SrCkuPLo+edgXXAggJlD25l4ngoOta7QNmronPHRc/rE8djBcreF53bM3p+SfT8tAJlvwv8CGhPWD12KfAh0D6nzP3AJqBn0t+3bpVz06gqyYIJ3vxRVUvzng8j/NHsaWaFrrEJ2M3M+hCWSu9MaNbK9xzhf/ottT+hqen7BboZ9ozuhxBqR/VeK3Cd1dF9x+h+v+j+2fyC7v673Odmdgch0RwN/LeZ7QwcA8xw9xXN+xjSFihxSFuT/8e9e3R/UHRryM7ATtHjtfkn3X2Lma1sRVzdCb+PlzURQ66NBcrU75NQn316RPfN2VL4dkLiGEtIUN8k9LdMbcZrpQ1R57i0dfVJ4Ep3t0Zu84H6xPD5/ItEo5F65h3O/yNeX3a7UU9RHG83EcOPW/H5Co0O62xmW/8GuPtrhJrTGDPrSEgcK4EHWvC+kmFKHNLWzYvuawqdNLMJZnZxtJ3wa8AaoCba6z3XPmxrHqq3KbrfMe/4Hg3E0dfMdi0Qw/FmdpWZ7VfgdU2ZH90fUODc9cC6aBvSerdH8Y4h7DB5t7sXqtlIG6bEIW2auy8BngRGmdlJuefM7EzgUuAYd9/k7nXAJMKw1Z/mlOsA/LrA5euXOjk+p2wVoTko32RCzeSG6Hr15T9PGKL7Mwo0kTXDHwmd/5eY2dYakZntQahRvBH9DOrdRUh4vyF0mquZSrajPg4RGA88BdxjZg8S5kV8hfAH/yPgezllfwEcDlxlZocR/kd/BNC3wHX/CFwJ/CT6Q/0GoeO5B2G0Vq7JwGjCnIn5ZvYw4ffzm4QmsIvd/fW4H8zdF0ad/lcAfzOz/yIkqFOBTsC38sqvjMp8g5BUno77npJ9qnFIm+fuiwijmn4P7AtcQBiNNBX4qrsvyCn7MWGS3A2E5qnvAu8QRh/lX/cDwtDZx4BRwHcI8ymGs62/pL6sAydF770O+DZwCrAAONHdC9Vomvv5royu9TZwFnAGYT7LCHd/vsBL7o7u72jpe0q2Wfj3KiKtYWbdCcngCXcfmXA4rWJmvyQ0jQ1098VJxyOVRzUOEdnKzPoC5wKPK2lIQ5Q4RAQzO93M5hLW6OpNmGkvUpASh4hAWAurH2Hm+Xfd/cmE45EKpj4OERGJRTUOERGJRYlDRERiUeIQEZFYlDhERCQWJQ4REYlFiUNERGL5/0YPDIPYXXaIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29096cf31d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r,c = filt_sig.shape\n",
    "\n",
    "# compute FFT - all trials at once\n",
    "sp = np.fft.rfft(filt_sig, axis=1)                    \n",
    "\n",
    "# compute amplitude\n",
    "amp = (np.abs(sp) / c) * 2\n",
    "\n",
    "mean_amp = np.mean(amp, axis=0)\n",
    "\n",
    "# get freq axis\n",
    "freq = np.fft.rfftfreq(c, 1/sr)   # Sample spacing is 2nd arg\n",
    "\n",
    "plt.plot(freq, mean_amp, 'b', linewidth=1)          # plot amplitude (power = np.absolute(sp)**2)\n",
    "\n",
    "plt.legend(['normalized fft'], loc='best', fancybox=1)  \n",
    "plt.xlabel('Frequency', **fig_font)\n",
    "plt.ylabel('Amplitude', **fig_font)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## You can also make a band-stop filter using a similar logic\n",
    "* Instead of convolving the lp and the hp filters, you can add them together\n",
    "* Use this to selectively get rid of the 60Hz line noise in the data by generated a band-stop filter over 50-70Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VNd99/HPbzZtCBAgFiNAxsZm8QYIr2A7XhLHSezYcfrYT9YmsdumbpumSeq+0peTJm3TOG3a52mcNI7TOEmbekuTED9OXBtv4BVhA2E1AoQRYJDYYSTNdp4/ZkaM9hkx0ujOfN+vFy9m7lzpHi6ar8787rnnmHMOEREpLr5CN0BERPJP4S4iUoQU7iIiRUjhLiJShBTuIiJFSOEuIlKEFO4iIkVI4S4iUoQU7iIiRShQqANPmjTJ1dfXF+rwIiKetGbNmjbnXO1g+xUs3Ovr62lsbCzU4UVEPMnMdmWzn8oyIiJFSOEuIlKEFO4iIkVI4S4iUoQU7iIiRUjhLiJShBTuIiJFSOEukiEaT/Do6t0kElp+UrxN4S6S4eXtB/nSz9fz5u4jhW6KyGlRuItkONERS/7dGStwS0ROj8JdJEM4kgz19ojCXbxN4S6SoT0aByAciRe4JSKnR+EukiEd6gp38bqswt3MbjCzrWbWZGb39PH6TDN7zszeNLP1ZnZj/psqMvzSod6ucBePGzTczcwP3A+8F5gP3GFm83vs9tfAo865hcDtwHfz3VCRkZCutavnLl6XTc/9YqDJObfDORcBHgZu7rGPA8amHo8D9uaviSIjp6ssE9UFVfG2bBbrmA7sznjeAlzSY5+vAv9jZn8CVAHX5aV1IiOsXWUZKRLZ9Nytj209b9+7A3jIOVcH3Aj81Mx6fW8zu8vMGs2ssbW1NffWigwzXVCVYpFNuLcAMzKe19G77PJp4FEA59wrQDkwqec3cs494JxrcM411NYOugSgyIgLR9Vzl+KQTbivBuaY2ZlmFiJ5wXR5j33eBq4FMLN5JMNdXXPxnFMXVFVzF28bNNydczHgbuApYDPJUTEbzexrZnZTare/AO40s3XAfwGfdM5p5iXxHJVlpFhkc0EV59yTwJM9tt2b8XgTcEV+myYy8rouqEYV7uJtukNVJIN67lIsFO4iGU6ma+6aFVI8TuEukqG96yYm9dzF2xTuIimRWIJYagUmlWXE6xTuIinpXvu4iiCRWIK4ltoTD1O4i6Sk55OZOCaUfK6x7uJhCneRlHQpZlJVGaC7VMXbFO4iKekwP9VzV7iLdyncRVLCCncpIgp3kZR0jX1iuiyjOd3FwxTuIinpsswk9dylCCjcRVJOlWXKuj0X8SKFu0hK+q7UiVXJnrtGy4iXKdxFUtJzuavnLsVA4S6S0lWWqdJNTOJ9CneRlPZInFDAx5jyQNdzEa9SuIukhCNxKkN+gn4fQb9pZkjxNIW7SEo4Eqcy6AegIuhXz108TeEuktIejVERSoZ7ZSigmrt4msJdJCVZlknW2ytDfo2WEU9TuIukhCPxrp57RUhlGfE2hbtISnskTlVXWUY9d/E2hbtISjgS6yrLVIQCGi0jnqZwF0lpzyjLVAb9XXesiniRwl0kJRxNjnMHlWXE+xTuIim6oCrFROEuAsQTjkgsQWUwWXOvKguo5y6epnAX4dQkYemyTEXQT3s0TiLhCtkskSFTuItwapKwioyaO0BHTL138SaFuwinpvut7BHuKs2IVyncRegd7hUhTfsr3qZwFyE5aRicCnX13MXrsgp3M7vBzLaaWZOZ3dPPPr9nZpvMbKOZ/Sy/zRQZXr177ulw141M4k2BwXYwMz9wP3A90AKsNrPlzrlNGfvMAf4KuMI5d9jMJg9Xg0WGQzrcK4Kn7lAFlWXEu7LpuV8MNDnndjjnIsDDwM099rkTuN85dxjAOXcgv80UGV7tvS6oJvs9KsuIV2UT7tOB3RnPW1LbMp0DnGNmL5nZq2Z2Q74aKDISTpVl0hOHpcoymjxMPGrQsgxgfWzreWdHAJgDXA3UASvN7Dzn3JFu38jsLuAugJkzZ+bcWJHhkq6t9xznrsnDxKuy6bm3ADMyntcBe/vY51fOuahzbiewlWTYd+Oce8A51+Cca6itrR1qm0XyrndZRqNlxNuyCffVwBwzO9PMQsDtwPIe+/wSeBeAmU0iWabZkc+GigyncDRO0G8E/cm3RIXCXTxu0HB3zsWAu4GngM3Ao865jWb2NTO7KbXbU8BBM9sEPAd80Tl3cLgaLZJv7ZF410gZgJDfh99nGi0jnpVNzR3n3JPAkz223Zvx2AGfT/0R8ZzMVZgAzIzKoOZ0F+/SHaoiJMsv6Tp7WkXI33XnqojXKNxF6L7EXppWYxIvU7iL0F/PXQt2iHcp3EVIjpapCHW/BFWppfbEwxTuIiRvVqoM9lWWUc1dvEnhLkI/ZRmNlhEPU7iLkAx3XVCVYqJwFyE9zl0XVKV4KNyl5CUSjo5oop8Lqqq5izcp3KXktUe7TxqWVhnyE47GSd6ALeItCncpeT2X2EurCPlxDjpjiUI0S+S0KNyl5LX3WGIvLT00UnV38SKFu5S8cGr+mMpeNff0Unuqu4v3KNyl5A1UlgEtki3epHCXktdVlunjgiqoLCPepHCXkjdYz13hLl6kcJeSl66p9x4Kmay5a0538SKFu5S8U4tj976JCdRzF29SuEvJ67cso6GQ4mEKdyl56TtU+7ugqtEy4kUKdyl54UgMv88I+bu/HU6Nc1e4i/co3KXkhSNxKoN+zKzb9vKgDzM0eZh4ksJdSl5fi2MDmJkW7BDPUrhLyetrFaa09MyQIl6jcJeSl1yFKdDnaxVaJFs8SuEuJa892nsVprTKYEATh4knKdyl5A1UlqnQOqriUQp3KXntkXivudzTKlWWEY9SuEvJG/SCqsJdPEjhLiVv4Auqga47WEW8ROEuJa89MtAFVb8uqIonKdylpDnnCEd1QVWKj8JdSlpnLIFzvScNS9MFVfGqrMLdzG4ws61m1mRm9wyw321m5sysIX9NFBk+XdP99jNapqosQCzhiMQSI9kskdM2aLibmR+4H3gvMB+4w8zm97FfNfCnwGv5bqTIcDm1ClM/F1SDmvZXvCmbnvvFQJNzbodzLgI8DNzcx35fB+4DOvLYPpFh1d/i2GldqzFpqT3xmGzCfTqwO+N5S2pbFzNbCMxwzj2Rx7aJDLv+VmFK0yLZ4lXZhLv1sc11vWjmA/4Z+ItBv5HZXWbWaGaNra2t2bdSZJiEB+25pxbJVriLx2QT7i3AjIzndcDejOfVwHnA82bWDFwKLO/roqpz7gHnXINzrqG2tnborRbJk/bowDV3LZItXpVNuK8G5pjZmWYWAm4HlqdfdM4ddc5Ncs7VO+fqgVeBm5xzjcPSYpE8yr4so5q7eMug4e6ciwF3A08Bm4FHnXMbzexrZnbTcDdQZDh1lWUGmDgMVJYR7+n7s2gPzrkngSd7bLu3n32vPv1miYyM9kF67pVBLZIt3qQ7VKWknSrL9L8SE6Cl9sRzFO5S0tojMcygPNj3W+FUWUY1d/EWhbuUtHBqoQ6zvkb8nqrFqywjXqNwl5I20IyQAD6fUR706YKqeI7CXUpaeyTe7w1MaZWhgHru4jkKdylp4Uisa0RMfyqCmtNdvEfhLiUtnFXP3d91J6uIVyjcpaS1D7A4dpoWyRYvUrhLSQtnEe5aak+8SOEuJS0ciVHRzw1MackLqirLiLco3KWkhSPxfpfYS1PPXbxI4S4lLauhkEEtki3eo3CXkuWcIxyNU1WmC6pSfBTuUrIi8QTxhOt30rC0ilBAPXfxHIW7lKz2QeZyT6sM+YnEE8TiiZFolkheKNylZA22ClNapab9FQ9SuEvJGmxx7LQKrcYkHqRwl5LVPshCHWlaJFu8SOEuJSt9Y9Kgd6h2LbWnG5nEOxTuUrLSNfRsJg4DlWXEWxTuUrIGWxw7TWUZ8SKFu5SsrtEyg83nrnAXD1K4S8lKL3qdzUpMgOZ0F09RuEvJynmcu3ru4iEKdylZ4SzvUNU4d/EihbuUrPZonPKgD5/PBtwvPSWweu7iJQp3KVnhSGzQG5gAAn4fIb9P4S6eonCXkhWOxActyaRVhPxdF2BFvEDhLiUrm8Wx0zSnu3iNwl1KVjaLY6dVhPyaFVI8ReEuJSubJfbSKkNaak+8ReEuJSscze6CKiTvYtXEYeIlCncpWeEceu4V6rmLx2QV7mZ2g5ltNbMmM7unj9c/b2abzGy9ma0ws1n5b6pIfrVH4l1j2AejC6riNYOGu5n5gfuB9wLzgTvMbH6P3d4EGpxzFwCPA/flu6Ei+ZbzBVWFu3hINj33i4Em59wO51wEeBi4OXMH59xzzrlw6umrQF1+mymSf8kLqlnW3EN+2jVaRjwkm3CfDuzOeN6S2tafTwO/6esFM7vLzBrNrLG1tTX7VorkWSyeIBJP5DDOXRdUxVuyCfe+Jt5wfe5o9lGgAfhWX6875x5wzjU45xpqa2uzb6VInqXHrGddlgn66YgmSCT6/NEXGXWyCfcWYEbG8zpgb8+dzOw64MvATc65zvw0T2R4pEe+5DLOHVBpRjwjm3BfDcwxszPNLATcDizP3MHMFgLfJxnsB/LfTJH8ynYu9zTN6S5eM2i4O+diwN3AU8Bm4FHn3EYz+5qZ3ZTa7VvAGOAxM1trZsv7+XYio0K6fl4xyBJ7aV2rMSncxSOy+sl2zj0JPNlj270Zj6/Lc7tEhlW2i2OndfXctdSeeITuUJWSlGtZRotki9co3KUkhXO+oKqyjHiLwl1KUnuqvJL1xGHquYvHKNylJA29LKOau3iDwl2KTiyeoO3EwLdaDHmc+yA99wPHO3Sjk4wKCncpOv/2wnbe9a3nBwzirp57trNCpoZMDlSWORKOcOV9z/HjV5qzbqvIcFG4S9F5Yv0+jnfGWLv7SL/7hCNxQn4fAX92b4GKLO5QbWw+TEc0wZO/25dbg0WGgcJdikrL4TBb3jkOQGPzoX73a4/EqCzLrtcOEAr4CPhswJr76l3J463ZdZhDJyNZf2+R4aBwl6KyYnNy9osJVSFW7zrc737hHBbqSBtsTvfG5sNMqAqRcPDcFs3CIYWlcJei8szm/cyeVMWN50/ljV2HifdzcTMczX6JvbSBFsnuiMZZ33KE2xbXMWVsGSu27M+57SL5pHCXonGiM8ZrOw5x7bzJLKmfwInOGFveOdbnvu2ReNZj3NOSc7r3He7rW44SjTsurp/ANXOn8OJbbURiiZz/DSL5onCXorHyrVYi8QTXzptCQ/0EIFkq6Us4Esu5514R7L8sszpV3188q4br5k1O/qLZeTCn7y+STwp3KRpPb97PuIogDbNqmD6+gjPGlXeFbk+5rJ+allwku+8Lqo3Nh5gzeQw1VSGuOHsS5UEfz2xSaUYKR+EuRSGecDy/tZV3nVvbNbyxoX4Cq5sP4VzvuvtQwr2/C6rxhKNx1+GuTwvlQT9Lz67lmc0H+jy2yEhQuEtRePPt5PDDa+dN6dq2pL6G/cc6aTnc3mv/9kg867nc0/q7oPrW/uMc74ixpL6ma9t18yaz50g7W/cfz+kYIvmicJei8PTm/QR8xlXnnlqbt6vuvqt3aSYciQ2hLBPocz739Hj6JanjAVwzdzKASjNSMAp3KQorNh/gktkTGFse7Np2zpRqqssDrO7joupQyzJ99dxXNx9mytgy6moqurZNHlvOhXXjeGazxrtLYSjcxfN2HTxJ04ETXDt3Srftfp+xeFZNrztV4wlHZyyR+zj3fkbLNDYfYkn9BMys2/br5k1hXcsRWo9rvXgZeQp38bx07/i6eVN6vbakfgJv7T/B4YzpANLzwwxltEx7NN7tIumeI+3sPdrRrSSTdu28KTjdrSoFonAXz3tm037OmTKGmRMre73WMCt5kXNNxlQEXYtj53gTU0UogHPQET11c1L6U0FDxsXUtHnTqjljXDlPb1bdXUaewl087Wh7lNXNh7qNksl04YzxBP3WNakXZCyOnePcMpV9LNixuvkQY8oCzJ06ttf+Zsa186awalsbHQPMJikyHBTuMiyOd0T56IOv8cjqt4f1OC+81Uos4bhu3uQ+Xy8P+jl/+rhud6rmugpTWl+LZDc2H2bRrBr8Puvza66bP4X2aJxXtg/v3ao/emknn3potdZ4lS4Kd8m7RMLxuYfXsqqpjf+7omlYVyZasXk/E6tCXDSjd1kkbcmZE1jfcqSr95zr4thplT3mdD8ajrJ1/3GWzOr/2JfOnkBVyM8zw1iaicYTfOfZJp7dcoC//Pl63TglgMJdhsG3n36LFVsOcNU5tew50s6rO4an1xqNJ3huywHeNXdyvz1ngCWzJhCNO9a3HAUyyjI5TxzWvef+xtuHce7UePq+lAX8LJtTy4phvFv1uS0HOHgywlXn1LJ83V6+/+KOYTmOeIvCXfLqifV7+c5zTdy+ZAbf/9hiqssDPLamZViO1dh8mGMdsX5LMmmLUz3r9Dwz6Zp5zmWZrqX2kl//evMhAj7johnjB/y6a+dN5p1jHWzc2/cMlafrsTUt1FaX8eAnGnj/BdP45m+3aISOKNwlfzbuPcoXH1vP4lk1/M3NCygP+vnAhWfwmw37ONYRzfvxVmzeT8jvY9mc2gH3q6kKMWfymK6RLemyypDLMqmee2PzIc6bPm7Q73PN3MmYMSylmbYTnTy35QC3LpxO0O/jvtsuYN7Usfzpw2+yvfVE3o8n3qFwl7w4eKKTu36yhnEVQb730UWUBZKBd9viuuS6ouvzu66oc45nNu/nsrMmUlU2eHmloX4CjbsOk0i4IV9QzSzLdETjrNt9tNt8Mv2ZOKaMRTNrhiXcf/nmHmIJx22L61JtDPDAxxcT9Pu48yeNw/JLVbxB4S6nLRpP8Nn/fIO2E5088PHFTK4u73pt4YzxnFVblffSzPbWkzQfDA9akklbUl/D8Y4Ybx04fircc5w4rCKj575hz1Ei8cSA9fZM186bzIY9x3jnaEdOxxyIc47H17Rw4YzxzJlS3bW9rqaS735kEW8fDPO5h9f2uxqVFDeFu5y2rz+xidd2HuIfPnQ+F9R1rz+bGR9umMGaXYfzWiZYkeoFX9PP+Pae0neQrm4+THvXTUy5TxwGyZp7er6ahgFGymRK3z2bz+X3Nuw5xpZ3jvPhVK8906WzJ/KVmxbw7JYD/NP/bM3bMcU7FO5yWv7r9bf5ySu7uOvK2dyysHfIANy6cDo+g8fz2HtfsfkA86eNZfr4isF3BupqKpgytozG5kOEI3ECPiMUyO3Hv6ssE43T2HyIs2qrmDimLKuvnTN5DDMnVHYt4J0Pj63ZTSjg4wMXntHn6x+9ZCZ3XDyT7z6/nV+v25u344o3KNxlyNbsOsy9v9rAsjmT+Msb5va73+Sx5Vx1Ti3//UZLXkoER8IRGncdyrokA8lPEA31E2hsPkw4kvvi2ABlAR9mEO6M07jrcJ/zyQx0/GvnTWZVU1tebjTqiMb51dq9vGfBVMZVBPvcx8z4m5sW0DCrhi8+vo6Ne4+e9nHFOxTuMiThSIw/f2QtU8eV8507Fg04zhzgww0z2H+sk5XbWk/72C81HSTh4Kpzsw93gCWzathzpJ3trSdyvpgKybCsDPpZ13KEo+3RrOvtaVefO5lILMHr/Sz9l4sVmw9wtD3aZ0kmUyjg43sfXcy4iiCff2QdnTHdwVoqsgp3M7vBzLaaWZOZ3dPH62Vm9kjq9dfMrD7fDZXR5b7fbuXtQ2H+8bYLGVfZd88x07XzJjO+MpiXC6urmlqpLg9wYd24nL4uHcav7zyU8w1MaRWhAK/vTC/OkV29Pe3i+gmE/D5W5eEX3GNrdjNtXDlXnD1p0H1rq8v4xq3ns3X/cf51RdNpH1u8YdCfcDPzA/cD1wMtwGozW+6c25Sx26eBw865s83sduCbwP8ajgbLwNpOdLJl33E27zvG5n3H2LTvGMfao/zxNWdzx5KZ+AbpYWfjtR0HeejlZj55eT2XzJ6Y1deUBfx88KLp/Oy1tzkSjjC+MjSkYzvnePGtNi6bPbFrrdRszZ1azZiyACc6Y1TkOGlYWmXIT9uJTmqry5g5ofcslAOpCPlpqK9h5ba2IR077Z2jHbz4ViufvfrsQT8xpV0zdwofWlTH917YznsWTOX8HH8x9iUWT/DjV3bxgxd3MHFMiLlTxzJvWjXzp41l3rSx1FQN7f9Y8iOb7svFQJNzbgeAmT0M3AxkhvvNwFdTjx8HvmNm5jTJxbA42Rljz5F2Wg6HaTnczu5DYbbuP8Hmfce6LQwxZWwZ86aNZUxZgC//YgOPr2nh7z54PvPP6D2DYbbCkRhf+vl6Zk6o5Es3nJvT1962uI6HXm5m+bq9fPyy+iEdv/lgmD1H2vnDq2bn/LUBv4+FM8ezclvbkMoycOqi6pL6ml6Lc2Rj6ZxJ3PfbrRw43tFtyGgu/vvNFhKOrrHt2br3/fNZ1dTKFx5bx6//ZGnOF5Qzrd19hC//4nds3HuMS2dPIOj38cJbrfz8jVOfzKaOLWfutGrOmVLNjJoK6moqqaupYHpNxZA/OUn2sjnD04HdGc9bgEv628c5FzOzo8BE4PS6KH14rHE3P1y1s9u29JvMup5nvpb6m1Mbo/EE0XiCWMIRizsi8QSxeIJo3A16wa8i5GdMWSD5pzxAdervMWUBKoJ+HOAcJJzDOYcj+TjhoDOaIByJcTISJ9wZ40RnjHAkTjgSoz0SpyzopyLopzLkp7IsQGXQT1WZn4pQgFg8kQr0dg5lLDwBybrq2bVjWDZnUlevad60sUxI9Zycc/xy7R7+9onNfOA7q/jUFfV87rpzsrr5p6dvPbWVXQfDPHzXpTm/QRecMZa5U6t5fE3LkMM9XdIY7K7U/iypn8DKbW1DuqAKp4ZPNszKrd6eduWcWu777VZeamrrd3TRQNJj25fU11A/qSqnrx1XGeQbt57Ppx5q5F+f3cZfvDu3X86QnGL5H5/ayn+8tovJ1WV87yOLuOG8qV3vwdbjnWx551jqk2PyE+TLTQeJxBPdvs+kMSGm11RSN74Cs+S9A+n3QjjjcSS1YlZVWYDKUICq1HtjTJmfylCAUMCHz8BnyXe4mSUfWzIPwtE4JzqS77UTHTGOd8Y40RnlREes27z8fQn4jIDfCPp9BP2+rscBX/JvgOQ7PCmzK5t+7Lqed8+Vu685m/df0Pcop3zJ5t3ZV/ekZwJmsw9mdhdwF8DMmTOzOHRvYyuCzJpY2cfJ633YnvsktzkCfh+hjP+soN8I+JL/gX4f/fbInHO0R+Oc7IxzvCP5Q7L/eAfbW08FtcGpHy6zUz94lixNVGb8cNZUVXb9sFYE/XTGUj/UnXHC0eQvgL1HorRH45glb045b/o46jJ6QXU1FUyqKhuw3GJm3LKwjmvOncI3n9rCD1bu5In1+/jqTQt49/wpWfdAX995iIdebuYTl83i0izLMT3b8eGGGXz9iU1sfec4506tHvyLeli5rY26mgpm9bEwRzbSi2qcfs99aOE+P/VLd+W2oYX7G28fYUfrSf7wyrOGdPxr5k7h1kXT+e7zyfLMedOzK8845/j1+n18/YlNHDzRyScvr+fz159DdXn36y211WXUVtd2++WbSDjaTnSy+/CpT5stqceb9x3DkTyvlSE/4ytDnDE+GdyVIT+hgI/2aLozlAz8o+1R9h1p52RnjEg80dWZSrhkO7s6VyTvTaguD1BVluyUTR9fQXV5sjxXHvQN+F6PJ0h1AhNEYo5YItGtM9izUwndO5bpV051ME/tM2YIHatcZXOEFmBGxvM6oOeg2fQ+LWYWAMYBvYYEOOceAB4AaGhoGFLJ5j0LpvKeBVOH8qUlb1xlkL+/5Xw+tKiOL//id/zBT9dw3bzJfOUDC5gxSP24PRLni4+vo66mgi8NMOxxMB+86Ay+8eRmHmvczV+/f35OXxuLJ3hl+0Hef+G0IZVEAC6aMZ6Az4Z+QTWYDJ1503L/xQTg8xmXnzWRVdvacM7l/O94fM1uKoJ+brxg2pCOD/CV9y9g1bY2vvDYOpbfPXh5ZmfbSe791QZWbmvjgrpx/OiTS7L+pQDJf/PkseVMHlveNYmbDL9sim6rgTlmdqaZhYDbgeU99lkOfCL1+DbgWdXbR6/Fs2r49Z8s5cs3zuOlpoNc++0X+Lv/t4kj4Ui/X5Mux9z3oQuHVM5JmzimjGvnTeaXa/cQjQ/8sbindS1HON4ZY+nZQyvJQLIn94X3nMvNFw3tI/GHG+r44nvOzflibqZlcyZx4Hgnb+3P7Y7d9kicX6/bx43nTzutnl/6l/yWd47znWe39btf24lOvvKrDVz/7RdY+/YRvnbzAn7x2StyCnYpnEF/QlI19LuBpwA/8O/OuY1m9jWg0Tm3HPgh8FMzayLZY799OBstpy/o93HnlbN53wXT+PbTb/Hgqp08sno3n33X2Xzy8nrKM0aTrG4+xI9e3snHL5vFZWflXo7p6bbFM3hq436e39rK9fOzmz4AkiUZM7j8NNvwh1cNraQB5OVT49JUyWLlttacSlNPbXyHE52xnC+k9uW6+VO4deF07n9+O+/uUZ452Rnjh6t28v0XttMRS3D7khn82XVzhnwBWArDCtXBbmhocI2NjQU5tvS25Z1jfPM3W3huayvTxpXz59efw4cW1RGJJXjv/3mRuHP89s+uPK1ee1o0nuCyb6xg0cwaHvh4Q9Zfd9v3XiYST7D87qWn3YZCu+Yfn2fmxEoe+v2Ls/6ajzz4Km8fCvPCF96VlyGtR8IRrv/nF5lYFWL53Usxg0dW7+ZfntlG24lOblgwlS/ecC5n1Y457WNJ/pjZGufcoG8cjUcSAOZOHcuPfv9iXtl+kH/4zWa+9Ph6frhyJ2dOqqL5YJif3XlJXoIdkp8ablk4nR+91Ezr8eSY8cEc74jy5u4jQxoCORotmzOJRxtb6IzFu6ZHHsjuQ2Fe3n6Qz117Tl6CHWB8ZYhv3HI+n/lJI59/dC2b9h5jR9tJltTX8P2PLVZ93OM0/YB0c9lZE/nlH1/B/f97EZ2xOL/d+A67ADAjAAAItklEQVQfu3QWl581+J2Qubj94pnEEo7/eHVXVvu/uuMQ8YQ7rXr7aLJ0Ti3t0Thv7DqS1f4/frkZnxkfbjj9kkym6+ZP4ZaF03li/T78PuPBjzfw6B9cpmAvAuq5Sy9mxvsumMa7F0zhpaa2vNTZezqrdgzXzJ3MT1/dxR9dfVa3Gn9fVm1rpSLoZ9GsgZe084pLZ0/A7zNWNbUOen6PdUR5ePVu3nf+NM7IchbMXPz9LefzwYXTueKs3O/6ldFL/5PSr6Dfx9XnTs6qbDAUn1l2JodORvjFm3sG3XfltjYumT1h2Noy0qrLgyycMT6rqQgeeX03Jzpj3LlseEpSFSE/V51Tq2AvMvrflIK5bPZEFpwxlgdX7iAxwJ3Be460s6PtJEuzmCTLS5bOmcTv9hzl8Mn+h6BG4wl+9NJOLjlzQl7mg5HSoXCXgjEz7lw2m+2tJ3n+rf4XsUhPOXDlOcVRb09bNqcW5+Dl7Qf73ec3G95h79GOYeu1S/FSuEtBve+CaUwbV84PXtzZ7z4rt7UxZWwZcyYX15C8C+vGUV0eYFVT31MAO+d4cOUOZtdWcc3c3OauF1G4S0EF/T4+eXk9r+w4yIY9vVcKSiQcL28/yBVnTxrylAOjVcDv47LZE1mZmoqgp9d3HmJ9y1E+vfTMvA1/lNKhcJeCu/3imVSF/L1m+wTYtO8Yh05GWDanuOrtacvmTKLlcDu7DoZ7vfaDlTupqQxy6xAmGBNRuEvBjasI8ntLZvDrdXvZd7S922svpurt2aw45EWZUxFk2tF6ghVb9vOxS2cNeXpiKW0KdxkVPnXFmSSc46GXm7ttX7WtjblTq4t2XpP6iZVMH1/Ra0jkD1ftJOj38bEhznsvonCXUWHGhEree940fvba25zojAHJWRAbmw8XbUkGkiOGrjxnEq9sP0gsNUvmoZMRHl/Twi0XTc9qagaRvijcZdT4zLIzOd4R49HVyYW/Xm8+RCSe6CpdFKulZ9dyvDPGupbkBeX/fHUXnbEEn1l2ZoFbJl6mcJdRY+HMGhpm1fDvL+0knnCs2tZKyO/j4iGueuQVl581EbNkCaojGufHr+zi6nNrmTNlaAuCiIDCXUaZzyybTcvhdp7a+A4rt7XRUF9T9BcUa6pCnD99HKuaWlm+di9tJzp105KcNoW7jCrXz5/CrImVfPvpt9jyznGWFnG9PdPSsyfxxttH+N4L25k7tfq0FyQRUbjLqOL3GZ+64kyaDiSXoFtWJFP8DmbpnEnEE46dbSe5c9nsorthS0aewl1GnQ831DGuIkhNZZAFZ4wtdHNGxOJZNVQE/UwZW8YHLhza+q4imTSfu4w6laEA/3Dr+bRH4yVz231ZwM+9H5jPlLFlhALqc8npU7jLqPTe86cVugkj7o6LZxa6CVJE1EUQESlCCncRkSKkcBcRKUIKdxGRIqRwFxEpQgp3EZEipHAXESlCCncRkSJkfS3MOyIHNmsFdg3xyycBbYPuNfLUrtyoXbkbrW1Tu3JzOu2a5ZwbdNKlgoX76TCzRudcQ6Hb0ZPalRu1K3ejtW1qV25Gol0qy4iIFCGFu4hIEfJquD9Q6Ab0Q+3KjdqVu9HaNrUrN8PeLk/W3EVEZGBe7bmLiMgAPBfuZnaDmW01syYzu6fQ7Ukzs2Yz+52ZrTWzxgK249/N7ICZbcjYNsHMnjazbam/a0ZJu75qZntS52ytmd1YgHbNMLPnzGyzmW00sz9LbS/oORugXQU9Z2ZWbmavm9m6VLv+JrX9TDN7LXW+HjGz0Chp10NmtjPjfF00ku3KaJ/fzN40sydSz4f/fDnnPPMH8APbgdlACFgHzC90u1JtawYmjYJ2XAksAjZkbLsPuCf1+B7gm6OkXV8FvlDg8zUNWJR6XA28Bcwv9DkboF0FPWeAAWNSj4PAa8ClwKPA7ant/wb80Shp10PAbYX8GUu16fPAz4AnUs+H/Xx5red+MdDknNvhnIsADwM3F7hNo4pz7kXgUI/NNwM/Tj3+MfDBEW0U/bar4Jxz+5xzb6QeHwc2A9Mp8DkboF0F5ZJOpJ4GU38ccA3weGp7Ic5Xf+0qODOrA94HPJh6bozA+fJauE8Hdmc8b2EU/MCnOOB/zGyNmd1V6Mb0MMU5tw+SoQFMLnB7Mt1tZutTZZsRLxdlMrN6YCHJXt+oOWc92gUFPmepEsNa4ADwNMlP00ecc7HULgV5X/Zsl3Mufb7+LnW+/tnMyka6XcC/AF8CEqnnExmB8+W1cO9rteRR8dsZuMI5twh4L/DHZnZloRvkAd8DzgIuAvYB/1SohpjZGODnwOecc8cK1Y6e+mhXwc+Zcy7unLsIqCP5aXpeX7uNbKt6t8vMzgP+CpgLLAEmAH85km0ys/cDB5xzazI397Fr3s+X18K9BZiR8bwO2FugtnTjnNub+vsA8AuSP/SjxX4zmwaQ+vtAgdsDgHNuf+oNmQB+QIHOmZkFSQbofzrn/ju1ueDnrK92jZZzlmrLEeB5krXt8WYWSL1U0PdlRrtuSJW3nHOuE/gRI3++rgBuMrNmkmXka0j25If9fHkt3FcDc1JXmkPA7cDyArcJM6sys+r0Y+DdwIaBv2pELQc+kXr8CeBXBWxLl3R4ptxCAc5Zqv75Q2Czc+7bGS8V9Jz1165CnzMzqzWz8anHFcB1JK8HPAfcltqtEOerr3ZtyfgFbSTr2iN6vpxzf+Wcq3PO1ZPMq2edcx9hJM5Xoa8iD+Gq840kRw5sB75c6Pak2jSb5MiddcDGQrYL+C+SH9ejJD/pfJpkjW8FsC3194RR0q6fAr8D1pMM02kFaNdSkh+J1wNrU39uLPQ5G6BdBT1nwAXAm6njbwDuTW2fDbwONAGPAWWjpF3Pps7XBuA/SI2oKcQf4GpOjZYZ9vOlO1RFRIqQ18oyIiKSBYW7iEgRUriLiBQhhbuISBFSuIuIFCGFu4hIEVK4i4gUIYW7iEgR+v8QRJU3ua6tOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29096b3e4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# figure out our target fc given a desired band pass range\n",
    "bstop_range = [50,70]\n",
    "# convert from Hz to fc\n",
    "fcs = bstop_range/sr  \n",
    "\n",
    "# build the low pass filter\n",
    "fc=fcs[0]\n",
    "b=0.1\n",
    "win = 'blackman'\n",
    "type_filt = 'lp'\n",
    "\n",
    "# make our function call to generate a lp filter\n",
    "n, sinc_func_lp = msf.define_filt(fc=fc, b=b, window=win, type_filt=type_filt)\n",
    "\n",
    "# now make the hp filter\n",
    "fc=fcs[1]\n",
    "b=0.1\n",
    "win = 'blackman'\n",
    "type_filt = 'hp'\n",
    "\n",
    "# make our function call to generate a hp filter\n",
    "n, sinc_func_hp = msf.define_filt(fc=fc, b=b, window=win, type_filt=type_filt)\n",
    "\n",
    "# now add the two filters\n",
    "bstop_filt = sinc_func_lp+sinc_func_hp\n",
    "plt.plot(bstop_filt)\n",
    "\n",
    "# filter the signal\n",
    "filt_sig = np.zeros((r,c))\n",
    "for i in np.arange(r):\n",
    "    filt_sig[i,:] = msf.apply_filt(input_sig=seg_data[i,:], input_filter=bstop_filt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then FFT and look at filtered signal in the frequency domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEYCAYAAACk+XocAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFNW5//HPM6wj+yKbIGgURY0YNO6aQbyKuUGJ/lQk8cZojImJ0STeuOtE/f2iRo3RXM2NC0KMoiYazCLGhdGYqHFDURFRZDWgwiAMMGzz/P441UxP0z109fT0Mnzfr1e9pqvqVPXTxdDPnHPqnDJ3R0REJI6KYgcgIiLlR8lDRERiU/IQEZHYlDxERCQ2JQ8REYlNyUNERGIrevIws05m9pKZvW5ms8zsqmj7JDObF21/zcz2LXasIiIStC92AO6+3sxGu/taM2sH/MPMpke7L3T3R4oZn4iIbK3oNQ8Ad18bvexESGgN0boVJyIREWlOSSQPM6sws9eBpcCT7v5ytOtaM5tpZjeZWYcihigiIkmslKYnMbPuwKPAecByd18WJY07gffd/dqiBigiIkAJ9Hkkc/dVZvYsMNbdb462bTSzScCP0x1jZqWT/UREyoi759w1UPRmKzPra2Y9oteVwNHAu2Y2INpmwHjgrUzncPeyXa666qqix7C9xl/OsSv+4i/lHn9LlULNYyAw2cwqCMnsQXf/q5k9bWZ9CZ3mM4HvFDNIERFpVPTk4e6zgFFpto8pQjgiIpKFojdbbe+qqqqKHUKLlHP85Rw7KP5iK/f4W6qk7rbKhZl5uX8GEZFCMzO8nDvMRaR1rV+/nssuu4wBAwZQUVGBmWlp40tFRQUDBgzgsssuY/369a3ye6Wah0gbN3bsWDp37szNN9/MzjvvTPv2Re/qlFa2adMmFi5cyI9+9CPq6+uZPn36VmVaWvNQ8hBp4zp06MCqVauorKwsdihSYOvWraN79+5s3Lhxq31KHkoeIs2KviSKHYYUSaZ//5YmD/V5iIhIbEoeIiIZjB49mnvuuQeA+++/n7Fjx+b1/AsWLKCiooKGhoa0+9977z1GjRpFjx49+NWvfkV9fT3jxo2jZ8+enHrqqXmNJS71nImIZGHixIlMnDgx7+c1y9xydMMNNzB69Ghee+01AO677z4++eQTamtrMTMmT57MXXfdxd///ve8x7UtqnmISJuwefPmYoeQdwsWLGDvvfdusj58+PAtCcfdm00+rUnJQ0SKapddduGmm25i5MiR9OrVi9NOO40NGzZs2X/nnXey++6707dvX8aPH8+///3vLfsqKiq4/fbbGT58OMOHD9+y7Y477mD48OH06NGDK6+8knnz5nHooYfSs2dPJkyYwKZNmwBYuXIl48aNo1+/fvTp04dx48axZMmStHFOnjyZI444AoCf//zndOvWje7du9O9e3c6duzImWeeCcCqVav41re+xaBBgxgyZAhXXHHFlg7rhoYGLrzwQnbccUd22203/vKXv2S8LmPGjGHGjBl873vfo3v37kycOJGrr76aqVOn0r17d26//Xa++93v8sILL9CtWzd69+7dgn+FHBR7Zsc8zAzpIpJZqf8fGTZsmB900EG+dOlSr62t9REjRvj//u//urv7008/7X379vWZM2f6hg0b/LzzzvMjjzxyy7Fm5sccc4zX1tZ6fX39lm0nnHCC19XV+TvvvOOdOnXyo48+2ufPn++rVq3yvfbay6dMmeLu7suXL/dHHnnE6+vrva6uzk855RQfP378lvNXVVX53Xff7e7u9957rx9xxBFbxb9o0SLfaaedfPr06e7ufsIJJ/h3v/tdX7dunX/yySd+0EEH+W9+8xt3d7/jjjt8xIgRvmTJEq+trfXRo0d7RUWFb968Oe21SX5/d/fq6mo//fTTt6xniilZpn//aHvO372qeYgIZvlZcnX++efTv39/evbsybhx45g5cyYQOqnPOussRo4cSYcOHfjZz37GCy+8wMKFC7cce+mll9KzZ086deq0ZdvFF19Mly5dGDFiBPvssw/HHHMMQ4cOpVu3bhx33HG8/vrrAPTu3ZuvfvWrdOrUiS5dunDJJZfw3HPPZR33unXrGD9+PBdccAHHHnssH3/8MdOnT+cXv/gFnTt3pm/fvlxwwQVMnToVgIcffpgLLriAQYMG0bNnTy655JLcL1qRqcNcRCj2MJD+/ftveb3DDjtsaZr66KOP2H///bfs69KlC3369GHJkiXsvPPOAAwePHir8/Xr12/L68rKyibnr6ysZNmyZUD48r/gggt44oknWLlyJe5OXV1d1n0JZ511FiNGjODCCy8EQp/Exo0bGThwINDYspOI9aOPPmLIkCFbjh86dOg236NUKXmISMkaNGgQCxYs2LK+Zs0ali9f3iRhtKTD+MYbb2Tu3Lm8/PLL7LjjjrzxxhuMGjUqq+Rx3XXXMXfuXP7xj39s2TZkyBA6d+7M8uXL0x4/cOBAFi1atGU9+bPlolid5aAOcxEpYRMnTmTSpEm8+eabrF+/nksvvZSDDz64yV/vLVFXV0dlZSXdu3dnxYoVVFdXZ3Xc448/zm233cYf//hHOnbsuGX7gAEDOOaYY/jhD3/I6tWrcXfmzZu3pSnslFNO4dZbb2XJkiXU1tZy/fXXtyj+/v37s3jx4rTTj7Q2JQ8RKarm/no+6qijuOaaazjxxBPZaaed+PDDD7f0H2Q6NnVbc+e/4IILWLt2LX379uXQQw/ly1/+clbHPvTQQ3z66aeMGDFiy11X5557LhDuytqwYQN77bUXvXv35uSTT2bp0qUAnH322Rx77LGMHDmSAw44gJNOOiljbNuKHcL12XvvvRkwYECTprpC0NxWIm2c5rbavmluKxERKRlKHiIiEpuSh4iIxKbkISIisSl5iLRx7du3Z926dcUOQ4pg3bp1rfbY4aInDzPrZGYvmdnrZjbLzK6Ktg8zsxfNbI6ZPWBmGtAokoMxY8YwYcIEPvjggy0TAkrbtmnTJj744AMmTJjAmDFjWuU9ip483H09MNrdvwDsBxxnZgcB1wM3ufsewErgrCKGKVK2pk2bxj777MPhhx9Op06dMDMtbXzp1KkThx9+OPvssw/Tpk1rld+rkhrnYWY7AM8B5wJ/Bga4e4OZHQxUu/tWj/EyjfMQEYnN2sI4DzOrMLPXgaXAk8AHwEp3TzybcTEwqFjxiYhIUyXRjxAliS+YWXfgUWBEumKZjk+ej6aqqoqqqqo8RygiUt5qamqoqanJ2/lKqtkKwMyuBNYCP6Fps9VV7n5cmvJqthIRiansm63MrK+Z9YheVwJHA+8AM4CTo2LfAFqn10dERGIres3DzD4PTCYksgrgQXf/v2a2CzAV6AW8Dnzd3bead1g1DxGR+Fpa8yh68mgpJQ8RkfjKvtlKRETKj5KHiIjEpuQhIiKxKXmIiEhsSh4iIhKbkoeIiMSm5CEiIrEpeYiISGxKHiIiEpuSh4iIxKbkISIisSl5iIhIbEoeIiISm5KHiIjEpuQhIiKxKXmIiEhsSh4iIhKbkoeIiMSm5CEiIrEpeYiISGxKHiIiEpuSh4iIxKbkISIisRU9eZjZYDN7xszeMbNZZnZetP0qM1tsZq9Fy9hixyoiIoG5e3EDMBsADHD3mWbWFXgVOAE4FVjt7jdv43gv9mcQESk3Zoa7W67Ht89nMLlw96XA0uh1nZnNBnaKduf8wUREpPUUvdkqmZkNA/YDXoo2fc/MZprZXWbWo2iBiYhIE0WveSRETVa/B86PaiC3A1e7u5vZtcDNwFnpjq2urt7yuqqqiqqqqtYPWESkjNTU1FBTU5O38xW9zwPAzNoDfwYed/dfptk/FPiTu++bZp/6PEREYmppn0epNFvdA7yTnDiijvSEE4G3Ch6ViIikVfSah5kdBjwHzAI8Wi4FJhL6PxqA+cA57r4szfGqeYiIxNTSmkfRk0dLKXmIiMRXlGYrM+tqZoeY2X9G671yDUBERMpPrORhZv3NbCqwAngemBbtOtfM3jezI/IdoIiIlJ6sk4eZ7Qi8AJxCGIfxGo2D+NYAQ4HHzezz+Q5SRERKS5yax0+BnYHj3f0Iwq21ALj7LcB/EMaNXJ7XCEVEpOTESR7HA4+4+5/T7XT3GuAR4JA8xCUiIiUsTvLoC8zbRpnFwI65hyMiIuUgTvJYDIzaRpmDonIiItKGxUkevwfGmNk56Xaa2Y+Aw4FH8xGYiIiUrqwHCZpZN+CfwF7A20A7YE/gt8D+0fb3gYPcfWWrRJs+Lg0SFBGJqaAjzM2sO/Az4HSga9Ku9cCDwIXu/mmuweRCyUNEJL6iTE9iZu2APYCeQB0wx93X5xpESyh5iIjEp7mtlDxERGJrtcfQmtmVOZ7T3f2aHI8VEZEykLHmYWYNaTYnCqfLVh5td3dvl5/wtk01DxGR+Fqt5gF8NWW9kvAoWAN+QbjzagWh4/yLwI+AzcCZuQYjIiLlIc6tujcC3wAOdPcP0+wfSJgs8TF3TzsWpDWo5iEiEl8hn+fxdcLcVlslDgB3/zdhbqv/k2swIiJSHuIkj84038wF0J30/SEiItKGxEkerwInmdmIdDvN7BBCreO5fAQmIiKlK06fxxHAU0A9cC8hmawGegCHAV8DNgKHuftbrRFshrjU5yEiElOhpycZC9wODKPxtl0ITVXvAGe6+79yDSYXSh4iIvEVfIS5mVUQpl4fCfQCaoFX3f3lXINoCSUPEZH4ND2JkoeISGytOUgw9Y1+kGVRd/fbYpx3MDAFGEAYZHinu99qZr0IM/UOBeYDp7j7Z9meV0REWk+cDvMGGqcgSZU8bUms6UnMbAAwwN1nmllXQkf8CcA3geXufoOZXQT0cveL0xyvmoeISEwFq3kQvszT2QHYjfCMj3eAH8YJwN2XAkuj13VmNhsYTEggX4qKTQZqgK2Sh4iIFF7e+jzMbCjwOnCtu9+c4zmGEZLEPsAid++VtG+5u/dJc4xqHiIiMRWy5tEsd19gZg8B3yFMoBhL1GT1e+D8qAaSdUaorq7e8rqqqoqqqqq4by8i0qbV1NRQU1OTt/Pl9W4rM7sFOMfdK2Me1x74M/C4u/8y2jYbqHL3ZVG/yAx332p0u2oeIiLxFXJixG0FMgKYSLgzKq57gHcSiSPyGHBG9PobwLRMByt3iIgUVpy7rV7LtIvwTI9dCMno/Ji36h5GmA9rFuGuLQcuBf4FPAQMARYCJ7v7yjTH+6ZNTruCPX5KRKT8FWyQYIYnCyZsAN4ljNH4n1yDyYWZ+YYNTocOhXxXEZHyVrAOc3fPWxNXvjU0l9ZERCTvsk4IZvZfZrbvNsocZmZXtDyseJQ8REQKK05t4l7CwL3mnAhcknM0Odq8udDvKCKyfcvYbGVmZwDHp2w+zcy+kOGQjsBoYHl+Qsueah4iIoWVscM8GlvxHuFOKsg8r1Wy9cC57j4pbxFug5l5ba3Ts2eh3lFEpPy1Woe5uy81s88R5q4yYB5wC/DLdMUJTxH81N035hpMrlTzEBEprGbvtnL3TxKvzeybwEx3X9DqUcWk5CEiUlhxbtWd3JqBtISSh4hIYTXXYb4C+Jm7/zxpPRuebvbb1qTkISJSWM3VPFYROsCT10tyFindqisiUljNdZgPa269lKjmISJSWCU75UgcSh4iIoXVXJ9H6gDBrLn7Y7kemwslDxGRwmquz+OPxO/jsOiYgk6QruQhIlJYzSWPqynRDvJUSh4iIoXVXId5dQHjaBElDxGRwsp6kGCCmRlwMDAS6AF8Crzk7m/lObas6VZdEZHCipU8zOxo4A5g18Sm6Keb2SvAt939jTzGlxXVPERECivr5GFmhwJ/Jtze+zvCM8aXAj2BQ4GJwDNmdqi7z2mFWDNS8hARKaw4NY+fApuA0e7+csq+u8zsLuBpQkf7qXmKLytKHiIihRVnkOAXgYfTJA4A3P2fwO+Bo/IRWBxKHiIihRUneWwAVm+jzKqY58wLJQ8RkcKK80X/APA1MxuWbmf05METgaktDyseJQ8RkcKK0+dxJ3AY8LqZ3Qo8BywBKglNWv9NGFn+bOrUJs1NV2JmdwNfAZa5+77RtquAs4GPo2KXuvv0TOfQrboiIoWV8RnmWxU0a6Dpc8xTD0y33QjP98g4XYmZHQ7UAVNSksdqd785i7j8+eedww7L6mOIiAit+AzzNFpluhJ3f97MhqbZlfWHUrOViEhhxXkMbXUrxpHO98zsdOAV4Mfu/lmmgkoeIiKFFXt6kgK5Hbja3d3MrgVuBs7KVPiee6qZMSO8rqqqoqqqqhAxioiUjZqaGmpqavJ2vqz7PADMbBxhJPkwoFOGYu7u+8cKIjRb/SnR55Htvmi/P/WUM2ZMnHcUEdm+FazPw8zOBW5j230RufSLWPJ5zWyAuy+NVk8Emp10UXdbiYgUVpxmq/MJM+ieTJhFd30+AjCz+4EqoI+ZLQSuAkab2X5AAzAfOKe5c6jPQ0SksOIkj8HAr939uXwG4O4T02yeFOccSh4iIoUVZ4T5bKBfawXSEkoeIiKFFSd5XA2cYmZjWyuYXCl5iIgUVpxxHo9F05L8xczmAB8C6fo93N1PyleA2VDyEBEprDh3W30N+BHhrqg9oyWdvI9C3xYlDxGRworTYX45sJaQQP4BrGmViHKgW3VFRAorTvIYBtzl7ne2Uiw5U81DRKSw4nSYfwh0bK1AWkLJQ0SksOIkj18AE6LBeyVFyUNEpLDiNFutAT4A/mVmfwfeJ32/h7v7j/MRXLaUPERECitO8rgv6fXoaEnHASUPEZE2LE7yyJQsUu2TSyAtoeQhIlJYcQYJPptpn5l1ASYA3yI8z/x/Wh5a9nSrrohIYbXoYVBmdhAhYZwKdCEMIFyZh7hiUc1DRKSwYicPM+sNnE5IGnsREkYD8DRhNtxH8hlgNpQ8REQKK870JGMICWM8YbxH4uFNM4Az3H1R/sPLjpKHiEhhNZs8zGwQ8E3gTMIIcwM+Bh4G7idMUzKnmIkDlDxERAotY/Iwsz8Bx0ZlPgOmAFOBJ929ISpTiBi3SclDRKSwmqt5/CdhEOA1wA35euxsa1DyEBEprOamJ3kC6AxUA0vM7AEzG29mJTe/VWvdqnvHHfDSS61zbhGRcpYxebj7cYTnlv8EWEK4HfcPwMdmdo+ZHVOYELettWoeTz0Ff/9765xbRKScNTsxorsvc/eb3H0kMBK4hfBMjzOAxwlTkYwys0NaO9DmtFbyqKuDhQtb59wiIuUs61l13X1WNOHhYODLwINAPXAg8LyZfWBmV5vZHq0TamZKHiIihRVnSnYA3L3B3ae7+0SgP2Hsx98Jt/JeDryd1wizoOQhIlJYsZNHMnevc/d73L0K2AW4Epgb5xxmdreZLTOzN5O29TKzv5nZHDN7wsx6NHcOJQ8RkcJqUfJI5u4L3f1adx8R89BJhPEkyS4GnnL3PYBngEuaO0FrJY81a2D5cli7tnXOLyJSrvKWPHLl7s8DtSmbTwAmR68nE6ZEyai1btWtq4MBA2BRUcfPi4iUnqInjwz6ufsyAHdfCuzYXOHWqHk0NIQax557qulKRCRVi6ZkLxXPPltNdXV4XVVVRVVVVYvPuXYtVFbCsGFKHiJS/mpqaqipqcnb+czd83aynIMwGwr8yd33jdZnA1XuvszMBgAzMvWlmJlfdJFz3XX5jWnpUhg5Er7znbD+05/m9/wiIsVkZrh7zhMUlkqzldE4xTvAY4SBiADfAKY1d3BrNFvV1UHXrtC3L6xYkf/zi4iUs6InDzO7H/gnMNzMFprZN4HrgP8wsznA0dF6Rq2ZPCorYd26/J9fRKScFb3PIxpsmM7R2Z6jNe62WrMmJI/OnaG+Pv/nFxEpZ0WveeRDa9Y8OndWzUNEJJWSRwZ1ddClS2i2Us1DRKQpJY8MVPMQEclMySOD5A5z1TxERJpS8shANQ8RkcyUPDJQzUNEJLM2kTxa41Zd1TxERDJrE8mjNWoeiXEeqnmIiGxNySODxK26GiQoIrK1NpE8WrPZStOTiIhsrU0kj9Wr83/ORPLo2BE2bmy9pxWKiJSjNpE8alOfQ5gHieRhpqYrEZFUSh4ZrF4dkgcoeYiIpGoTyaM1nrdRWwu9e4fX6vcQEWmqTSSPfNc83ENC6tUrrKvmISLSVJtIHg0N+f1yX7MGOnQISQNU8xARSdUmkkevXvmtfaxY0dhkBap5iIikajPJI5/9HitWQJ8+jeuqeYiINNUmkkfv3vmteSxfrpqHiEhz2kTyaO1mK9U8RESaajPJI9/NVqp5iIhk1iaSR76brVTzEBFpXptIHrrbSkSksNoXO4DmmNl84DOgAdjo7gemK9erF7z/fv7ed8UK2HPPxnXVPEREmirp5EFIGlXu3my9QjUPEZHCKvVmKyOLGNXnISJSWKWePBx4wsxeNrOzMxVSzUNEpLBKvdnqUHdfamY7Ak+a2Wx3fz610AMPVDN7NlRXQ1VVFVVVVS1609RBgqp5iEi5q6mpoaamJm/nM3fP28lak5ldBax295tTtvvGjU6PHrBsWeMzOHLlHpJFbW34CfCb38Arr4SfIiJtgZnh7pbr8SXbbGVmO5hZ1+h1F+AY4K10Zdu3h733hjfeaPn7rl4dzpdIHKCah4hIqpJNHkB/4Hkzex14EfiTu/8tU+H99oOZM1v+pu+/D5/7XNNt6vMQEWmqZPs83P1DYL9sy3/hC/Dqqy1/3/ffh913b7pNNQ8RkaZKueYRS75qHnPnwm67Nd2mmoeISFNtJnnsuy+88w5s2NCy87z//tbJo2/f0BkvIiJBm0keXbrA5z8Pz291I288c+du3Wy1227wwQfhcbciItKGkgfAccfB44+37Bzpah5du4ZxH4sWtezcIiJtRZtKHl/+Mvz1r7kfv2pVuFV34MCt9w0fDu+9l/u5RUTakjaVPA44AD75BD78MLfj584Nt+lWpLkqe+yh5CEiktCmkkdFBZxxBlx3XW7Hz5gBhx2Wft/w4TBnTs6hiYi0KW0qeQBccgk8+ii8/Xb8Y6dPh7Fj0+9Ts5WISKM2lzx69QoJ5KKL4h1XVwcvvQRHHZV+/x57qOYhIpLQ5pIHwLnnhjEfM2Zkf8zTT8MXvwjduqXfv8susGZN7v0pIiJtSZtMHp06wU03wTnnhDuotmX9erjsMvj2tzOXadcOTjoJHnoof3GKiJSrNpk8AL76VRg9Gr71rTC4b/NmOOss+O1vty577bVhYOCppzZ/zlNPhalTWydeEZFyUjbP88jEzDzTZ1i3LvRh7LdfmF5kwQJYujQMBExMuT5/Puy/P7z5Juy0U/PvtXkzDB0Kf/gDHHRQfj+HiEghtdnneeRDZSVMmxYe7LT//vDss+FL/+bocVIbNsB558EPfrDtxAGh6er660Pz1saNrRu7iEgpa9M1j3Tmz4cjjoAJE8I8WDvuCA8+2PThT81xh698JUxhcsstYEl5e84c+Pe/YdQo6N493ucQESmkltY8trvkAfDuu3DbbXDIITBxYvoR5c2prQ0JqKoKrroqJKBp00L/ym67hf1PPglDhsQ7r4hIoSh55JA88uHjj6G6Gu67DwYPDv0rDz4IBx4YmsVuvDHUTEaPDslFRKSUKHkUKXkkrFoFs2aFvpT2Sc9lfO45uPzysK99ezj6aDj88DBD78aNYdm8OcwEnPrYWxGR1qbkUeTksS3u8NFH8Oc/wxtvhFl7O3YMy6ZNYSqV738frrgidMiLiBSCkkeJJ49tWbo0jB/ZYQf43e/Cc0NERFqbbtUtcwMGwFNPwYgRsM8+YV6ut98OSWXdumJHJ4VUX1/sCESyp5pHCZk1C6ZMCR3v9fVhLq1+/UJzVrt24a6wXr3CaPj99w+P3R00KNzhpSav8vbss3D22eF2b8v5b0GR7LXpZiszGwvcQqgh3e3u16cp02aSR6r6+tBfkpheZdOmcBvwu+/Cq6+GGsqSJaGWsuOOYVLHHj2gT5+w9OsHK1eGub523jmMjt9558YnJXbsGMp37tx4/g4d4t+6XK7Wrg3NhaXg/PPh1lvh9dfDjAgira3NJg8zqwDeA8YAHwEvAxPc/d2UcmWdPGpqaqiqqmrROVavhk8/DT8/+wyWLw/rH38cksP69bBwYZieZeHCMFULhBH2K1eGpOEeEkcigXTunH6prGy6vnJlDcOGVVFZCT17hgTVrRt06RLuLOvSJXxBd+sW+nM6d275NcuH2bPhgANqWLy4il69ihuLO+y6a6hBHn54GDuUjXz87hST4i+uliaP9tsuUjQHAnPdfQGAmU0FTgDebfaoMpOPX8Bu3TJPJb8t7qHW0a5daC5xD0mlvj79sm5d0/WpU2vYb78q6utD0nrmmZDE1qwJy9q14efq1bBiRbhtOVEz6t278XW6JbG/V6/814buvRfWravhoYeqOOec/J57WxoaQnJfvz5c+9mzw/bLLw+PE9h//7A9UeNMXSoqQoJ/9NEaVqyoYtWqULvs3r3xLr5Nm0LZ9u1Dwu7UKfz71teHMoltHTqE7alLRUXT1+3bh7Lt2zc2oyYvifLZeuUVuPvuGo48sqpsa7rlnjxaqpSTx07AoqT1xYSEInlk1nR8iln4UunUKdRatuXdd+E738nuvdxDIlm+vHFZsaLx9YIF8NprTfcvXx4ST3JzXDbJp0uX8GXXsePWTXGbN4fBnWPGwOTJpE0emzeHRJm8JJJnXV2IafXqxtfZ/tywIZy/R49wjdu1CzW0Cy8Mj0A+8ED49a9DvIkv5dQv6oaGME5o1qzws3v3cN5Vq0JCSv6S37gxbKuvD9e/c+dQdv368FkStc50S0ND489Nm8K5EkkpdXHfOtaKisakkvyzQ4ewf8OG0LTatWtjzIlklnjvDh3CdUr8jiYSWuJ1um2ZXidvy5QsM5VJfJbkz/Tmm+F3NvmzpZZJXp8wITwzqK0o5eSR7u+Y8m2fEszCF0XXrqF5K1ubNoXmtdSkklgWLdo6Ga1dG76cNmwIX3rt2jUmEvfw1/3BB4ck0rNnY3wNDY1fqpWVTZdEs13Xro21vcTrxGdKXk9XpmPH8D4dO6bqa1qBAAAK7UlEQVT/rJMmZX9dqqvDUgoSNdjkJTn5JH42NITENXBgeBTCOec0/hslklNDQ+MNIBs3hsSXOEeihTr5ZzavU7clny8RW3KZTMk0eVm9Go48sum21HLJ69nOn1cuSrnP42Cg2t3HRusXA57aaW5mpfkBRERKXFvtMG8HzCF0mP8b+BdwmrvPLmpgIiJSus1W7r7ZzL4P/I3GW3WVOERESkDJ1jxERKR0lelNcmEAoZm9a2bvmdlFxY5nW8xssJk9Y2bvmNksM/tBtL2Xmf3NzOaY2RNmlsU9TsVjZhVm9pqZPRatDzOzF6P4HzCzkq3NmlkPM3vYzGab2dtmdlA5XX8z+6GZvWVmb5rZ78ysYylffzO728yWmdmbSdsyXm8zu9XM5prZTDMr+lDJDPHfEP3+zDSzP5hZ96R9l0TxzzazY4oT9ZZYtoo9ad+FZtZgZr2TtsW+9mWZPKIBhL8CjgX2Bk4zsz2LG9U2bQJ+5O57AYcA34tivhh4yt33AJ4BLilijNk4H3gnaf164KYo/pXAWUWJKju/BP7q7iOAkYQxQ2Vx/c1sEHAeMMrd9yU0OZ9GaV//SYT/o8nSXm8zOw74nLvvDpwD/LqQgWaQLv6/AXu7+37AXBrj3ws4BRgBHAfcblbUiWbSxY6ZDQaOBhYkbcvp2pdl8iBpAKG7bwQSAwhLlrsvdfeZ0es6YDYwmBD35KjYZGB8cSLctugX78vAXUmbjwL+EL2eDHy10HFlw8y6AUe4+yQAd9/k7p9RRtcfaAd0iWoXlYSZF0ZTotff3Z8HalM2p17vE5K2T4mOewnoYWb9CxFnJunid/en3D26sZcXCf+HAY4Hpka/V/MJiaVo49IyXHuAXwD/nbItp2tfrskj3QDCnYoUS2xmNgzYj/DL19/dl0FIMEApP3cw8YvnAGbWB6hN+s+0GBhUpNi2ZVfgUzObFDW7/cbMdqBMrr+7fwTcBCwElgCfAa8BK8vk+if0S7ne/aLtqf+nl1D6/6fPBP4avS75+M1sHLDI3Wel7Mop9nJNHmU7gNDMugK/B86PaiDlEvd/Asui2lPi+htb/1uU6udpD4wC/sfdRwFrCE0opRpvE2bWk/AX4lBCguhCaB5JVRafJ42y+j9tZpcBG939gcSmNMVKJn4zqwQuA9LNnJZT7OWaPBYDOyetDyZU4Uta1Nzwe+C37j4t2rwsUUU0swHAx8WKbxsOA443s3nAA4TmqlsIVdzE71Ep/zssJvzV9Uq0/gdCMimX6380MM/dV7j7ZuBR4FCgZ5lc/4RM13sxMCSpXMl+FjP7BqH5dmLS5lKP/3PAMOANM/uQEN9rZtaPHGMv1+TxMrCbmQ01s47ABOCxIseUjXuAd9z9l0nbHgPOiF5/A5iWelApcPdL3X1nd9+VcL2fcfevAzOAk6NipRz/MmCRmQ2PNo0B3qZMrj+huepgM+scdcQm4i/1659aO02+3mfQGO9jwH/BltklViaat4qsSfwWHhPxE+B4d1+fVO4xYEJ0B9wuwG6Egc3FtCV2d3/L3Qe4+67uvgshYXzB3T8m12vv7mW5AGMJI9DnAhcXO54s4j0M2AzMBF4ntFePBXoDT0Wf5UmgZ7FjzeKzfAl4LHq9C/ASYfr8B4EOxY6vmbhHEv7wmAk8AvQop+tPaHKYDbxJ6GzuUMrXH7if8BfsekLy+ybQK9P1JtxB+T7wBuGuslKMfy7hTqXXouX2pPKXRPHPBo4ptdhT9s8Derfk2muQoIiIxFauzVYiIlJESh4iIhKbkoeIiMSm5CEiIrEpeYiISGxKHiIiEpuSh5QFM/tGNI10c8tmMzu+2LGKbA9KZu5/kSzVREsm7xYmDJHtm5KHlJsad7+62EGIbO/UbCUiIrEpeUibY2bVUR/IUWb2kpnVR48G3SHa39/MbjezRdG+eWZ2XTRdfuq5hpjZlOiRnqvM7JFo23wzeyap3L3Re+6b5hwNZvZamu3fMbNXzWytma0ws2mpjwA1sy9Fx/+XmZ1p4RG066LYfx5NtZ163hPMbIaZ1ZrZJ2b2pJkdEe2rjD7H/AzX7p7o/T6XxaWW7ZiSh7RFiQnbfkd4bsethOautWa2M/AK8O3o582EfpKfADXJX8ZmNoQwM+rXgH8CdwB7EWay3SHNe2Y9UZyZTQFuJ0xueAfwEHAE8E8zq0pzyHlR+VmEx+muA34cvU4+7yWE6dr3BB4mTJ//BeBpMzva3dcRHgswxMwOTzm2E3Ai8E93/yDbzyLbJ/V5SLkZ3cyzoSe5+8Kk9fnuflRKmTuAgcA4d388sdHMvk9IMlcRHhIFcAPhSXdnufu9UbkrCclj11w/gJmdDHwduA84w6MnAZrZz4BXgSlmtqu7b0o6bCRwuLv/Kyr7/wgzvH7NzM5393VmtjvwU8Iz5ke7+ydR2V8SZhK+kfAEyymE6dBPA55Peo/jge7Ab3P9bLIdKfa0x1q0ZLMQnlXRQJjWPtNyZFT2qmj9opRzDIi2P5bm/EaYantZtN6ZUGt5K03ZQ6JYnknaNik6975pyjcAryWtT4/K9klT9ppo33HR+pei459MU/bRqOzwaP2yaP3UNGXPBs4H2kXrHxIexNQuqcw0oB7oVex/by2lv6jmIeWm2t2vybLs/JT1UYQk0cfMUh/HacAGYLCZDSQ856OS8PyPVC8Smo1yNYrwJf39NJWoPaNY9gMeT9r+XprzfBb97BT9TPS3vJha0N3vTNl0H3ApcAzwuJn1Ao4F/uLutdl9DNmeKXlIucnUZJVO6hd8z+jnwdGSjhMSR49ovW6rAu5uZi35gu0JtAOubCaGXinb1mcoB43XJHHMqiximEKoqZxGSFKnEvpf1GQlWVHykO1JIhFc4+7VzRU0sxHRy4EZivRJWU/9Ik+cZ6u7oaI4Vrn7sOZiyEHi83UDmiQ3M+sMrHd3B3D3uWb2InCCmXUATomO+UueY5I2SndbyfbkzejnAel2mtlPzewiM2tPaCZaDRxgZu1Syn2exqaihA3Rz9TbfdPd8vomoXmsX5oYvmJm10TvEdes6OeBafbdBqw1s6FJ26ZE8Y4nPCb5QXffmMP7ynZIyUO2G+4+H3gOOM7MTkreZ2anA1cAx7r7JnffTOgE35lwG2+iXAfg+jSnf5dQ6/hKUlkDLk9T9l7C/71fRedLlB9IuBvsYtI0l2XhfkIN6DIz65103s8BJwPvu/uCpPJTgY2Eu8raE/pBRLKiZitpqzL1jXybkEAeNrPHgbeAPQhf+p8C300qewVwFHCtmR1F+Mt+DLATW4/puJ9wp9R/m9luwDxCZ3RPIPn2Ydz93mgCx5OAfc3sCcL/xVMI/S0XufuHcT+wu88xs2rC7bpvmNmfouswgVBT+mZK+ZVRmZOAD9z9hbjvKdsv1TyknMQZiJe2nLu/B+wP/Ab4PPADwl1Kk4ED3X1OUtnVwOHAr4C9CYlnETCWlOTk7h8Tbqt9Otr/LeBtQnNQbZp4TiLcOrsGOIuQON4Cxrv7jbl+bne/ltD5vRA4nTDA8UXCbcyvpDnkoeinOsolFov6z0QkS2bWg5AQanzrQYhlJRpseBGwu7vPK3Y8Uj5U8xDZTpnZYOBMwmBHJQ6JRclDZDtjZhPN7HVCJ39foLq4EUk5UvIQyU2siRBLzGJgMKHp7dvu/o8ixyNlSH0eIiISm2oeIiISm5KHiIjEpuQhIiKxKXmIiEhsSh4iIhKbkoeIiMT2/wHeE4X8+U42/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x17a532fd400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r,c = filt_sig.shape\n",
    "\n",
    "# compute FFT - all trials at once\n",
    "sp = np.fft.rfft(filt_sig, axis=1)                    \n",
    "\n",
    "# compute amplitude\n",
    "amp = (np.abs(sp) / c) * 2\n",
    "\n",
    "mean_amp = np.mean(amp, axis=0)\n",
    "\n",
    "# get freq axis\n",
    "freq = np.fft.rfftfreq(c, 1/sr)   # Sample spacing is 2nd arg\n",
    "\n",
    "plt.plot(freq, mean_amp, 'b', linewidth=1)          # plot amplitude (power = np.absolute(sp)**2)\n",
    "\n",
    "plt.legend(['normalized fft'], loc='best', fancybox=1)  \n",
    "plt.xlabel('Frequency', **fig_font)\n",
    "plt.ylabel('Amplitude', **fig_font)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}