Changeset 487 for trunk

Timestamp:

07/15/11 18:52:03 (10 months ago)

Author:

mpereira

bzr:base-revision:

michael.fsp@gmail.com-20110714163140-v3qk0357xh8xn2zv

bzr:committer:

Michael Pereira <michael.fsp@gmail.com>

bzr:file-ids:

src/analysis.py 106@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fanalysis.py
src/signals/spikes.py 375@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fsignals%2Fspikes.py

bzr:mapping-version:

v4

bzr:merge:

michael.fsp@gmail.com-20110715165055-30cyagkrbyn2q7yv

bzr:repository-uuid:

400bb3d0-8d2e-0410-a2c0-f692ac833539

bzr:revision-id:

michael.fsp@gmail.com-20110715165126-dd4oub61g4lous2j

bzr:revno:

431

bzr:revprop:branch-nick:

trunk

bzr:root:

trunk

bzr:text-revisions:

src/analysis.py michael.fsp@gmail.com-20110715165055-30cyagkrbyn2q7yv
src/signals/spikes.py michael.fsp@gmail.com-20110715165055-30cyagkrbyn2q7yv

bzr:timestamp:

2011-07-15 18:51:26.930000067 +0200

bzr:user-agent:

bzr2.3.1+bzr-svn1.0.5dev0

svn:original-date:

2011-07-15T16:51:26.930000Z

Message:

added new method SpikeList.averaged_instantaneous_rate

Location:

trunk/src

Files:

• analysis.py (modified) (3 diffs)
• signals/spikes.py (modified) (4 diffs)

trunk/src/analysis.py

@@ -294 +294 @@
                                 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
+        time_stamp_resolution - temporal resolution of input and output in ms
         direction             - Asymmetric kernels have two possible directions.
                                 The values are -1 or 1, default is 1. The
@@ -323 +321 @@
     
     To obtain single trial rate function of trial one should use
-        r = norm * numpy.convolve(sua, kernel)
+        r = norm * scipy.signal.fftconvolve(sua, kernel)
     To obtain trial-averaged spike train one should use
-        r_avg = norm * numpy.convolve(sua, np.mean(X,1))
+        r_avg = norm * scipy.signal.fftconvolve(sua, numpy.mean(X,1))
         
         where X is an array of shape (l,n), where n is the number of trials and
         l the length of each trial
-
-    Note that filter assumes a causal filtering. To represent the resulting
-    firing rate function at center of gravity (acausal) use
-
-        r=r(length(kernel):end);    -> discard first length(kernel)-1
-                                       elements which considered zero entries
-                                       (border effect)
-        t=((1:length(r))+m_idx)*TimeStampResolution -> in seconds
-        plot(t,r)   -> in second time scale
+    
+    Note that the output of scipy.signal.fftconvolve needs to trimmed acordingly
+    before being displayed. For more information see the source of the method
+    SpikeTrain.instantaneous_rate()
+    
+    See also:
+        SpikeTrain.instantaneous_rate, SpikeList.averaged_instantaneous_rate
     """
     assert form.upper() in ('BOX','TRI','GAU','EPA','EXP','ALP'), "form must \
@@ -343 +339 @@
     
     assert direction in (1,-1), "direction must be either 1 or -1"
+    
+    time_stamp_resolution = time_stamp_resolution / 1000. #convert to SI units
     
     norm = 1./time_stamp_resolution

trunk/src/signals/spikes.py

@@ -573 +573 @@
         return hist
     
-    def instantaneous_rate(self, resolution, kernel, norm, t_start=None,
-                           t_stop=None, acausal=True, m_idx=None,
-                           trim=False):
+    def instantaneous_rate(self, resolution, kernel, norm, m_idx=None,
+                           t_start=None, t_stop=None, acausal=True, trim=False):
         """
         Estimate instantaneous firing rate by kernel convolution.
@@ -603 +602 @@
                           NOTE: if True and an assymetrical kernel is provided
                           the output will not be aligned with [t_start, t_stop]
