SpikingNeurons.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:739d6da00ffd7048e2cb308701f5eb7b7b50fa393903ad0c955435fe4ac44f8e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Sugandha Sharma (ID: 20401693)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      " 1) Generating a random input signal"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1.1) Gaussian white noise"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "a), b)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "#import matplotlib.pyplot as plt\n",
      " \n",
      "\n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample\n",
      "        \n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) and X(w) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    # plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\n",
      "    \n",
      "    return x, X\n",
      "    \n",
      "   \n",
      "\n",
      "# Same as generate_signal but without plots\n",
      "def gen_sig(T, dt, rms, limit, seed):\n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))   # Fourier Transform in switch to frequency domain\n",
      "    \n",
      "    return x, X\n",
      " \n",
      "\n",
      "# calculate and plot the average X omega value for 100 signals   \n",
      "def plot_avgXomega(T, dt, rms, limit): \n",
      "    seed = 3\n",
      "    sum_X = 0\n",
      "    count = 100             # number of signals to compute the average\n",
      "    for num in xrange(count):\n",
      "        seed = seed + 1\n",
      "        x, X = gen_sig(T, dt, rms, limit, seed)\n",
      "        sum_X += np.abs(X)\n",
      "    \n",
      "    avg_X = sum_X / count\n",
      "    freq = np.arange(int(T/dt)) - (T/dt)/2      # spread the frequencies evenly from -500 to 500 (enforcing symmetry to get a real signal)\n",
      "    omega = 2*numpy.pi*freq\n",
      "\n",
      "\n",
      "    figure()\n",
      "    plot(omega, avg_X)\n",
      "    xlabel('$\\omega$')\n",
      "    ylabel('$|X(\\omega)|$')\n",
      "    #xlim(-100,100)\n",
      "\n",
      "    \n",
      "def main():\n",
      "    \n",
      "    # Generating plots of x(t)\n",
      "    T = 1\n",
      "    dt = 0.001\n",
      "    seed = 3\n",
      "    rms = 0.5\n",
      "    generate_signal(T, dt, rms, 5, seed)   \n",
      "    generate_signal(T, dt, rms, 10, seed)\n",
      "    generate_signal(T, dt, rms, 20, seed)\n",
      "    \n",
      "    \n",
      "    # Generate plot of average X(omega)\n",
      "    plot_avgXomega(T, dt, rms, 10)\n",
      "    \n",
      "   \n",
      "            \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1.2) Gaussian power spectrum noise"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "a) "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "#import matplotlib.pyplot as plt\n",
      " \n",
      "\n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, bandwidth, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.bandwidth = bandwidth\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"    \n",
      "        \n",
      "        # num of points(coefficients) is halved since  rfftfreq adn irfft methods used later return n+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "        \n",
      "        # generate the matrix containing the points \n",
      "        y_values = np.zeros((1, num_points + 1))   \n",
      "        x_values = np.zeros((1, num_points + 1))  \n",
      "             \n",
      "        # Create sample frequencies and compute the coefficients depending upon omega from a gaussian distribution\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "        sigmas = np.exp( (-1 * sample_freqs**2) / (2 * self.bandwidth**2) )\n",
      "        \n",
      "        i = 0\n",
      "        for count in xrange(num_points + 1):\n",
      "            if sigmas[i] > 0 :\n",
      "                 std_dev = sigmas[i] \n",
      "            else:\n",
      "                std_dev = 0.0001\n",
      "            y_values[:,i] = rng.normal(0.0, std_dev, 1 )\n",
      "            x_values[:,i] = rng.normal(0.0, std_dev, 1 )\n",
      "            i = i + 1\n",
      "        \n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                            # set the DC component to zero \n",
      "        \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)    # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "       \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample\n",
      "\n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    # plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\n",
      "    \n",
      "    return x, X\n",
      "    \n",
      "   \n",
      "\n",
      "# Same as generate_signal but without plots\n",
      "def gen_sig(T, dt, rms, limit, seed):\n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))   # Fourier Transform in switch to frequency domain\n",
      "    \n",
      "    return x, X\n",
      " \n",
      "\n",
      "# calculate and plot the average X omega value for 100 signals   \n",
      "def plot_avgXomega(T, dt, rms, limit): \n",
      "    seed = 3\n",
      "    sum_X = 0\n",
      "    count = 100             # number of signals to compute the average\n",
      "    for num in xrange(count):\n",
      "        seed = seed + 1\n",
      "        x, X = gen_sig(T, dt, rms, limit, seed)\n",
      "        sum_X += np.abs(X)\n",
      "    \n",
      "    avg_X = sum_X / count\n",
      "    freq = np.arange(int(T/dt)) - (T/dt)/2      # spread the frequencies evenly from -500 to 500 (enforcing symmetry to get a real signal)\n",
      "    omega = 2*numpy.pi*freq\n",
      "\n",
      "\n",
      "    figure()\n",
      "    plot(omega, avg_X)\n",
      "    xlabel('$\\omega$')\n",
      "    ylabel('$|X(\\omega)|$')\n",
      "    #xlim(-500,500)\n",
      "\n",
      "\n",
      "    \n",
      "def main():\n",
      "    \n",
      "    # Generating plots of x(t)\n",
      "    T = 1\n",
      "    dt = 0.001\n",
      "    seed = 3\n",
      "    rms = 0.5\n",
      "    generate_signal(T, dt, rms, 5, seed)   \n",
      "    generate_signal(T, dt, rms, 10, seed)\n",
      "    generate_signal(T, dt, rms, 20, seed)\n",
      "    \n",
      "    # Generate plot of average X(omega)\n",
      "    plot_avgXomega(T, dt, rms, 10)\n",
      "        \n",
      "            \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "2) Simulating a Spiking Neuron"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "a), b) "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Spiky Leaky Integrate and Fire Neuron Model\n",
      "\n",
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "    \n",
      "# subclass of LeakyIntegNeuron class    \n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        try:\n",
      "            bias = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 40.0)) / self.tau_rc))\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')    \n",
      "        return bias\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        try:\n",
      "            x = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 150.0)) / self.tau_rc))\n",
      "            gain = x - bias\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')\n",
      "        return gain\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        spike_time = list()\n",
      "        \n",
      "        num_spikes = 0\n",
      "        test = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time):  \n",
      "            if t > t_rest:\n",
      "                V[i] = V[i-1] + ((J*R - V[i-1]) / self.tau_rc) * dt\n",
      "                test = test + 1\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                spike_time.append(t)\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0\n",
      "\n",
      "        print \"Number of Spikes = \" , num_spikes  #, test\n",
      "        ## plot membrane potential trace\n",
      "        figure()\n",
      "        vlines(spike_time, 0, 1.5, linewidth=1)\n",
      "        #plot(time, V)\n",
      "        title('Leaky Integrate and Fire Neuron')\n",
      "        ylabel('Membrane Potential (V)')\n",
      "        xlabel('Time (sec)')\n",
      "        ylim([0,2])\n",
      "        show()\n",
      "    \n",
      "\n",
      "def main():\n",
      "    \n",
      "    neuron = spikyLIF(0.02, 0.002) \n",
      "    encoder = 1\n",
      "    bias = neuron.cal_bias()\n",
      "    gain = neuron.cal_gain(bias)\n",
      "    print \"gain = \" , gain\n",
      "    print \"Jbias = \", bias\n",
      "    \n",
      "    # Demo of spike generation\n",
      "    neuron2 = spikyLIF(10, 4)\n",
      "    neuron2.generate_spikes(T=50, dt=0.001, J=1.4, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # 2a) spike output for a constant input of x=0 \n",
      "    x = 0\n",
      "    J = neuron.cal_current(x, gain, encoder, bias)\n",
      "    neuron.generate_spikes(T=1, dt=0.001, J=J, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # 2b) spike output for a constant input of x=1 \n",
      "    x = 1\n",
      "    J = neuron.cal_current(x, gain, encoder, bias)\n",
      "    neuron.generate_spikes(T=1, dt=0.001, J=J, Vth=1, V_spike=0.5, R=1)\n",
      "   \n",
      "\n",
      "    \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "c), d)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "import matplotlib.patches as mpatches       # for legend on the graphs\n",
      " \n",
      "# subclass of LeakyIntegNeuron class    \n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        try:\n",
      "            bias = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 40.0)) / self.tau_rc))\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')    \n",
      "        return bias\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        try:\n",
      "            x = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 150.0)) / self.tau_rc))\n",
      "            gain = x - bias\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')\n",
      "        return gain\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over timet\n",
      "        spike_time = list()\n",
      " \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                V[i] = V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                spike_time.append(t)\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0\n",
      "                \n",
      "\n",
      "        #figure()\n",
      "        #vlines(spike_time, 0, 1.5, linewidth=3)\n",
      "        print \"Number of Spikes = \" , num_spikes\n",
      "        ## plot membrane potential trace\n",
      "        #figure()\n",
      "        plot(time, V, 'g')\n",
      "        #vlines(spike_time, 0, 1.5, linewidth=1)\n",
      "        title('Leaky Integrate and Fire Neuron')\n",
      "        ylabel('Membrane Potential (V)')\n",
      "        xlabel('Time (sec)')\n",
      "        #xlim(0,0.2)\n",
      "        show()\n",
      "        \n",
      "\n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    # plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\n",
      "    #xlim(0,0.2)\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "  \n",
      "def main():\n",
      "    \n",
      "    # spike output for a time varying signal as an input  \n",
      "    print \"test\"\n",
      "    # parameters for x(t)\n",
      "    T = 1\n",
      "    dt = 0.001\n",
      "    seed = 3\n",
      "    rms = 0.5\n",
      "    limit = 30\n",
      "       \n",
      "    # Create an LIF neuron\n",
      "    neuron = spikyLIF(0.02, 0.002) \n",
      "    encoder = 1.0\n",
      "    bias = neuron.cal_bias()\n",
      "    gain = neuron.cal_gain(bias)\n",
      "    \n",
      "    # 2c plot x(t) and spikes on the same plot\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "    J = neuron.cal_current(x, gain, encoder, bias)\n",
      "    neuron.generate_spikes(T=1, dt=0.001, J=J, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    \n",
      "    # 2d plot voltage and spikes for 0.2s\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed) \n",
      "    J = neuron.cal_current(x, gain, encoder, bias)\n",
      "    neuron.generate_spikes(T=0.2, dt=0.001, J=J, Vth=1, V_spike=0.5, R=1)\n",
      "\n",
      "        \n",
      "            \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Bonus Question"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# TESTTING BONUS\n",
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "import matplotlib.patches as mpatches       # for legend on the graphs\n",
      " \n",
      "# subclass of LeakyIntegNeuron class    \n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        try:\n",
      "            bias = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 40.0)) / self.tau_rc))\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')    \n",
      "        return bias\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        try:\n",
      "            x = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 150.0)) / self.tau_rc))\n",
      "            gain = x - bias\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')\n",
      "        return gain\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        spike_time = list()\n",
      " \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                V[i] = J[t*1000]*R*(1-exp(-t/self.tau_rc))      #V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                spike_time.append(t)\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0\n",
      "                \n",
      "\n",
      "        print \"Number of Spikes = \" , num_spikes\n",
      "        ## plot membrane potential trace\n",
      "        #figure()\n",
      "        #plot(time, V, 'g')\n",
      "        vlines(spike_time, 0, 1.5, linewidth=1)\n",
      "        title('Leaky Integrate and Fire Neuron')\n",
      "        ylabel('Membrane Potential (V)')\n",
      "        xlabel('Time (sec)')\n",
      "        #xlim(0,0.2)\n",
      "        show()\n",
      "        \n",
      "\n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    # plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "  \n",
      "def main():\n",
      "    \n",
      "    # spike output for a time varying signal as an input  \n",
      "    print \"test\"\n",
      "    # parameters for x(t)\n",
      "    T = 1\n",
      "    dt = 0.001\n",
      "    seed = 3\n",
      "    rms = 0.5\n",
      "    limit = 30\n",
      "       \n",
      "    # Create an LIF neuron\n",
      "    neuron = spikyLIF(0.02, 0.002) \n",
      "    encoder = 1.0\n",
      "    bias = neuron.cal_bias()\n",
      "    gain = neuron.cal_gain(bias)\n",
      "    \n",
      "    # 2c plot x(t) and spikes on the same plot\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "    J = neuron.cal_current(x, gain, encoder, bias)\n",
      "    neuron.generate_spikes(T=1, dt=0.001, J=J, Vth=1, V_spike=0.5, R=1)\n",
      "\n",
      "        \n",
      "            \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "3) Simulating Two Spiking Neurons"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "a), b), c), d)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "import matplotlib.patches as mpatches       # for legend on the graphs\n",
      "    \n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        try:\n",
      "            bias = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 40.0)) / self.tau_rc))\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')    \n",
      "        return bias\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        try:\n",
      "            x = 1.0 / (1 - np.exp((self.tau_ref - (1.0 / 150.0)) / self.tau_rc))\n",
      "            gain = x - bias\n",
      "        finally:\n",
      "            np.seterr(divide='ignore')\n",
      "        return gain\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1, flag=0):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        spike_time = list()\n",
      "        \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                if numpy.isscalar(J):\n",
      "                    V[i] = V[i-1] + ((J*R - V[i-1]) / self.tau_rc) * dt\n",
      "                else:\n",
      "                    V[i] = V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                spike_time.append(t)\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0\n",
      "\n",
      "        print \"Number of Spikes = \" , num_spikes\n",
      "        ## plot membrane potential trace\n",
      "        #figure()\n",
      "        if flag == 1:\n",
      "            #plot(time, V-1.5, 'g')\n",
      "            vlines(spike_time, -1.5, 0, linewidth=1)\n",
      "        else:    \n",
      "            #plot(time, V, 'g')    \n",
      "            vlines(spike_time, 0, 1.5, linewidth=1)    \n",
      "        title('Leaky Integrate and Fire Neuron')\n",
      "        ylabel('Membrane Potential (V)')\n",
      "        xlabel('Time (sec)')\n",
      "        #ylim([-0.5,2])\n",
      "        #show()\n",
      "  \n",
      "        \n",
      "\n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    # plot the time domain signal\n",
      "    #figure()\n",
      "    plot(samples_array*dt, x, 'r')          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "  \n",
      "def main():\n",
      "     \n",
      "       \n",
      "    # Create an LIF neuron (1)\n",
      "    neuron1 = spikyLIF(0.02, 0.002) \n",
      "    encoder1 = 1.0\n",
      "    bias1 = neuron1.cal_bias()\n",
      "    gain1 = neuron1.cal_gain(bias1)\n",
      "    \n",
      "    \n",
      "    # Create another LIF neuron (2)\n",
      "    neuron2 = spikyLIF(0.02, 0.002) \n",
      "    encoder2 = -1.0\n",
      "    bias2 = neuron2.cal_bias()\n",
      "    gain2 = neuron2.cal_gain(bias2)\n",
      "    \n",
      "    time = np.arange(0, 1, 0.001) \n",
      "    x_values = np.zeros(time.shape)\n",
      "   \n",
      "    # 3a Plot x(t) and the spiking output for x(t)=0 \n",
      "    x = 0\n",
      "    figure()\n",
      "    plot(time, x_values, 'r')\n",
      "    J1 = neuron1.cal_current(x, gain1, encoder1, bias1)\n",
      "    neuron1.generate_spikes(T=1, dt=0.001, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    figure()\n",
      "    plot(time, x_values, 'r')\n",
      "    J2 = neuron2.cal_current(x, gain2, encoder2, bias2)\n",
      "    neuron2.generate_spikes(T=1, dt=0.001, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "\n",
      "    \n",
      "    # 3b Plot x(t) and the spiking output for x(t)=1 \n",
      "    x_values = np.ones(time.shape)\n",
      "    x = 1\n",
      "    figure()\n",
      "    plot(time, x_values, 'r')\n",
      "    J1 = neuron1.cal_current(x, gain1, encoder1, bias1)\n",
      "    neuron1.generate_spikes(T=1, dt=0.001, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    figure()\n",
      "    plot(time, x_values, 'r')\n",
      "    J2 = neuron2.cal_current(x, gain2, encoder2, bias2)\n",
      "    neuron2.generate_spikes(T=1, dt=0.001, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # 3c Plot x(t) and the spiking output for x(t)=1/2sin(10\u03c0t)\n",
      "    x_values = 0.5 * np.sin(10*np.pi*time) \n",
      "    figure()\n",
      "    plot(time, x_values, 'r')\n",
      "    J1 = neuron1.cal_current(x_values, gain1, encoder1, bias1)\n",
      "    neuron1.generate_spikes(T=1, dt=0.001, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    #figure()\n",
      "    #plot(time, x_values, 'r')\n",
      "    J2 = neuron2.cal_current(x_values, gain2, encoder2, bias2)\n",
      "    neuron2.generate_spikes(T=1, dt=0.001, J=J2, Vth=1, V_spike=0.5, R=1, flag=1)\n",
      "    \n",
      "    # 3d Plot x(t) and the spiking output for a random signal\n",
      "    \n",
      "    # parameters for x(t)\n",
      "    T = 2\n",
      "    dt = 0.001\n",
      "    seed = 3\n",
      "    rms = 0.5\n",
      "    limit = 5\n",
      "    \n",
      "    figure()\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "    J1 = neuron1.cal_current(x, gain1, encoder1, bias1)\n",
      "    neuron1.generate_spikes(T=2, dt=0.001, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "    J2 = neuron2.cal_current(x, gain2, encoder2, bias2)\n",
      "    neuron2.generate_spikes(T=2, dt=0.001, J=J2, Vth=1, V_spike=0.5, R=1, flag=1)\n",
      "    \n",
      "    \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "4) Computing an Optimal Filter"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "a)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Optimal Filter for decoding pairs of spikes\n",
      "import numpy\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "\n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        pass\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        pass\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        S = zeros(len(time))\n",
      " \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                if numpy.isscalar(J):\n",
      "                    V[i] = V[i-1] + ((J*R - V[i-1]) / self.tau_rc) * dt\n",
      "                else:\n",
      "                    V[i] = V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                S[i] = 1.5\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0  \n",
      "         \n",
      "        \n",
      "        #print \"Number of Spikes = \" , num_spikes\n",
      "        return S\n",
      "    \n",
      "        \n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    \"\"\"# plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\"\"\"\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "\n",
      "\n",
      "def main():\n",
      "    T = 2.0         # length of signal in seconds\n",
      "    dt = 0.001      # time step size\n",
      "    \n",
      "    # parameters for x(t)\n",
      "    seed = 3     # the random number seed to use (so we can regenerate the same signal again)\n",
      "    rms = 0.5    # the root mean square power level of the signa\n",
      "    limit = 5    # the maximum frequency for the signal x(t) (in Hz)\n",
      "\n",
      "    # Generate bandlimited white noise (use your own function from part 1.1)\n",
      "    #x, X = syde556.generate_signal(T, dt, rms=0.5, limit=5, seed=3)\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "\n",
      "    Nt = len(x)                # Count the number of sample points in the time varying signal\n",
      "    t = numpy.arange(Nt) * dt  # Generate evenly spaced (spaced by 1) Nt values between the interval  [0,Nt) and then multipy  \n",
      "                               # by dt to change the spacing between two adjacent points to dt\n",
      "\n",
      "    # Neuron parameters\n",
      "    tau_ref = 0.002          # Absolute Refractory period of the neuron\n",
      "    tau_rc = 0.02            # The RC time constant of the LIF neuron \n",
      "    x0 = 0.0                 # firing rate at x=x0 is a0\n",
      "    a0 = 40.0\n",
      "    x1 = 1.0                 # firing rate at x=x1 is a1\n",
      "    a1 = 150.0\n",
      "\n",
      "    # Calculate the gain, bias and the threshold of the LIF neuron\n",
      "    eps = tau_rc/tau_ref\n",
      "    r1 = 1.0 / (tau_ref * a0)\n",
      "    r2 = 1.0 / (tau_ref * a1)\n",
      "    f1 = (r1 - 1) / eps\n",
      "    f2 = (r2 - 1) / eps\n",
      "    alpha = (1.0/(numpy.exp(f2)-1) - 1.0/(numpy.exp(f1)-1))/(x1-x0) \n",
      "    x_threshold = x0-1/(alpha*(numpy.exp(f1)-1))              \n",
      "    Jbias = 1-alpha*x_threshold;   \n",
      "\n",
      "    # Simulate the two neurons (use your own function from part 3)\n",
      "    #spikes = syde556.two_neurons(x, dt, alpha, Jbias, tau_rc, tau_ref)\n",
      "    \n",
      "\n",
      "    # Create an LIF neuron (1)\n",
      "    neuron1 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder1 = 1.0\n",
      "    \n",
      "    # Create another LIF neuron (2)\n",
      "    neuron2 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder2 = -1.0\n",
      "    \n",
      "    \n",
      "    # generate the spikes for the first neuron\n",
      "    J1 = neuron1.cal_current(x, alpha, encoder1, Jbias)\n",
      "    V1 = neuron1.generate_spikes(T, dt, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # generate the spikes for the second neuron\n",
      "    J2 = neuron2.cal_current(x, alpha, encoder2, Jbias)\n",
      "    V2 = neuron2.generate_spikes(T, dt, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # create and array of spikes of both neurons\n",
      "    spikes = np.array([V1, V2])\n",
      "    \n",
      "    freq = numpy.arange(Nt)/T - Nt/(2.0*T)   # spread the frequencies evenly from -500 to 500 (enforcing symmetry to get a real signal)\n",
      "    omega = freq*2*numpy.pi                  # Calculate the corresponding angular frequencies in randians\n",
      "\n",
      "    r = spikes[0] - spikes[1]                # Compute the resultant response of the two balanced neurons combined together\n",
      "    R = numpy.fft.fftshift(numpy.fft.fft(r)) # Compute the discrete fourier transform of the response and \n",
      "                                             # then shift the zero frequency component to the centre of the spectrum  \n",
      "\n",
      "    \n",
      "    sigma_t = 0.025                          # The filter should decay towards zero. \n",
      "                                             # A spike shouldn't affect values far away from when it happens\n",
      "    W2 = numpy.exp(-omega**2*sigma_t**2)     # Therefore, create a gaussian window to accomplish this\n",
      "    W2 = W2 / sum(W2)                        # Scale the window appropriately\n",
      "\n",
      "                                          # Computing H(\u03c9)= [(X(\u03c9)R(\u03c9)*)\u2217W(\u03c9)] / [|R(\u03c9)|^2\u2217W(\u03c9)]\n",
      "    CP = X*R.conjugate()                  # Compute  X(\u03c9)R(\u03c9)*\n",
      "    WCP = numpy.convolve(CP, W2, 'same')  # (X(\u03c9)R(\u03c9)*)\u2217 W(\u03c9) Convolve with the gausian window\n",
      "    RP = R*R.conjugate()                  # Compute |R(\u03c9)|^2 -> Square of magnitude of R\n",
      "    WRP = numpy.convolve(RP, W2, 'same')  # Compute the convolution with the gausian window |R(\u03c9)|^2\u2217W(\u03c9)\n",
      "   \n",
      "    XP = X*X.conjugate()                  # Compute |X(\u03c9)|^2 -> Square of magnitude of X\n",
      "    WXP = numpy.convolve(XP, W2, 'same')  # Compute the convolution with the gaussian window |X(\u03c9)|^2\u2217W(\u03c9)\n",
      "\n",
      "    H = WCP / WRP                         # Compute the filter H(\u03c9)= [(X(\u03c9)R(\u03c9)*)\u2217W(\u03c9)] / [|R(\u03c9)|^2\u2217W(\u03c9)] \n",
      "\n",
      "    h = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(H))).real     # Compute the filter in time domain\n",
      "                                                                            # and shift the zero frequency component to the centre\n",
      "                                                                            # Only get the real part of the filter\n",
      "\n",
      "    XHAT = H*R                            # Compute a decoded estimation of x in frequency domain ( XHAT(\u03c9)=R(\u03c9)H(\u03c9) )\n",
      "    \n",
      "\n",
      "    xhat = numpy.fft.ifft(numpy.fft.ifftshift(XHAT)).real  # Compute xhat(t) by taking the inverse fourier transform of XHAT(\u03c9)\n",
      "\n",
      "\n",
      "    import pylab\n",
      "\n",
      "    pylab.figure(0)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.sqrt(XP), label='$|X(\\omega)|$')  # The original signal plotted in frequency domain \n",
      "    pylab.legend()                                           # magnitude of X -> [X(\u03c9)]  versus frequency (Hz) \n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|X(\\omega)|$')\n",
      "    xlim(-50, 50)\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "\n",
      "    figure(1)\n",
      "    pylab.subplot(1,2,2)\n",
      "    pylab.plot(freq, numpy.sqrt(RP), label='$|R(\\omega)|$')  # The response r(t)=a1(t)\u2212a2(t) plotted in frequency domain\n",
      "    pylab.legend()                                           # magnitude of R -> |R(\u03c9)| versus frequency (Hz)\n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|R(\\omega)|$')\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "\n",
      "    pylab.figure(2)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, H.real)            # The Optimal filter plotted in frequency domain \n",
      "    pylab.xlabel('freq (Hz)')           # real part of H versus frequency (Hz)\n",
      "    pylab.title('H (real) versus frequency')\n",
      "    pylab.xlim(-50, 50)\n",
      "\n",
      "    #figure(3)\n",
      "    pylab.subplot(1,2,2)\n",
      "    pylab.plot(t-T/2, h)        # The plot of the optimal filter in time domain (only the real part hence t-T/2)\n",
      "    pylab.title('$h(t)$ versus time')\n",
      "    pylab.xlabel('t')\n",
      "    pylab.xlim(-0.5, 0.5)\n",
      "\n",
      "\n",
      "    pylab.figure(4)\n",
      "    pylab.plot(t, r, color='k', label='spikes', alpha=0.2)                # Plot response versus time [ r(t)=a1(t)\u2212a2(t) ]\n",
      "    pylab.plot(t, x, linewidth=2, label='input signal x(t)')              # Plot the original input signal x(t) versus time\n",
      "    pylab.plot(t, xhat, label='decoded estimate xhat(t)')                     # Plot the decoded estimate xhat versus time\n",
      "    pylab.title('Reconstruction of the signal using an optimal filter')\n",
      "    pylab.legend(loc='best')\n",
      "    pylab.xlabel('time')\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "    \n",
      "    \n",
      "    pylab.figure(5)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.abs(XHAT), label='$|XHAT(\\omega)|$')  # The estimated signal plotted in frequency domain \n",
      "    pylab.legend()                                           # magnitude of XHAT -> [XHAT(\u03c9)]  versus frequency (Hz) \n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|XHAT(\\omega)|$')\n",
      "    xlim(-50, 50)\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "    \n",
      "    pylab.figure(6)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.abs(H))           # The Optimal filter plotted in frequency domain\n",
      "    pylab.xlabel('freq (Hz)')           # magnitude of H -> |H(\u03c9)| versus frequency (Hz)\n",
      "    pylab.title('$|H(\\omega)|$ versus frequency')\n",
      "    pylab.xlim(-50, 50)\n",
      "    \n",
      "    pylab.show()\n",
      "\n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Plotting graphs for comparison"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Optimal Filter for decoding pairs of spikes\n",
      "import numpy\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "import pylab\n",
      "\n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        pass\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        pass\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        S = zeros(len(time))\n",
      " \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                if numpy.isscalar(J):\n",
      "                    V[i] = V[i-1] + ((J*R - V[i-1]) / self.tau_rc) * dt\n",
      "                else:\n",
      "                    V[i] = V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                S[i] = 1.5\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0  \n",
      "         \n",
      "        \n",
      "        #print \"Number of Spikes = \" , num_spikes\n",
      "        return S\n",
      "    \n",
      "        \n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    \"\"\"# plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\"\"\"\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "\n",
      "def comp(duration, lim, colour):\n",
      "    T = duration        # length of signal in seconds\n",
      "    dt = 0.001      # time step size\n",
      "    \n",
      "    # parameters for x(t)\n",
      "    seed = 3     # the random number seed to use (so we can regenerate the same signal again)\n",
      "    rms = 0.5    # the root mean square power level of the signa\n",
      "    limit = lim    # the maximum frequency for the signal x(t) (in Hz)\n",
      "\n",
      "    # Generate bandlimited white noise (use your own function from part 1.1)\n",
      "    #x, X = syde556.generate_signal(T, dt, rms=0.5, limit=5, seed=3)\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "\n",
      "    Nt = len(x)                # Count the number of sample points in the time varying signal\n",
      "    t = numpy.arange(Nt) * dt  # Generate evenly spaced (spaced by 1) Nt values between the interval  [0,Nt) and then multipy  \n",
      "                               # by dt to change the spacing between two adjacent points to dt\n",
      "\n",
      "    # Neuron parameters\n",
      "    tau_ref = 0.002          # Absolute Refractory period of the neuron\n",
      "    tau_rc = 0.02            # The RC time constant of the LIF neuron \n",
      "    x0 = 0.0                 # firing rate at x=x0 is a0\n",
      "    a0 = 40.0\n",
      "    x1 = 1.0                 # firing rate at x=x1 is a1\n",
      "    a1 = 150.0\n",
      "\n",
      "    # Calculate the gain, bias and the threshold of the LIF neuron\n",
      "    eps = tau_rc/tau_ref\n",
      "    r1 = 1.0 / (tau_ref * a0)\n",
      "    r2 = 1.0 / (tau_ref * a1)\n",
      "    f1 = (r1 - 1) / eps\n",
      "    f2 = (r2 - 1) / eps\n",
      "    alpha = (1.0/(numpy.exp(f2)-1) - 1.0/(numpy.exp(f1)-1))/(x1-x0) \n",
      "    x_threshold = x0-1/(alpha*(numpy.exp(f1)-1))              \n",
      "    Jbias = 1-alpha*x_threshold;   \n",
      "\n",
      "    # Simulate the two neurons (use your own function from part 3)\n",
      "    #spikes = syde556.two_neurons(x, dt, alpha, Jbias, tau_rc, tau_ref)\n",
      "    \n",
      "\n",
      "    # Create an LIF neuron (1)\n",
      "    neuron1 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder1 = 1.0\n",
      "    \n",
      "    # Create another LIF neuron (2)\n",
      "    neuron2 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder2 = -1.0\n",
      "    \n",
      "    \n",
      "    # generate the spikes for the first neuron\n",
      "    J1 = neuron1.cal_current(x, alpha, encoder1, Jbias)\n",
      "    V1 = neuron1.generate_spikes(T, dt, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # generate the spikes for the second neuron\n",
      "    J2 = neuron2.cal_current(x, alpha, encoder2, Jbias)\n",
      "    V2 = neuron2.generate_spikes(T, dt, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # create and array of spikes of both neurons\n",
      "    spikes = np.array([V1, V2])\n",
      "    \n",
      "    freq = numpy.arange(Nt)/T - Nt/(2.0*T)   # spread the frequencies evenly from -500 to 500 (enforcing symmetry to get a real signal)\n",
      "    omega = freq*2*numpy.pi                  # Calculate the corresponding angular frequencies in randians\n",
      "\n",
      "    r = spikes[0] - spikes[1]                # Compute the resultant response of the two balanced neurons combined together\n",
      "    R = numpy.fft.fftshift(numpy.fft.fft(r)) # Compute the discrete fourier transform of the response and \n",
      "                                             # then shift the zero frequency component to the centre of the spectrum  \n",
      "\n",
      "    \n",
      "    sigma_t = 0.025                          # The filter should decay towards zero. \n",
      "                                             # A spike shouldn't affect values far away from when it happens\n",
      "    W2 = numpy.exp(-omega**2*sigma_t**2)     # Therefore, create a gaussian window to accomplish this\n",
      "    W2 = W2 / sum(W2)                        # Scale the window appropriately\n",
      "\n",
      "                                          # Computing H(\u03c9)= [(X(\u03c9)R(\u03c9)*)\u2217W(\u03c9)] / [|R(\u03c9)|^2\u2217W(\u03c9)]\n",
      "    CP = X*R.conjugate()                  # Compute  X(\u03c9)R(\u03c9)*\n",
      "    WCP = numpy.convolve(CP, W2, 'same')  # (X(\u03c9)R(\u03c9)*)\u2217 W(\u03c9) Convolve with the gausian window\n",
      "    RP = R*R.conjugate()                  # Compute |R(\u03c9)|^2 -> Square of magnitude of R\n",
      "    WRP = numpy.convolve(RP, W2, 'same')  # Compute the convolution with the gausian window |R(\u03c9)|^2\u2217W(\u03c9)\n",
      "   \n",
      "    XP = X*X.conjugate()                  # Compute |X(\u03c9)|^2 -> Square of magnitude of X\n",
      "    WXP = numpy.convolve(XP, W2, 'same')  # Compute the convolution with the gaussian window |X(\u03c9)|^2\u2217W(\u03c9)\n",
      "\n",
      "    H = WCP / WRP                         # Compute the filter H(\u03c9)= [(X(\u03c9)R(\u03c9)*)\u2217W(\u03c9)] / [|R(\u03c9)|^2\u2217W(\u03c9)] \n",
      "\n",
      "    h = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(H))).real     # Compute the filter in time domain\n",
      "                                                                            # and shift the zero frequency component to the centre\n",
      "                                                                            # Only get the real part of the filter\n",
      "\n",
      "    XHAT = H*R                            # Compute a decoded estimation of x in frequency domain ( XHAT(\u03c9)=R(\u03c9)H(\u03c9) )\n",
      "    \n",
      "\n",
      "    xhat = numpy.fft.ifft(numpy.fft.ifftshift(XHAT)).real  # Compute xhat(t) by taking the inverse fourier transform of XHAT(\u03c9)\n",
      "    \n",
      "\n",
      "    \n",
      "\n",
      "    \"\"\"pylab.figure(0)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.sqrt(XP), label='$|X(\\omega)|$')  # The original signal plotted in frequency domain \n",
      "    pylab.legend()                                           # magnitude of X -> [X(\u03c9)]  versus frequency (Hz) \n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|X(\\omega)|$')\n",
      "    xlim(-50, 50)\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "\n",
      "    figure(1)\n",
      "    pylab.subplot(1,2,2)\n",
      "    pylab.plot(freq, numpy.sqrt(RP), label='$|R(\\omega)|$')  # The response r(t)=a1(t)\u2212a2(t) plotted in frequency domain\n",
      "    pylab.legend()                                           # magnitude of R -> |R(\u03c9)| versus frequency (Hz)\n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|R(\\omega)|$')\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "\n",
      "    pylab.figure(2)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, H.real)            # The Optimal filter plotted in frequency domain \n",
      "    pylab.xlabel('freq (Hz)')           # real part of H versus frequency (Hz)\n",
      "    pylab.title('H (real) versus frequency')\n",
      "    pylab.xlim(-50, 50)\n",
      "\n",
      "    #figure(3)\n",
      "    #pylab.subplot(1,2,2)\n",
      "    pylab.plot(t-T/2, h, colour)        # The plot of the optimal filter in time domain (only the real part hence t-T/2)\n",
      "    pylab.title('$h(t)$ versus time')\n",
      "    pylab.xlabel('t')\n",
      "    pylab.xlim(-0.5, 0.5)\n",
      "\n",
      "\n",
      "    pylab.figure(4)\n",
      "    pylab.plot(t, r, color='k', label='spikes', alpha=0.2)                # Plot response versus time [ r(t)=a1(t)\u2212a2(t) ]\n",
      "    pylab.plot(t, x, linewidth=2, label='input signal x(t)')              # Plot the original input signal x(t) versus time\n",
      "    pylab.plot(t, xhat, label='decoded estimate xhat(t)')                     # Plot the decoded estimate xhat versus time\n",
      "    pylab.title('Reconstruction of the signal using an optimal filter')\n",
      "    pylab.legend(loc='best')\n",
      "    pylab.xlabel('time')\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
      "    \n",
      "    \n",
      "    pylab.figure(5)\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.abs(XHAT), label='$|XHAT(\\omega)|$')  # The estimated signal plotted in frequency domain \n",
      "    pylab.legend()                                           # magnitude of XHAT -> [XHAT(\u03c9)]  versus frequency (Hz) \n",
      "    pylab.xlabel('freq (Hz)')\n",
      "    pylab.ylabel('$|XHAT(\\omega)|$')\n",
      "    xlim(-50, 50)\n",
      "    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\"\"\"\n",
      "    \n",
      "    #pylab.figure()\n",
      "    pylab.subplot(1,2,1)\n",
      "    pylab.plot(freq, numpy.abs(H), colour)           # The Optimal filter plotted in frequency domain \n",
      "    pylab.xlabel('freq (Hz)')           # magnitude of H -> |H(\u03c9)| versus frequency (Hz)\n",
      "    pylab.title('$|H(\\omega)|$ versus frequency')\n",
      "    pylab.xlim(-50, 50)\n",
      "    #pylab.plot(t, x, linewidth=2, label='input signal x(t)')\n",
      "\n",
      "    #pylab.show()\n",
      "    \n",
      "    \n",
      "def main():\n",
      "    \n",
      "    #comp(1.0, 5, 'r')\n",
      "    #comp(4.0, 5, 'g')\n",
      "    #comp(10.0, 5, 'b')\n",
      "    \n",
      "    comp(2.0, 2.0, 'r')\n",
      "    comp(2.0, 10.0, 'g')\n",
      "    comp(2.0, 30.0, 'b')\n",
      "   \n",
      "\n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "5) Using Post-Synaptic Currents as a Filter"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Optimal Filter for decoding pairs of spikes\n",
      "#import syde556\n",
      "import numpy\n",
      "import random\n",
      "import numpy as np\n",
      "import functools\n",
      "\n",
      "class spikyLIF():\n",
      "    \n",
      "    def __init__(self, tau_rc=0.02, tau_ref=0.002):\n",
      "        self.tau_rc = tau_rc\n",
      "        self.tau_ref = tau_ref \n",
      "        \n",
      "    def cal_bias(self):\n",
      "        pass\n",
      "\n",
      "    def cal_gain(self, bias): \n",
      "        pass\n",
      "    \n",
      "    def cal_current(self, x, gain, encoder, bias):\n",
      "        J = (gain * x * encoder) + bias  \n",
      "        return J\n",
      "    \n",
      "    def generate_spikes(self, T, dt, J, Vth=1, V_spike=0.5, R=1):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        ----------\n",
      "        T - total simulation time (sec)    \n",
      "        dt - simulation time step (sec)  \n",
      "        Vth - spike threshold voltage (V)\n",
      "        V_spike - spike delta voltage above Vth (V)\n",
      "        R  - resistance (Ohm)\n",
      "        \"\"\"\n",
      "        time = np.arange(0, T, dt)       # time array\n",
      "        t_rest = 0                      # initial refractory time\n",
      "        \n",
      "        V = zeros(len(time))    # potential (V) trace over time\n",
      "        S = zeros(len(time))\n",
      " \n",
      "        num_spikes = 0\n",
      "        # iterate over each time step\n",
      "        for i, t in enumerate(time): \n",
      "            if t > t_rest:\n",
      "                if numpy.isscalar(J):\n",
      "                    V[i] = V[i-1] + ((J*R - V[i-1]) / self.tau_rc) * dt\n",
      "                else:\n",
      "                    V[i] = V[i-1] + ((J[t*1000]*R - V[i-1]) / self.tau_rc) * dt\n",
      "            if V[i] >= Vth:\n",
      "                num_spikes += 1\n",
      "                V[i] += V_spike\n",
      "                t_rest = t + self.tau_ref\n",
      "                S[i] = 1.5\n",
      "            if V[i] < 0:  \n",
      "                V[i] = 0  \n",
      "         \n",
      "        \n",
      "        #print \"Number of Spikes = \" , num_spikes\n",
      "        return S\n",
      "    \n",
      "        \n",
      "class GaussianWhiteNoise():\n",
      "    \n",
      "    def __init__(self, T, limit, rms):\n",
      "        \"\"\"\n",
      "        Parameters:\n",
      "        -----------\n",
      "        T: the length of the signal in seconds\n",
      "        limit: the maximum frequency for the signal (in Hz)\n",
      "        rms: the root mean square power level of the signal.\n",
      "        \"\"\"\n",
      "        self.T = T\n",
      "        self.limit = limit\n",
      "        self.rms = rms\n",
      "              \n",
      "    \n",
      "    def create_sample(self, dt, rng):\n",
      "        \"\"\"\n",
      "        Generate random fourier coefficients and convert them back to \n",
      "        time domain to get a signal with a limit on its max frequecy\n",
      "        Parameters:\n",
      "        -----------\n",
      "        dt: the time step in seconds\n",
      "        rng: A random number generatorn\n",
      "        \"\"\"\n",
      "                    \n",
      "        # set value for the std. deviation\n",
      "        std_dev = 1\n",
      "        \n",
      "        # num of points(coefficients) is halved since irfft method used later returns n output points if input has n//2+1 points\n",
      "        num_points = int(np.ceil(( self.T / dt) / 2.0 ))\n",
      "      \n",
      "        # generate the matrix containing the points\n",
      "        y_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        x_values = rng.normal(0.0, std_dev, (1, num_points + 1) )\n",
      "        points = x_values + 1j * y_values\n",
      "        points[:, 0] = 0.0                  # set the DC component to zero \n",
      "       \n",
      "             \n",
      "        # Check that the max frequency is below the limit and adjust the power accordingly\n",
      "        sample_freqs = np.fft.rfftfreq(num_points * 2, dt)      # compute discrete fourier transform sample frequencies  \n",
      "                                                                # returns 501 frequency samples (n//2 + 1)       \n",
      "        zero_points = sample_freqs > self.limit\n",
      "        points[:, zero_points] = 0.0                            # set the coefficients to zero for frequencies higher than the limit\n",
      " \n",
      "        a = np.array([0])\n",
      "        signal = np.fft.irfft(points, axis=1)          # compute the inverse of the Discrete Fourier transform to get the real signal\n",
      "        \n",
      "        # compute the rms of the signal\n",
      "        rms_orig = matplotlib.mlab.rms_flat(signal)\n",
      "        \n",
      "        # scale the signal so that it has the rms specified\n",
      "        signal = signal * self.rms / rms_orig\n",
      "          \n",
      "        return functools.partial(self.sample_func, a, signal)\n",
      "\n",
      "\n",
      "    @staticmethod\n",
      "    def sample_func(a, signal):\n",
      "        size =  signal.shape[1] - 1\n",
      "        b = a % size\n",
      "        sample = signal[:,np.array(b)[0]]\n",
      "        a += 1\n",
      "        return sample        \n",
      "\n",
      "\n",
      "# A function to generate a randomly varying x(t) signal\n",
      "def generate_signal(T, dt, rms, limit, seed):\n",
      "    \n",
      "    obj = GaussianWhiteNoise(T, limit, rms)\n",
      "    sample = obj.create_sample(dt, rng=np.random.RandomState(seed))\n",
      "\n",
      "    samples_array = np.arange(T/dt)             # generate evenly spaced samples\n",
      "    x = np.array([sample() for i in samples_array])[:,0]\n",
      "    \n",
      "    X = np.fft.fftshift(np.fft.fft(x))          # compute the discrete fourier transform and \n",
      "                                                      # shift the zero frequency component to the center of the spectrum\n",
      "    \"\"\"# plot the time domain signal\n",
      "    figure()\n",
      "    plot(samples_array*dt, x)          # samples_array is multiplied by dt to limit the duration of the signal to 'T'\n",
      "    xlabel('$t$')\n",
      "    ylabel('$x(t)$')\"\"\"\n",
      "    \n",
      "    return x, X \n",
      "    \n",
      "\n",
      "def compute_filter(dt, tau, n, T, colour):\n",
      "    \n",
      "    t = np.arange(T/dt)*dt\n",
      "    t = t - T/2\n",
      "    h = (t**n) * np.exp(-t/tau) \n",
      "    h[np.where(t<0)]=0\n",
      "    h = h/norm(h,1)\n",
      "\n",
      "    #figure()\n",
      "    plot(t, h, colour)\n",
      "    xlabel('t')\n",
      "    ylabel('h(t)')\n",
      "    xlim(-0.05, 0.08)\n",
      "    \n",
      "    return h\n",
      "\n",
      "\n",
      "def cal_decoders(activity_matrix, x, S):\n",
      "    \n",
      "    gamma = (( np.dot(activity_matrix, activity_matrix.transpose()) ) / S )            # Y = AA^T/S         \n",
      "    upsilon = ( np.dot(activity_matrix, x) ) / S         # v = A.x/S\n",
      "    decoders = np.dot(np.linalg.pinv(gamma), upsilon)        # compute the decoders (using psedo inverse to deal with singular matrices)\n",
      "    return decoders \n",
      "\n",
      "    \n",
      "\n",
      "def main():\n",
      "    T = 2.0         # length of signal in seconds\n",
      "    dt = 0.001      # time step size\n",
      "    \n",
      "    # parameters for x(t)\n",
      "    seed = 3     # the random number seed to use (so we can regenerate the same signal again)\n",
      "    rms = 0.5    # the root mean square power level of the signa\n",
      "    limit = 5    # the maximum frequency for the signal x(t) (in Hz)\n",
      "\n",
      "    # Generate bandlimited white noise (use your own function from part 1.1)\n",
      "    #x, X = syde556.generate_signal(T, dt, rms=0.5, limit=5, seed=3)\n",
      "    x, X = generate_signal(T, dt, rms, limit, seed)\n",
      "\n",
      "    Nt = len(x)                # Count the number of sample points in the time varying signal\n",
      "    t = numpy.arange(Nt) * dt  # Generate evenly spaced (spaced by 1) Nt values between the interval  [0,Nt) and then multipy  \n",
      "                               # by dt to change the spacing between two adjacent points to dt\n",
      "        \n",
      "    # ---------------------------------------------------------------------------------------------------------------------------------\n",
      "    # Part a)\n",
      "    \n",
      "    tau = 0.007\n",
      "    n = 0.0\n",
      "    h = compute_filter(dt, tau, n, T, 'r')\n",
      "\n",
      "    n = 1.0\n",
      "    h = compute_filter(dt, tau, n, T, 'g')\n",
      "    \n",
      "    n = 2.0\n",
      "    h = compute_filter(dt, tau, n, T, 'b')\n",
      "    \n",
      "   # --------------------------------------------------------------------------------------------------------------------------------- \n",
      "    # Part b)\n",
      "    \n",
      "    n = 0.0\n",
      "    \n",
      "    figure()\n",
      "    tau = 0.002\n",
      "    h = compute_filter(dt, tau, n, T, 'r')\n",
      "\n",
      "    tau = 0.005\n",
      "    h = compute_filter(dt, tau, n, T, 'g')\n",
      "    \n",
      "    tau = 0.01\n",
      "    h = compute_filter(dt, tau, n, T, 'b')\n",
      "    \n",
      "    tau = 0.02\n",
      "    h = compute_filter(dt, tau, n, T, 'm')\n",
      "    \n",
      "    # ---------------------------------------------------------------------------------------------------------------------------------\n",
      "    # Part c)\n",
      "    \n",
      "    # Neuron parameters\n",
      "    tau_ref = 0.002          # Absolute Refractory period of the neuron\n",
      "    tau_rc = 0.02            # The RC time constant of the LIF neuron \n",
      "    x0 = 0.0                 # firing rate at x=x0 is a0\n",
      "    a0 = 40.0\n",
      "    x1 = 1.0                 # firing rate at x=x1 is a1\n",
      "    a1 = 150.0\n",
      "\n",
      "    # Calculate the gain, bias and the threshold of the LIF neuron\n",
      "    eps = tau_rc/tau_ref\n",
      "    r1 = 1.0 / (tau_ref * a0)\n",
      "    r2 = 1.0 / (tau_ref * a1)\n",
      "    f1 = (r1 - 1) / eps\n",
      "    f2 = (r2 - 1) / eps\n",
      "    alpha = (1.0/(numpy.exp(f2)-1) - 1.0/(numpy.exp(f1)-1))/(x1-x0) \n",
      "    x_threshold = x0-1/(alpha*(numpy.exp(f1)-1))              \n",
      "    Jbias = 1-alpha*x_threshold;   \n",
      "\n",
      "    # Simulate the two neurons (use your own function from part 3)\n",
      "    #spikes = syde556.two_neurons(x, dt, alpha, Jbias, tau_rc, tau_ref)\n",
      "    \n",
      "    # Create an LIF neuron (1)\n",
      "    neuron1 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder1 = 1.0\n",
      "    \n",
      "    # Create another LIF neuron (2)\n",
      "    neuron2 = spikyLIF(tau_rc, tau_ref) \n",
      "    encoder2 = -1.0\n",
      "    \n",
      "    \n",
      "    # generate the spikes for the first neuron\n",
      "    J1 = neuron1.cal_current(x, alpha, encoder1, Jbias)\n",
      "    V1 = neuron1.generate_spikes(T, dt, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # generate the spikes for the second neuron\n",
      "    J2 = neuron2.cal_current(x, alpha, encoder2, Jbias)\n",
      "    V2 = neuron2.generate_spikes(T, dt, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # create and array of spikes of both neurons\n",
      "    #spikes = np.array([V1, V2])\n",
      " \n",
      "    freq = numpy.arange(Nt)/T - Nt/(2.0*T)   # spread the frequencies evenly from -500 to 500 (enforcing symmetry to get a real signal)\n",
      "    omega = freq*2*numpy.pi                  # Calculate the corresponding angular frequencies in randians\n",
      "    \n",
      "    \n",
      "    n = 0\n",
      "    tau = 0.007\n",
      "    figure()\n",
      "    h = compute_filter(dt, tau, n, T, 'r')\n",
      "    \n",
      "    fspikes1 = np.convolve(V1, h, mode='same')\n",
      "    fspikes2 = np.convolve(V2, h, mode='same')\n",
      "\n",
      "    A = np.array([fspikes1, fspikes2])\n",
      "    decoders = cal_decoders(A, x, Nt)\n",
      "    xhat = np.dot(A.T, decoders)\n",
      "    \n",
      "    # ---------------------------------------------------------------------------------------------------------------------------------\n",
      "    # PLOTS\n",
      "    \n",
      "    # frequency plot for h(t)\n",
      "    H = np.fft.fftshift(np.fft.fft(h))\n",
      "    figure()\n",
      "    plot(freq,np.abs(H))\n",
      "    xlabel('freq(Hz)')\n",
      "    ylabel('|H(w)|')\n",
      "    \n",
      "    # Plot x(t), spikes and decoded xhat value\n",
      "    figure()\n",
      "    plot(t,V1, 'y')\n",
      "    plot(t,V2-1.5, 'm')\n",
      "    plot(t,xhat, 'cyan')\n",
      "    plot(t,x,'r')\n",
      "    \n",
      "    rmse = np.sqrt(np.average((x-xhat)**2)) \n",
      "    print \"RMSE of the original signal:\", rmse\n",
      "       \n",
      "    # ---------------------------------------------------------------------------------------------------------------------------------\n",
      "    # Part d) \n",
      "    limit = 5.0\n",
      "    seed = 10\n",
      "    x2, X2 = generate_signal(T, dt, rms, limit, seed)\n",
      "    \n",
      "    # generate the spikes for the first neuron\n",
      "    J1 = neuron1.cal_current(x2, alpha, encoder1, Jbias)\n",
      "    V1 = neuron1.generate_spikes(T, dt, J=J1, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    # generate the spikes for the second neuron\n",
      "    J2 = neuron2.cal_current(x2, alpha, encoder2, Jbias)\n",
      "    V2 = neuron2.generate_spikes(T, dt, J=J2, Vth=1, V_spike=0.5, R=1)\n",
      "    \n",
      "    \n",
      "    fspikes1_2 = np.convolve(V1, h, mode='same')\n",
      "    fspikes2_2 = np.convolve(V2, h, mode='same')\n",
      "\n",
      "    A2 = np.array([fspikes1_2, fspikes2_2])\n",
      "    xhat2 = np.dot(A2.T, decoders)\n",
      "    \n",
      "    # PLOTS\n",
      "    \n",
      "    # Plot x(t), spikes and decoded xhat value\n",
      "    figure()\n",
      "    plot(t,V1, 'y')\n",
      "    plot(t,V2-1.5, 'm')\n",
      "    plot(t,xhat2, 'cyan')\n",
      "    plot(t,x2,'r')\n",
      "    \n",
      "    rmse2 = np.sqrt(np.average((x2-xhat2)**2)) \n",
      "    print \"RMSE of the new signal:\", rmse2\n",
      "    \n",
      "if __name__ == \"__main__\": main()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}