stuff for plotting

code/analysis_rs.py

@@ -1,59 +1,54 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from read_baseline_data import *
-from IPython import embed
 from NixFrame import *
 from utility import *
+from IPython import embed
 
+# plot and data values
 inch_factor = 2.54
 data_dir = '../data'
 dataset = '2018-11-09-ad-invivo-1'
+
+# read eod and time of baseline
 time, eod = read_baseline_eod(os.path.join(data_dir, dataset))
 
-#fig = plt.figure(figsize=(12/inch_factor, 8/inch_factor))
+fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
-#ax = fig.add_subplot(111)
+ax.plot(time[:1000], eod[:1000])
-#ax.plot(time[:1000], eod[:1000])
+ax.set_xlabel('time [ms]', fontsize=12)
-#ax.set_xlabel('time [ms]', fontsize=12)
+ax.set_ylabel('voltage [mV]', fontsize=12)
-#ax.set_ylabel('voltage [mV]', fontsize=12)
+plt.xticks(fontsize=8)
-#plt.xticks(fontsize = 8)
+plt.yticks(fontsize=8)
-#plt.yticks(fontsize = 8)
+fig.tight_layout()
-#fig.tight_layout()
 #plt.savefig('eod.pdf')
+plt.show()
 
-#interspikeintervalhistogram, windowsize = 1 ms
+# read spikes during baseline activity
-#plt.hist
-#coefficient of variation
-#embed()
-#exit()
-
 spikes = read_baseline_spikes(os.path.join(data_dir, dataset))
+# calculate interpike intervals and plot them
 interspikeintervals = np.diff(spikes)
-#fig = plt.figure()
-#plt.hist(interspikeintervals, bins=np.arange(0, np.max(interspikeintervals), 0.0001))
-#plt.show()
 
+fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
+plt.hist(interspikeintervals, bins=np.arange(0, np.max(interspikeintervals), 0.0001))
+plt.show()
+
+# calculate coefficient of variation
 mu = np.mean(interspikeintervals)
 sigma = np.std(interspikeintervals)
 cv = sigma/mu
-#print(cv)
+print(cv)
-
-# calculate zero crossings of the eod
-# plot mean of eod circles
-# plot std of eod circles
-# plot psth into the same plot
-# calculate vector strength
 
+# calculate eod times and indices by zero crossings
 threshold = 0
 shift_eod = np.roll(eod, 1)
 eod_times = time[(eod >= threshold) & (shift_eod < threshold)]
 sampling_rate = 40000.0
 eod_idx = eod_times*sampling_rate
 
-
+# align eods and spikes to eods
 max_cut = int(np.max(np.diff(eod_idx)))
 eod_cuts = np.zeros([len(eod_idx)-1, max_cut])
-# eods 15 + 16 are to short
+spike_times = []
-relative_times = []
 eod_durations = []
 
 for i, idx in enumerate(eod_idx[:-1]):
@@ -62,25 +57,23 @@ for i, idx in enumerate(eod_idx[:-1]):
 	eod_cuts[i, len(eod_cut):] = np.nan
 	time_cut = time[int(idx):int(eod_idx[i+1])]
 	spike_cut = spikes[(spikes > time_cut[0]) & (spikes < time_cut[-1])]
-	relative_time  = spike_cut - time_cut[0]
+	spike_time = spike_cut - time_cut[0]
-	if len(relative_time) > 0:
+	if len(spike_time) > 0:
-		relative_times.append(relative_time[:][0]*1000)
+		spike_times.append(spike_time[:][0]*1000)
-		eod_durations.append(len(eod_cut))
+		eod_durations.append(len(eod_cut)/sampling_rate*1000)
 
+# calculate vector strength
+vs = vector_strength(spike_times, eod_durations)
 
+# determine means and stds of eod for plot
+# determine time axis
 mu_eod = np.nanmean(eod_cuts, axis=0)
 std_eod = np.nanstd(eod_cuts, axis=0)*3
+time_axis = np.arange(max_cut)/sampling_rate*1000
 
-vs = vector_strength(relative_times, eod_durations)
+# plot eod form and spike histogram
-embed()
-exit()
-
-#time_axis = np.arange(max_cut)/sampling_rate*1000
-
-#fig = plt.figure(figsize=(12/inch_factor, 8/inch_factor))
-'''
 fig, ax1 = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
-ax1.hist(relative_times, color='crimson')
+ax1.hist(spike_times, color='crimson')
 ax1.set_xlabel('time [ms]', fontsize=12)
 ax1.set_ylabel('number', fontsize=12)
 ax1.tick_params(axis='y', labelcolor='crimson')
@@ -88,7 +81,6 @@ plt.yticks(fontsize = 8)
 ax1.spines['top'].set_visible(False)
 
 ax2 = ax1.twinx()
-
 ax2.fill_between(time_axis, mu_eod+std_eod, mu_eod-std_eod, color='dodgerblue', alpha=0.5)
 ax2.plot(time_axis, mu_eod, color='black', lw=2)
 ax2.set_ylabel('voltage [mV]', fontsize=12)
@@ -98,5 +90,3 @@ plt.xticks(fontsize = 8)
 plt.yticks(fontsize=8)
 fig.tight_layout()
 plt.show()
-'''
-

code/utility.py

@@ -12,9 +12,8 @@ def zero_crossing(eod,time):
 
 def vector_strength(spike_times, eod_durations):
 	n = len(spike_times)
-	phase_times = np.zeros(n)
+	phase_times = np.zeros(len(spike_times))
 	for i, idx in enumerate(spike_times):
-		phase_times[i]= spike_times[i]/eod_durations[i]
+		phase_times[i] = (spike_times[i] / eod_durations[i]) * 2 * np.pi
-
 	vs = np.sqrt((1/n*sum(np.cos(phase_times)))**2 + (1/n*sum(np.sin(phase_times)))**2)
 	return vs