beat

code/response_beat.py

@@ -1,30 +1,72 @@
 import matplotlib.pyplot as plt
 import numpy as np
-import scipy.stats as ss
 from read_chirp_data import *
 from utility import *
 from IPython import embed
 
-# define sampling rate and data path
+# define data path and important parameters
-sampling_rate = 40 #kHz
 data_dir = "../data"
-cut_window = 20
+sampling_rate = 40 #kHz
+cut_window = 40
+cut_range = np.arange(-cut_window * sampling_rate, 0, 1)
+window = 1
 
-data = ["2018-11-13-aa-invivo-1", "2018-11-13-ac-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1",
+'''
+# norm: -150, 150, 300
+data = ["2018-11-13-aa-invivo-1", "2018-11-13-ac-invivo-1","2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1",
         "2018-11-13-ai-invivo-1", "2018-11-13-aj-invivo-1", "2018-11-13-ak-invivo-1", "2018-11-13-al-invivo-1"]
 
+# norm: -50
+data = ["2018-11-20-aa-invivo-1", "2018-11-20-ab-invivo-1", "2018-11-20-ac-invivo-1","2018-11-20-ad-invivo-1", 
+        "2018-11-20-ae-invivo-1", "2018-11-20-af-invivo-1", "2018-11-20-ag-invivo-1", "2018-11-20-ah-invivo-1",
+        "2018-11-20-ai-invivo-1"]
+
+'''
+data = ["2018-11-14-aa-invivo-1", "2018-11-14-ac-invivo-1", "2018-11-14-ad-invivo-1", "2018-11-14-af-invivo-1",
+        "2018-11-14-ag-invivo-1", "2018-11-14-ah-invivo-1", "2018-11-14-ai-invivo-1", "2018-11-14-ak-invivo-1",
+        "2018-11-14-al-invivo-1", "2018-11-14-am-invivo-1", "2018-11-14-an-invivo-1"]
+
+rates = {}
+
 for dataset in data:
     spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
     df_map = map_keys(spikes)
     print(dataset)
     for df in df_map.keys():
-        beat_duration = 1/df
+        '''
+        if df == 50:
+            pass
+        else:
+            continue
+        '''
+        
+        print(df)
+        rep_rates = []
+        beat_duration = int(abs(1 / df) * 1000)
         beat_window = 0
-        while beat_window + beat_duration <= cut_window:
+        while beat_window + beat_duration <= cut_window/2:
             beat_window = beat_window + beat_duration
         for rep in df_map[df]:
             for phase in spikes[rep]:
-                response = spikes[rep][phase]
+                # get spikes 40 ms before the chirp first chirp
+                spikes_to_cut = np.asarray(spikes[rep][phase])
+                spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window) & (spikes_to_cut < 0)]
+                spikes_idx = np.round(spikes_cut * sampling_rate)
+                # also save as binary, 0 no spike, 1 spike
+                binary_spikes = np.isin(cut_range, spikes_idx) * 1
+                smoothed_data = smooth(binary_spikes, window, 1 / sampling_rate)
+                train = smoothed_data[window:beat_window+window]
+                rep_rates.append(np.std(train))
                 break
-            #cut = response[response[]]
+        df_rate = np.mean(rep_rates)
+        if df in rates.keys():
+            rates[df].append(df_rate)
+        else:
+            rates[df] = [df_rate]
 
+fig, ax = plt.subplots()
+for i, k in enumerate(sorted(rates.keys())):
+    ax.plot(np.ones(len(rates[k]))*k, rates[k], 'o')
+#ax.legend(sorted(rates.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
+fig.tight_layout()
+plt.show()

code/spikes_analysis.py

@@ -71,7 +71,7 @@ for dataset in data:
                     # check the phase
                     if phase[1] > phase_vec[idx] and phase[1] < phase_vec[idx+1]:
 
-                        # get spikes between 50 ms before and after the chirp
+                        # get spikes between 40 ms before and after the chirp
                         spikes_to_cut = np.asarray(spikes[rep][phase])
                         spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window*2) & (spikes_to_cut < cut_window*2)]
                         spikes_idx = np.round(spikes_cut*sampling_rate)