+                          
+        See also:
+            analysis.make_kernel
         """
         
@@ -1038 +1040 @@
             return self.id_list
         elif type(sub_list) == int:
+            import ipdb; ipdb.set_trace()
             return numpy.random.permutation(self.id_list)[0:sub_list]
         else:
@@ -1786 +1789 @@
             return result
 
+    def averaged_instantaneous_rate(self, resolution, kernel, norm, m_idx=None,
+                                    t_start=None, t_stop=None, acausal=True,
+                                    trim=False):
+        """
+        Estimate the instantaneous firing rate averaged across neurons in the
+        SpikeList, by kernel convolution.
+        
+        Inputs:
+            resolution  - time stamp resolution of the spike times (ms). the
+                          same resolution will be assumed for the kernel
+            kernel      - kernel function used to convolve with
+            norm        - normalization factor associated with kernel function
+                          (see analysis.make_kernel for details)
+            t_start     - start time of the interval used to compute the firing
+                          rate
+            t_stop      - end time of the interval used to compute the firing
+                          rate (included)
+            acausal     - if True, acausal filtering is used, i.e., the gravity
+                          center of the filter function is aligned with the
+                          spike to convolve
+            m_idx       - index of the value in the kernel function vector that
+                          corresponds to its gravity center. this parameter is
+                          not mandatory for symmetrical kernels but it is
+                          required when assymmetrical kernels are to be aligned
+                          at their gravity center with the event times
+            trim        - if True, only the 'valid' region of the convolved
+                          signal are returned, i.e., the points where there
+                          isn't complete overlap between kernel and spike train
+                          are discarded
+                          NOTE: if True and an assymetrical kernel is provided
+                          the output will not be aligned with [t_start, t_stop]
+                          
+        See also:
+            analysis.make_kernel, SpikeTrain.instantaneous_rate
+        """
+        
+        if t_start is None:
+            t_start = self.t_start
+            
+        if t_stop is None:
+            t_stop = self.t_stop
+        
+        if m_idx is None:
+            m_idx = kernel.size / 2
+            
+        spikes_slice = []
+        for i in self:
+            train_slice = i.spike_times[(i.spike_times >= t_start) & (
+            i.spike_times <= t_stop)]
+            spikes_slice += train_slice.tolist()
+        
+        time_vector = numpy.zeros((t_stop - t_start)/resolution + 1)
+        
+        for spike in spikes_slice:
+            index = (spike - t_start) / resolution
+            time_vector[index] += 1.
+            
+        avg_time_vector = time_vector / float(self.id_list.size)
+        
+        r = norm * scipy.signal.fftconvolve(avg_time_vector, kernel, 'full')
+
+        if acausal is True:
+            if trim is False:
+                r = r[m_idx:-(kernel.size - m_idx)]
+                t_axis = numpy.linspace(t_start, t_stop, r.size)
+                return t_axis, r
+            
+            elif trim is True:
+                r = r[2 * m_idx:-2*(kernel.size - m_idx)]
+                t_start = t_start + m_idx * resolution
+                t_stop = t_stop - ((kernel.size) - m_idx) * resolution
+                t_axis = numpy.linspace(t_start, t_stop, r.size)
+                return t_axis, r
+            
+        if acausal is False:
+            if trim is False:
+                r = r[m_idx:-(kernel.size - m_idx)]
+                t_axis = (numpy.linspace(t_start, t_stop, r.size) +
+                          m_idx * resolution)
+                return t_axis, r
+            
+            elif trim is True:
+                r = r[2 * m_idx:-2*(kernel.size - m_idx)]
+                t_start = t_start + m_idx * resolution
+                t_stop = t_stop - ((kernel.size) - m_idx) * resolution
+                t_axis = (numpy.linspace(t_start, t_stop, r.size) +
+                          m_idx * resolution)
+                return t_axis, r
+    
     def fano_factor(self, time_bin):
         """