Changeset 185

Timestamp:

07/24/08 17:23:32 (4 months ago)

Author:

apdavison

Message:

Added a method to calculate the frequency spectrum of a SpikeTrain. This is a very simple method with no corrections for the finite duration of the spike train, etc., so feel free to improve it.

Files:

• trunk/src/analysis.py (modified) (1 diff)
• trunk/src/signals.py (modified) (5 diffs)
• trunk/test/test_analysis.py (added)

trunk/src/analysis.py

@@ -96 +96 @@
-        if D[k] == max_val:
-            return k
+        
+def simple_frequency_spectrum(x):
+    """
+    Very simple calculation of frequency spectrum with no detrending,
+    windowing, etc. Just the first half (positive frequency components) of
+    abs(fft(x))
+    """
+    spec = numpy.absolute(numpy.fft.fft(x))
+    spec = spec[:len(x)/2] # take positive frequency components
+    spec /= len(x)         # normalize
+    spec *= 2.0            # to get amplitudes of sine components, need to multiply by 2
+    spec[0] /= 2.0         # except for the dc component
+    return spec
-
-class TuningCurve(object):

trunk/src/signals.py

@@ -10 +10 @@
-except Exception:
-    print "Warning: pylab not present"
+from NeuroTools import analysis
-
-
@@ -120 +121 @@
-        return len(self.spike_times)
-
+    def duration(self):
+        return self.t_stop - self.t_start
+
-    def format(self, relative=False, quantized=False):
-        """
@@ -308 +312 @@
-        return tuning_curve
-
+    def frequency_spectrum(self, time_bin):
+        """
+        Returns the frequency spectrum of the time histogram together with the
+        frequency axis.
+        """
+        hist = self.time_histogram(time_bin, normalized=False)
+        freq_spect = analysis.simple_frequency_spectrum(hist)
+        freq_bin = 1000.0/self.duration() # Hz
+        freq_axis = numpy.arange(len(freq_spect)) * freq_bin    
+        return freq_spect, freq_axis    
+
+    def f1f0_ratio(self, time_bin, f_stim):
+        """
+        Returns the F1/F0 amplitude ratio where the input stimulus frequency is
+        f_stim.
+        """
+        freq_spect = self.frequency_spectrum(time_bin)[0]
+        F0 = freq_spect[0]
+        freq_bin = 1000.0/self.duration()
+        try:
+            F1 = freq_spect[int(round(f_stim/freq_bin))]
+        except IndexError:
+            errmsg = "time_bin (%f) too large to estimate F1/F0 ratio for an input frequency of %f" % (time_bin, f_stim)
+            errmsg += "\nFrequency_spectrum: %s" % freq_spect
+            raise Exception(errmsg)
+        return F1/F0
+    
-
-class SpikeList(object):
@@ -412 +443 @@
-            self.spiketrains[id] = spiketrain
-            self.id_list.append(id)
-            self.N += 1
-        self.__calc_startstop()
-
@@ -915 +945 @@
-    def as_pyNN_SpikeArray(self):
-        """
-o
+u
-
-        >>> pop = Population((10,),SpikeSourceArray, {'spike_times': sl.as_pyNN_SpikeArray() })