peak finding comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "from math import pi, log\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import fft, ifft\n",
    "from scipy.optimize import curve_fit\n",
    " \n",
    "def _datacheck_peakdetect(x_axis, y_axis):\n",
    "    if x_axis is None:\n",
    "        x_axis = range(len(y_axis))\n",
    "    \n",
    "    if len(y_axis) != len(x_axis):\n",
    "        raise (ValueError, \n",
    "                'Input vectors y_axis and x_axis must have same length')\n",
    "    \n",
    "    #needs to be a numpy array\n",
    "    y_axis = np.array(y_axis)\n",
    "    x_axis = np.array(x_axis)\n",
    "    return x_axis, y_axis\n",
    "    \n",
    "def _peakdetect_parabole_fitter(raw_peaks, x_axis, y_axis, points):\n",
    "    \"\"\"\n",
    "    Performs the actual parabole fitting for the peakdetect_parabole function.\n",
    "    \n",
    "    keyword arguments:\n",
    "    raw_peaks -- A list of either the maximium or the minimum peaks, as given\n",
    "        by the peakdetect_zero_crossing function, with index used as x-axis\n",
    "    x_axis -- A numpy list of all the x values\n",
    "    y_axis -- A numpy list of all the y values\n",
    "    points -- How many points around the peak should be used during curve\n",
    "        fitting, must be odd.\n",
    "    \n",
    "    return -- A list giving all the peaks and the fitted waveform, format:\n",
    "        [[x, y, [fitted_x, fitted_y]]]\n",
    "        \n",
    "    \"\"\"\n",
    "    func = lambda x, k, tau, m: k * ((x - tau) ** 2) + m\n",
    "    fitted_peaks = []\n",
    "    for peak in raw_peaks:\n",
    "        index = peak[0]\n",
    "        x_data = x_axis[index - points // 2: index + points // 2 + 1]\n",
    "        y_data = y_axis[index - points // 2: index + points // 2 + 1]\n",
    "        # get a first approximation of tau (peak position in time)\n",
    "        tau = x_axis[index]\n",
    "        # get a first approximation of peak amplitude\n",
    "        m = peak[1]\n",
    "        \n",
    "        # build list of approximations\n",
    "        # k = -m as first approximation?\n",
    "        p0 = (-m, tau, m)\n",
    "        popt, pcov = curve_fit(func, x_data, y_data, p0)\n",
    "        # retrieve tau and m i.e x and y value of peak\n",
    "        x, y = popt[1:3]\n",
    "        \n",
    "        # create a high resolution data set for the fitted waveform\n",
    "        x2 = np.linspace(x_data[0], x_data[-1], points * 10)\n",
    "        y2 = func(x2, *popt)\n",
    "        \n",
    "        fitted_peaks.append([x, y, [x2, y2]])\n",
    "        \n",
    "    return fitted_peaks\n",
    "    \n",
    "    \n",
    "def peakdetect(y_axis, x_axis = None, lookahead = 300, delta=0):\n",
    "    \"\"\"\n",
    "    Converted from/based on a MATLAB script at: \n",
    "    http://billauer.co.il/peakdet.html\n",
    "    \n",
    "    function for detecting local maximas and minmias in a signal.\n",
    "    Discovers peaks by searching for values which are surrounded by lower\n",
    "    or larger values for maximas and minimas respectively\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- (optional) A x-axis whose values correspond to the y_axis list\n",
    "        and is used in the return to specify the postion of the peaks. If\n",
    "        omitted an index of the y_axis is used. (default: None)\n",
    "    lookahead -- (optional) distance to look ahead from a peak candidate to\n",
    "        determine if it is the actual peak (default: 200) \n",
    "        '(sample / period) / f' where '4 >= f >= 1.25' might be a good value\n",
    "    delta -- (optional) this specifies a minimum difference between a peak and\n",
    "        the following points, before a peak may be considered a peak. Useful\n",
    "        to hinder the function from picking up false peaks towards to end of\n",
    "        the signal. To work well delta should be set to delta >= RMSnoise * 5.\n",
    "        (default: 0)\n",
    "            delta function causes a 20% decrease in speed, when omitted\n",
    "            Correctly used it can double the speed of the function\n",
    "    \n",
    "    return -- two lists [max_peaks, min_peaks] containing the positive and\n",
    "        negative peaks respectively. Each cell of the lists contains a tupple\n",
    "        of: (position, peak_value) \n",
    "        to get the average peak value do: np.mean(max_peaks, 0)[1] on the\n",
    "        results to unpack one of the lists into x, y coordinates do: \n",
    "        x, y = zip(*tab)\n",
    "    \"\"\"\n",
    "    max_peaks = []\n",
    "    min_peaks = []\n",
    "    dump = []   #Used to pop the first hit which almost always is false\n",
    "       \n",
    "    # check input data\n",
    "    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)\n",
    "    # store data length for later use\n",
    "    length = len(y_axis)\n",
    "    \n",
    "    \n",
    "    #perform some checks\n",
    "    if lookahead < 1:\n",
    "        raise ValueError, \"Lookahead must be '1' or above in value\"\n",
    "    if not (np.isscalar(delta) and delta >= 0):\n",
    "        raise ValueError, \"delta must be a positive number\"\n",
    "    \n",
    "    #maxima and minima candidates are temporarily stored in\n",
    "    #mx and mn respectively\n",
    "    mn, mx = np.Inf, -np.Inf\n",
    "    \n",
    "    #Only detect peak if there is 'lookahead' amount of points after it\n",
    "    for index, (x, y) in enumerate(zip(x_axis[:-lookahead], \n",
    "                                        y_axis[:-lookahead])):\n",
    "        if y > mx:\n",
    "            mx = y\n",
    "            mxpos = x\n",
    "        if y < mn:\n",
    "            mn = y\n",
    "            mnpos = x\n",
    "        \n",
    "        ####look for max####\n",
    "        if y < mx-delta and mx != np.Inf:\n",
    "            #Maxima peak candidate found\n",
    "            #look ahead in signal to ensure that this is a peak and not jitter\n",
    "            if y_axis[index:index+lookahead].max() < mx:\n",
    "                max_peaks.append([mxpos, mx])\n",
    "                dump.append(True)\n",
    "                #set algorithm to only find minima now\n",
    "                mx = np.Inf\n",
    "                mn = np.Inf\n",
    "                if index+lookahead >= length:\n",
    "                    #end is within lookahead no more peaks can be found\n",
    "                    break\n",
    "                continue\n",
    "            #else:  #slows shit down this does\n",
    "            #    mx = ahead\n",
    "            #    mxpos = x_axis[np.where(y_axis[index:index+lookahead]==mx)]\n",
    "        \n",
    "        ####look for min####\n",
    "        if y > mn+delta and mn != -np.Inf:\n",
    "            #Minima peak candidate found \n",
    "            #look ahead in signal to ensure that this is a peak and not jitter\n",
    "            if y_axis[index:index+lookahead].min() > mn:\n",
    "                min_peaks.append([mnpos, mn])\n",
    "                dump.append(False)\n",
    "                #set algorithm to only find maxima now\n",
    "                mn = -np.Inf\n",
    "                mx = -np.Inf\n",
    "                if index+lookahead >= length:\n",
    "                    #end is within lookahead no more peaks can be found\n",
    "                    break\n",
    "            #else:  #slows shit down this does\n",
    "            #    mn = ahead\n",
    "            #    mnpos = x_axis[np.where(y_axis[index:index+lookahead]==mn)]\n",
    "    \n",
    "    \n",
    "    #Remove the false hit on the first value of the y_axis\n",
    "    try:\n",
    "        if dump[0]:\n",
    "            max_peaks.pop(0)\n",
    "        else:\n",
    "            min_peaks.pop(0)\n",
    "        del dump\n",
    "    except IndexError:\n",
    "        #no peaks were found, should the function return empty lists?\n",
    "        pass\n",
    "        \n",
    "    return [max_peaks, min_peaks]\n",
    "    \n",
    "    \n",
    "def peakdetect_fft(y_axis, x_axis, pad_len = 5):\n",
    "    \"\"\"\n",
    "    Performs a FFT calculation on the data and zero-pads the results to\n",
    "    increase the time domain resolution after performing the inverse fft and\n",
    "    send the data to the 'peakdetect' function for peak \n",
    "    detection.\n",
    "    \n",
    "    Omitting the x_axis is forbidden as it would make the resulting x_axis\n",
    "    value silly if it was returned as the index 50.234 or similar.\n",
    "    \n",
    "    Will find at least 1 less peak then the 'peakdetect_zero_crossing'\n",
    "    function, but should result in a more precise value of the peak as\n",
    "    resolution has been increased. Some peaks are lost in an attempt to\n",
    "    minimize spectral leakage by calculating the fft between two zero\n",
    "    crossings for n amount of signal periods.\n",
    "    \n",
    "    The biggest time eater in this function is the ifft and thereafter it's\n",
    "    the 'peakdetect' function which takes only half the time of the ifft.\n",
    "    Speed improvementd could include to check if 2**n points could be used for\n",
    "    fft and ifft or change the 'peakdetect' to the 'peakdetect_zero_crossing',\n",
    "    which is maybe 10 times faster than 'peakdetct'. The pro of 'peakdetect'\n",
    "    is that it resutls in one less lost peak. It should also be noted that the\n",
    "    time used by the ifft function can change greatly depending on the input.\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- A x-axis whose values correspond to the y_axis list and is used\n",
    "        in the return to specify the postion of the peaks.\n",
    "    pad_len -- (optional) By how many times the time resolution should be\n",
    "        increased by, e.g. 1 doubles the resolution. The amount is rounded up\n",
    "        to the nearest 2 ** n amount (default: 5)\n",
    "    \n",
    "    return -- two lists [max_peaks, min_peaks] containing the positive and\n",
    "        negative peaks respectively. Each cell of the lists contains a tupple\n",
    "        of: (position, peak_value) \n",
    "        to get the average peak value do: np.mean(max_peaks, 0)[1] on the\n",
    "        results to unpack one of the lists into x, y coordinates do: \n",
    "        x, y = zip(*tab)\n",
    "    \"\"\"\n",
    "    # check input data\n",
    "    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)\n",
    "    zero_indices = zero_crossings(y_axis, window = 11)\n",
    "    #select a n amount of periods\n",
    "    last_indice = - 1 - (1 - len(zero_indices) & 1)\n",
    "    # Calculate the fft between the first and last zero crossing\n",
    "    # this method could be ignored if the begining and the end of the signal\n",
    "    # are discardable as any errors induced from not using whole periods\n",
    "    # should mainly manifest in the beginning and the end of the signal, but\n",
    "    # not in the rest of the signal\n",
    "    fft_data = fft(y_axis[zero_indices[0]:zero_indices[last_indice]])\n",
    "    padd = lambda x, c: x[:len(x) // 2] + [0] * c + x[len(x) // 2:]\n",
    "    n = lambda x: int(log(x)/log(2)) + 1\n",
    "    # padds to 2**n amount of samples\n",
    "    fft_padded = padd(list(fft_data), 2 ** \n",
    "                n(len(fft_data) * pad_len) - len(fft_data))\n",
    "    \n",
    "    # There is amplitude decrease directly proportional to the sample increase\n",
    "    sf = len(fft_padded) / float(len(fft_data))\n",
    "    # There might be a leakage giving the result an imaginary component\n",
    "    # Return only the real component\n",
    "    y_axis_ifft = ifft(fft_padded).real * sf #(pad_len + 1)\n",
    "    x_axis_ifft = np.linspace(\n",
    "                x_axis[zero_indices[0]], x_axis[zero_indices[last_indice]],\n",
    "                len(y_axis_ifft))\n",
    "    # get the peaks to the interpolated waveform\n",
    "    max_peaks, min_peaks = peakdetect(y_axis_ifft, x_axis_ifft, 500,\n",
    "                                    delta = abs(np.diff(y_axis).max() * 2))\n",
    "    #max_peaks, min_peaks = peakdetect_zero_crossing(y_axis_ifft, x_axis_ifft)\n",
    "    \n",
    "    # store one 20th of a period as waveform data\n",
    "    data_len = int(np.diff(zero_indices).mean()) / 10\n",
    "    data_len += 1 - data_len & 1\n",
    "    \n",
    "    \n",
    "    fitted_wave = []\n",
    "    for peaks in [max_peaks, min_peaks]:\n",
    "        peak_fit_tmp = []\n",
    "        index = 0\n",
    "        for peak in peaks:\n",
    "            index = np.where(x_axis_ifft[index:]==peak[0])[0][0] + index\n",
    "            x_fit_lim = x_axis_ifft[index - data_len // 2:\n",
    "                                    index + data_len // 2 + 1]\n",
    "            y_fit_lim = y_axis_ifft[index - data_len // 2:\n",
    "                                    index + data_len // 2 + 1]\n",
    "            \n",
    "            peak_fit_tmp.append([x_fit_lim, y_fit_lim])\n",
    "        fitted_wave.append(peak_fit_tmp)\n",
    "    \n",
    "    #pylab.plot(range(len(fft_data)), fft_data)\n",
    "    #pylab.show()\n",
    "    \n",
    "    pylab.plot(x_axis, y_axis)\n",
    "    pylab.hold(True)\n",
    "    pylab.plot(x_axis_ifft, y_axis_ifft)\n",
    "    #for max_p in max_peaks:\n",
    "    #    pylab.plot(max_p[0], max_p[1], 'xr')\n",
    "    pylab.show()\n",
    "    return [max_peaks, min_peaks]\n",
    "    \n",
    "    \n",
    "def peakdetect_parabole(y_axis, x_axis, points = 9):\n",
    "    \"\"\"\n",
    "    Function for detecting local maximas and minmias in a signal.\n",
    "    Discovers peaks by fitting the model function: y = k (x - tau) ** 2 + m\n",
    "    to the peaks. The amount of points used in the fitting is set by the\n",
    "    points argument.\n",
    "    \n",
    "    Omitting the x_axis is forbidden as it would make the resulting x_axis\n",
    "    value silly if it was returned as index 50.234 or similar.\n",
    "    \n",
    "    will find the same amount of peaks as the 'peakdetect_zero_crossing'\n",
    "    function, but might result in a more precise value of the peak.\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- A x-axis whose values correspond to the y_axis list and is used\n",
    "        in the return to specify the postion of the peaks.\n",
    "    points -- (optional) How many points around the peak should be used during\n",
    "        curve fitting, must be odd (default: 9)\n",
    "    \n",
    "    return -- two lists [max_peaks, min_peaks] containing the positive and\n",
    "        negative peaks respectively. Each cell of the lists contains a list\n",
    "        of: (position, peak_value) \n",
    "        to get the average peak value do: np.mean(max_peaks, 0)[1] on the\n",
    "        results to unpack one of the lists into x, y coordinates do: \n",
    "        x, y = zip(*max_peaks)\n",
    "    \"\"\"\n",
    "    # check input data\n",
    "    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)\n",
    "    # make the points argument odd\n",
    "    points += 1 - points % 2\n",
    "    #points += 1 - int(points) & 1 slower when int conversion needed\n",
    "    \n",
    "    # get raw peaks\n",
    "    max_raw, min_raw = peakdetect_zero_crossing(y_axis)\n",
    "    \n",
    "    # define output variable\n",
    "    max_peaks = []\n",
    "    min_peaks = []\n",
    "    \n",
    "    max_ = _peakdetect_parabole_fitter(max_raw, x_axis, y_axis, points)\n",
    "    min_ = _peakdetect_parabole_fitter(min_raw, x_axis, y_axis, points)\n",
    "    \n",
    "    max_peaks = map(lambda x: [x[0], x[1]], max_)\n",
    "    max_fitted = map(lambda x: x[-1], max_)\n",
    "    min_peaks = map(lambda x: [x[0], x[1]], min_)\n",
    "    min_fitted = map(lambda x: x[-1], min_)\n",
    "    \n",
    "    \n",
    "    #pylab.plot(x_axis, y_axis)\n",
    "    #pylab.hold(True)\n",
    "    #for max_p, max_f in zip(max_peaks, max_fitted):\n",
    "    #    pylab.plot(max_p[0], max_p[1], 'x')\n",
    "    #    pylab.plot(max_f[0], max_f[1], 'o', markersize = 2)\n",
    "    #for min_p, min_f in zip(min_peaks, min_fitted):\n",
    "    #    pylab.plot(min_p[0], min_p[1], 'x')\n",
    "    #    pylab.plot(min_f[0], min_f[1], 'o', markersize = 2)\n",
    "    #pylab.show()\n",
    "    \n",
    "    return [max_peaks, min_peaks]\n",
    "    \n",
    " \n",
    "def peakdetect_sine(y_axis, x_axis, points = 9, lock_frequency = False):\n",
    "    \"\"\"\n",
    "    Function for detecting local maximas and minmias in a signal.\n",
    "    Discovers peaks by fitting the model function:\n",
    "    y = A * sin(2 * pi * f * x - tau) to the peaks. The amount of points used\n",
    "    in the fitting is set by the points argument.\n",
    "    \n",
    "    Omitting the x_axis is forbidden as it would make the resulting x_axis\n",
    "    value silly if it was returned as index 50.234 or similar.\n",
    "    \n",
    "    will find the same amount of peaks as the 'peakdetect_zero_crossing'\n",
    "    function, but might result in a more precise value of the peak.\n",
    "    \n",
    "    The function might have some problems if the sine wave has a\n",
    "    non-negligible total angle i.e. a k*x component, as this messes with the\n",
    "    internal offset calculation of the peaks, might be fixed by fitting a \n",
    "    k * x + m function to the peaks for offset calculation.\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- A x-axis whose values correspond to the y_axis list and is used\n",
    "        in the return to specify the postion of the peaks.\n",
    "    points -- (optional) How many points around the peak should be used during\n",
    "        curve fitting, must be odd (default: 9)\n",
    "    lock_frequency -- (optional) Specifies if the frequency argument of the\n",
    "        model function should be locked to the value calculated from the raw\n",
    "        peaks or if optimization process may tinker with it. (default: False)\n",
    "    \n",
    "    return -- two lists [max_peaks, min_peaks] containing the positive and\n",
    "        negative peaks respectively. Each cell of the lists contains a tupple\n",
    "        of: (position, peak_value) \n",
    "        to get the average peak value do: np.mean(max_peaks, 0)[1] on the\n",
    "        results to unpack one of the lists into x, y coordinates do: \n",
    "        x, y = zip(*tab)\n",
    "    \"\"\"\n",
    "    # check input data\n",
    "    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)\n",
    "    # make the points argument odd\n",
    "    points += 1 - points % 2\n",
    "    #points += 1 - int(points) & 1 slower when int conversion needed\n",
    "    \n",
    "    # get raw peaks\n",
    "    max_raw, min_raw = peakdetect_zero_crossing(y_axis)\n",
    "    \n",
    "    # define output variable\n",
    "    max_peaks = []\n",
    "    min_peaks = []\n",
    "    \n",
    "    # get global offset\n",
    "    offset = np.mean([np.mean(max_raw, 0)[1], np.mean(min_raw, 0)[1]])\n",
    "    # fitting a k * x + m function to the peaks might be better\n",
    "    #offset_func = lambda x, k, m: k * x + m\n",
    "    \n",
    "    # calculate an approximate frequenzy of the signal\n",
    "    Hz = []\n",
    "    for raw in [max_raw, min_raw]:\n",
    "        if len(raw) > 1:\n",
    "            peak_pos = [x_axis[index] for index in zip(*raw)[0]]\n",
    "            Hz.append(np.mean(np.diff(peak_pos)))\n",
    "    Hz = 1 / np.mean(Hz)\n",
    "    \n",
    "    # model function\n",
    "    # if cosine is used then tau could equal the x position of the peak\n",
    "    # if sine were to be used then tau would be the first zero crossing\n",
    "    if lock_frequency:\n",
    "        func = lambda x, A, tau: A * np.sin(2 * pi * Hz * (x - tau) + pi / 2)\n",
    "    else:\n",
    "        func = lambda x, A, Hz, tau: A * np.sin(2 * pi * Hz * (x - tau) + \n",
    "                                                pi / 2)\n",
    "    #func = lambda x, A, Hz, tau: A * np.cos(2 * pi * Hz * (x - tau))\n",
    "    \n",
    "    \n",
    "    #get peaks\n",
    "    fitted_peaks = []\n",
    "    for raw_peaks in [max_raw, min_raw]:\n",
    "        peak_data = []\n",
    "        for peak in raw_peaks:\n",
    "            index = peak[0]\n",
    "            x_data = x_axis[index - points // 2: index + points // 2 + 1]\n",
    "            y_data = y_axis[index - points // 2: index + points // 2 + 1]\n",
    "            # get a first approximation of tau (peak position in time)\n",
    "            tau = x_axis[index]\n",
    "            # get a first approximation of peak amplitude\n",
    "            A = peak[1]\n",
    "            \n",
    "            # build list of approximations\n",
    "            if lock_frequency:\n",
    "                p0 = (A, tau)\n",
    "            else:\n",
    "                p0 = (A, Hz, tau)\n",
    "            \n",
    "            # subtract offset from waveshape\n",
    "            y_data -= offset\n",
    "            popt, pcov = curve_fit(func, x_data, y_data, p0)\n",
    "            # retrieve tau and A i.e x and y value of peak\n",
    "            x = popt[-1]\n",
    "            y = popt[0]\n",
    "            \n",
    "            # create a high resolution data set for the fitted waveform\n",
    "            x2 = np.linspace(x_data[0], x_data[-1], points * 10)\n",
    "            y2 = func(x2, *popt)\n",
    "            \n",
    "            # add the offset to the results\n",
    "            y += offset\n",
    "            y2 += offset\n",
    "            y_data += offset\n",
    "            \n",
    "            peak_data.append([x, y, [x2, y2]])\n",
    "       \n",
    "        fitted_peaks.append(peak_data)\n",
    "    \n",
    "    # structure date for output\n",
    "    max_peaks = map(lambda x: [x[0], x[1]], fitted_peaks[0])\n",
    "    max_fitted = map(lambda x: x[-1], fitted_peaks[0])\n",
    "    min_peaks = map(lambda x: [x[0], x[1]], fitted_peaks[1])\n",
    "    min_fitted = map(lambda x: x[-1], fitted_peaks[1])\n",
    "    \n",
    "    \n",
    "    #pylab.plot(x_axis, y_axis)\n",
    "    #pylab.hold(True)\n",
    "    #for max_p, max_f in zip(max_peaks, max_fitted):\n",
    "    #    pylab.plot(max_p[0], max_p[1], 'x')\n",
    "    #    pylab.plot(max_f[0], max_f[1], 'o', markersize = 2)\n",
    "    #for min_p, min_f in zip(min_peaks, min_fitted):\n",
    "    #    pylab.plot(min_p[0], min_p[1], 'x')\n",
    "    #    pylab.plot(min_f[0], min_f[1], 'o', markersize = 2)\n",
    "    #pylab.show()\n",
    "    \n",
    "    return [max_peaks, min_peaks]\n",
    " \n",
    "    \n",
    "def peakdetect_sine_locked(y_axis, x_axis, points = 9):\n",
    "    \"\"\"\n",
    "    Convinience function for calling the 'peakdetect_sine' function with\n",
    "    the lock_frequency argument as True.\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- A x-axis whose values correspond to the y_axis list and is used\n",
    "        in the return to specify the postion of the peaks.\n",
    "    points -- (optional) How many points around the peak should be used during\n",
    "        curve fitting, must be odd (default: 9)\n",
    "        \n",
    "    return -- see 'peakdetect_sine'\n",
    "    \"\"\"\n",
    "    return peakdetect_sine(y_axis, x_axis, points, True)\n",
    "    \n",
    "    \n",
    "def peakdetect_zero_crossing(y_axis, x_axis = None, window = 11):\n",
    "    \"\"\"\n",
    "    Function for detecting local maximas and minmias in a signal.\n",
    "    Discovers peaks by dividing the signal into bins and retrieving the\n",
    "    maximum and minimum value of each the even and odd bins respectively.\n",
    "    Division into bins is performed by smoothing the curve and finding the\n",
    "    zero crossings.\n",
    "    \n",
    "    Suitable for repeatable signals, where some noise is tolerated. Excecutes\n",
    "    faster than 'peakdetect', although this function will break if the offset\n",
    "    of the signal is too large. It should also be noted that the first and\n",
    "    last peak will probably not be found, as this function only can find peaks\n",
    "    between the first and last zero crossing.\n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find peaks\n",
    "    x_axis -- (optional) A x-axis whose values correspond to the y_axis list\n",
    "        and is used in the return to specify the postion of the peaks. If\n",
    "        omitted an index of the y_axis is used. (default: None)\n",
    "    window -- the dimension of the smoothing window; should be an odd integer\n",
    "        (default: 11)\n",
    "    \n",
    "    return -- two lists [max_peaks, min_peaks] containing the positive and\n",
    "        negative peaks respectively. Each cell of the lists contains a tupple\n",
    "        of: (position, peak_value) \n",
    "        to get the average peak value do: np.mean(max_peaks, 0)[1] on the\n",
    "        results to unpack one of the lists into x, y coordinates do: \n",
    "        x, y = zip(*tab)\n",
    "    \"\"\"\n",
    "    # check input data\n",
    "    x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)\n",
    "    \n",
    "    zero_indices = zero_crossings(y_axis, window = window)\n",
    "    period_lengths = np.diff(zero_indices)\n",
    "            \n",
    "    bins_y = [y_axis[index:index + diff] for index, diff in \n",
    "        zip(zero_indices, period_lengths)]\n",
    "    bins_x = [x_axis[index:index + diff] for index, diff in \n",
    "        zip(zero_indices, period_lengths)]\n",
    "        \n",
    "    even_bins_y = bins_y[::2]\n",
    "    odd_bins_y = bins_y[1::2]\n",
    "    even_bins_x = bins_x[::2]\n",
    "    odd_bins_x = bins_x[1::2]\n",
    "    hi_peaks_x = []\n",
    "    lo_peaks_x = []\n",
    "    \n",
    "    #check if even bin contains maxima\n",
    "    if abs(even_bins_y[0].max()) > abs(even_bins_y[0].min()):\n",
    "        hi_peaks = [bin.max() for bin in even_bins_y]\n",
    "        lo_peaks = [bin.min() for bin in odd_bins_y]\n",
    "        # get x values for peak\n",
    "        for bin_x, bin_y, peak in zip(even_bins_x, even_bins_y, hi_peaks):\n",
    "            hi_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])\n",
    "        for bin_x, bin_y, peak in zip(odd_bins_x, odd_bins_y, lo_peaks):\n",
    "            lo_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])\n",
    "    else:\n",
    "        hi_peaks = [bin.max() for bin in odd_bins_y]\n",
    "        lo_peaks = [bin.min() for bin in even_bins_y]\n",
    "        # get x values for peak\n",
    "        for bin_x, bin_y, peak in zip(odd_bins_x, odd_bins_y, hi_peaks):\n",
    "            hi_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])\n",
    "        for bin_x, bin_y, peak in zip(even_bins_x, even_bins_y, lo_peaks):\n",
    "            lo_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])\n",
    "    \n",
    "    max_peaks = [[x, y] for x,y in zip(hi_peaks_x, hi_peaks)]\n",
    "    min_peaks = [[x, y] for x,y in zip(lo_peaks_x, lo_peaks)]\n",
    "    \n",
    "    return [max_peaks, min_peaks]\n",
    "        \n",
    "    \n",
    "def _smooth(x, window_len=11, window='hanning'):\n",
    "    \"\"\"\n",
    "    smooth the data using a window of the requested size.\n",
    "    \n",
    "    This method is based on the convolution of a scaled window on the signal.\n",
    "    The signal is prepared by introducing reflected copies of the signal \n",
    "    (with the window size) in both ends so that transient parts are minimized\n",
    "    in the begining and end part of the output signal.\n",
    "    \n",
    "    input:\n",
    "        x: the input signal \n",
    "        window_len: the dimension of the smoothing window; should be an odd \n",
    "            integer\n",
    "        window: the type of window from 'flat', 'hanning', 'hamming', \n",
    "            'bartlett', 'blackman'\n",
    "            flat window will produce a moving average smoothing.\n",
    " \n",
    "    output:\n",
    "        the smoothed signal\n",
    "        \n",
    "    example:\n",
    " \n",
    "    t = linspace(-2,2,0.1)\n",
    "    x = sin(t)+randn(len(t))*0.1\n",
    "    y = _smooth(x)\n",
    "    \n",
    "    see also: \n",
    "    \n",
    "    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, \n",
    "    numpy.convolve, scipy.signal.lfilter\n",
    " \n",
    "    TODO: the window parameter could be the window itself if a list instead of\n",
    "    a string   \n",
    "    \"\"\"\n",
    "    if x.ndim != 1:\n",
    "        raise ValueError, \"smooth only accepts 1 dimension arrays.\"\n",
    " \n",
    "    if x.size < window_len:\n",
    "        raise ValueError, \"Input vector needs to be bigger than window size.\"\n",
    "    \n",
    "    if window_len<3:\n",
    "        return x\n",
    "    \n",
    "    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:\n",
    "        raise(ValueError,\n",
    "            \"Window is not one of '{0}', '{1}', '{2}', '{3}', '{4}'\".format(\n",
    "            *('flat', 'hanning', 'hamming', 'bartlett', 'blackman')))\n",
    "    \n",
    "    s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]\n",
    "    #print(len(s))\n",
    "    if window == 'flat': #moving average\n",
    "        w = np.ones(window_len,'d')\n",
    "    else:\n",
    "        w = eval('np.' + window + '(window_len)')\n",
    " \n",
    "    y = np.convolve(w / w.sum(), s, mode = 'valid')\n",
    "    return y\n",
    "    \n",
    "    \n",
    "def zero_crossings(y_axis, window = 11):\n",
    "    \"\"\"\n",
    "    Algorithm to find zero crossings. Smoothens the curve and finds the\n",
    "    zero-crossings by looking for a sign change.\n",
    "    \n",
    "    \n",
    "    keyword arguments:\n",
    "    y_axis -- A list containg the signal over which to find zero-crossings\n",
    "    window -- the dimension of the smoothing window; should be an odd integer\n",
    "        (default: 11)\n",
    "    \n",
    "    return -- the index for each zero-crossing\n",
    "    \"\"\"\n",
    "    # smooth the curve\n",
    "    length = len(y_axis)\n",
    "    x_axis = np.asarray(range(length), int)\n",
    "    \n",
    "    # discard tail of smoothed signal\n",
    "    y_axis = _smooth(y_axis, window)[:length]\n",
    "    zero_crossings = np.where(np.diff(np.sign(y_axis)))[0]\n",
    "    indices = [x_axis[index] for index in zero_crossings]\n",
    "    \n",
    "    # check if zero-crossings are valid\n",
    "    diff = np.diff(indices)\n",
    "    if diff.std() / diff.mean() > 0.2:\n",
    "        print diff.std() / diff.mean()\n",
    "        print np.diff(indices)\n",
    "        raise(ValueError, \n",
    "            \"False zero-crossings found, indicates problem {0} or {1}\".format(\n",
    "            \"with smoothing window\", \"problem with offset\"))\n",
    "    # check if any zero crossings were found\n",
    "    if len(zero_crossings) < 1:\n",
    "        raise(ValueError, \"No zero crossings found\")\n",
    "    \n",
    "    return indices \n",
    "    # used this to test the fft function's sensitivity to spectral leakage\n",
    "    #return indices + np.asarray(30 * np.random.randn(len(indices)), int)\n",
    "    \n",
    "############################Frequency calculation#############################\n",
    "#    diff = np.diff(indices)\n",
    "#    time_p_period = diff.mean()\n",
    "#    \n",
    "#    if diff.std() / time_p_period > 0.1:\n",
    "#        raise ValueError, \n",
    "#            \"smoothing window too small, false zero-crossing found\"\n",
    "#    \n",
    "#    #return frequency\n",
    "#    return 1.0 / time_p_period\n",
    "##############################################################################\n",
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[315.0, 369.5, 378.5, 444.0, 459.0, 477.0, 485.5, 490.5]\n",
      "[ 295.   315.5  327.5  368.5  379.5  409.   428.   444.   458.   477.   486.\n",
      "  494.5]\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAA6YAAAFwCAYAAACrembwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcHWWZ9//PRSAbCWEPCgGUACIRAQXEDVyB+HNfcBtX\n",
       "GJeBMOPMz+1xBJ95nM3HBXDHdVxAHTdUCKNoEBFEQFZBWUQSUNYEEkhIAtfzR1WTTuf0fs6pOnU+\n",
       "79fLF+nq01VXjun+9nXXXfcdmYkkSZIkSVXZrOoCJEmSJEn9zcZUkiRJklQpG1NJkiRJUqVsTCVJ\n",
       "kiRJlbIxlSRJkiRVysZUkiRJklSpURvTiPhSRNweEVcN8/nXRcQVEXFlRFwQEfu1v0xJkjTAbJYk\n",
       "Nc1Y7ph+GThyhM/fBDwzM/cD/gX4fDsKkyRJwzKbJUmNMmpjmpnnA8tH+PyFmXlv+eFvgF3aVJsk\n",
       "SWrBbJYkNU27nzF9K3BWm88pSZImzmyWJNXe5u06UUQ8C3gL8LR2nVOSJE2c2SxJ6hVtaUzLRRVO\n",
       "A47MzJZTiyIi23EtSZIGZGZUXUNdmc2SpCpMNJsn3ZhGxK7A94DXZ+YNI73WXyA2iIiTMvOkquuo\n",
       "A9+Ljfl+bOB7sTHfj43ZVA3PbJ4Yv8c28L3YmO/HBr4XG/P92NhksnnUxjQiTgcOA7aPiKXAicAW\n",
       "AJn5OeCDwDbAZyICYF1mHjzRgiRJ0sjMZklS04zamGbma0b5/DHAMW2rSJIkjchsliQ1TbtX5dXY\n",
       "Lam6gBpZUnUBNbOk6gJqZEnVBdTMkqoLkBpuSdUF1MiSqguomSVVF1AjS6ouoGaWVF1AU0Rmdx7R\n",
       "iYj0ORZJUruYK5PneyhJaqfJ5Ip3TCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJ\n",
       "kiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJ\n",
       "kiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZU\n",
       "kiRJklQpG1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlbIx\n",
       "lSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVs\n",
       "TCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQp\n",
       "G1NJkiRJUqVsTCVJkiRJlbIxlSRJkiRVysZUkiRJklQpG1NJkiRJUqVsTCVJkiRJlRq1MY2IL0XE\n",
       "7RFx1QivOSUiro+IKyLigPaWKEmSBjObJUlNM5Y7pl8GjhzukxGxEJifmXsCfwt8pk21SZKk1sxm\n",
       "SVKjbD7aCzLz/IjYfYSXvAj4avna30TE1hExNzNvb0+JUvctiFg4DxbNhukrYc1SOOXqzLOqrkuS\n",
       "wGyWJDXPqI3pGOwMLB308TJgF8DwU09aELHwUDj5NJg/cOxY2GNBBDanknqE2SxJ6intWvwohnyc\n",
       "bTqv1HXzYNHgphTgNJi/KxxfVU2SNAFmsxongv+I4PFV1yGp/dpxx/RWYN6gj3cpj20iIk4a9OGS\n",
       "zFzShutLbTUbprc6PgtmdLsWSRtExOHA4RWX0SvMZjXVu4GtgHdUXYik9mZzOxrTM4HjgDMi4inA\n",
       "iuGeYcnMk9pwPamjVsKaVsdXwepu1yJpg7JhWjLwcUScWFkx9Wc2q8l2rroASYV2ZvOojWlEnA4c\n",
       "BmwfEUuBE4EtykI+l5lnRcTCiLgBuB9480SLkepgKZxyLOwxeDrvMXDjLXBqlXVJ0gCzWX3u0VUX\n",
       "IKn9IrM7j5xERGbm0OddpFpaELFwBkf8ZBar1s/ggp/dAqe68JFUL+bK5PkeqtdEkMBtmd41lepo\n",
       "MrliYyoNI4IHganAlEwerroeSRszVybP91C9pmxMATbLdEEvqW4mkyvtWpVXaqI1wApgz6oLkSRJ\n",
       "MIszHzyI565/KVteeFTE4gURC6uuSVJ7tGPxI6mpZgA/Bg4B/lBxLZIk9bXHx5SFC3nMtG9xIxTZ\n",
       "7D7jUoN4x1RqIYItKL4/zgcOrrgcSZL63q7EP5RN6SPcZ1xqDhtTqbUZwAPAMmCnimuRJKnvbclm\n",
       "LfcTd59xqRlsTKXWZlDsW3ovMKfiWiRJ6nsr2WJ9q+PuMy41g42p1NpMijumNqaSJNXAH3nKN17D\n",
       "rusGH3Ofcak5XPxIam2gMb0PG1NJkir3Z8797XK+tey5vP22ZK8DpnHxL91nXGoOG1OpNafySpJU\n",
       "LzPu4+g7zuXo1wG/yuSoqguS1D5O5ZVaGzyVd6uKa5EkSWaz1Gg2plJrA3dMHwC2KLePkSRJ1RnI\n",
       "5pXArAh/j5WaxG9oqbWZwAOZJD5nKklSHQxk80PA/cDsiuuR1EY2plJrA9OFwMZUkqQ6MJulBrMx\n",
       "lVrbkg37orkAkiRJ1TObpQazMZVa2we4ofyziyxIklS9fYDryz+bzVLD2JhKrR0E/Lb8s6OykiRV\n",
       "z2yWGszGVBoigs2B/YFLykP3U0wfkiRJFYhgKrAAuKw8ZDZLDWNjKm1qJ+C+TO4tPzb8JEmq1qOB\n",
       "uzNZVX5sNksNY2MqbWoOPNKUguEnSVLV5gArBn1sNksNY2MqbapVYzqrolokSZLZLDWejam0Ke+Y\n",
       "SpJUL2az1HA2ptKmnC4kSVK92JhKDWdjKm1qaPitwvCTJKlKZrPUcDam0qYclZUkqV7MZqnhbEyl\n",
       "TRl+kiTVi9ksNZyNqbQpw0+SpHpx/Qep4WxMpU25JL0kSfViNksNZ2Mqbco7ppIk1YvZLDWcjam0\n",
       "qTnAfYM+duU/SZKqZTZLDWdjKm1qJsVI7ABHZSVJqpbZLDWcjam0qRnA6kEfG36SJFXLbJYazsZU\n",
       "2tTQ8FsNTI9gSkX1SJLU71o2phFERfVIajMbU2lTMxkUfpk8XH48s7KKJEnqb0OzeR3wMDCtsook\n",
       "tZWNqbSpGcADQ445ZUiSpOqYzVLD2ZhKmxo6XQgMP0mSKlFO1zWbpYazMZUGKZ8j3QJ4cMinXJZe\n",
       "kqRqTAUeyuShIcfNZqlBbEyljU0H1mSSQ447KitJUjVa3S0Fs1lqFBtTaWMbLa4wiOEnSVI1bEyl\n",
       "PmBjKm2s1eIKYPhJklQVs1nqAzam0sZGGpWd1eVaJEnSyLOZzGapIWxMpY05XUiSpHoxm6U+YGMq\n",
       "bWy4UVlX/pMkqRrDNaZms9QgNqbSxhyVlSSpXsxmqQ/YmEobc4EFSZLqZSZms9R4NqbSxhyVlSSp\n",
       "XsxmqQ/YmEobc1VeSZLqxWyW+oCNqbSxkZakd1RWkqTuM5ulPjBqYxoRR0bEdRFxfUS8p8Xnt4+I\n",
       "xRFxeURcHRFv6kilUnfMpAi6oQw/SbVhNqvPmM1SHxixMY2IKcAngSOBxwOviYh9hrzsOOB3mbk/\n",
       "cDjw0YjYvAO1St2wDbC8xXGXpJdUC2az+pDZLPWB0e6YHgzckJk3Z+Y64AzgxUNe8xdgq/LPWwF3\n",
       "Z+b69pYpdc02wIoWxx2VlVQXZrP6jdks9YHRRk93BpYO+ngZcMiQ15wG/DwibgNmA69qX3lS1w03\n",
       "Kmv4SaoLs1n9xmyW+sBod0xzDOd4P3B5Zj4a2B/4VETMnnRlUjW2xvCTVG9ms/qN2Sz1gdHumN4K\n",
       "zBv08TyKkdnBngp8GCAzb4yIPwF7A5cMPVlEnDTowyWZuWSc9UqdNtKorEvSSxWKiMMpnpfsd2az\n",
       "+s1IU3nNZqlC7czmyBx+4LVcKOEPwHOA24CLgddk5rWDXvMx4N7M/FBEzAUuBfbLzHuGnCszM9pR\n",
       "tNQpEVwLvCKTa4Yc3xK4K5MZ1VQmaah+zRWzWf0mghuAhZn8ccjxrYFbMh95nlpSxSaTKyPeMc3M\n",
       "9RFxHHAOMAX4YmZeGxFvKz//OeBfgS9HxBUUU4PfPTT4pB4y3HSh1cC0CKZk8lCXa5KkR5jN6kMj\n",
       "TuWNIDLHNMVdUo2NeMe0rRdyVFY9III1wDaZm27kHcEq4FGZrOx+ZZKGMlcmz/dQdRdBAOuAGZms\n",
       "a/H5B4E5mazpenGSNjGZXBlt8SOpb0QU03RbNaUlF1mQJKm7ZgNrWjWlJbNZaggbU2mD4RZXGGD4\n",
       "SZLUXcMtSjjAbJYawsZU2mC4Z1gGGH6SJHXX1jhoLPUFG1Npg7GMyrosvSRJ3WM2S33CxlTawKm8\n",
       "kiTVi1N5pT5hYyptMFr4rcLwkySpm8xmqU/YmEob+IypJEn14jOmUp+wMZU2cLqQJEn1YjZLfcLG\n",
       "VNrAZ0wlSaoXG1OpT9iYShuMZSqvK/9JktQ9Y5nKazZLDWBjKm3gqKwkSfViNkt9wsZU2sCV/yRJ\n",
       "qhezWeoTNqbSBjsCd47weUdlJUnqLrNZ6hM2ptIGOwO3jvB5w0+SpC6JIDCbpb5hYyoBEcwGNgfu\n",
       "HeFlhp8kSd2zNbA2k1UjvMZslhrCxlQqPBq4NZMc4TWGnyRJ3TPa3VIwm6XGsDGVCmMNP5eklySp\n",
       "O8xmqY/YmEoFR2UlSaoXs1nqIzamUmEn4C+jvMYl6SVJ6h6zWeojNqZSYQ4jL3wEjspKktRNW2E2\n",
       "S33DxlQqzAJWjvIaw0+SpO6ZDSOuyAtms9QYNqZSYSzhtxqYFsGULtQjSVK/m80YB43LPU8l9TAb\n",
       "U6kwaviVW8k8gCOzkiR1w1iyeT2wHpjelYokdYyNqVQYy6gsOGVIkqRuGctjNmA2S41gYyoVxtqY\n",
       "uvqfJEndMZbHbMBslhrBxlQqjDX8HJWVJKk7nM0k9REbU6ngdCFJkurFxlTqIzamUsHwkySpXmbh\n",
       "bCapb9iYSgUbU0mSaqLc/sVslvqIjan6XgSbA1Mp9ikdzf0UI7iSJKlzpgEPZ7J2DK81m6UGsDGV\n",
       "yoWPyn1KR+OorCRJnTfWu6VgNkuNYGMqjX3hI3BJekmSumGsz5eC2Sw1go2pVITZ/WN8raOykiR1\n",
       "ntks9RkbUwlmAg+M8bWGnyRJnTcTG1Opr9iYSjamkiTVzUzGtighmM1SI9iYSuNvTF35T5KkzjKb\n",
       "pT5jYyrBDByVlSSpTmbgbCapr9iYSk7llSSpbpzKK/UZG1Np/I3pjA7WIkmSzGap79iYSuObyvsA\n",
       "RVhKkqTOGc9UXrNZagAbU2l8o7KGnyRJnWc2S33GxlQy/CRJqpvxPGNqNksNYGMqOZVXkqS6cSqv\n",
       "1GdsTCXvmEqSVDdms9RnbEwlw0+SpLpxKq/UZ2xMpfFN5V0HbBbBFh2sR5KkfjeeqbyrgRkRRAfr\n",
       "kdRhNqbSOO6YZpI4MitJUqeNJ5sfBh7EvUylnjZqYxoRR0bEdRFxfUS8Z5jXHB4Rv4uIqyNiSdur\n",
       "lDprPFN5wcZUUsXMZvUBs1nqM5uP9MmImAJ8EngucCvw24g4MzOvHfSarYFPAUdk5rKI2L6TBUsd\n",
       "MJ6pvGD4SaqQ2aw+YTZLfWa0O6YHAzdk5s2ZuQ44A3jxkNe8FvhuZi4DyMy72l+m1FHjWWABDD9J\n",
       "1TKb1Q/MZqnPjNaY7gwsHfTxsvLYYHsC20bELyLikoj4m3YWKHXB1sCKcbze8JNUJbNZ/cBslvrM\n",
       "iFN5gRzDObYADgSeQ/ED4cKIuCgzrx/6wog4adCHSzJzyRjrlDppe+DOcbze8JMqEBGHA4dXXEYd\n",
       "mM3qBztgNku1185sHq0xvRWYN+jjeRQjs4MtBe7KzNXA6oj4JfBEYJPwy8yTJl6q1H4RTKMIsnvH\n",
       "8WWGn1SBsmFaMvBxRJxYWTHVMpvVaBHMBAIXP5Jqr53ZPNpU3kuAPSNi94iYChwNnDnkNT8Enh4R\n",
       "UyJiJnAI8PuJFiR12XbA3eU2MGNl+EmqktmsptsBuMtslvrLiHdMM3N9RBwHnANMAb6YmddGxNvK\n",
       "z38uM6+LiMXAlcDDwGmZafipV4x3Gi/AKmB2B2qRpFGZzeoDZrPUhyJzPINRk7hQRGZmdOVi0hhF\n",
       "8BzgA5k8axxf83FgaSYf61xlkkZjrkye76HqKIIjgXdl8vxxfM0ngT9kcmrnKpM0msnkymhTeaWm\n",
       "m8io7HJg2w7UIkmSzGapL9mYqt/tAIx3f797gG06UIskSSoaU7NZ6jM2pup3uwB/GefXLMfwkySp\n",
       "U+ZhNkt9x8ZU/e4Q4Lfj/BqnC0mS1DkTyeZ7MJulnmZjqr4VwRbAk4GLxvmlTheSJKkDyv3FDwAu\n",
       "HueXesdU6nE2pupnewG3ZbJinF9n+EmS1Bn7ADdnsnKcX+egsdTjbEzVz3YCbpvA1zmVV5KkzpgL\n",
       "3DqBrzObpR5nY6p+tiNwxwS+bjmwdQTu/SdJUntNJpu9Yyr1MBtT9bMdmED4ZbIOeBCY1faKJEnq\n",
       "bzsw/j1MyWQ1EBFMb39JkrrBxlT9bEcmEH6lldiYSpLUbhO9Ywpms9TTbEzVlyKYQbFP2kTDbxWG\n",
       "nyRJbVNm8y5MfNDYbJZ6mI2p+tWJwBuwMZUkqS7+D/A6zGapL21edQFSNyyIWDgPFs2G6SthzRy+\n",
       "vM29vAngrgme0vCTJKm99i3/6x1TqQ/ZmKrxFkQsPBROPg3mDxx7DSeu+zE7sIoXLJvgaQ0/SZLa\n",
       "a6/yv2az1IecyqvGmweLToP5P5kKRzwaDt8N7nn0LVs8Ztr/d1kmN03wtIafJEmTFMGXI5gdwTTg\n",
       "0cDUTJZO8HRms9TDvGOqxpsN038yFU7YE2585YbjM/6bfWJuLMzb86wJnNbwkyRpEiLYAXgT8B3g\n",
       "dIqmdN0kTmk2Sz3MO6ZqvJWw5pTtN25KAVa/ghlsxfETPO0qYMvJ1iZJUh97YvnfLwAXAYdM8nxm\n",
       "s9TDvGOqxlsKp9wxdcqz4aEtNvnkVGZM8LSOykqSNDlPBJYDjwKen8nVkzyf2Sz1MO+YqvGuzjzr\n",
       "nvW7rGj5ybWsnuBpDT9JkiZnH+DDwJFtaErBbJZ6mo2p+sJDd514A9+f9deNDp7JjdzHqRM8peEn\n",
       "SdLkzAWuz+ScNp3PbJZ6mFN51R/WvDm56ecf4+tffzZTmcFaVnMfp05w4SMw/CRJmqwdgTvaeL5V\n",
       "wB5tPJ+kLrIxVb/YkZVf+2He97WPtOl8NqaSJE3OXOD2Np7PbJZ6mFN51S/m0v5RWcNPkqQJiCDo\n",
       "zB1Ts1nqUTamarwIZgDTgHvbeFrDT5KkCVgQsfAItvjp4Txz2pHEdxdELGzTqc1mqYc5lVf9YAfg\n",
       "jkyyjec0/CRJGqcFEQsPhZNPY/18OB/giGNhjwURXJ0TXvdhgNks9TDvmKoftHsaLxh+kiSN2zxY\n",
       "dBrMH3zsNJi/KxzfhtObzVIPszFVPzgQuLbN5zT8JEkap9kwvdXxWTCjDac3m6UeZmOqfvAC4Cdt\n",
       "PqfhJ0nSOK2ENa2Or4LVbTi92Sz1MBtT9YNnAOe2+Zz3A7PKVQUlSdIYLIVTjoUbBh87Bm68BU5t\n",
       "w+ltTKUeFpntXA9mhAtFZGb6S7y6KoKpFEE1rc2LHxHBGmDrzNajv5I6y1yZPN9DVeFxseNLt2L/\n",
       "/96dn56/ClbfAqe2YeGjgS1o1gHTM1k/+UoljddkcsVVedV02wN3t7spLQ2MzNqYSpI0Rn/gjkuB\n",
       "Wy9ODm/neTPJCFYBW9LeLeIkdYFTedV02wN3dejcThmSJGn8dsBsljSEjamarpON6f0YfpIkjZfZ\n",
       "LGkTNqZqOu+YSpJUL2azpE3YmKrpnC4kSVK92JhK2oSNqZpue+DODp3b8JMkafzMZkmbsDFV020P\n",
       "3N2hcxt+kiSNn9ksaRM2pmq6OcCKDp17FTC7Q+eWJKmpzGZJm7AxVdPNolihrxNWYvhJkjReZrOk\n",
       "TdiYqum2pBg97YT7MPwkSRqvTmazjanUo2xM1XRb4qisJEl1YjZL2oSNqZrO6UKSJNWL2SxpEzam\n",
       "ajqnC0mSVC8+ZiNpEzamajqnC0mSVC+dzuatOnRuSR1kY6qmc7qQJEn1YjZL2oSNqRorgsA7ppIk\n",
       "1UYEmwHTgNUduoTZLPWoURvTiDgyIq6LiOsj4j0jvO6giFgfES9rb4nShE0D1meyrkPnN/wkVcJs\n",
       "Vg+bCazO5OEOnd9slnrUiI1pREwBPgkcCTweeE1E7DPM6/4DWAxEB+qUJqKTU4XABRYkVcBsVo/r\n",
       "dDbbmEo9arQ7pgcDN2TmzZm5DjgDeHGL1x0P/DdwZ5vrkyajk6v+geEnqRpms3qZ2SyppdEa052B\n",
       "pYM+XlYee0RE7EwRiJ8pD2XbqpMmp5PPlwI8CEyJYGoHryFJQ5nN6mWdzub7gRkRTOngNSR1wGiN\n",
       "6ViC7BPAezMzKaYKOV1IddHR6UKZJMV0Xpell9RNZrN6Waez+eHy/N41lXrM5qN8/lZg3qCP51GM\n",
       "zA72JOCMiADYHjgqItZl5plDTxYRJw36cElmLhlvwdI4dHq6EMAKYA5wV4evI/W9iDgcOLziMurA\n",
       "bFYv62Y2r+jwdaS+185sjmIwddgLbQ78AXgOcBtwMfCazLx2mNd/GfhRZn6vxecyMx2xVddE8ELg\n",
       "7Zm8oIPXuBT420wu7dQ1JLXWr7liNquXRfBS4I2ZvKSD17gS+JtMrujUNSS1NplcGfGOaWauj4jj\n",
       "gHOAKcAXM/PaiHhb+fnPTeSiUpfMpphq20nLgW06fA1JeoTZrB5nNktqabSpvGTm2cDZQ461DL3M\n",
       "fHOb6pLaYSs6H34rgK07fA1J2ojZrB5mNktqabTFj6ReNpti2fhOWoGjspIkjZXZLKklG1M1WTdG\n",
       "ZZfjqKwkSWNlNktqycZUTdatUVnDT5KksTGbJbVkY6om69ZzLE4XkiRpbMxmSS3ZmKrJujEq63Qh\n",
       "SZLGzmyW1JKNqZrMUVlJkuplK1z8SFILNqZqsm6Myt4DbNvha0iS1BTd2MfUbJZ6kI2pmqwbd0zv\n",
       "Arbv8DUkSWqKbtwxNZulHmRjqibrxh1Tw0+SpLHrxh1Ts1nqQTamarJu3DG9F5gZwdQOX0eSpCbo\n",
       "xh3T5cCcCDbv8HUktZGNqRopgqALd0wzSeBuYLtOXkeSpF4XwWbATGBVJ6+TyXqKgWMXQJJ6iI2p\n",
       "mmoGsC6TdV24llOGJEka3SzggUwe7sK1zGapx9iYqqm68XzpAMNPkqTRdeP50gFms9RjbEzVVIaf\n",
       "JEn10o3nSweYzVKPsTFVUxl+kiTVi4PGkoZlY6qm6mb4LQe27tK1JEnqVd0cNDabpR5jY6qm6mb4\n",
       "rcDwkyRpNN0cNDabpR5jY6qmMvwkSaoXB40lDcvGVE1l+EmSVC8OGksalo2pmsrwkySpXhw0ljQs\n",
       "G1M1leEnSVK9dHOPcbNZ6jE2pmoq75hKklQvW2E2SxqGjamaqtt3TLfp0rUkSepV3jGVNCwbUzVV\n",
       "N++Y3gtsHUF06XqSJPUi75hKGpaNqZqqa3dMM1kDPATM6Mb1JEnqUd28Y3o/MC2CqV26nqRJ2rzq\n",
       "AqQO6eYdUyhGZrcFHujiNdVFCyIWzoNFs2H6SlizFE65OvOsquuSpB7StTummWQEyyketbm9G9dU\n",
       "95nNzWJjqqbq5jOmAHcAOwLLunhNdcmCiIWHwsmnwfyBY8fCHgsiMAAlacy6eccUNmSzjWkDmc3N\n",
       "41ReNVW375jeDszt4vXURfNg0eDgAzgN5u8Kxw/3NQsiFh4VsfhVEUuOili8IGJh5yuVpFrr5jOm\n",
       "YDY32kSyWfXmHVM1VbfvmBp+DTYbprc6PmuY54odxZWklrp9x9RsbrDxZrPqzzumapxydVzDT22z\n",
       "Eta0Or4KVrc67iiuJG0sgs2AmcCqLl7WbG6w8WYzwGPj0Dc9NZ78l1fGZuc5m6l+bEzVRNOBdZms\n",
       "6+I1Db8GWwqnHAs3DD52DNx4C5za6vWO4krSJmYBD2TycBevaTY32HizuZjN9OdP/ppLd/oO+cyz\n",
       "4YhD4WSb0/pwKq+aqNvTeKEIv/27fE11ydWZZy2I4GCe9d05LH9gNblqBVcsGm5a7kRGcSWp4br9\n",
       "fCkU2bx3l6+pLtmQzc/+/hzuWbWah1eu4Mphs3keLPoGf9ly8LHTYP7CYjaTj9nUgI2pmqjbCx+B\n",
       "o7KNdw15DjAF2Ae4HrhouNfezDafPZptnv8tboqBY69iD27hxpajuJLUB7r9iA2YzY13DfkzINiQ\n",
       "zRcM99otmbJlse38xpzNVB82pmqiKu6Y3gns0OVrqrvmAvdkckcE/wO8HDit1Quv454/3ca3b1nI\n",
       "0ddOZc4Of+GJe/+ed89axQsWu+eapD5VxR1Ts7n5HgXcXmbzL4CXAV9u9cIVzJkC92xyfKWzmWrD\n",
       "xlRNVMUd07uB7bp8TXXXq4Bfl38+C3g2LRrTcvGtE+/jVeeela9666Djf96dw958KLzX1Xol9aEq\n",
       "7piazc03NJufSYvGNILN5vDReD3vue/r3LHVwPGj2f2hP7LTt7pTqkbj4kdqoirumN6F4dd0bwP+\n",
       "o/zzVcAThnndnsBBwDuHHL9yd248rtVqvbuw2QkRvKCt1UpSvVSRzTamzTfWbN7nXt6006U8+IaF\n",
       "sPhVcN5CWHwOH7r8Ji68GyCCfSO4NIKfRzA/ghkRuDBSF3nHVE20DbC8y9d8AIgIZmbyQJevrQ4r\n",
       "74LuDlxdHvo9sHcEW7RY/flQ4IJMHhxy/LJpbLlfq/NvwQ47Az+KYN9Mrm1j6ZJUF1Vk873AjAim\n",
       "ZrK2y9dWh5VbEM0DrikPXQ3sE8GUzE0eJj0U+NW1ueKHwA8HneNfKQaTfwS8EJgDXEIxI+pTs/jx\n",
       "t54Xr71gG1Y+7OM3necdUzXRNrR6iKCDMkkcmW2y7YHVmdwPUA4+LKP1ao+HAhe2OH7BvWw3u9XJ\n",
       "b2fBnsCtwAfbU64k1U5V2bwc2Lab11XXzAVWDAwEZ7ISuIMhM5NKw2Yz8LTyz4cD/wj8HbB+Nj/4\n",
       "wgs4/uGfsvIZ34bD3F6m82xM1UTb0uXwK9mYNtc8YOmQYz+HllN8ngOc1+L4RVfz3unHEDcOPvgW\n",
       "Nv/ztZwA8A7g2RHs24Z6JaluzGa125iyuZz1NFw2XwgcHMEM4KnA+ZmsB/52bz691RncvNHs0tNg\n",
       "/q7F9jKQ++SiAAAgAElEQVTqABtTNZHhp3ZrFX7fplh04RER7AlsCVyxyRl2jKev2vkj67+yx3Yz\n",
       "5jx664cOmcZVC2Hxxax/5ypeuB/wE+ATwPs68jeQpGpVlc2uAdFcY8pmiq1kNqN4DGcjmdwD3AK8\n",
       "Bbix/JhM/rQLF/621UXdXqZzfMZUTbQt3X+OBWxMm2wXiqm7g50HzItgfiY3lMeOAhaX08ceEXNj\n",
       "IbtwMi/61eyHYPZ9wMVnMoNlnJq3b3hWJYLPArdEMD2TNR38+0hStzlorHZrlc0/B74ewW6Z/Lk8\n",
       "dhRw9tBsHuQC4J+B0wcfXMOqlr9LrnJ7mY7xjqmayFFZtdtewE2DD2Synumf/i277nN27BtLYs9Y\n",
       "zLSv/g1w9iZfvRWLeNGQZ15exHy22ng6UCbLKRZveGq7/wKSVDEbU7Vbq2xeB3wfeO2gw0fRKps3\n",
       "uJjiedWvDT64FE45lkcGngE4Bm68BU6dTNEanndM1URVhd9fKTZ6VvM8Azhu8IGYGwvZZ+oTePHa\n",
       "XRlYaOE7H3qYv5wE/Gnjr57K9JZnndpyOtBPgSMpRn0lqSnMZrXbM4GvtDj+CeC8CL4GPAgcArx0\n",
       "hPN8A7g8k8sGH7w686wFxX4xxz/MXo+/ny3WLeeaRa7K2zk2pmqirq/8V7oFOKyC66qDItiaYm/S\n",
       "Szb6xFYsKpvSDV75p834Om8BvrPR8bXDTMtd23I60PeB70XwnhGmHUlSr6kymw+q4LrqoAi2B3YD\n",
       "fjf0c5n8PoIvAh8FZgJfKFfsbSmT1QzN+FLZhJ4VwUEUDexId141SU7lVRNVNSr7Z2DXUV+lXvME\n",
       "4KpN9sAbz13Q+ziFMzeeDsSZ3Mh9LacDXQ6sB/afWLmSVEtms9rpCRR3OdcP8/mPUNyA+z3wgTZc\n",
       "75LyfAvacC4NwzumapRys+U5wIoKLn8LxeidmmV3Npmby7jugubteVbMDfg6xzOVGaxlNfdtvPDR\n",
       "I69NMoILgCfTYiRYknpNBFtQ3Lka9q5VB5nNzbQ7rbK5lMndwMvbdbEymy+kyOar2nVebczGVE2z\n",
       "LXDvCCNonbQU2DmCKZk8VMH11Rm7AzdvcrS4C7rHRosaDX8XlLIJHetzKVcATxxvoZJUU9sDd2Xy\n",
       "cAXXvgXYNYLw8YhG2Z1W2dxZl2M2d9SYpvJGxJERcV1EXB8R72nx+ddFxBURcWVEXBAR+7W/VGlM\n",
       "5gJ3VHHhcnuP5cBOVVxfHbM7PLLk/CPy9jyLZZzA11nMtzmPr7OYZSxqdRd0Ai7HqbwahdmsHlJl\n",
       "Nq8CHgB2qOL66pjd6X5jegVmc0eNesc0IqYAnwSeC9wK/DYizszMawe97CbgmZl5b0QcCXweeEon\n",
       "CpZGsSMVhV9p4FmWWyusQe21O/CtVp8Y513Q8bgMeHwEu2ZySwfOrx5nNqvH1CWbq6xB7bU73W9M\n",
       "LwH2i+DRmdzW5Wv3hbHcMT0YuCEzb87MdcAZwIsHvyAzL8zMe8sPf0Ox4a1UharDz2dZmucxdDn8\n",
       "MrkX+AKwqJvXVU8xm9VLzGa1WxXZfA/wVczmjhlLY7ozxbNzA5aVx4bzVjpzB0Eai6rDz9X/GiSC\n",
       "mRT73w27wEIHfQ94TgXXVW8wm9VLzGa1TQSzKZ5brmJGkdncQWNZ/GjMD4pHxLOAtwBPm3BF0uRU\n",
       "HX63AHtXeH211+OAGzJZV8G1LwPmR7B1ZiWrTKvezGb1kjpks3dMm2Mf4LqKFpr8LbBPBLNH2htV\n",
       "EzOWxvRWYN6gj+dRjMxupFxU4TTgyMxc3upEEXHSoA+XZOaSMVcqjc2OVLvFxp+B51V4fbXXvsA1\n",
       "VVw4k7UR/IaimfhJFTXUTUQcDhxecRl1YTarl+wIXF/h9f+MAzNNUmU2r4ngUuCpwDlV1FA37czm\n",
       "sTSmlwB7RsTuwG3A0cBrhhS0K8Wt7ddn5g1DTzAgM0+aaKHSGDkqq3Z6PMXm3FX5JfBMbEwBKBum\n",
       "JQMfR8SJlRVTPbNZvcRsVjvVJZttTGlvNo/6jGlmrgeOo3jzfw98KzOvjYi3RcTbypd9ENgG+ExE\n",
       "/C4iLp5oQdIk7QjcXuH1fY6lWeZT7Sj/QPhJGzGb1WOqbkzN5mYxmxsqMruz13BEZGZGVy6mvhXB\n",
       "9cDCzGp+YEUQwEpg53JlVfWwcrrOOzKp5Bf6cvGlu4A5FT3nWmvmyuT5HqobIrgZeFZmJQvJEcFm\n",
       "wP3A9pncX0UNap8IrgTemFnNo1vl4ku3A7Mres611iaTK2NZlVfqJZVt4g2QSeLIbCOUgwx7ADdW\n",
       "VUMmD1CsvLpXVTVI0mSUP0urzuaHKX6Wzhvttaq38t/TYyn2aa5EuejRXyju3KqNbEzVGBHMAKYB\n",
       "91Vcis+yNMO2QAD3VFzHlcB+FdcgSRM1C3i4BncqzeZmmAusqcGstKswm9vOxlRNsgNwR3nXskre\n",
       "MW2GPSm2iqn639NVwBMrrkGSJqrq50sHmM3NsCcw7GJuXeSgcQfYmKpJ6hJ+f6KYAqredgjFfmVV\n",
       "+zVukSKpd1W9KOEAs7kZzOYGszFVk9SlMf0DsHfVRWjSnkoRPFX7JcVm3jtUXYgkTYDZrHaqSzYv\n",
       "AfaLYNuqC2kSG1M1ieGndnoKcFHVRWSyFjgXOKrqWiRpAsxmtVNdsnkNcB5wRNW1NImNqZqkLuF3\n",
       "IzAvgqlVF6KJKZeC34EKV+Qd4ifAC6ouQpImoC7ZfD3w2Ag2r7oQTUwEWwOzgZsrLmWA2dxmNqZq\n",
       "klo8x1Le4VpKsZy5etNewPU12p/sLOD5EWxRdSGSNE61aEwzWQ38FVfm7WV7A3+owaKEA84Cjoxg\n",
       "StWFNIWNqZqkFuFX+iNOGeple1NM+6qFTP5Ccff2aVXXIknjZDarXR5HvbJ5KXArxfRitYGNqZqk\n",
       "TuHnsyy97XHAdVUXMYRThiT1IrNZ7bI3ZnOj2ZiqSeZi+Kk99gOurrqIIQw/Sb3IbFa7mM0NZ2Oq\n",
       "JnFUVu1yIHBZ1UUMcQmwXQR7VV2IJI2D2ax2qWM2XwTsHOG6Iu1gY6pGiCAoVlG9s+paSlcDCyL8\n",
       "Hus1EWwPzAFuqrqWwTJ5GPhPYEkEM6uuR5JGUy4Ksy1wV9W1lK4CnlD+zqAeEsGjgGnALVXXMli5\n",
       "SOL/BX4ZwfSq6+l1/tKspngUsDyTB6suBCCTu4B7cWXeXnQA8LuyEayVTD4KXIv7pknqDbsAt2ey\n",
       "vupCADL5K7AW2LXqWjRuBwCX1WhF3kdk8u8Ug9nPrbqWXmdjqqbYi2KPsjq5jGLaiXpLHacKDfYd\n",
       "4H0RbFd1IZI0ij0pVsKtE7O5Nx1AvbP528D/imDbqgvpZTamaoo9qV9jeilwcNVFaNwOAH5XdREj\n",
       "+ALFL3rvqroQSRpFHQeNzebedCD1zubPAX8GFlVdSC+zMVVT1HFU9nzgGVUXoXF7EjUelS2nxP07\n",
       "8CY39ZZUc3UcNDabe1Pds3kdRTa/2fVFJs43Tk1Rx/D7DbBvBLOrLkRjE8FjKBY+qts+aRvJ5Grg\n",
       "btzUW1K91TGbLwT2dxG53lGuRj8VuKHqWkaSyeXA/cCTq66lV9mYqilqN10okzUUU4aeWnUtGrMX\n",
       "Aj8uV9mruzOBl1RdhCSNoHazmTK5H7gSB/Z6yQuBH9VxUcIWzOZJsDFVzyunTDyWeo6k/RI4rOoi\n",
       "NGYvpAiVXvBVium886suRJKGimBzYDdqtvVW6TzM5l7SS9n8FeDYcgaWxsnGVE0wD7inHAWtG8Ov\n",
       "R0QwBzgE+GnVtYxFJtcDHwX+repaJKmF3YC/lrOH6sZs7hHlCvQHAD+vupaxyOQ64BTgw1XX0os2\n",
       "r7oATc6CiIXzYNFsmL4S1iyFU67OPKvqurqsjs+wDHjkWZZMHqi6mCZrw/fCi4HzajrAMZxTgVsj\n",
       "2DGTO6ouRpIGqXM2/xp4cgTT6rL/uYb1EuDcTFZXXcg4fJwim7fL5O6qi+klNqY9bEHEwkPh5NPY\n",
       "MJXvWNhjQQR91pzuQ00Xq8nk/ohHnmXpidG+XtSm74U3A5/qTIWdUf77uhg4CPhJ1fVI2sCB41pn\n",
       "830RXEuxbcz5VdfTRG389/9m4CPtrq+TMlkVwWUUiyCdU3U9vcSpvD1sHiwa/Is4wGkwf1c4vqqa\n",
       "KnIAcHnVRYzgPOBZVRfRZJP9XohgG4oA+VEn6uuw3+KefFKtDAyWnQ1HfBsOOxuOOBROXhCxsOra\n",
       "ushs7lOt/v0/mdlf3D2e/7fjOU8E2wNPAHpxQOdizOZxszHtYbNheqvjwbzHRvD4btdToQOo96bL\n",
       "P8YV2jpqOrO3bXU82G2PMe4n9jTgNz06pWsJ8M4IDq26EEmFnZjz/laDZY9iq/dWVVMFzOY+1Wqw\n",
       "+Cus3GlHHjo1gl9F8I4IYgynejrw63KP0F5zHnBchM3peNiY9rCVtF5Q4B523g74fgSHRtDyF/am\n",
       "iGAq8DjgqqprGcGvge0ieFzVhTRRBIfeykFPbPW5u9llB+DM8t/JSJ5JsYJyz8nkp8AngDdWXYuk\n",
       "wmbM27PV8YfZ50llNm/T7Zq6KYLpwB7ANVXXMoLzgV1c2bz9hrtx8mh+dw3wLxTTc38QwRajnOow\n",
       "ejebzwI+A7yh6lp6iY1pj4ogruZV57yax2y0p9MxcONKLnoDxVYS3wM+F8F2ERxRBkXTHAhcV+eF\n",
       "hcp9t74KvLPqWpomgrnAd5ex178eO2S7oGPgxvu49A3AeuC75aq7rc6xGfAKenOq0IAfAc+uugip\n",
       "30UQETx1BTu2bDzXce+tFNn86Qi2jeDICKa1eu2CiIVHRSx+VcSSoyIW99g04CcDV9d5Fkq5X/VX\n",
       "gXdUXUvTrGBGtjq+luW3Z3IOxSylKcB3Ipjd6rVlNr8cs7mvRGbLfzvtv1BEZuZYbttrFOX0h28C\n",
       "z5vLu798IB9ZMAtmrILVt8CpAw+XRzALuBGYCcwCLqDYVHo28BZgK+C1wN3AD3tsNVIAIvhH4DGZ\n",
       "HFd1LSOJYBeK9/4xmdxbdT1NEcG7gb0yOWZBxMJd4fih3wvlgMxXgVvYMX7BVixiKtNZyxru4xTu\n",
       "yLUUdxyfkEl3fiC2WRngNwHvyuR7VdfTLebK5Pketk/5ffgt4PC5/NOXXshHXzZ4OuMxcONFsOga\n",
       "8lcU2TydIpt/xcbZPGcH/veHj+L/vuSrrJw78PXHwg0Xwgm9sIBSBO8F5mbyD1XXMpIIdgMuA3bP\n",
       "ZGXV9TTFDvHvXzmCj7/069yx1cCxgX//g35HnQF8Dbg+k/cNPUcERwL/nsn+XSu8zSKYAtwM/F1m\n",
       "z+zDOmmTyRUb0x5Uzlf/JrBgtP3BypHYJ1GMTn2AYoRqbfm/mRSBOAXYHjgik7s6WHrbRXAmcHom\n",
       "p1ddy2giOB24OJOPV11LE5Q/8K8C3p458lSfCB7F1B9cy+Nfu4qXrN75kU/8cOot/P6/1vPg0f8r\n",
       "kzM6XHJHRfA0ip8Lu/dqgz1e5srk+R62T/k9+EXgiZk8ONxgWfna6RR3FZ8OvAeYCjxIkc0znsIh\n",
       "D17IxTsMvcZCWHxW5lFd+0tNUARnAV/K5L+rrmU0Efw3sCSTT1ZdSxNEsDlwzS685vNP4Izntvr3\n",
       "P+i1u1A8h/zSTH5VHpsF7AJ8Bzgpk+92++/QThEcBnwJmG82j+FrbUx7TwSnAPdkctI4v242sBNw\n",
       "X/nfqzJ5uLwD+wngb4ATgc9lsra9VbdfuenyjRS/iK+oup7RRPAUisZhz3IKkSYhgrdSPKfyjLH8\n",
       "sI/H7vRb3nj7kzf5xBf3X8bS3+3a64FRfh//HnhLJhdWXU83mCuT53vYPhF8Brg1k/8zzq8byOaV\n",
       "wI7AVa8kfvHt4vm6jbwKzvt25uHtqLdTItgR+CMwrxfuQkbwDIoBhceVj95oEiJ4O/Bq4Fljyubi\n",
       "zuhXgP8sD72bYqDmW8A7G5LNfwRek8klVdfTDZPJFZ8x7TERPAZ4HfCF8X5tJiszuT6T2zO5YuAH\n",
       "cCaZyQnAM4CXAouHex6vZl4NnNULTSlAJhcBdwAvrLqWXlf+oP974J/HHFozbm89VX325Tf2evBB\n",
       "8X1MEe79tl2UVLkI9gReSdHgjMugbP5rJldmksMtbngfm/fC6qSvBc7shaa09CuKQYHa34muu4lk\n",
       "cyaLKRbv2xc4BHg+sF0m72hQNn8Vs3lMbEx7z0eBj2WyrN0nzuQa4HkUq+j9LoKXtfsabfYmil/E\n",
       "e8nJwLuHW+xCY3YQMINiq5SxWTvMtPe1rG5LRfXwaeA5EexbdSFSvyh/Gf8Y8J+Z/KUd51wKpwxd\n",
       "0O117Hjfrzl93whe3I5rdNCb6KFsLhuHk4H3jGEFd43sqeV/fzWeL8rknEzemsnRA4MzHaitSqcC\n",
       "CyPYu+pC6s6pvD2ivIN5EvBc4KDRni1tw/WeTfFQ+qszOb+T15qICBYAi4HdemlabLk0+neBOzI5\n",
       "pup6elUEnwVuyeRfx/w1c2Mhu3AyLxq0NcCZ3MgyFuXt9V9MZKwi+P+Bp2by0qpr6TRzZfJ8Dycn\n",
       "gq2B/02x5dTB7XwMptUzqteQ64D/Al5Wxyn7EewP/AB4bC9Niy0b0u9T5Iqr9E5QBF8E/pD5yLRc\n",
       "lSJ4P7B/Jq+qupZO8xnThitXLvsZsAw4IZO/dum6RwP/Bjw/c+OR2yqVU6Y+DlyZyfurrme8ItgB\n",
       "uB7YL5Nbqq6n15SDNDdRvH+3jutr58ZCtuJ4pjKDtazmPk5tUlMKj/y8uBL4h0x+XHU9nWSuTJ7v\n",
       "4cRFsCXwc4qfRydkckeXrvs6ioHqIzK5qRvXHIvybtAnKBb5O7HqesYrgp2AaylWaG/7rLSmK/fm\n",
       "vQl4fLtmDjRJ+fPiSornZs+pup5OmkyubN7uYtRe5epmZ1B8s7+xmyOQmXyrbAJ+GcGTavSD5kMU\n",
       "y+x/tupCJiKTOyP4JHBFBI/L5Paqa+oxbwX+Z7xNKUDZhDaqER0qk9URvAv4cARn9dJdC6lXlLNf\n",
       "vs2GBce6NvUwk2+U2Xx+BAfWKEP+hWKV/9OqLmQiMvlrBJ+jyOa9e22Xgho4FvhRjX5XrJVM7i9n\n",
       "NP1rBP/TwOnKbeEd0xorV+r7IjAHeGFVK+VG8GHgUODoTO6sooZBtQysxLtHJndXWctkRfBpij3u\n",
       "3l6nUe9uWBCxcB4smg3TV8KapXDKWPbmK5/lupriPavdFPO6KN+nXwP/lclnqq6nU8yVyfM9HL+y\n",
       "KfwSxQDpSzKpZEGiCP4TOIDikZtK83DQSryPyWR5lbVMVgSfB+ZR5Myfq66nmyaZzddR3EC5qPOV\n",
       "9qZyr+PfUOx+Me5FTHuFq/I2UNmADawU9+KKt2/5EHA5xejsthXWAfA+4IyqQ7hN3g9cAvwsgudW\n",
       "XUy3LIhYeCicfDYc8W047Gw44lA4eUHEwjF8+f4U++9e0OEye1o5Evtm4F8i2KvqeqSmiGB7ip8/\n",
       "d1E851nlKrn/THHH9vzyWdcqvR/4Wq83paX3UOyRfW653kZfmGQ2HwSPNF0aRjmD6c3Av0WwR9X1\n",
       "1JGNaQ2VK+7dBfwVOKbTCx2NJpO1mbwLOBs4s6pvpgieTrHX6oequH67ZbKifEb2H4HPRnBOBC8r\n",
       "Rx4baxem/MNpDFqACDgN5u86tqXU3wl8wempo8vkOuC9FL+07lZ1PVKvi+DlwJ3Anynupj1YZT2Z\n",
       "PFhu9fZz4IfldnJdF8HhwNEUU3l7XibLM3k3xUD4FyI4O4KXmM0jGshmp6eOIpOrKQaVLohgl6rr\n",
       "qRsb05optxE5BTiSYopQnb7J303xfN653f5miuCJFCvmvaFpzy9k8n2K/btOBz4ArIjgggjeUC5k\n",
       "0xgRPGcdBz291edmMfLfNYKDKPbZ7cnnl6pQThX6BvC2qmuRelkE0ym2FHke8IqaZfO7KBZIPDeC\n",
       "R3fzwhEcCHwHeH23Fn/qlky+A+xD8SzxicDyCH4VwesamM3PXzvxbD4UeAE0d2pqu2XyWYrvm7+t\n",
       "upa6sTGtkQhmUuz99ftyT6da7a+Yybpye45PAeeU0407LoKdKbaueV9TVzIrR76/ksmBwB4U2w/8\n",
       "HbAsgs9GsEe5EFbPimA+8I3klmtbff5Odh12Klo5Uv154Lim/fLTBZ8Bjo14ZH85SeMQwSyKLVp+\n",
       "l8nPapjNazP5F4rGYHG5OmrHRTCPIpv/KZNzu3HNbiuz+cuZHEBxN/HDwAnA0gg+E8FjhmbzgoiF\n",
       "R0UsflXEkqMiFo9xKmxlytWUvwY3X9Pq83cyb84IX7sZRTa/oyGPWHXTp4C3R3BI1YXUiY1pvbwe\n",
       "2Al4Q9WFjCSTjwA/As6L4GmdvFb5DfsbiruJX+zkteoik7vKgYlDgCcDdwC/BR6M4BcRvCmCAyPY\n",
       "tdpKxyaCiOA44CLgA3dw2weGbhz/ZqYtvZRPzIvgnyKImBsLY89YHPvGktgzFrPlez9C8fPqW1X8\n",
       "HXpZJtdT/Gz5QQSPr7oeqQe9EdiG4tmwOvs34H8osvnQTl6oHOj6DfBViqa98cpsPjuTg4GDgbuB\n",
       "yyiy+dwI3rBrvHzRIWz+qQk+p9lVZTafQLFQ3nvu4K8f3DSbpy+7lJN3jeDvh5nK/GZgLcX+7BqH\n",
       "8nGbN1E8Ird3xeXUhqvy1kQ5LeQS4B8zWVx1PaMpf0C9DvgoReN0aiafH9PXFntJLmIq01nLGu7j\n",
       "lFZ7SUbwUuBzwFsz+VE76+815fs9FXghxS9JTwJWUQwQ/BT4ecULZD0igudQbEW1LfAYikUR9qSY\n",
       "mn4DDLtx/NXAj5lx6p/YZ9HjedGgZ12+s/tDLNt3Ua748ae7/zdqhgjeCvwD8MxM7qm6nnYwVybP\n",
       "93Bk5Uymyyhma/ys6npGU2bFG4CPALcDn8hs76BuBK8APg28KbPZ22+Npny/p1Fm80E8+7kX84tp\n",
       "Q1+3EBaflXlUO689nhV0yzqfQ7Gdz7bAYyma692Blw7sDDBMNl8D/ARYDLw7k4cj2B/4e+AI4DmZ\n",
       "/L6df7d+EsHbgXcAhzdk8bBJ5YqNaQ2UPzD+i+IHxutq9uzKiMpnP+cDHwMuBL5MsUH1na2mO8Xc\n",
       "WMgunLxR03EmN7CMEwaa0wieBLyGogE7MpNLO/836T3lD7O5FMEwn+KO8mLgl1X8Gyp/gfsg8Grg\n",
       "emAFsAz4E/DVTO4dwzm2Zpen3MAxv9l0mvjXWZzXtzfY+0n5c+Y/KJ6Re1MmV1Rc0qSZK5PXT+/h\n",
       "eLfCKL9nTqe4I/TGHsvmAykGBj9BscL/l4A/AHdMdEHF8jn/V1M0vs/L5PI2ldsYr4xY8h04bOjx\n",
       "F7Pt5Wdy94Ht+jc0sILu4MWKjoUbLoQThv6bjmBL4CTg5RTb7a0AlgI3A1/OZOVo1yt3ZPguRUN7\n",
       "M8Vg86nAZ5vSTFWl/DnzMeAZFHsiX1lxSZM2mVzp6WfW6qz8h7Ygk6vG8PJ/AhYAT+ul4AMof7m9\n",
       "IoLLgBdRTCXaAdgygguAhyhWMPwzsD277vkyXnT9Rqu+8SLmc9rTvxnBN4EdgadTzL03+EZQPjwP\n",
       "8KEI9qVo5j8D7BbBDyiakGsyeWjo147jrvXjgOcDB1IsgHA6xejwFIr9dWcC9wBPAF5BMbXrqZnc\n",
       "NsG/04rY9zdX0yLYmdqsxSa6LZOM4D3AbRRbFDnoo74xzC/yeyyIYITm9H0Uz/w/swez+TLgsgh+\n",
       "B7yYIg92AGZG8GtgPRuyeUeKnQDWAbdS/J0Po1iU707gYuBRFPuZf4riDlnP//LcCato3fT/lSfs\n",
       "CayK4LsUd7N/3yqbx+pRzHn/ady7yQq6h7Pbv0WwO8UeuyvYkM2/psjmv07kepncU26d8ziKAY+L\n",
       "mjLzpmplNr+LYhGzcyP6+3dfG9MOKPcNPAV4VgSXUrzPPwKuAH4y8MOo3HfssxS/9D8vkwcqKnnS\n",
       "MvkTxYqFJwNEsB9FqD1IMVVkN+BBNp/aeoGbWdfcRDGaewnF1N1R765pg0yuAT4QwQeB7YDjgO8B\n",
       "O0WwguLf3+cyubzlXesfbrZnbPWGBaz8r1UUz5h8gKIB3R34McVeecspfnDeQxG+qynuJGxD8W/7\n",
       "tZn8atJ/mbXDjOavrdeCI72o/OX6ExH8BTgjgidlcl/VdUmdNg8WtdoK45k87hMRTAHOGpTN21As\n",
       "6LIAeH7dFjsaj3KK5sfL/1FOwdyH4mf3bhQ/49cA+1EMJD+fIgO+STHQuBPFWge/pZhp4c+LESyF\n",
       "U46FPQb/WzsGblzJBSdQvIfHAz8Adozg3vLPn2/V6EewE0XOHkDx/9lNFNuMTH0W+z6x6DU3NpWt\n",
       "HkXx/9d6isHjK4GjM1u8eJzK/Li2/J/aqHxvPzoom5+cyaqq66qCU3nbqNxD7K0Uc8W/SXH38CCK\n",
       "H/rPA54G7EKx3cUNFMuPn0vxXGnPBt94xJ6xmNdzxCafcJpmR5SDHztR/Ls8GtiBnZ+5nmPPn7XJ\n",
       "i087fDW3/uJ0YG+KpvY64Nxu79U3zHTvG1n2/9q7/1C76zqO48/X5swfi5YIc7WBEjrMoIlDh784\n",
       "RppbokUgDjQDjaJ0ktAi/0hBLBukpYJkrkgnylImE+e0mpoE/qotf2yWKwfdtK2gsjJyY+/++HxO\n",
       "53zvPefcu7tz7vnxeT3gcs/O/Z677/fNe+e1z/l+v58Pq1qd1bXpkbgLOJK0zMNQnQ2qKyFXeq2E\n",
       "Gkp86HxO2PIYr0+YMO58Fo89zmtvks4I3kW6TPF60rrdq/u9jrgNn1b3aba4vPb9pGy+ErgYOJo0\n",
       "8NxL+oBApA8N3gO8BTxNyuYHgd9+gtlf2cz+c8f/3b24l9VmlsQPAUUM/GRrbfke0z7LM12eSTrL\n",
       "9AywptX9W/ny3pNI64EuAG4tbeIADzr6J/ffXBa/bwsr3146YYP1PB2vRm3Gd6yFfKnx1RzK4bzL\n",
       "f3ib290f3ZXvCX4BuKXbk6PMlFHOlZkyyjWU+AhwOvCN0zll1i/59YLx26yAzY8RK0iXPK4mXdZ6\n",
       "yzBMQmij4f/ZnC6RnZ2/IN0POgf4awT7m1/T6tL0K+H3z8KqTvdN2+DL9wT/CvhmxHDOeN3Tgamk\n",
       "82hy63YAAAdpSURBVEk3z88G7o6Ib7fY5jZgOfAO8LmI2NrNnRxk+Q1lC+lG8JsiuLPPuzTwPOjo\n",
       "L5+1tjqJE4GnSJOuDfyMo+ONaq5MhbO5s5zNvyCddbrxJDTm/8jbKJnKmVkbTvlDtS2ky7Cf7Pf+\n",
       "HKieDUwlzSbd9/dx0g3xLwArI2JH0zYrgKsiYoWk04DvRcSybu7koMqXYjwAzCNNXLRv6q9VLSKe\n",
       "6tW+DRPXoqrX9Rims9bujape1EPiLGADcCNw2zBd1juKuTIVzubO8gyi60kTtp1dv3d0apdY+j2n\n",
       "zrWocj0aXIuqHmXzOcBPgBsiuKObv7vXDiZXZk3y81OBnRGxKyL2kgZhF43b5kLSAstExHPAPEnz\n",
       "p7MzQ+iTpJvAD2hQmtW6vztDq9bvHRgwtV7+8tgdmxjjGtaxmfU8zTo2D+KgNKv1ewcGTK3bvzCC\n",
       "Z0jv9VcAj0pc2GYhdRsczubOLgK20jQoBXglYtOmiOXrI2qbIpa3ObtUm6mdHAK1fu/AgKn1ewcG\n",
       "SK3fOzBgat3+hflM6TLgixKPSFxQQjZPNivvB0lrHdWNAadNYZuFpIWdR1oE64B1/d4PswOVB6GD\n",
       "OBC1PojgD3mNwkuBm4DbJZ4CXia9v/8R2EOaZXtsmM6qjihncwcR/Kjf+2BmdrAi2CmxFLiMtOTT\n",
       "HTmbX2FEs3mygelUD3D8CH7oC2NmVpI8+/LaPCPgCaRPgBeTBjyLyGsgAvMk3iUtR7AX+DfwTxqz\n",
       "STZ/7W96/N8IPjWDhzTKnM1mZgXIM4P/QOJuUibX8vdlpGw+mmo27yXl81Sz+Z0IPjODh9TRZAPT\n",
       "P5EOum4R6VPXTtsszM9NIMmh2ETS9f3eh0HhWlS5Hg2uRdWA1OOw6bxI8sCoS5zNPTQg/8YGgmtR\n",
       "5Xo0uBZVA1KPoc/myQamLwLHSzoWeJO0DuLKcdtsBK4CHpC0DPh7REy4VGjUJlcwMzPrE2ezmZmN\n",
       "nI4D04jYJ+kq4HHSlPRrI2KHpC/kn38/IjZJWiFpJ+m08dAuCGtmZjbonM1mZjaKJl3H1MzMzMzM\n",
       "zKyXJlsuZkokLZL0pKRXJb0iaVV+/gZJY5K25q/lTa/5uqTXJb0m6bxu7MegkHSYpOckbZO0XdK3\n",
       "8vNHSfqppN9JekLSvKbXjGQ9OtSiyN6okzQ7H/cj+c/F9UazFvUosj8k7ZL0Uj7m5/NzxfZGm3oU\n",
       "2RvT4WyucjY3OJtbczZXOZsTZ3NVT7M5Ig76CzgGWJIfzyUt/H0icD1wbYvtPwxsA+YAxwI7gVnd\n",
       "2JdB+QKOyN8PAZ4FzgTWAKvz818Dbi6hHm1qUWxv5OO8FrgP2Jj/XGRvdKhHkf0BvAEcNe65Ynuj\n",
       "TT2K7I1p1s/ZPPEYnc2da1Fsb+TjdDZ3rkeR/eFsnlI9utIbXTljGhF/joht+fG/gB2kNdRg4nT1\n",
       "kBa/vj8i9kbErryTp3ZjXwZFRLyTHx5KugfobzQteJ6/15dOGOl6tKkFFNobkhYCK4C7adSgyN6A\n",
       "tvUQhfYHE4+72N7IWvVBqb1xQJzNEzmbG5zNVc7mKmfzBM7mqp5kc1cGppU9SrMEnkz69A3gakm/\n",
       "kbS26TT3B6hObT9GIyxHgqRZkraRFjN/MiJeBeZHY1bE3cD8/Hik69GmFlBobwC3Al8lrSNVV2Rv\n",
       "ZK3qEZTZHwH8TNKLkj6fnyu5N1rVA8rsjYPibE6czQ3O5gmczVXO5gZnc1XPsrmrA1NJc4EHgWvy\n",
       "p7N3AscBS4C3gO90ePlIzcIUEfsjYglp7bizJZ0z7udB52MemXq0qEWNQntD0gXAnojYSutPlorq\n",
       "jQ71KLI/gDMi4mRgOfBlSWc1/7Ck3sha1aPU3pg2Z3ODs7nB2dzgbK5yNk/gbK7qWTZ3bWAqaQ7w\n",
       "ELAuIh4GiIg9kZEuBaifup3ywt/DLiL+ATwKnALslnQMgKQFwJ68WRH1aKrF0oJ743TgQklvAPcD\n",
       "H5N0L+X2Rqt63FNqf0TEW/n7X4ANpOMutTda1qPU3pguZ3NrzuYGZzPgbB7P2dzE2VzVy2zu1qy8\n",
       "AtYC2yPiu03PL2ja7NPAy/nxRuASSYdKOg44Hni+G/syCCQdXT+FLelw4FxgK+m4L8+bXQ48nB+P\n",
       "bD3a1aL+jzkrpjci4rqIWBQRxwGXAFsi4jIK7A1oW4/PlvjeIekISe/Nj48EziMdd5G90a4epb53\n",
       "TIezucrZ3OBsrnI2VzmbG5zNVb3O5kO6tJ9nAJcCL0namp+7DlgpaQnplO0bQH3x7+2S1gPbgX3A\n",
       "l/IIe1QsAH4saRZp8H9vRPw812a9pCuAXcDFMPL1aFeLewrtjfHqx3Yz5fXGeKJRjzWSPkpZ/TEf\n",
       "2JDGEhwC3BcRT0h6kTJ7o109/N4xdc7mKmdzg7O5M2dzg7PZ2dysp9ms0aqVmZmZmZmZDZuuz8pr\n",
       "ZmZmZmZmdiA8MDUzMzMzM7O+8sDUzMzMzMzM+soDUzMzMzMzM+srD0zNzMzMzMysrzwwNTMzMzMz\n",
       "s77ywNTMzMzMzMz6ygNTMzMzMzMz66v/AcKYoCqaO7LnAAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106ec51d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "path=\"/Users/stanpol/Python/\"\n",
    "fig = plt.figure(figsize=[16,6])\n",
    "ax = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "\n",
    "with open(path+\"sample4nf.dat\") as f:\n",
    "\tdata = np.loadtxt(f, unpack=True, skiprows=2)\n",
    "data[-1] = data[-1]/data[-1].max()\n",
    "a1, = ax.plot(data[0], data[-1], 'b')\n",
    "a2, = ax2.plot(data[0], data[-1], 'b')\n",
    "ax.set_ylim([0,1.2])\n",
    "ax2.set_ylim([0,1.2])\n",
    "\n",
    "_max, _min = peakdetect(data[-1], data[0], 1, 0.01)\n",
    "xm = [p[0] for p in _max]\n",
    "ym = [p[1] for p in _max]\n",
    "xn = [p[0] for p in _min]\n",
    "yn = [p[1] for p in _min]\n",
    "print xm\n",
    "ax.plot(xm, ym, 'ro')\n",
    "ax.plot(xn, yn, 'go')\n",
    "\n",
    "peakind = signal.find_peaks_cwt(data[-1], np.arange(1,25)) \n",
    "ax2.plot(data[0][peakind], data[-1][peakind], 'ro')\n",
    "print data[0][peakind]\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}