Median Filter for Noise Reduction in Atmospheric Readings.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Median Filter For Noise Reduction in Atmospheric Readings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline Ideas\n",
    "* Problem \n",
    "    * Separate noise from signal \n",
    "    * Possible lack of communication \n",
    "* Possible Solution \n",
    "* Evaluation of solution "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Statement\n",
    "Earth's atmosphere is an incredibly dynamic and hard-to-predict environment $-$ temperatures, wind speeds, and UV-radiation acting on an object will vary wildly as it ascends to ever-greater heights. We face this challenging environment with the goal of accurately sensing the altitude of a small payload tethered to a high-altitude weather balloon. As the payload must execute certain mission-critical operations at a specific altitude, an investigation into noise reduction in our atmospheric sensing systems has been deemed necessary. In addition to the unavoidable noise and variation in atmospheric conditions, there is also the possibility of the loss of communication between the ground station and payload, further necessitating reliance on the payload's automatic decision-making in the potential absence of manually-sent commands. Our primary objective, therefore, is to develop a method to effectively separate noise from the signal and accurately assess the altitude of the payload."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Possible Solutions\n",
    "A number of potential solutions have been proposed: \n",
    "* Bandpass filtering \n",
    "* Arithmetic Mean filtering\n",
    "* Median filtering\n",
    "Here, each of these solutions will be detailed\n",
    "\n",
    "### Bandpass Filtering\n",
    "\n",
    "### Arithmetic Mean Filtering\n",
    "\n",
    "### Median filtering\n",
    "Median filtering is the process of replacing each entry in a signal with the median of itself and surrounding data points. We will define a median filter to be of algorithmic complexity $\\mathcal{O}(n)$ if it calculates the median for a given data point with the $\\it{n}$ nearest data points. In contrast to both bandpass and arithmetic mean filtering, a median filter has the advantage of effectively filtering even the most extreme outliers from the data set provided that there are a sufficient number of valid data points. \n",
    "\n",
    "Simplicity is another advantage of median filter over more sophisticated approaches. A median filter can be implemented by maintaining a list of last $\\it{n}$ measured values, sorting the list of values, and selecting the midpoint of the newly-sorted list. Due to the nature of a sorted list, outlying values will be pushed to the extreme ends of either side of the list and will therefore be ignored. By selecting the median of this list, if the range of surrounding values is of sufficient size, the signal will be effectively separated from the noise inherent to the system. \n",
    "\n",
    "However, median filtering has some weaknesses. For our purposes, the method is still limited by the constant change in altitude and pressure. If a median filter is employed over too large a range of values, there will be an unnecessary loss of resolution in the signal. The median strategy is also somewhat dependent on noise being evenly distributed on both sides of the true signal. If the system produces noise that is disproportionally high or low, the median will be likewise skewed, and further mathematical adjustments may be required. Additionally, although the median filter is effective at filtering out the occasional outlier, if too many extreme values are received $-$ especially sequentially $-$ it will render the filter ineffective at reducing the data to a useful signal.\n",
    "\n",
    "## // I hope, further analysis of the data is needed to confirm this hypothesis - Jack \n",
    "## // It seems like the acutal most likely scenario is just a hiccup in recording. At least from the two available data sets we have --Nick\n",
    "\n",
    "The below plots show an example of applying a median filter to a subset of barometric pressure data. Figure 2 is a particularly noisy set of data that has two abnormalities that could potentially be experienced during a flight: inexplicably reading zero for a duration of time and a hiccup in responses from the sensor.\n",
    "\n",
    "Figure 2 contains 76 readings $-$ roughly amounting to a minute's worth of data. In this test-data set, we have introduced two groups of sequential readings of zero in our set $-$ for a total of fifteen zero-readings $-$ as well as a gap of ten readings to simulate a lack of response from the sensor for roughly 10-12 seconds.\n",
    "\n",
    "As shown in Figure 3, the filter effectively removes outlying data points when applied within a sufficiently large order, despite a nearly 20% loss in data. It does not, however, interpolate between data points very well. \n",
    "\n",
    "# //Should we maybe apply an arithmetic mean filter after our median filter? --Nick? \n",
    "\n",
    "### Failure of the Median Filter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named 'matplotlib'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-a2405b2a2501>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmagic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcsv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mmagic\u001b[0;34m(self, arg_s)\u001b[0m\n\u001b[1;32m   2156\u001b[0m         \u001b[0mmagic_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marg_s\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpartition\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m' '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2157\u001b[0m         \u001b[0mmagic_name\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmagic_name\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprefilter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mESC_MAGIC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2158\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2159\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2160\u001b[0m     \u001b[0;31m#-------------------------------------------------------------------------\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[0;34m(self, magic_name, line)\u001b[0m\n\u001b[1;32m   2077\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'local_ns'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getframe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_locals\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2078\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2079\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2080\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2081\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-104>\u001b[0m in \u001b[0;36mmatplotlib\u001b[0;34m(self, line)\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    187\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 188\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/magics/pylab.py\u001b[0m in \u001b[0;36mmatplotlib\u001b[0;34m(self, line)\u001b[0m\n\u001b[1;32m     98\u001b[0m             \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Available matplotlib backends: %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mbackends_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     99\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 100\u001b[0;31m             \u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshell\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menable_matplotlib\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    101\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_show_matplotlib_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    102\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36menable_matplotlib\u001b[0;34m(self, gui)\u001b[0m\n\u001b[1;32m   2933\u001b[0m         \"\"\"\n\u001b[1;32m   2934\u001b[0m         \u001b[0;32mfrom\u001b[0m \u001b[0mIPython\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcore\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpylabtools\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mpt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2935\u001b[0;31m         \u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfind_gui_and_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgui\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpylab_gui_select\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2936\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2937\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mgui\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36mfind_gui_and_backend\u001b[0;34m(gui, gui_select)\u001b[0m\n\u001b[1;32m    257\u001b[0m     \"\"\"\n\u001b[1;32m    258\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 259\u001b[0;31m     \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    260\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    261\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mgui\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mgui\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m'auto'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: No module named 'matplotlib'"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "import time\n",
    "\n",
    "\n",
    "# Import entire data set from ERRNO_3 flightE\n",
    "csvfile = open(\"baro.csv\")\n",
    "reader = csv.reader(csvfile)\n",
    "t = []\n",
    "alt = []\n",
    "for row in reader: \n",
    "    t.append(row[0])\n",
    "    alt.append(row[1])\n",
    "\n",
    "\n",
    "# Subset of barometer data, I deleted sections and added zeros to show the kind of noise we might expect \n",
    "# The most likely kind of noise is no reading, and that has historically happened when it gets really cold. \n",
    "\n",
    "subset = open(\"subset.csv\")\n",
    "subset_reader = csv.reader(subset)\n",
    "pressure = []\n",
    "for row in subset_reader: \n",
    "    pressure.append(row[0])\n",
    "\n",
    "def medianFilter(data, n):\n",
    "    \"\"\"Applies a simple median filter to a list of data\"\"\"\n",
    "    filtered = []\n",
    "    for i,p in enumerate(pressure):\n",
    "        local = data[i:i+n]\n",
    "        local.sort()\n",
    "        mid = int(len(local) / 2)\n",
    "        filtered.append(local[mid])\n",
    "    return filtered\n",
    "\n",
    "\n",
    "\n",
    "## Plot entire data set, Figure 1\n",
    "plt.plot(alt)\n",
    "plt.ylabel('Pressure')\n",
    "plt.xlabel('Time')\n",
    "plt.title('Entire Data Set (Figure 1)')\n",
    "plt.show()\n",
    "\n",
    "## Plot noisy subset of pressure data, Figure 2 \n",
    "plt.plot(pressure)\n",
    "plt.ylabel('Pressure')\n",
    "plt.xlabel('Time')\n",
    "plt.title('Noisy Subset (Figure 2)')\n",
    "plt.axis([0, 80, 0, 200])\n",
    "plt.show()\n",
    "\n",
    "## Apply median filter to noisy section of data, Figure 3\n",
    "plt.plot(medianFilter(pressure, 25))\n",
    "plt.ylabel('Pressure')\n",
    "plt.xlabel('Time')\n",
    "plt.title('Median-filtered Subset (Figure 3)')\n",
    "plt.axis([0, 80, 0, 200])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}