diff --git a/code/analysis_rs.py b/code/analysis_rs.py index 356d461..950b572 100644 --- a/code/analysis_rs.py +++ b/code/analysis_rs.py @@ -1,58 +1,55 @@ 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)) -#ax = fig.add_subplot(111) -#ax.plot(time[:1000], eod[:1000]) -#ax.set_xlabel('time [ms]', fontsize=12) -#ax.set_ylabel('voltage [mV]', fontsize=12) -#plt.xticks(fontsize = 8) -#plt.yticks(fontsize = 8) -#fig.tight_layout() +fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor)) +ax.plot(time[:1000], eod[:1000]) +ax.set_xlabel('time [ms]', fontsize=12) +ax.set_ylabel('voltage [mV]', fontsize=12) +plt.xticks(fontsize=8) +plt.yticks(fontsize=8) +fig.tight_layout() #plt.savefig('eod.pdf') +plt.show() -#interspikeintervalhistogram, windowsize = 1 ms -#plt.hist -#coefficient of variation -#embed() -#exit() - +# read spikes during baseline activity 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) - -# 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 +print(cv) -threshold = 0; +# 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 -relative_times = [] +spike_times = [] +eod_durations = [] for i, idx in enumerate(eod_idx[:-1]): eod_cut = eod[int(idx):int(eod_idx[i+1])] @@ -60,34 +57,38 @@ 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] - if len(relative_time) > 0: - relative_times.append(relative_time[:][0]*1000) + spike_time = spike_cut - time_cut[0] + if len(spike_time) > 0: + spike_times.append(spike_time[:][0]*1000) + 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 -#fig = plt.figure(figsize=(12/inch_factor, 8/inch_factor)) +# plot eod form and spike histogram 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') -plt.yticks(fontsize = 8) +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) ax2.tick_params(axis='y', labelcolor='dodgerblue') -plt.xticks(fontsize = 8) -plt.yticks(fontsize = 8) +plt.xticks(fontsize=8) +plt.yticks(fontsize=8) fig.tight_layout() plt.show() +#NixToFrame(data_dir) \ No newline at end of file diff --git a/code/utility.py b/code/utility.py index 86f2407..1ad9601 100644 --- a/code/utility.py +++ b/code/utility.py @@ -1,8 +1,8 @@ import numpy as np -def zero_crossing(eod,time): - threshold = 0; +def zero_crossing(eod, time): + threshold = 0 shift_eod = np.roll(eod, 1) eod_times = time[(eod >= threshold) & (shift_eod < threshold)] sampling_rate = 40000.0 @@ -10,9 +10,10 @@ def zero_crossing(eod,time): return eod_idx -def vector_strength(spike_times, eod_durations) - alphas = spike_times/ eod_durations - cs = (1/len(spike_times))*np.sum(np.cos(alphas))^2 - sn = (1/len(spike_times))*np.sum(np.sin(alphas))^2 - vs = np.sprt(cs+sn) +def vector_strength(spike_times, eod_durations): + n = len(spike_times) + phase_times = np.zeros(len(spike_times)) + for i, idx in enumerate(spike_times): + 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