Changeset 451 for trunk/src/analysis.py

Timestamp:

08/23/10 13:47:57 (21 months ago)

Author:

mpereira

Message:

Added make_kernel() to analysis.py

Files:

• trunk/src/analysis.py (modified) (1 diff)

trunk/src/analysis.py

@@ -92 +92 @@
             return k
         
+def make_kernel(form, sigma, time_stamp_resolution, direction=1, kernel=None,
+                m_idx=None):
+    """
+    Constructs a numeric linear convolution kernel of basic shape to be used
+    for data smoothing (linear low pass filtering) and firing rate estimation
+    from single trial or trial-averaged spike trains.
+    
+    Exponential and alpha kernels may also be used to represent postynaptic
+    currents / potentials in a linear (current-based) model.
+    
+    Adapted from original script written by Martin P. Nawrot for the FIND MATLAB
+    toolbox.
+    See FIND - a unified framework for neural data analysis,
+        Meier R, Egert U, Aertsen A, Nawrot MP; Neural Netw. 2008 Oct;
+        21(8):1085-93. 
+        
+        Nawrot M, Aertsen A, Rotter S (1999)
+        Single-trial estimation of neuronal firing rates - From single neuron
+        spike trains to population activity.
+        J Neurosci Meth 94: 81-92
+        
+    Inputs:
+        form                  - kernel form (string) Currently implemented forms
+                                are BOX - boxcar, TRI - triangle, GAU - gaussian
+                                , EPA - epanechnikov, EXP - exponential, ALP -
+                                alpha function.
+                                EXP and ALP are aymmetric kernel forms and
+                                assume optional parameter 'direction'
+        sigma                 - standard deviation of the distribution
+                                associated with kernel shape. This parameter
+                                defines the time resolution of the kernel
+                                estimate and makes different kernels comparable
+                                (cf. [1] for symetric kernels). This is used
+                                here as an alternative definition to the cut-off
+                                frequency of the associated linear filter.
+        time_stamp_resolution - temporal resolution of input and output in SI
+                                units. E.g. value 0.001 means 1ms time
+                                resolution
+        direction             - Asymmetric kernels have two possible directions.
+                                The values are -1 or 1, default is 1. The
+                                definition here is that for direction = 1 the
+                                kernel represents the impulse response function
+                                of the linear filter. Default value is 1.
+    
+    Outputs:
+        kernel  - array of kernel. The length of this array is always an odd
+                  number to represent symmetric kernels such that the center bin
+                  coincides with the median of the numeric array, i.e for a
+                  triangle, the maximum will be at the center bin with equal
+                  number of bins to the right and to the left.
+        norm    - for rate estimates. The kernel vector is normalized such that
+                  the sum of all entries equals unity sum(kernel)=1. When
+                  estimating rate functions from discrete spike data (0/1) the
+                  additional parameter 'norm' allows for the normalization to
+                  rate in spikes per second.
+                  Use: rate = norm * scipy.signal.lfilter(kernel, 1, spike_data)
+        m_idx   - index of the numerically determined median (center of gravity)
+                  of the kernel function
+    """
+    assert form.upper() in ('BOX','TRI','GAU','EPA','EXP','ALP'), "form must \
+    be one either 'BOX','TRI','GAU','EPA','EXP' or 'ALP'!"
+    
+    assert direction in (1,-1), "direction must be either 1 or -1"
+    
+    norm = 1./time_stamp_resolution
+
+    if form.upper() == 'BOX':
+        w = 2.0 * sigma * numpy.sqrt(3)
+        width = 2 * numpy.floor(w / 2.0 / time_stamp_resolution) + 1  # always odd number of bins
+        height = 1. / width
+        kernel = numpy.ones((1, width)) * height  # area = 1
+        
+    elif form.upper() == 'TRI':
+        w = 2 * sigma * numpy.sqrt(6)
+        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution)
+        trileft = numpy.arange(1, halfwidth + 2)
+        triright = numpy.arange(halfwidth, 0, -1)  # odd number of bins
+        triangle = numpy.append(trileft, triright)
+        kernel = triangle / triangle.sum()  # area = 1
+        
+    elif form.upper() == 'EPA':
+        w = 2.0 * sigma * numpy.sqrt(5)
+        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution)
+        base = numpy.arange(-halfwidth, halfwidth + 1)
+        parabula = base**2
+        epanech = parabula.max() - parabula  # inverse parabula
+        kernel = epanech / epanech.sum()  # area = 1
+        
+    elif form.upper() == 'GAU':
+        SI_sigma = sigma / 1000.0
+        w = 2.0 * sigma * 2.7  # > 99% of distribution weight
+        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution)  # always odd 
+        base = numpy.arange(-halfwidth, halfwidth + 1) / 1000.0 * (
+            time_stamp_resolution)
+        g = numpy.exp(-(base**2) / 2.0 / SI_sigma**2) / SI_sigma / numpy.sqrt(
+            2.0 * numpy.pi)
+        kernel = g / g.sum()
+        
+    elif form.upper() == 'ALP':
+        SI_sigma = sigma / 1000.0
+        w = 5.0 * sigma
+        alpha = numpy.arange(1, (2.0 * numpy.floor(
+            (w / time_stamp_resolution / 2.0)) + 1) + 1) / 1000.0 * \
+        time_stamp_resolution
+        alpha = (2.0 / SI_sigma**2) * alpha * numpy.exp(-alpha * numpy.sqrt(2) \
+                                                        / SI_sigma)
+        kernel = alpha / alpha.sum()  # normalization
+        if direction == -1:
+            kernel = numpy.flipud(kernel)
+            
+    elif form.upper() == 'EXP':
+        SI_sigma = sigma / 1000.0
+        w = 5.0 * sigma
+        expo = numpy.arange(1, (2.0 * numpy.floor(w / time_stamp_resolution / (
+            2.0)) + 1) + 1) / 1000.0 * time_stamp_resolution
+        expo = numpy.exp(-expo / SI_sigma)
+        kernel = expo / expo.sum()
+        if direction == -1:
+            kernel = numpy.flipud(kernel)
+    
+    kernel = kernel.ravel()
+    m_idx = numpy.nonzero(kernel.cumsum() >= 0.5)[0].min()
+    
+    return kernel, norm, m_idx
+
 def simple_frequency_spectrum(x):
     """