diff --git a/differentcells_filter.py b/differentcells_filter.py
new file mode 100644
index 0000000..aeb12e8
--- /dev/null
+++ b/differentcells_filter.py
@@ -0,0 +1,209 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+#from axes import label_axes, labelaxes_params
+from myfunctions import auto_rows
+from myfunctions import default_settings
+from myfunctions import remove_tick_marks
+from myfunctions import remove_tick_ymarks
+import matplotlib.gridspec as gridspec
+from rotated import ps_df
+
+def plot_single_cells(ax, data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']):
+    colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+    # labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
+    #baseline = pd.read_pickle('data_baseline.pkl')
+    #data_beat = pd.read_pickle('data_beat.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    #data = np.unique(data_all['dataset'])[0]
+    #data = np.unique(data_all['dataset'])
+    end = ['original', '005','05', '2' ]
+
+    end = ['original','005']
+    y_sum = [[]]*len(data)*len(end)
+    counter = 0
+    for dd,set in enumerate(data):
+        for ee, e in enumerate(end):
+            d = data_all[data_all['dataset'] == set]
+            #x = d['delta_f'] / d['eodf'] + 1
+            #embed()
+            #y = d['result_frequency_' + e]
+            y = d['result_amplitude_max_' + e]
+            #y2 = d['result_amplitude_max_' + e]
+            y_sum[counter] = np.nanmax(y)
+            counter += 1
+            #print(np.nanmax(y))
+    #embed()
+    lim = np.max(y_sum)
+    hd = 0.3
+    ws = 0.35
+    rows = 1
+    cols = 3
+    #embed()
+    main_grid = gridspec.GridSpec(1, 2,hspace=0.4, width_ratios = [1,7])
+    filter_grid = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=main_grid[0], hspace=0.4)
+    upper_filter_g = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=filter_grid[0], wspace=ws, hspace=0.4)  # ,
+    grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=upper_filter_g[0], wspace=ws, hspace=0.4)#,
+    wish_df = 150
+    fc = 'lightgrey'
+    ec = 'grey'
+    sampling_rate = 40000
+    data_beat = pd.read_pickle('data_beat.pkl')
+    df, p, f, db = ps_df(data_beat, d=set, wish_df=wish_df, window='no', sampling_rate=sampling_rate)
+    ax = {}
+    ax_nr = 0
+    colors = ['brown']
+    #embed()
+    ax[ax_nr] = plt.subplot(grid1[0])
+    ax[ax_nr].spines['right'].set_visible(False)
+    ax[ax_nr].spines['left'].set_visible(False)
+    ax[ax_nr].spines['top'].set_visible(False)
+    ax[ax_nr].spines['bottom'].set_visible(False)
+    ax[ax_nr].set_yticks([])
+    ax[ax_nr].set_xticks([])
+    ax[ax_nr].plot(p, f, color=colors[0])
+    eodf = d['eodf'].iloc[0]
+    mult = 1.1
+    ax[ax_nr].fill_between([np.min(p), np.max(p)], [eodf / 2, eodf / 2], color=fc, edgecolor=ec)
+    ax[ax_nr].set_ylim([0, eodf *0.6])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    upper_filter_g = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=filter_grid[1], wspace=ws, hspace=0.4)  # ,
+    grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=upper_filter_g[0], wspace=ws, hspace=0.4)#,
+    i = 0
+    sigma = 0.0005
+    df, p, f, db = ps_df(data_beat, d=set, wish_df=wish_df, window=sigma * sampling_rate,
+                                                     sampling_rate=sampling_rate)
+    ax_nr = 1
+    sigmaf = 1 / (2 * np.pi * sigma)
+    gauss = np.exp(-(f ** 2 / (2 * sigmaf ** 2)))
+    stepsize = np.abs(f[0]-f[1])
+    wide = 2
+    scale = 1
+    prev_height = np.max((p[int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p[int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    ax[ax_nr] = plt.subplot(grid1[1])
+    ax[ax_nr].spines['right'].set_visible(False)
+    ax[ax_nr].spines['left'].set_visible(False)
+    ax[ax_nr].spines['top'].set_visible(False)
+    ax[ax_nr].set_yticks([])
+    ax[ax_nr].set_xticks([])
+    ax[ax_nr].spines['bottom'].set_visible(False)
+    ax[ax_nr].fill_between(max(p) * gauss ** 2, f, facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    grid = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=main_grid[1], hspace=0.4)
+
+
+    #grid = gridspec.GridSpecFromSubplotSpec(2,1,subplot_spec=grid[1], hspace = 0.4)
+
+    grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=grid[0], wspace=ws, hspace=0.4)#,
+
+
+    end = ['original']
+
+    color_modul = 'steelblue'
+    color_mpf = 'red'
+
+    y_sum1 = plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,title = True,xlabel = False, label =False,remove =True)
+
+    grid2 = gridspec.GridSpecFromSubplotSpec(rows, cols, subplot_spec=grid[1], wspace=ws, hspace=0.4)#,
+
+    end = ['05']
+    #y_sum = [[]] * len(data)
+    color_modul = 'steelblue'
+    color_mpf = 'red'
+    y_sum2 = plot_single(lim, data, end, data_all, grid2, color_mpf, color_modul,title = False,  label = True,remove =False)
+
+    #for dd, d in enumerate(data):
+        # embed()
+        #embed()
+        #ax[1].set_ylim([0, np.nanmax([y_sum1,y_sum2])])
+        #ax[1, dd].set_ylim([0, 350])
+    plt.subplots_adjust(wspace = 0.4,left = 0.17, right = 0.96,bottom = 0.2)
+
+def plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,label = True, xlabel = True,remove =True, title = True):
+    y_sum = [[]] * len(data)
+    ax = {}
+    for dd,set in enumerate(data):
+        for ee, e in enumerate(end):
+            if title == True:
+                plt.title('Cell '+str(dd))
+            d = data_all[data_all['dataset'] == set]
+            x = d['delta_f'] / d['eodf'] + 1
+            #embed()
+            y = d['result_frequency_' + e]
+            y2 = d['result_amplitude_max_' + e]
+            y_sum[dd] = np.nanmax(y)
+            ff = d['delta_f'] / d['eodf'] + 1
+            fe = d['beat_corr']
+            #fig.suptitle(set)
+            grid2 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=grid1[dd], wspace=0.02, hspace=0.2)  # ,
+            ax[0] = plt.subplot(grid2[0])
+            ax[0].plot(ff, fe, color='grey', linestyle='--')
+            ax[0].plot(x, y, color=color_mpf)
+
+            ax[set+e+str(1)] = plt.subplot(grid2[1])
+            if (set == data[0]) and (label == True):
+                ax[set+e+str(1)].set_ylabel('Modulation ')
+                ax[0].set_ylabel('MPF [EODf]')
+            if (set == data[0]) and (xlabel == True):
+                ax[set + e + str(1)].set_xlabel('EOD multiples')
+            #ax[0, dd].set_title(e + ' ms')
+            ax[0].set_xlim([0, 4])
+
+            ax[set+e+str(1)].plot(x, y2, color=color_modul)
+            # ax[1,0].set_ylabel('modulation depth [Hz]')
+            #ax[2, ee].plot(x, y2, color=colors[0])
+            #ax[2, 0].set_ylabel('         modulation depth [Hz]')
+            # ax[1,ee].annotate("", xy=(0.53, 16.83), xytext=(0.53, 17.33), arrowprops=dict(arrowstyle="->"))
+            # ax[1,ee].annotate("", xy=(1.51, 16.83), xytext=(1.51, 17.33), arrowprops=dict(arrowstyle="->"))
+
+            #ax[1, 0].set_xlabel('stimulus frequency [EODf]')
+
+            #ax[1, 2].set_xlabel('stimulus frequency [EODf]')
+            #ax[2, 3].set_xlabel('stimulus frequency [EODf]')
+
+            ax[0].spines['right'].set_visible(False)
+            ax[0].spines['top'].set_visible(False)
+            ax[set+e+str(1)].spines['right'].set_visible(False)
+            ax[set+e+str(1)].spines['top'].set_visible(False)
+            #ax[2, ee].spines['right'].set_visible(False)
+            #ax[2, ee].spines['top'].set_visible(False)
+            ax[0].set_xlim([0, 5])
+            #plt.tight_layout()
+            # fig.label_axes()
+            ax[set+e+str(1)].set_ylim([0, lim])
+            #ax[0].set_ylim([0, 240])
+            #ax[set+e+str(1)] = remove_tick_marks(ax[0])
+            ax[0] = remove_tick_marks(ax[0])
+            if set != data[0]:
+                ax[0] = remove_tick_ymarks(ax[0])
+                ax[set + e + str(1)] = remove_tick_ymarks(ax[set+e+str(1)] )
+            if remove == True:
+                ax[set+e+str(1)] = remove_tick_marks(ax[set+e+str(1)])
+            #ax[0].set_ylim([0, lim])
+    return y_sum
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']
+    data = ['2019-10-21-aa-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
+    default_settings(data,intermediate_width = 6.7,intermediate_length = 5)
+    #fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
+    ax = {}
+    plot_single_cells(ax, data = data)
+    #fig.savefig()
+    plt.savefig('differentcells_filter.pdf')
+    plt.savefig('../highbeats_pdf/differentcells_filter.pdf')
+    # plt.subplots_adjust(left = 0.25)
+    plt.show()
+    #plt.close()
diff --git a/differentcells_trans.py b/differentcells_trans.py
new file mode 100644
index 0000000..2c3c186
--- /dev/null
+++ b/differentcells_trans.py
@@ -0,0 +1,191 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+#from axes import label_axes, labelaxes_params
+from myfunctions import auto_rows
+from myfunctions import default_settings
+from myfunctions import remove_tick_marks
+import matplotlib.gridspec as gridspec
+import string
+
+def plot_single_cells(ax, data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']):
+    colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+    # labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
+    #baseline = pd.read_pickle('data_baseline.pkl')
+    #data_beat = pd.read_pickle('data_beat.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    #data = np.unique(data_all['dataset'])[0]
+    #data = np.unique(data_all['dataset'])
+    end = ['original', '005','05', '2' ]
+
+    end = ['original','005']
+    y_sum = [[]]*len(data)*len(end)
+    counter = 0
+    for dd,set in enumerate(data):
+        for ee, e in enumerate(end):
+            d = data_all[data_all['dataset'] == set]
+            #x = d['delta_f'] / d['eodf'] + 1
+            #embed()
+            #y = d['result_frequency_' + e]
+            y = d['result_amplitude_max_' + e]
+            #y2 = d['result_amplitude_max_' + e]
+            y_sum[counter] = np.nanmax(y)
+            counter += 1
+            #print(np.nanmax(y))
+    #embed()
+    lim = np.max(y_sum)
+
+    #embed()
+    grid = gridspec.GridSpec(1,3,width_ratios = [0.2,4,4],)
+    grid0 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=grid[0], wspace=0.02, hspace=0.1)  # ,
+
+    label1 = plt.subplot(grid0[0])
+    label2 = plt.subplot(grid0[1])
+    label3 = plt.subplot(grid0[2])
+    labels = [label1,label2,label3]
+    titels = ['Cell 1','Cell 2','Cell 3']# < 0.5 EODf,with 0.5 ms wide
+    fs_big = 11
+    weight = 'bold'
+    for ll,l in enumerate(labels):
+        #embed()
+        l.spines['right'].set_visible(False)
+        l.spines['left'].set_visible(False)
+        l.spines['top'].set_visible(False)
+        l.spines['bottom'].set_visible(False)
+        l.set_yticks([])
+        l.set_xticks([])
+        l.set_ylabel(titels[ll],labelpad = 15, fontsize = fs_big, fontweight=weight, rotation = 0)
+    hd = 0.3
+    grid1 = gridspec.GridSpecFromSubplotSpec(3,1,subplot_spec=grid[1], wspace=0.02, hspace=0.4)#,
+
+
+    end = ['original']
+
+    color_modul = 'steelblue'
+    color_mpf = 'red'
+
+    y_sum1,moduls = plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,arrows = False,nrs = [0,1,2], title = 'Binary spike trains',xlabel = True,yticks = False)
+
+    grid2 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=grid[2], wspace=0.02, hspace=0.4)#,
+
+    end = ['05']
+    end = ['2']
+    #y_sum = [[]] * len(data)
+    color_modul = 'steelblue'
+    color_mpf = 'red'
+    y_sum2, moduls = plot_single(lim, data, end, data_all, grid2, color_mpf, color_modul,arrows = True, mods = moduls, nrs = [3,4,5], title = '2 ms Gaussian',xlabel = False)
+
+    #for dd, d in enumerate(data):
+        # embed()
+        #embed()
+        #ax[1].set_ylim([0, np.nanmax([y_sum1,y_sum2])])
+        #ax[1, dd].set_ylim([0, 350])
+    plt.subplots_adjust(wspace = 0.4,left = 0.1, right = 0.96,bottom = 0.1)
+
+def plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,mods = [], arrows = True, nr_size = 12, nrs = [0,1,2], xlabel = True,title = '',yticks =True):
+    y_sum = [[]] * len(data)
+    ax = {}
+    moduls = [[]]*len(data)
+    for dd,set in enumerate(data):
+        for ee, e in enumerate(end):
+            d = data_all[data_all['dataset'] == set]
+            x = d['delta_f'] / d['eodf'] + 1
+            #embed()
+            y = d['result_frequency_' + e]
+            y2 = d['result_amplitude_max_' + e]
+            y_sum[dd] = np.nanmax(y)
+            ff = d['delta_f'] / d['eodf'] + 1
+            fe = d['beat_corr']
+            #fig.suptitle(set)
+            grid2 = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=grid1[dd], wspace=0.02, hspace=0.2)  # ,
+            ax[0] = plt.subplot(grid2[0])
+            ax[0].plot(ff, fe, color='grey', linestyle='--')
+            ax[0].plot(x, y, color=color_mpf)
+            ax[0].text(-0.1, 1.1, string.ascii_uppercase[nrs[dd]], transform=ax[0].transAxes,
+                       size=nr_size, weight='bold')
+
+            if dd == 0:
+                plt.title(title)
+            ax[set+e+str(1)] = plt.subplot(grid2[1])
+            #ax[0, dd].set_title(e + ' ms')
+            ax[0].set_xlim([0, 4])
+            #if (e == 2) and xlabel == True:
+            ax[set+e+str(1)].plot(x, y2, color=color_modul,zorder = 2)
+            if arrows == True:
+                plt.fill_between(x, y2,mods[dd], color = 'gainsboro', edgecolor= 'grey',zorder = 1)
+                array = [0.65, 1.65, 2.65]
+                small_arrows = False
+                if small_arrows == True:
+                    for a in range(len(array)):
+                        #embed()
+                        pos = np.argmin(np.abs(np.array(x)-array[a]))
+                        x_present = np.array(x)[pos]
+                        y2 = np.array(y2)
+                        #embed()
+                        pos_change = 2
+                        nr = 25
+                        #embed()
+                        if (np.array(mods[dd])[pos]-y2[pos])>nr:
+                            plt.plot([x_present, x_present],
+                                     [y2[pos] + nr, np.max(np.array(mods[dd])[pos - pos_change:pos + pos_change])],
+                                     color='black')
+                            plt.scatter([x_present],[y2[pos]+nr], marker = 'v',s = 10, color='black')
+            moduls[dd] = y2
+            # ax[1,0].set_ylabel('modulation depth [Hz]')
+            #ax[2, ee].plot(x, y2, color=colors[0])
+            #ax[2, 0].set_ylabel('         modulation depth [Hz]')
+            # ax[1,ee].annotate("", xy=(0.53, 16.83), xytext=(0.53, 17.33), arrowprops=dict(arrowstyle="->"))
+            # ax[1,ee].annotate("", xy=(1.51, 16.83), xytext=(1.51, 17.33), arrowprops=dict(arrowstyle="->"))
+
+            #ax[1, 0].set_xlabel('stimulus frequency [EODf]')
+            if (dd == 2) and xlabel == True:
+                ax[set+e+str(1)].set_xlabel('stimulus frequency [EODf]')
+                ax[0].set_ylabel('MPF [EODf]')
+                ax[set + e + str(1)].set_ylabel('Modulation ')
+            #ax[1, 2].set_xlabel('stimulus frequency [EODf]')
+            #ax[2, 3].set_xlabel('stimulus frequency [EODf]')
+
+            ax[0].spines['right'].set_visible(False)
+            ax[0].spines['top'].set_visible(False)
+            ax[set+e+str(1)].spines['right'].set_visible(False)
+            ax[set+e+str(1)].spines['top'].set_visible(False)
+            #ax[2, ee].spines['right'].set_visible(False)
+            #ax[2, ee].spines['top'].set_visible(False)
+            ax[0].set_xlim([0, 5])
+            ax[set+e+str(1)].set_xlim([0, 5])
+            #plt.tight_layout()
+            # fig.label_axes()
+            ax[set+e+str(1)].set_ylim([0, lim])
+            #ax[0].set_ylim([0, 240])
+            ax[0] = remove_tick_marks(ax[0])
+            if set != data[-1]:
+                ax[set+e+str(1)] = remove_tick_marks(ax[set+e+str(1)])
+            if yticks == True:
+                remove_tick_ymarks(ax[set+e+str(1)])
+                remove_tick_ymarks(ax[0])
+            #ax[0].set_ylim([0, lim])
+    return y_sum, moduls
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']
+    data = ['2019-10-21-aa-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
+    data = ['2019-09-23-ad-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
+    default_settings(data,intermediate_width = 6.7,intermediate_length = 7.5)
+    #fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
+    ax = {}
+    plot_single_cells(ax, data = data)
+    #fig.savefig()
+    plt.savefig('differentcells_trans.pdf')
+    plt.savefig('../highbeats_pdf/differentcells_trans.pdf')
+    # plt.subplots_adjust(left = 0.25)
+    plt.show()
+    #plt.close()
diff --git a/examples.py b/examples.py
new file mode 100644
index 0000000..59f6cea
--- /dev/null
+++ b/examples.py
@@ -0,0 +1,428 @@
+
+
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+#from axes import label_axes, labelaxes_params
+from myfunctions import auto_rows
+#from differentcells import default_settings
+#from differentcells import plot_single_cells
+import matplotlib.gridspec as gridspec
+from functionssimulation import single_stim
+import math
+from functionssimulation import find_times
+from functionssimulation import rectify
+from functionssimulation import global_maxima
+from functionssimulation import integrate_chirp
+from functionssimulation import find_periods
+from myfunctions import default_settings
+from axes import labelaxes_params,label_axes
+from mpl_toolkits.axes_grid1 import host_subplot
+import mpl_toolkits.axisartist as AA
+import string
+
+
+def plot_single_cells(ax,colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1'], var = '05'):
+
+    # labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    end = ['original', '005','05', '2' ]
+    end = [var]
+    y_sum = [[]] * len(data)
+    axis = {}
+    for dd,set in enumerate(data):
+        for ee, e in enumerate(end):
+            d = data_all[data_all['dataset'] == set]
+            eod = d['eodf'].iloc[0]
+            x = d['delta_f'] / d['eodf'] + 1
+            xx = d['delta_f']
+            y = d['result_frequency_' + e]
+            y2 = d['result_amplitude_max_' + e]
+            y_sum[dd] = np.nanmax(y)
+            axis[1] = plt.subplot(ax[0])
+            axis[1].plot(x, y, zorder = 1,color=colors[0])
+            axis[1].set_ylabel('AF [Hz]')
+            axis[1].set_xlim([0, 4])
+            labels = [item.get_text() for item in axis[1].get_xticklabels()]
+            empty_string_labels = [''] * len(labels)
+            axis[1].set_xticklabels(empty_string_labels)
+
+            axis[2] = host_subplot(ax[1], axes_class=AA.Axes)
+            #axis[2] = plt.subplot(ax[1])
+            #host = host_subplot(ax[1], axes_class=AA.Axes)
+            #host.spines['right'].set_visible(False)
+            #host.spines['top'].set_visible(False)
+            #axis[2] = host.twiny()
+            axis[2].plot(xx, y2,  label="Beats [Hz]", zorder = 2,color=colors[0])
+            axis[2].set_ylabel('Modulation ')
+            axis[1].spines['right'].set_visible(False)
+            axis[1].spines['top'].set_visible(False)
+            axis[2].spines['right'].set_visible(False)
+            axis[2].spines['top'].set_visible(False)
+            axis[1].set_xlim([0, np.max(x)])
+            axis[2].set_xlim([-eod, np.max(xx)])
+
+            nr_size = 10
+            axis[2].text(-0.02, 1.1, string.ascii_uppercase[4],
+                                              transform=axis[2].transAxes,
+                                              size=nr_size, weight='bold')
+            axis[1].text(-0.02, 1.1, string.ascii_uppercase[3],
+                                              transform=axis[1].transAxes,
+                                              size=nr_size, weight='bold')
+            axis[3] = axis[2].twiny()
+            axis[3].set_xlabel('EOD multiples')
+            offset = -40
+            new_fixed_axis = axis[3].get_grid_helper().new_fixed_axis
+            axis[3].axis["bottom"] = new_fixed_axis(loc="bottom", axes=axis[3],
+                                                offset=(0,offset))
+            axis[3].spines['right'].set_visible(False)
+            axis[3].spines['top'].set_visible(False)
+            axis[3].axis["bottom"].toggle(all=True)
+            axis[2].set_xlabel("Difference frequency [Hz]")
+            #par2.set_xlim([np.min(xx), np.max(xx)])
+            axis[3].set_xlim([0, np.max(x)])
+            #p1, = host.plot([0, 1, 2], [0, 1, 2], label="Density")
+            #p2, = par1.plot([0, 1, 2], [0, 3, 2], label="Temperature")
+            p3, = axis[3].plot(x, y2,color = 'grey',zorder = 1)
+            #embed()
+            axis[2].set_xticks(np.arange(-eod,np.max(xx),eod/2))
+            #ax['corr'].set_yticks(np.arange(eod_fe[0] - eod_fr, eod_fe[-1] - eod_fr, eod_fr / 2))
+
+    axis[2].set_ylim([0, np.nanmax(y_sum)])
+    plt.subplots_adjust(wspace = 0.4,left = 0.17, right = 0.96,bottom = 0.2)
+    return axis,y2
+
+def plot_beat_corr(ax,lower,beat_corr_col = 'red',df_col = 'pink',ax_nr = 3,multiple = 3):
+    eod_fr = 500
+    eod_fe = np.arange(0, eod_fr * multiple, 5)
+    beats = eod_fe - eod_fr
+    beat_corr = eod_fe % eod_fr
+    beat_corr[beat_corr > eod_fr / 2] = eod_fr - beat_corr[beat_corr > eod_fr / 2]
+    #gs0 = gridspec.GridSpec(3, 1, height_ratios=[4, 1, 1], hspace=0.7)
+
+    #plt.figure(figsize=(4.5, 6))
+    style = 'dotted'
+    color_v = 'black'
+    color_b = 'silver'
+    # plt.subplot(3,1,1)
+    ax['corr'] = plt.subplot(lower[ax_nr])
+    np.max(beats) / eod_fr
+
+    ax['corr'].set_xticks(np.arange((eod_fe[0]-eod_fr)/eod_fr+1, (eod_fe[-1]-eod_fr)/eod_fr+1,(eod_fr/2)/eod_fr+1))
+    ax['corr'].set_yticks(np.arange((eod_fe[0]-eod_fr)/eod_fr+1, (eod_fe[-1]-eod_fr)/eod_fr+1,(eod_fr/2)/eod_fr+1))
+    ax['corr'].set_xticks(np.arange(0,10,0.5))
+    ax['corr'].set_yticks(np.arange(0, 10, 0.5))
+    # plt.axvline(x = -250, Linestyle = style,color = color_v)
+    # plt.axvline(x = 250, Linestyle = style,color = color_v)
+    # plt.axvline(x = 750, Linestyle = style,color = color_v)
+    # plt.axvline(x = 1500, Linestyle = style)
+    # plt.subplot(3,1,2)
+    plt.xlabel('Beats [Hz]')
+    plt.ylabel('Difference frequency [Hz]')
+    #plt.subplot(gs0[1])
+    if beat_corr_col != 'no':
+        plt.plot(beats/eod_fr+1, beat_corr/(eod_fr+1), color=beat_corr_col, alpha = 0.7)
+        plt.ylim([0,np.max(beat_corr/(eod_fr+1))*1.4])
+        plt.xlim([(beats/eod_fr+1)[0],(beats/eod_fr+1)[-1]])
+    if df_col != 'no':
+        plt.plot(beats/eod_fr+1, np.abs(beats)/(eod_fr+1), color=df_col, alpha = 0.7)
+    #plt.axvline(x=-250, Linestyle=style, color=color_v)
+    #plt.axvline(x=250, Linestyle=style, color=color_v)
+    #plt.axvline(x=750, Linestyle=style, color=color_v)
+    plt.xlabel('EOD adjusted beat [Hz]')
+    ax['corr'] .spines['right'].set_visible(False)
+    ax['corr'] .spines['top'].set_visible(False)
+    ax['corr'] .spines['left'].set_visible(True)
+    ax['corr'] .spines['bottom'].set_visible(True)
+    # plt.axvline(x = 1250, Linestyle = style,color = color_v)
+    # plt.axvline(x = 1500, Linestyle = style,color = color_v)
+    mult = np.array(beats) / eod_fr + 1
+    # plt.subplot(3,1,3)
+    plt.xlabel('EOD multiples')
+    plt.ylabel('EOD adj. beat [Hz]', fontsize = 10)
+    plt.grid()
+    #plt.subplot(gs0[2])
+    #plt.plot(mult, beat_corr, color=color_b)
+    # plt.axvline(x = 0, Linestyle = style)
+    #plt.axvline(x=0.5, Linestyle=style, color=color_v)
+    # plt.axvline(x = 1, Linestyle = style)
+    #plt.axvline(x=1.5, Linestyle=style, color=color_v)
+    #plt.axvline(x=2.5, Linestyle=style, color=color_v)
+    #plt.xlabel('EOD multiples')
+    #plt.ylabel('EOD adj. beat [Hz]', fontsize = 10)
+    return ax
+
+def try_resort_automatically():
+    diffs = np.diff(dfs)
+    fast_sampling = dfs[np.concatenate([np.array([True]),diffs <21])]
+    second_derivative = np.diff(np.diff(fast_sampling))
+    first_index = np.concatenate([np.array([False]),second_derivative <0])
+    second_index = np.concatenate([second_derivative > 0,np.array([False])])
+    remaining = fast_sampling[np.concatenate([np.array([True]),second_derivative == 0, np.array([True])])]
+    first = np.arange(0,len(first_index),1)[first_index]
+    second = np.arange(0, len(second_index), 1)[second_index]-1
+    residual = []
+    indeces = []
+    for i in range(len(first)):
+        index = np.arange(first[i],second[i],2)
+        index2 = np.arange(first[i], second[i], 1)
+        indeces.append(index2)
+        residual.append(fast_sampling[index])#first[i]:second[i]:2
+    indeces = np.concatenate(indeces)
+    remaining = fast_sampling[~indeces]
+    residual = np.concatenate(residual)
+    new_dfs = np.sort(np.concatenate([residual, remaining]))
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    data = ['2019-09-23-ad-invivo-1']
+    labelaxes_params(xoffs=1, yoffs=0, labels='A', font=dict(fontweight='bold'))
+    labelaxes_params(xoffs=-6, yoffs=1, labels='A', font=dict(fontweight='bold'))
+
+
+    default_settings(data,intermediate_width = 6.29,intermediate_length = 7.5, ts = 6, ls = 8, fs = 9)
+    fig = plt.figure()
+    #fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
+    ax = {}
+    #ax = plt.subplot(grid[2])
+    data_all = pd.read_pickle('data_beat.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+    eod = d['eodf'].iloc[0]
+    dfs = np.unique(d['df'])
+    #embed()
+    grid = gridspec.GridSpec(2, 4, wspace=0.0, height_ratios=[6, 2], width_ratios=[1,1,0.3,3], hspace=0.2)
+
+    low_nr = 60
+    from_middle = 45 #20
+    example_df = [low_nr- eod,eod / 2 - from_middle - eod,eod - low_nr - eod,low_nr,eod / 2 - from_middle, low_nr + eod]
+    #example_df = [1, eod / 2 - 20 - eod, eod - low_nr - eod, low_nr, eod / 2 - 20, low_nr + eod]
+
+    rows = len(example_df)
+    cols = 1
+    power_raster = gridspec.GridSpecFromSubplotSpec(rows, cols,
+                                                    subplot_spec=grid[0, 0],wspace = 0.05, hspace=0.3)
+    max_p = [[]]*len(example_df)
+    for i in range(len(example_df)):
+
+        power = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1,1.7],hspace = 0.2, wspace = 0.2, subplot_spec = power_raster[i])
+
+        first = ['crimson', 'lightcoral', 'darkviolet']
+        second = ['hotpink', 'deeppink', 'mediumvioletred']
+        third = ['khaki', 'yellow', 'gold']
+        third = ['orange', 'orangered', 'darkred']
+        fourth = ['DarkGreen', 'LimeGreen', 'YellowGreen']
+        fith = ['SkyBlue', 'DeepSkyBlue', 'Blue']
+        colors = np.concatenate([fourth, third, first])
+
+        ax_nr = 0
+        ax['scatter_small'+str(i)] = plt.subplot(power[ax_nr])
+        eod_fr = eod
+        eod_fe = [example_df[i] + eod]
+        e = 0
+        factor = 200
+        sampling = 500 * factor
+        minus_bef = -250
+        plus_bef = -200
+        #minus_bef =  -2100
+        #plus_bef = -100
+        f_max, lims, _ = single_stim(ax, [colors[i]], 1, 1, eod_fr, eod_fe, e, power,delta_t = 0.001, add = 'no',minus_bef =minus_bef, plus_bef = plus_bef,sampling = sampling,
+                                     col_basic = 'silver',col_hline = 'no',labels = False,a_fr=1, ax_nr=ax_nr , phase_zero=[0], shift_phase=0,df_col = 'no',beat_corr_col='no', size=[120], a_fe=0.8)
+
+
+
+    ax['between'] = plt.subplot(grid[0, 2])
+    ax['between'].spines['right'].set_visible(False)
+    ax['between'].spines['top'].set_visible(False)
+    ax['between'].spines['left'].set_visible(False)
+    ax['between'].spines['bottom'].set_visible(False)
+    ax['between'].set_ylim([np.min(dfs), np.max(dfs)])
+    ax['between'].set_xlim([-0.5,30])
+    ax['between'].set_xticks([])
+    ax['between'].set_yticks([])
+    ax['between'].set_ylim(ax['between'].get_ylim()[::-1])
+    nr_size = 10
+
+    ax['scatter'] = plt.subplot(grid[0,1])
+    ax['scatter'].spines['right'].set_visible(False)
+    ax['scatter'].spines['top'].set_visible(False)
+    counter = 0
+    new_dfs = np.concatenate([dfs[0:25], dfs[25:40:2], dfs[40:53:2], dfs[54:-1]])
+    for i in range(len(new_dfs)):
+        spikes = d[d['df'] == new_dfs[i]]['spike_times']
+        counter += 1
+        ll = 0.1
+        ul = 0.3
+        transformed_spikes = spikes.iloc[0]-spikes.iloc[0][0]
+        used_spikes = transformed_spikes[transformed_spikes>ll]
+        used_spikes = used_spikes[used_spikes<ul]*1000
+        ax['scatter'].scatter(used_spikes,np.ones(len(used_spikes))*new_dfs[i],s = 0.2,color = 'silver')
+
+    #plt.gca().invert_yaxis()
+    #ax = plt.gca()
+    ax['scatter'].set_ylim([np.min(dfs),np.max(dfs)])
+    ax['scatter'].set_xlim([ll*1000,ul*1000])
+    ax['scatter'].set_ylabel('Difference frequency [Hz]')
+    ax['scatter'].set_xlabel('Time [ms]')
+    ax['scatter'].set_ylim(ax['scatter'].get_ylim()[::-1])
+    ax['scatter'].text(-0.1, 1.025, string.ascii_uppercase[0], transform=ax['scatter'].transAxes,
+            size= nr_size, weight='bold')
+
+
+
+    #embed()
+    axis = gridspec.GridSpecFromSubplotSpec(2, 1,
+                                             subplot_spec=grid[1,:], wspace=0, hspace=0.5)
+    x = d['df'] / d['eodf'] + 1
+    main_color = 'darkgrey'
+    var = '05'
+    var = 'original'
+    axis,y2 = plot_single_cells(axis,colors = [main_color], data = data,var = var)
+
+
+
+
+    new_dfs = np.concatenate([dfs[0:25], dfs[25:40:2], dfs[40:53:2], dfs[54:-1]])
+    low_nr = 60
+    low = [low_nr- eod, low_nr, low_nr + eod, low_nr + eod * 2, low_nr + eod * 3]
+    high_nr = eod / 2 - 20
+    high = [high_nr - eod, high_nr, high_nr + eod, high_nr + eod * 2, high_nr + eod * 3]
+    high_nr = eod - low_nr
+    low_mirrowed = [high_nr - eod, high_nr, high_nr + eod, high_nr + eod * 2, high_nr + eod * 3]
+    first = ['crimson','lightcoral','darkviolet']
+    second = ['hotpink','deeppink','mediumvioletred']
+    third = ['khaki','yellow','gold']
+    third = ['orange','orangered','darkred']
+    fourth = ['DarkGreen','LimeGreen','YellowGreen']
+    fith = ['SkyBlue','DeepSkyBlue','Blue']
+    colors = np.concatenate([fourth, third,first ])
+    example_df = np.concatenate([first, high, low_mirrowed])
+    #embed()
+    new = np.transpose([low, high, low_mirrowed])
+    example_df = np.concatenate([new[0], new[1]])#new[2]new[2]
+    rows = len(example_df)
+    cols = 1
+    power_raster = gridspec.GridSpecFromSubplotSpec(rows, cols,
+                                                    subplot_spec=grid[0, 3],wspace = 0.05, hspace=0.3)
+    #plt.tight_layout(power_raster)
+    #embed()
+    low_nr = 60
+    from_middle = 45 #20
+    example_df = [low_nr- eod,eod / 2 - from_middle - eod,eod - low_nr - eod,low_nr,eod / 2 - from_middle, low_nr + eod]
+    #example_df = [1, eod / 2 - 20 - eod, eod - low_nr - eod, low_nr, eod / 2 - 20, low_nr + eod]
+
+    max_p = [[]]*len(example_df)
+    for i in range(len(example_df)):
+
+        power = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1,1.7],hspace = 0.2, wspace = 0.2, subplot_spec = power_raster[i])
+
+
+        ax_nr = 0
+        ax['scatter_small'+str(i)] = plt.subplot(power[ax_nr])
+        eod_fr = eod
+        eod_fe = [example_df[i] + eod]
+        e = 0
+        factor = 200
+        sampling = 500 * factor
+        minus_bef = -250
+        plus_bef = -200
+        #minus_bef =  -2100
+        #plus_bef = -100
+        f_max, lims, _ = single_stim(ax, [colors[i]], 1, 1, eod_fr, eod_fe, e, power,delta_t = 0.001, add = 'no',minus_bef =minus_bef, plus_bef = plus_bef,sampling = sampling,
+                                     col_basic = 'silver',col_hline = 'no',labels = False,a_fr=1, ax_nr=ax_nr , phase_zero=[0], shift_phase=0,df_col = 'no',beat_corr_col='no', size=[120], a_fe=0.8)
+
+        ax['between'].scatter(0.12,example_df[i],zorder=2,  s=25,marker = '<',color=colors[i])
+        ax['between'].scatter(0.12,example_df[i], zorder=2, s=25,marker = '<',color=colors[i])
+        ll = np.abs(plus_bef)
+        ul = np.abs(minus_bef)
+        df = new_dfs[np.argmin(np.abs(new_dfs - example_df[i]))]
+        spikes = d[d['df'] == df]['spike_times']
+
+        tranformed_spikes = spikes.iloc[0]*1000-spikes.iloc[0][0]*1000
+        used_spikes = tranformed_spikes[tranformed_spikes>ll]
+        used_spikes = used_spikes[used_spikes<ul]
+        used_spikes = used_spikes-used_spikes[0]
+        ax['scatter_small'+str(i)].scatter((used_spikes),np.ones(len(used_spikes))*-2,zorder=2,s = 2,marker = '|',color = 'black')#color = colors[i]
+        ax['scatter_small'+str(i)].set_ylim([-2.5,2.2])
+
+        ax['power'+str(i)] = plt.subplot(power[1])
+        nfft = 4096  #
+        sampling_rate = 40000
+        # embed()
+
+        pp = [[]]*len(spikes)
+        for s in range(len(spikes)):
+            new_spikes = list(map(int, (spikes.iloc[s] - spikes.iloc[s][0]) * sampling_rate))
+            array = np.zeros(new_spikes[-1] + 2)
+            array[new_spikes] = 1
+            array = array*sampling_rate
+            if var == '05':
+                window05 = 0.0005 * sampling_rate
+                array = gaussian_filter(array, sigma=window05) * sampling_rate
+            pp[s], f = ml.psd(array - np.mean(array), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #embed()
+        p = np.mean(pp, axis = 0)
+        #embed()
+        diff = d['eodf'].iloc[0] * (df / d['eodf'].iloc[0] - int(df / d['eodf'].iloc[0]))
+        if diff > d['eodf'].iloc[0] * 0.5:
+            diff = diff - d['eodf'].iloc[0] * 0.5
+        plt.plot(f, p, zorder=1 ,color=main_color)
+        max_p[i] = np.max(p)
+        ax['power'+str(i)].scatter(f[np.argmax(p[f < 0.5 * eod])],max(p[f < 0.5 * eod]),zorder=2,color = colors[i], s = 25)
+        ax['power' + str(i)].scatter(f[f == f[np.argmin(np.abs(f-eod))]], p[f == f[np.argmin(np.abs(f-eod))]]*0.90, zorder=2,
+                                     s=25, color = 'darkgrey',edgecolor = 'black')
+        ax['power' + str(i)].axvline(x = eod/2, color = 'black', linestyle = 'dashed', lw = 0.5)
+        plt.xlim([-40, 1600])
+        axis[3].scatter(example_df[i]/(eod)+1, np.sqrt(np.max(p[f < 0.5 * eod])*np.abs(f[0]-f[1])),zorder=3, s=20,marker = 'o',color=colors[i])
+        axis[1].scatter(example_df[i]/(eod)+1,f[np.argmax(p[f < 0.5 * eod])],zorder=2, s=20,marker = 'o',color=colors[i])
+
+
+        if i != rows-1:
+            #ax['power'+str(i)].set_xticks([])
+            #ax['scatter_small'].set_xticks([])
+            labels = [item.get_text() for item in ax['scatter_small'+str(i)].get_xticklabels()]
+            empty_string_labels = [''] * len(labels)
+            ax['scatter_small'+str(i)].set_xticklabels(empty_string_labels)
+            labels = [item.get_text() for item in ax['power'+str(i)].get_xticklabels()]
+            empty_string_labels = [''] * len(labels)
+            ax['power'+str(i)].set_xticklabels(empty_string_labels)
+        else:
+            ax['power'+str(i)].set_xlabel('Frequency [Hz]')
+            ax['scatter_small'+str(i)].set_xlabel('Time [ms]')
+        ax['power' + str(i)].set_yticks([])
+        ax['power'+str(i)].spines['left'].set_visible(False)
+        ax['scatter_small'+str(i)].spines['left'].set_visible(False)
+        ax['scatter_small'+str(i)].set_yticks([])
+        ax['power'+str(i)].spines['right'].set_visible(False)
+        ax['power'+str(i)].spines['top'].set_visible(False)
+        ax['scatter_small'+str(i)].spines['right'].set_visible(False)
+        ax['scatter_small'+str(i)].spines['top'].set_visible(False)
+    for i in range(len(example_df)):
+        ax['power'+str(i)].set_ylim([0,np.max(max_p)])
+    ax['power'+str(0)].text(-0.1, 1.1, string.ascii_uppercase[2], transform=ax['power'+str(0)].transAxes,
+            size= nr_size, weight='bold')
+    ax['scatter_small'+str(0)].text(-0.1, 1.1, string.ascii_uppercase[1], transform=ax['scatter_small'+str(0)].transAxes,
+            size= nr_size, weight='bold')
+    plt.subplots_adjust(left = 0.11, bottom = 0.18, top = 0.94)
+    #fig.label_axes()
+    #fig.label_axes()
+    #embed()
+    #grid.format(
+    #    xlabel='xlabel', ylabel='ylabel', suptitle=titles[mode],
+    #    abc=True, abcloc='ul',
+    #    grid=False, xticks=25, yticks=5)
+    plt.savefig('singlecellexample5.pdf')
+    plt.savefig('../highbeats_pdf/singlecellexample5.pdf')
+    # plt.subplots_adjust(left = 0.25)
+    plt.show()
+    #plt.close()
diff --git a/localmaxima.py b/localmaxima.py
new file mode 100644
index 0000000..4b0fcdd
--- /dev/null
+++ b/localmaxima.py
@@ -0,0 +1,185 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed
+import matplotlib as matplotlib
+import math
+import scipy.integrate as integrate
+from scipy import signal
+from scipy.interpolate import interp1d
+from scipy.interpolate import CubicSpline
+import scipy as sp
+import pickle
+from scipy.spatial import distance
+from myfunctions import *
+import time
+from matplotlib import gridspec
+from matplotlib_scalebar.scalebar import ScaleBar
+import matplotlib.mlab as ml
+import scipy.integrate as si
+import pandas as pd
+from functionssimulation import find_times
+from functionssimulation import find_periods
+from functionssimulation import integrate_chirp
+from functionssimulation import rectify, find_beats,find_dev
+from functionssimulation import global_maxima, find_lm, conv
+
+def snip(left_c,right_c,e,g,sampling, deviation_s,d,eod_fr, a_fr, eod_fe,phase_zero,p, size,s, sigma,a_fe,deviation,beat_corr):
+    time, time_cut, cut = find_times(left_c[g], right_c[g], sampling, deviation_s[d])
+    eod_fish_r, period_fish_r, period_fish_e = find_periods(time, eod_fr, a_fr, eod_fe, e)
+    #embed()
+    eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
+    eod_rec_down, eod_rec_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify
+    eod_overlayed_chirp = (eod_fish_r + eod_fe_chirp)[cut:-cut]
+
+    maxima_values, maxima_index, maxima_interp = global_maxima(period_fish_e, period_fish_r,
+                                                               eod_rec_up[cut:-cut])  # global maxima
+    index_peaks, value_peaks, peaks_interp = find_lm(eod_rec_up[cut:-cut])  # local maxima
+    middle_conv, eod_conv_down, eod_conv_up, eod_conv_downsampled = conv(eod_fr,sampling, cut, deviation[d], eod_rec_up,
+                                                                         eod_rec_down)  # convolve
+    eod_fish_both = integrate_chirp(a_fe, time, eod_fe[e] - eod_fr, phase_zero[p], size[s], sigma)
+    am_corr_full = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
+                                   sigma)  # indirect am calculation
+    _, time_fish, cut_f = find_times(left_c[g], right_c[g], eod_fr, deviation_s[d])  # downsampled through fish EOD
+    am_corr_ds = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
+    am_df_ds = integrate_chirp(a_fe, time_fish, eod_fe[e] - eod_fr, phase_zero[p], size[s],
+                              sigma)  # indirect am calculation
+    return time_cut, eod_conv_up, am_corr_full, peaks_interp, maxima_interp, am_corr_ds,am_df_ds,eod_fish_both,eod_overlayed_chirp
+
+
+
+def power_func(bef_c, aft_c, win, deviation_s, sigma, sampling, d_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation, show_figure = False, plot_dist = False, save = False):
+    results = [[]]*7
+    for d in range(len(deviation)):
+        bef_c = 0.3
+        aft_c = -0.1
+        left_c, right_c, left_b, right_b, period_distance_c, period_distance_b, _, period, to_cut,exclude,consp_needed,deli,interval = period_calc(
+            [float('Inf')]*len(beat_corr), 50, win, deviation_s[d], sampling, beat_corr, bef_c, aft_c, 'stim')
+        save_n = win
+        for s in range(len(size)):
+            for p in range(len(phase_zero)):
+                beats = eod_fe - eod_fr
+                for e in range(len(eod_fe)):
+                    left_b = [-0.3*sampling]*len(beat_corr)
+                    right_b = [-0.1 * sampling]*len(beat_corr)
+                    time_b, conv_b, am_corr_b, peaks_interp_b, maxima_interp_b,am_corr_ds_b, am_df_ds_b, am_df_b,eod_overlayed_chirp = snip(left_b, right_b,e,e,
+                                                                                                sampling, deviation_s,
+                                                                                                d, eod_fr, a_fr, eod_fe,
+                                                                                                phase_zero, p, size, s,
+                                                                                                sigma, a_fe, deviation,
+                                                                                                beat_corr)
+                    #time_c, conv_c, am_corr_c, peaks_interp_c, maxima_interp_c,am_corr_ds_c, am_df_ds_c, am_df_c = snip(left_c, right_c,e,e,
+                    #                                                                            sampling, deviation_s,
+                    #                                                                            d, eod_fr, a_fr, eod_fe,
+                    #                                                                            phase_zero, p, size, s,
+                    #                                                                            sigma, a_fe, deviation,
+                    #                                                                            beat_corr)
+
+                    #embed()
+                    nfft = 4096
+                    name = ['conv','df','df ds','corr','corr ds','global max','local max']
+                    var = [conv_b,am_df_b,am_df_ds_b, am_corr_b, am_corr_ds_b,maxima_interp_b,peaks_interp_b ]
+                    samp = [sampling,sampling,eod_fr,sampling,eod_fr,sampling,sampling]
+                    #pp, f = ml.psd(eod_overlayed_chirp - np.mean(eod_overlayed_chirp), Fs=sampling, NFFT=nfft, noverlap=nfft / 2)
+                    for i in range(len(results)):
+
+
+                        plot = False
+                        pp, f = ml.psd(var[i] - np.mean(var[i]), Fs=samp[i], NFFT=nfft,
+                                       noverlap=nfft / 2)
+
+                        if plot == True:
+                            plt.figure()
+                            plt.subplot(1,2,1)
+                            plt.plot(var[i])
+                            plt.title(name[i])
+                            plt.subplot(1, 2, 2)
+                            plt.plot(f,pp)
+                            #plt.xlim([0,2000])
+                            plt.show()
+                        #print(results)
+                        #embed()
+                        if type(results[i]) != dict:
+                            results[i] = {}
+                            results[i]['type'] = name[i]
+                            #embed()
+                            results[i]['f'] = list([f[np.argmax(pp[f < 0.5 * eod_fr])]])
+                            results[i]['amp'] = list([np.sqrt(si.trapz(pp, f, np.abs(f[1]-f[0])))])
+                            results[i]['max'] = list([np.sqrt(np.max(pp[f < 0.5 * eod_fr])*np.abs(f[1]-f[0]))])
+                        else:
+                            results[i]['f'].extend(list([f[np.argmax(pp[f < 0.5 *eod_fr])]]))
+                            #embed()
+                            results[i]['amp'].extend(list([np.sqrt(si.trapz(pp, f, np.abs(f[1]-f[0])))]))
+                            results[i]['max'].extend(list([np.sqrt(np.max(pp[f < 0.5 * eod_fr]) * np.abs(f[1] - f[0]))]))
+                #if save:
+                #    results = pd.DataFrame(results)
+                #    results.to_pickle('../data/power_simulation.pkl')
+                #    np.save('../data/Ramona/power_simulation.npy', results)
+    return results
+
+
+def plot_power(results):
+    plt.rcParams['figure.figsize'] = (3, 3)
+    plt.rcParams["legend.frameon"] = False
+    colors = ['black', 'magenta', 'pink', 'orange', 'moccasin', 'red', 'green', 'silver']
+    colors = ['red','pink']
+    results = [results[5]]
+    fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True)
+    all_max = [[]] * len(results)
+    #embed()
+    for i in range(len(results)):
+        #embed()
+        #ax[0].set_ylabel(results[i]['type'], rotation=0, labelpad=40, color=colors[i])
+        ax[0].plot(beats / eod_fr + 1, np.array(results[i]['f']) / eod_fr + 1, color=colors[i])
+        # plt.title(results['type'][i])
+        ax[1].plot(beats / eod_fr + 1, np.array(results[i]['amp']), color=colors[0])
+        ax[1].plot(beats / eod_fr + 1, np.array(results[i]['max']), color=colors[1])
+        #ax[2].plot(beats / eod_fr + 1, np.array(results[i]['amp']), color=colors[i])
+        all_max[i] = np.max(np.array(results[i]['amp']))
+    #for i in range(len(results)):
+    #    ax[2].set_ylim([0, np.max(all_max)])
+    plt.subplots_adjust(left=0.25)
+    #ii, jj = np.shape(ax)
+    ax[0].set_ylabel('MPF')
+    ax[1].set_ylabel('Modulation depth')
+    #ax[0, 2].set_title('Modulation depth (same scale)')
+    for i in range(len(ax)):
+        ax[1].set_xlabel('EOD multiples')
+        ax[i].spines['right'].set_visible(False)
+        ax[i].spines['top'].set_visible(False)
+    plt.subplots_adjust(bottom = 0.2)
+    plt.savefig('localmaxima.pdf')
+    plt.savefig('../highbeats_pdf/localmaxima.pdf')
+    plt.show()
+
+delta_t = 0.014  # ms
+interest_interval = delta_t * 1.2
+bef_c = interest_interval / 2
+aft_c = interest_interval / 2
+sigma = delta_t / math.sqrt((2 * math.log(10)))  # width of the chirp
+size = [120]  # maximal frequency excursion during chirp / 60 or 100 here
+phase_zero = [0]  # phase when the chirp occured (vary later) / zero at the peak o a beat cycle
+phase_zero = np.arange(0,2*np.pi,2*np.pi/10)
+eod_fr = 500  # eod fish reciever
+a_fr = 1  # amplitude fish reciever
+amplitude = a_fe = 0.2  # amplitude fish emitter
+factor = 200
+sampling = eod_fr * factor
+sampling_fish = 500
+#start = 0
+#end = 2500
+#step = 10
+start = 510
+end = 3500
+step = 500
+win = 'w2'
+d = 1
+x = [ 1.5, 2.5,0.5,]
+x = [ 1.5]
+time_range = 200 * delta_t
+deviation_ms, deviation_s, deviation_dp = find_dev(x, sampling)
+start = 5
+end = 2500
+step = 25
+eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
+results = power_func( bef_c, aft_c, 'w2', deviation_s, sigma, sampling, deviation_ms, beat_corr, size, [phase_zero[0]], delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp, show_figure = True, plot_dist = False, save = True)
+plot_power(results)
\ No newline at end of file
diff --git a/rotated_singlethree.py b/rotated_singlethree.py
new file mode 100644
index 0000000..afee2e1
--- /dev/null
+++ b/rotated_singlethree.py
@@ -0,0 +1,522 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+from myfunctions import auto_rows
+from functionssimulation import  default_settings
+import matplotlib.gridspec as gridspec
+from myfunctions import remove_tick_marks
+
+def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000):
+
+    #nfft = 4096
+    #trial_cut = 0.1
+    #freq_step = sampling_rate / nfft
+    data_cell = data[data['dataset'] == d]#
+    dfs = np.unique(data_cell['df'])
+    df_here = dfs[np.argmin(np.abs(dfs - wish_df))]
+    dfs310 = data_cell[data_cell['df'] == df_here]
+    #pp = [[]]*len(dfs310)
+    pp = []
+    ppp = []
+    trial_cut = 0.1
+    for i in range(len(dfs310)):
+        duration = dfs310.iloc[i]['durations']
+        #cut_vec = np.arange(0, duration, trial_cut)
+        cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        #for j, cut in enumerate(cut_vec):
+        #    # print(j)
+        #    spike_times = dfs310.iloc[i]['spike_times']
+        #    spikes = spike_times - spike_times[0]
+        #    spikes_cut = spikes[(spikes > cut) & (spikes < cut_vec[j + 1])]
+        #    if cut == cut_vec[-2]:
+        #        #counter_cut += 1
+        #        break
+        #    if len(spikes_cut) < 10:
+        #        #counter_spikes += 1
+        #        break
+        #    spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #    spikes_idx = np.round((spikes_cut - trial_cut * j) * sampling_rate)
+        #    for spike in spikes_idx:
+        #        spikes_mat[int(spike)] = 1#
+        #
+        #    #spikes_mat = np.zeros(int(spikes[-1]* sampling_rate + 5))
+        #    #spikes_idx = np.round((spikes) * sampling_rate)
+        #    #for spike in spikes_idx:
+        #    #    spikes_mat[int(spike)] = 1
+        #    spikes_mat = spikes_mat * sampling_rate
+        #    if type(window) != str:
+        #        spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+        #        # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+        #        # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        #    else:
+        #        smoothened = spikes_mat * 1
+        #    nfft = 4096
+        #    p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #    pp.append(p)
+        spike_times = dfs310.iloc[i]['spike_times']
+        if len(spike_times) < 3:
+            counter_spikes += 1
+            break
+
+        spikes = spike_times - spike_times[0]
+        spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        if len(spikes_cut) < 3:
+            counter_cut += 1
+            break
+        spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5))
+        spikes_idx = np.round((spikes) * sampling_rate)
+        for spike in spikes_idx:
+            spikes_mat[int(spike)] = 1
+        spikes_mat = spikes_mat * sampling_rate
+        if type(window) != str:
+            spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+            # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+            # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        else:
+            spikes_mat = spikes_mat*1
+        nfft = 4096
+        p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        ppp.append(p)
+        #spike_times = data_cell.iloc[i]['spike_times']#
+
+        #if len(spike_times) < 3:
+        #    counter_spikes += 1
+        #    break
+
+        #spikes = spike_times - spike_times[0]
+
+        # cut trial into snippets of 100 ms
+        #cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+
+        #if len(spikes_cut) < 3:
+        #    counter_cut += 1
+        #    break
+        #spikes_new = spikes_cut - spikes_cut[0]
+
+        #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
+        # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #spikes_idx = np.round((spikes_new) * sampling_rate)
+        #for spike in spikes_idx:
+        #    spikes_mat[int(spike)] = 1
+        #spikes_mat = spikes_mat * sampling_rate
+
+        #nfft = 4096
+        #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #ppp.append(p)
+    #p_mean = np.mean(pp,axis = 0)
+    p_mean2 = np.mean(ppp, axis=0)
+    #ref = (np.max(p_mean2))
+    #
+    db = 10 * np.log10(p_mean2 / np.max(p_mean2))
+    #ref = (np.max(p_mean2))
+    #db2 = 10 * np.log10(p_mean2 / ref)
+    #embed()
+    return df_here,p_mean2,f,db
+
+
+
+def plot_example_ps(grid,colors = ['brown'],fc = 'lightgrey',line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+
+
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl')
+    #iter = np.unique(data['dataset'])
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+
+
+        # fig.suptitle(d, labelpad = 25)
+        #print(d)
+        ax = {}
+
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        #ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        #ax[0].legend( loc=(0,1),
+        #                 ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+
+        #ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 0,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0] = remove_tick_marks(ax[0])
+        #ax[0].set_ylim([0, 2000])
+
+        wide = 2
+        #embed()
+        nr = 1
+        for i in range(len(sigma)):
+            ax[i+nr] = plt.subplot(grid[i+nr])
+            plot_filter(ax, i+nr, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            ax[i+nr].set_ylim([0, eodf*1.5])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))-1].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(len(df)):
+            ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot(p[nr], f[nr], color=colors[0])
+        ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
+        ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(0, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    #plt.plot([np.min(p),np.max(p)],[eodf,eodf], color = 'red')
+    #embed()
+    #ax[nr].plot([0]*5)
+    #ax[nr].plot([1000]*5)
+    # ax[0].fill_between( [max(p[0])]*len(f[1]),f[0], facecolor='lightgrey', edgecolor='grey')
+
+
+    ax[ax_nr].set_ylim([0, eodf * 1.5])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_filter(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot( p4[array_nr],f, color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    ax[ax_nr].scatter( now_height+440,  np.abs(df), marker = '>', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f,  facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df),
+               color=color_df, marker='o', linestyle='None')
+    ax[ax_nr].plot(
+               np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf,
+               color=color_stim, marker='o', linestyle='None')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+
+def plot_amp(ax, mean1, dev,name = 'amp',nr = 1):
+    np.unique(mean1['type'])
+    all_means = mean1[mean1['type'] == name +' mean']
+    original = all_means[all_means['dev'] == 'original']
+    #m005 = all_means[all_means['dev'] == '005']
+    m05 = all_means[all_means['dev'] == '05']
+    m2 = all_means[all_means['dev'] == '2']
+    # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True)
+    versions = [original, m05, m2] #m005,
+    for i in range(len(versions)):
+        keys = [k for k in versions[i]][2::]
+        try:
+            data = np.array(versions[i][keys])[0]
+        except:
+            break
+        axis = np.arange(0, len(data), 1)
+        axis_new = axis * 1
+        similarity = [keys, data]
+        sim = np.argsort(similarity[0])
+        # similarity[sim]
+        all_means = mean1[mean1['type'] == name+' std']
+        std = all_means[all_means['dev'] == dev[i]]
+        std = np.array(std[keys])[0]
+        #ax[1, 1].set_ylabel('Modulation depth')
+        #ax[nr,i].set_title(dev[i] + ' ms')
+        all_means = mean1[mean1['type'] == name+' 95']
+        std95 = all_means[all_means['dev'] == dev[i]]
+        std95 = np.array(std95[keys])[0]
+        all_means = mean1[mean1['type'] == name+' 05']
+        std05 = all_means[all_means['dev'] == dev[i]]
+        std05 = np.array(std05[keys])[0]
+        ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]),
+                              color='gainsboro')
+        ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]),
+                              color='darkgrey')
+
+        # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf')
+        ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black')
+        # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %')
+    #embed()
+    return ax
+
+def create_beat_corr(hz_range, eod_fr):
+    beat_corr = hz_range%eod_fr
+    beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
+    return beat_corr
+
+
+def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    half_page_width = 7.9 / inch_factor
+    intermediate_width = 12 / inch_factor
+    whole_page_width = 16 * 2 / inch_factor
+
+    small_length = 6 / inch_factor
+    intermediate_length = 12 * 1.5 / inch_factor
+    max_length = 25 / inch_factor
+    whole_page_width = 6.7
+    intermediate_length = 3.7
+
+    #plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    # load data for plot
+
+    # data1 = pd.read_csv('ma_allcells_unsmoothed.csv')
+    # data2 = pd.read_csv('ma_allcells_05.csv')
+    # data3 = pd.read_csv('ma_allcells_2.csv')
+    # data_tob = pd.read_csv('ma_toblerone.csv')
+
+    # smothed = df_beat[df_beat['dev'] == 'original']
+    # data1 = smothed[smothed['type'] == 'amp']
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,2,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+    #fig.suptitle(set)
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
+        ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
+        #embed()
+        eod = d['eodf'].iloc[0]
+        ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        #embed()
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+        #ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
+
+        #ax[0].set_ylabel('MPF [EODf]')
+        #ax[0].set_ylabel('Modulation ')
+        #ax[0, dd].set_title(e + ' ms')
+        ax[0].set_xlim([0, 4])
+
+        ax[1] = plt.subplot(plots[1])
+        if e != sigma[-1]:
+            ax[1] = remove_tick_marks(ax[1])
+        else:
+            ax[1].set_ylabel('Modulation ')
+        ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        if e == 'original':
+            whole_height = height
+            whole_height = np.max(y2)
+        if (e != 'whole') and (e!= 'original'):
+            ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            color = 'black')
+        #ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        #if i == len(sigma)-1:
+        #    ax[1].set_xlabel('stimulus frequency [EODf]')
+        ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        ax[1].spines['top'].set_visible(False)
+        ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        #ax[1].spines['top'].set_visible(False)
+        ax[0].set_xlim([0, 5])
+        #embed()
+        ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+    #fig.tight_layout()
+    # fig.label_axes()
+
+
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    #fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
+    trans = False
+    sigma = [0.0005, 0.002]  # 0.00005,0.00025,
+    if trans == True:
+        col = 1
+        row = 2
+        col_small = len(sigma)+2
+        row_small = 1
+        l = 6
+        t = 'horizontal'
+        wd = [1]
+        hd = [1,2.5]
+        left = 1
+        right = 0
+    else:
+        col = 2
+        row = 1
+        row_small = len(sigma)+1
+        col_small = 1
+        t = 'vertical'
+        l = 9
+        wd = [2, 4]
+        hd = [1]
+        left = 0
+        right = 1
+    default_settings(data, intermediate_width=7.4, intermediate_length=9, ts=6, ls=10, fs=10)
+    grid = gridspec.GridSpec(row, col, right = 0.95,wspace=0.02,top = 0.95,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
+    wish_df = 240#324
+    color_eod = 'darkgreen'#'orange'
+    color_stim = 'navy'
+    color_df = 'orange'#'green'
+    colors = ['brown']
+    colors_mpf = ['red']
+    colors_mod = ['steelblue']
+    line_col = 'black'
+    plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],fc = 'lightgrey',line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
+    #plt.show()
+    #embed()
+    #fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
+    #embed()
+    plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025','whole',
+    plt.savefig('rotatedps_singlethree.pdf')
+    plt.savefig('../highbeats_pdf/rotatedps_singlethree.pdf')
+    #plt.savefig('.pdf')
+    plt.show()
+    # plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
+    # plt.show()
+    # plt.close()
+    # embed()
+    # plot_single_tublerones() # original beat_activity
\ No newline at end of file
diff --git a/rotatedps_single.py b/rotatedps_single.py
new file mode 100644
index 0000000..e68e40c
--- /dev/null
+++ b/rotatedps_single.py
@@ -0,0 +1,521 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+from myfunctions import auto_rows
+from functionssimulation import  default_settings
+import matplotlib.gridspec as gridspec
+from myfunctions import remove_tick_marks
+
+def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000):
+
+    #nfft = 4096
+    #trial_cut = 0.1
+    #freq_step = sampling_rate / nfft
+    data_cell = data[data['dataset'] == d]#
+    dfs = np.unique(data_cell['df'])
+    df_here = dfs[np.argmin(np.abs(dfs - wish_df))]
+    dfs310 = data_cell[data_cell['df'] == df_here]
+    #pp = [[]]*len(dfs310)
+    pp = []
+    ppp = []
+    trial_cut = 0.1
+    for i in range(len(dfs310)):
+        duration = dfs310.iloc[i]['durations']
+        #cut_vec = np.arange(0, duration, trial_cut)
+        cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        #for j, cut in enumerate(cut_vec):
+        #    # print(j)
+        #    spike_times = dfs310.iloc[i]['spike_times']
+        #    spikes = spike_times - spike_times[0]
+        #    spikes_cut = spikes[(spikes > cut) & (spikes < cut_vec[j + 1])]
+        #    if cut == cut_vec[-2]:
+        #        #counter_cut += 1
+        #        break
+        #    if len(spikes_cut) < 10:
+        #        #counter_spikes += 1
+        #        break
+        #    spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #    spikes_idx = np.round((spikes_cut - trial_cut * j) * sampling_rate)
+        #    for spike in spikes_idx:
+        #        spikes_mat[int(spike)] = 1#
+        #
+        #    #spikes_mat = np.zeros(int(spikes[-1]* sampling_rate + 5))
+        #    #spikes_idx = np.round((spikes) * sampling_rate)
+        #    #for spike in spikes_idx:
+        #    #    spikes_mat[int(spike)] = 1
+        #    spikes_mat = spikes_mat * sampling_rate
+        #    if type(window) != str:
+        #        spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+        #        # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+        #        # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        #    else:
+        #        smoothened = spikes_mat * 1
+        #    nfft = 4096
+        #    p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #    pp.append(p)
+        spike_times = dfs310.iloc[i]['spike_times']
+        if len(spike_times) < 3:
+            counter_spikes += 1
+            break
+
+        spikes = spike_times - spike_times[0]
+        spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        if len(spikes_cut) < 3:
+            counter_cut += 1
+            break
+        spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5))
+        spikes_idx = np.round((spikes) * sampling_rate)
+        for spike in spikes_idx:
+            spikes_mat[int(spike)] = 1
+        spikes_mat = spikes_mat * sampling_rate
+        if type(window) != str:
+            spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+            # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+            # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        else:
+            spikes_mat = spikes_mat*1
+        nfft = 4096
+        p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        ppp.append(p)
+        #spike_times = data_cell.iloc[i]['spike_times']#
+
+        #if len(spike_times) < 3:
+        #    counter_spikes += 1
+        #    break
+
+        #spikes = spike_times - spike_times[0]
+
+        # cut trial into snippets of 100 ms
+        #cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+
+        #if len(spikes_cut) < 3:
+        #    counter_cut += 1
+        #    break
+        #spikes_new = spikes_cut - spikes_cut[0]
+
+        #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
+        # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #spikes_idx = np.round((spikes_new) * sampling_rate)
+        #for spike in spikes_idx:
+        #    spikes_mat[int(spike)] = 1
+        #spikes_mat = spikes_mat * sampling_rate
+
+        #nfft = 4096
+        #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #ppp.append(p)
+    #p_mean = np.mean(pp,axis = 0)
+    p_mean2 = np.mean(ppp, axis=0)
+    #ref = (np.max(p_mean2))
+    #
+    db = 10 * np.log10(p_mean2 / np.max(p_mean2))
+    #ref = (np.max(p_mean2))
+    #db2 = 10 * np.log10(p_mean2 / ref)
+    #embed()
+    return df_here,p_mean2,f,db
+
+
+
+def plot_example_ps(grid,colors = ['brown'],fc = 'lightgrey',line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+
+
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl')
+    #iter = np.unique(data['dataset'])
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+
+
+        # fig.suptitle(d, labelpad = 25)
+        #print(d)
+        ax = {}
+
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        #ax[0].set_ylim([0, 2000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            ax[i+2].set_ylim([0, eodf*1.5])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(len(df)+1):
+            ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot(p[nr], f[nr], color=colors[0])
+        ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
+        ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(0, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    #plt.plot([np.min(p),np.max(p)],[eodf,eodf], color = 'red')
+    #embed()
+    #ax[nr].plot([0]*5)
+    #ax[nr].plot([1000]*5)
+    # ax[0].fill_between( [max(p[0])]*len(f[1]),f[0], facecolor='lightgrey', edgecolor='grey')
+
+
+    ax[ax_nr].set_ylim([0, eodf * 1.5])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_filter(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot( p4[array_nr],f, color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    ax[ax_nr].scatter( now_height+440,  np.abs(df), marker = '>', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f,  facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df),
+               color=color_df, marker='o', linestyle='None')
+    ax[ax_nr].plot(
+               np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf,
+               color=color_stim, marker='o', linestyle='None')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+
+def plot_amp(ax, mean1, dev,name = 'amp',nr = 1):
+    np.unique(mean1['type'])
+    all_means = mean1[mean1['type'] == name +' mean']
+    original = all_means[all_means['dev'] == 'original']
+    #m005 = all_means[all_means['dev'] == '005']
+    m05 = all_means[all_means['dev'] == '05']
+    m2 = all_means[all_means['dev'] == '2']
+    # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True)
+    versions = [original, m05, m2] #m005,
+    for i in range(len(versions)):
+        keys = [k for k in versions[i]][2::]
+        try:
+            data = np.array(versions[i][keys])[0]
+        except:
+            break
+        axis = np.arange(0, len(data), 1)
+        axis_new = axis * 1
+        similarity = [keys, data]
+        sim = np.argsort(similarity[0])
+        # similarity[sim]
+        all_means = mean1[mean1['type'] == name+' std']
+        std = all_means[all_means['dev'] == dev[i]]
+        std = np.array(std[keys])[0]
+        #ax[1, 1].set_ylabel('Modulation depth')
+        #ax[nr,i].set_title(dev[i] + ' ms')
+        all_means = mean1[mean1['type'] == name+' 95']
+        std95 = all_means[all_means['dev'] == dev[i]]
+        std95 = np.array(std95[keys])[0]
+        all_means = mean1[mean1['type'] == name+' 05']
+        std05 = all_means[all_means['dev'] == dev[i]]
+        std05 = np.array(std05[keys])[0]
+        ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]),
+                              color='gainsboro')
+        ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]),
+                              color='darkgrey')
+
+        # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf')
+        ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black')
+        # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %')
+    #embed()
+    return ax
+
+def create_beat_corr(hz_range, eod_fr):
+    beat_corr = hz_range%eod_fr
+    beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
+    return beat_corr
+
+
+def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    half_page_width = 7.9 / inch_factor
+    intermediate_width = 12 / inch_factor
+    whole_page_width = 16 * 2 / inch_factor
+
+    small_length = 6 / inch_factor
+    intermediate_length = 12 * 1.5 / inch_factor
+    max_length = 25 / inch_factor
+    whole_page_width = 6.7
+    intermediate_length = 3.7
+
+    #plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    # load data for plot
+
+    # data1 = pd.read_csv('ma_allcells_unsmoothed.csv')
+    # data2 = pd.read_csv('ma_allcells_05.csv')
+    # data3 = pd.read_csv('ma_allcells_2.csv')
+    # data_tob = pd.read_csv('ma_toblerone.csv')
+
+    # smothed = df_beat[df_beat['dev'] == 'original']
+    # data1 = smothed[smothed['type'] == 'amp']
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,2,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+    #fig.suptitle(set)
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
+        ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
+        #embed()
+        eod = d['eodf'].iloc[0]
+        ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        #embed()
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+        #ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
+
+        #ax[0].set_ylabel('MPF [EODf]')
+        #ax[0].set_ylabel('Modulation ')
+        #ax[0, dd].set_title(e + ' ms')
+        ax[0].set_xlim([0, 4])
+
+        ax[1] = plt.subplot(plots[1])
+        if e != sigma[-1]:
+            ax[1] = remove_tick_marks(ax[1])
+        else:
+            ax[1].set_ylabel('Modulation ')
+        ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        if e == 'original':
+            whole_height = height
+            whole_height = np.max(y2)
+        if (e != 'whole') and (e!= 'original'):
+            ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            color = 'black')
+        #ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        #if i == len(sigma)-1:
+        #    ax[1].set_xlabel('stimulus frequency [EODf]')
+        ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        ax[1].spines['top'].set_visible(False)
+        ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        #ax[1].spines['top'].set_visible(False)
+        ax[0].set_xlim([0, 5])
+        #embed()
+        ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+    #fig.tight_layout()
+    # fig.label_axes()
+
+
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    #fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
+    trans = False
+    sigma = [0.0005, 0.002]  # 0.00005,0.00025,
+    if trans == True:
+        col = 1
+        row = 2
+        col_small = len(sigma)+2
+        row_small = 1
+        l = 6
+        t = 'horizontal'
+        wd = [1]
+        hd = [1,2.5]
+        left = 1
+        right = 0
+    else:
+        col = 2
+        row = 1
+        row_small = len(sigma)+2
+        col_small = 1
+        t = 'vertical'
+        l = 9
+        wd = [2, 4]
+        hd = [1]
+        left = 0
+        right = 1
+    default_settings(data, intermediate_width=7.4, intermediate_length=9, ts=6, ls=10, fs=10)
+    grid = gridspec.GridSpec(row, col, right = 0.95,wspace=0.02,top = 0.97,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
+    wish_df = 240#324
+    color_eod = 'darkgreen'#'orange'
+    color_stim = 'navy'
+    color_df = 'orange'#'green'
+    colors = ['brown']
+    colors_mpf = ['red']
+    colors_mod = ['steelblue']
+    line_col = 'black'
+    plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],fc = 'lightgrey',line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
+    #plt.show()
+    #embed()
+    #fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
+    #embed()
+    plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
+    plt.savefig('rotatedps_single.pdf')
+    plt.savefig('../highbeats_pdf/rotatedps_single.pdf')
+    #plt.savefig('.pdf')
+    plt.show()
+    # plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
+    # plt.show()
+    # plt.close()
+    # embed()
+    # plot_single_tublerones() # original beat_activity
\ No newline at end of file
diff --git a/rotatedps_singleall.py b/rotatedps_singleall.py
new file mode 100644
index 0000000..6607686
--- /dev/null
+++ b/rotatedps_singleall.py
@@ -0,0 +1,744 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+from myfunctions import auto_rows
+from functionssimulation import  default_settings
+import matplotlib.gridspec as gridspec
+from myfunctions import remove_tick_marks
+
+def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000):
+
+    #nfft = 4096
+    #trial_cut = 0.1
+    #freq_step = sampling_rate / nfft
+    data_cell = data[data['dataset'] == d]#
+    dfs = np.unique(data_cell['df'])
+    df_here = dfs[np.argmin(np.abs(dfs - wish_df))]
+    dfs310 = data_cell[data_cell['df'] == df_here]
+    #pp = [[]]*len(dfs310)
+    pp = []
+    ppp = []
+    trial_cut = 0.1
+    for i in range(len(dfs310)):
+        duration = dfs310.iloc[i]['durations']
+        #cut_vec = np.arange(0, duration, trial_cut)
+        cut_vec = np.arange(0, duration, trial_cut)
+        spike_times = dfs310.iloc[i]['spike_times']
+        if len(spike_times) < 3:
+            counter_spikes += 1
+            break
+
+        spikes = spike_times - spike_times[0]
+        spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        if len(spikes_cut) < 3:
+            counter_cut += 1
+            break
+        spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5))
+        spikes_idx = np.round((spikes) * sampling_rate)
+        for spike in spikes_idx:
+            spikes_mat[int(spike)] = 1
+        spikes_mat = spikes_mat * sampling_rate
+        if type(window) != str:
+            spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+            # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+            # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        else:
+            spikes_mat = spikes_mat*1
+        nfft = 4096
+        p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        ppp.append(p)
+        #spike_times = data_cell.iloc[i]['spike_times']#
+
+        #if len(spike_times) < 3:
+        #    counter_spikes += 1
+        #    break
+
+        #spikes = spike_times - spike_times[0]
+
+        # cut trial into snippets of 100 ms
+        #cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+
+        #if len(spikes_cut) < 3:
+        #    counter_cut += 1
+        #    break
+        #spikes_new = spikes_cut - spikes_cut[0]
+
+        #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
+        # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #spikes_idx = np.round((spikes_new) * sampling_rate)
+        #for spike in spikes_idx:
+        #    spikes_mat[int(spike)] = 1
+        #spikes_mat = spikes_mat * sampling_rate
+
+        #nfft = 4096
+        #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #ppp.append(p)
+    #p_mean = np.mean(pp,axis = 0)
+    p_mean2 = np.mean(ppp, axis=0)
+    #ref = (np.max(p_mean2))
+    #
+    db = 10 * np.log10(p_mean2 / np.max(p_mean2))
+    #ref = (np.max(p_mean2))
+    #db2 = 10 * np.log10(p_mean2 / ref)
+    #embed()
+    return df_here,p_mean2,f,db
+
+def plot_example_ps_trans(grid,colors = ['brown'],line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+        ax = {}
+        fc = 'lightgrey'
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps_trans(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+        ax[0].set_xlim([0, 1000])
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        ax[1].set_xlim([0, 1000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter_trans(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            #ax[i+2].set_ylim([0, eodf*1.5])
+            ax[i + 2].set_xlim([0, 1000])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(1,len(df)+1):
+            ax[i].axvline(x = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_example_ps(grid,colors = ['brown'],line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+
+
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl')
+    #iter = np.unique(data['dataset'])
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+
+
+        # fig.suptitle(d, labelpad = 25)
+        #print(d)
+        ax = {}
+        fc = 'lightgrey'
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        #ax[0].set_ylim([0, 2000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            ax[i+2].set_ylim([0, eodf*1.5])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(1,len(df)+1):
+            ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot( f[nr],p[nr], color=colors[0])
+        ax[ax_nr].fill_between( [f[0][-1],f[0][-1]],[np.min(p), np.max(p)], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(df[nr] + eodf,p[nr][int((df[nr] + eodf) / stepsize) + 1],  color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(eodf - 1,np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale,  color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(f[nr][f[nr]<eodf/2],p[nr][f[nr]<eodf/2],  color=colors[0])
+        ax[ax_nr].plot(f[nr][f[nr] > eodf / 2],np.zeros(len(f[nr][f[nr] > eodf / 2])),  color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between( [eodf/2,eodf/2],[np.min(p),np.max(p)], color=fc,edgecolor=ec)
+        ax[ax_nr].plot( df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(df[nr] + eodf,0,  color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(eodf - 1,0,  color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    #ax[ax_nr].set_ylim([0, eodf * 1.5])
+    #ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot(p[nr], f[nr], color=colors[0])
+        ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
+        ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(0, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    ax[ax_nr].set_ylim([0, eodf * 1.5])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+def plot_filter_trans(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot(f,  p4[array_nr],color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    ax[ax_nr].plot([np.abs(df), np.abs(df)],[prev_height, now_height+440], color = 'black')
+    ax[ax_nr].scatter( np.abs(df), now_height+440,  marker = 'v', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(f, max(p4[0]) * gauss3 ** 2, facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot( eodf, np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale,color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( abs(df),np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,
+               color=color_df, marker='o', linestyle='None')
+    #ax[ax_nr].plot(df + eodf,
+    #           np.max(df + eodf,p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,
+    #           color=color_stim, marker='o', linestyle='None')
+    #ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_filter(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot( p4[array_nr],f, color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    #ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    #ax[ax_nr].scatter( now_height+440,  np.abs(df), marker = '>', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f,  facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df),
+               color=color_df, marker='o', linestyle='None')
+    ax[ax_nr].plot(
+               np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf,
+               color=color_stim, marker='o', linestyle='None')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+
+def plot_amp(ax, mean1, dev,name = 'amp',nr = 1):
+    np.unique(mean1['type'])
+    all_means = mean1[mean1['type'] == name +' mean']
+    original = all_means[all_means['dev'] == 'original']
+    #m005 = all_means[all_means['dev'] == '005']
+    m05 = all_means[all_means['dev'] == '05']
+    m2 = all_means[all_means['dev'] == '2']
+    # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True)
+    versions = [original, m05, m2] #m005,
+    for i in range(len(versions)):
+        keys = [k for k in versions[i]][2::]
+        try:
+            data = np.array(versions[i][keys])[0]
+        except:
+            break
+        axis = np.arange(0, len(data), 1)
+        axis_new = axis * 1
+        similarity = [keys, data]
+        sim = np.argsort(similarity[0])
+        # similarity[sim]
+        all_means = mean1[mean1['type'] == name+' std']
+        std = all_means[all_means['dev'] == dev[i]]
+        std = np.array(std[keys])[0]
+        #ax[1, 1].set_ylabel('Modulation depth')
+        #ax[nr,i].set_title(dev[i] + ' ms')
+        all_means = mean1[mean1['type'] == name+' 95']
+        std95 = all_means[all_means['dev'] == dev[i]]
+        std95 = np.array(std95[keys])[0]
+        all_means = mean1[mean1['type'] == name+' 05']
+        std05 = all_means[all_means['dev'] == dev[i]]
+        std05 = np.array(std05[keys])[0]
+        ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]),
+                              color='gainsboro')
+        ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]),
+                              color='darkgrey')
+
+        # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf')
+        ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black')
+        # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %')
+    #embed()
+    return ax
+
+def create_beat_corr(hz_range, eod_fr):
+    beat_corr = hz_range%eod_fr
+    beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
+    return beat_corr
+
+def plot_mean_cells_modul( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,1,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+        ax[0] = plt.subplot(plots[0])
+        eod = d['eodf'].iloc[0]
+
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        ax[0].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        if e == 'original':
+            whole_height = height
+            whole_height = np.max(y2)
+        if (e != 'whole') and (e!= 'original'):
+            ax[0].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            ax[0].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            color = 'black')
+        ax[0].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        if i == len(sigma)-1:
+            ax[0].set_xlabel('stimulus frequency [EODf]')
+        ax[0].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        #ax[1].spines['top'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        #ax[0].set_xlim([0, 5])
+        #embed()
+        #ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+
+
+def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    half_page_width = 7.9 / inch_factor
+    intermediate_width = 12 / inch_factor
+    whole_page_width = 16 * 2 / inch_factor
+
+    small_length = 6 / inch_factor
+    intermediate_length = 12 * 1.5 / inch_factor
+    max_length = 25 / inch_factor
+    whole_page_width = 6.7
+    intermediate_length = 3.7
+
+    #plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    # load data for plot
+
+    # data1 = pd.read_csv('ma_allcells_unsmoothed.csv')
+    # data2 = pd.read_csv('ma_allcells_05.csv')
+    # data3 = pd.read_csv('ma_allcells_2.csv')
+    # data_tob = pd.read_csv('ma_toblerone.csv')
+
+    # smothed = df_beat[df_beat['dev'] == 'original']
+    # data1 = smothed[smothed['type'] == 'amp']
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,1,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+    #fig.suptitle(set)
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        if e != 'whole':
+            ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
+            ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
+        ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
+        #embed()
+        eod = d['eodf'].iloc[0]
+
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        #embed()
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+        #ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
+
+        #ax[0].set_ylabel('MPF [EODf]')
+        #ax[0].set_ylabel('Modulation ')
+        #ax[0, dd].set_title(e + ' ms')
+        ax[0].set_xlim([0, 4])
+
+        #ax[1] = plt.subplot(plots[1])
+        #if e != sigma[-1]:
+        #    ax[1] = remove_tick_marks(ax[1])
+        #ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        #height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        #if e == 'original':
+        #    whole_height = height
+        #    whole_height = np.max(y2)
+        #if (e != 'whole') and (e!= 'original'):
+            #ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            #ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            #color = 'black')
+        #ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        #if i == len(sigma)-1:
+        #    ax[1].set_xlabel('stimulus frequency [EODf]')
+        #ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        #ax[1].spines['top'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        #ax[1].spines['top'].set_visible(False)
+        ax[0].set_xlim([0, 5])
+        #embed()
+        ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+    #fig.tight_layout()
+    # fig.label_axes()
+
+
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    #fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
+    trans = False
+    sigma = [0.0005, 0.002]  # 0.00005,0.00025,
+    if trans == True:
+        col = 1
+        row = 2
+        col_small = len(sigma)+2
+        row_small = 1
+        l = 6
+        t = 'horizontal'
+        wd = [1]
+        hd = [1,2.5]
+        left = 1
+        right = 0
+    else:
+        col = 2
+        row = 1
+        row_small = len(sigma)+2
+        col_small = 1
+        t = 'vertical'
+        l = 9
+        wd = [3, 4]
+        hd = [1]
+        left = 0
+        right = 1
+    default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
+    grid1 = gridspec.GridSpec(1, 2, left=0.15, right=0.95, wspace=0.25, height_ratios=hd, width_ratios=[1,1],
+                             hspace=0.2)  # ,
+
+    grid2 = gridspec.GridSpecFromSubplotSpec(row, col, subplot_spec=grid1[0], wspace=0.02,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid2[0,0], wspace=0.15, hspace=hs)
+    wish_df = 240#324
+    color_eod = 'darkgreen'#'orange'
+    color_stim = 'navy'
+    color_df = 'orange'#'green'
+    colors = ['brown']
+    colors_mpf = ['brown']
+    colors_mod = ['steelblue']
+    line_col = 'black'
+    plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
+    #plt.show()
+    #embed()
+    #fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid2[left,right], wspace=0.25, hspace=hs)
+    #embed()
+    plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
+
+
+    grid3 = gridspec.GridSpecFromSubplotSpec(row, col,subplot_spec=grid1[1],  wspace=0.02, height_ratios=hd,
+                                             width_ratios=wd, hspace=0.2)  # ,
+
+
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
+                                            subplot_spec=grid3[0, 0], wspace=0.15, hspace=hs)
+    wish_df = 240  # 324
+    color_eod = 'darkgreen'  # 'orange'
+    color_stim = 'navy'
+    color_df = 'orange'  # 'green'
+    colors = ['brown']
+    colors_mpf = ['brown']
+    colors_mod = ['steelblue']
+    line_col = 'black'
+    plot_example_ps_trans(axis, input=['2019-10-21-aa-invivo-1'], line_col=line_col, colors=colors, sigma=sigma,
+                    wish_df=wish_df, color_eod=color_eod, color_stim=color_stim, color_df=color_df)
+    # plt.show()
+    # embed()
+    # fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
+                                            subplot_spec=grid3[left, right], wspace=0.25, hspace=hs)
+    # embed()
+    plot_mean_cells_modul(axis, data=['2019-10-21-aa-invivo-1'], line_col=line_col, lw=1.23, size=22,
+                    sigma=['whole', 'original', '05', '2'], colors=colors_mpf, wish_df=wish_df, color_eod=color_eod,
+                    color_df=color_df, color_modul=colors_mod)  # '005','025',
+
+    plt.savefig('rotatedps_singleall.pdf')
+    plt.savefig('../highbeats_pdf/rotatedps_singleall.pdf')
+    #plt.savefig('.pdf')
+    plt.show()
+    # plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
+    # plt.show()
+    # plt.close()
+    # embed()
+    # plot_single_tublerones() # original beat_activity
\ No newline at end of file
diff --git a/rotatedps_singleamodul.py b/rotatedps_singleamodul.py
new file mode 100644
index 0000000..9b7230b
--- /dev/null
+++ b/rotatedps_singleamodul.py
@@ -0,0 +1,722 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+from myfunctions import auto_rows
+from functionssimulation import  default_settings
+import matplotlib.gridspec as gridspec
+from myfunctions import remove_tick_marks
+
+def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000):
+
+    #nfft = 4096
+    #trial_cut = 0.1
+    #freq_step = sampling_rate / nfft
+    data_cell = data[data['dataset'] == d]#
+    dfs = np.unique(data_cell['df'])
+    df_here = dfs[np.argmin(np.abs(dfs - wish_df))]
+    dfs310 = data_cell[data_cell['df'] == df_here]
+    #pp = [[]]*len(dfs310)
+    pp = []
+    ppp = []
+    trial_cut = 0.1
+    for i in range(len(dfs310)):
+        duration = dfs310.iloc[i]['durations']
+        #cut_vec = np.arange(0, duration, trial_cut)
+        cut_vec = np.arange(0, duration, trial_cut)
+        spike_times = dfs310.iloc[i]['spike_times']
+        if len(spike_times) < 3:
+            counter_spikes += 1
+            break
+
+        spikes = spike_times - spike_times[0]
+        spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        if len(spikes_cut) < 3:
+            counter_cut += 1
+            break
+        spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5))
+        spikes_idx = np.round((spikes) * sampling_rate)
+        for spike in spikes_idx:
+            spikes_mat[int(spike)] = 1
+        spikes_mat = spikes_mat * sampling_rate
+        if type(window) != str:
+            spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+            # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+            # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        else:
+            spikes_mat = spikes_mat*1
+        nfft = 4096
+        p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        ppp.append(p)
+        #spike_times = data_cell.iloc[i]['spike_times']#
+
+        #if len(spike_times) < 3:
+        #    counter_spikes += 1
+        #    break
+
+        #spikes = spike_times - spike_times[0]
+
+        # cut trial into snippets of 100 ms
+        #cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+
+        #if len(spikes_cut) < 3:
+        #    counter_cut += 1
+        #    break
+        #spikes_new = spikes_cut - spikes_cut[0]
+
+        #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
+        # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #spikes_idx = np.round((spikes_new) * sampling_rate)
+        #for spike in spikes_idx:
+        #    spikes_mat[int(spike)] = 1
+        #spikes_mat = spikes_mat * sampling_rate
+
+        #nfft = 4096
+        #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #ppp.append(p)
+    #p_mean = np.mean(pp,axis = 0)
+    p_mean2 = np.mean(ppp, axis=0)
+    #ref = (np.max(p_mean2))
+    #
+    db = 10 * np.log10(p_mean2 / np.max(p_mean2))
+    #ref = (np.max(p_mean2))
+    #db2 = 10 * np.log10(p_mean2 / ref)
+    #embed()
+    return df_here,p_mean2,f,db
+
+def plot_example_ps_trans(grid,colors = ['brown'],line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+        ax = {}
+        fc = 'lightgrey'
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps_trans(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+        ax[0].set_xlim([0, 1000])
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        ax[1].set_xlim([0, 1000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter_trans(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            #ax[i+2].set_ylim([0, eodf*1.5])
+            ax[i + 2].set_xlim([0, 1000])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))].set_xlabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_ylabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(1,len(df)+1):
+            ax[i].axvline(x = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_example_ps(grid,colors = ['brown'],line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+
+
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl')
+    #iter = np.unique(data['dataset'])
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+
+
+        # fig.suptitle(d, labelpad = 25)
+        #print(d)
+        ax = {}
+        fc = 'lightgrey'
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        #ax[0].set_ylim([0, 2000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            ax[i+2].set_ylim([0, eodf*1.5])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(1,len(df)+1):
+            ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot( f[nr],p[nr], color=colors[0])
+        ax[ax_nr].fill_between( [f[0][-1],f[0][-1]],[np.min(p), np.max(p)], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(df[nr] + eodf,p[nr][int((df[nr] + eodf) / stepsize) + 1],  color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(eodf - 1,np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale,  color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(f[nr][f[nr]<eodf/2],p[nr][f[nr]<eodf/2],  color=colors[0])
+        ax[ax_nr].plot(f[nr][f[nr] > eodf / 2],np.zeros(len(f[nr][f[nr] > eodf / 2])),  color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between( [eodf/2,eodf/2],[np.min(p),np.max(p)], color=fc,edgecolor=ec)
+        ax[ax_nr].plot( df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(df[nr] + eodf,0,  color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(eodf - 1,0,  color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    #ax[ax_nr].set_ylim([0, eodf * 1.5])
+    #ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot(p[nr], f[nr], color=colors[0])
+        ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
+        ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(0, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    ax[ax_nr].set_ylim([0, eodf * 1.5])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+def plot_filter_trans(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot(f,  p4[array_nr],color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    ax[ax_nr].plot([np.abs(df), np.abs(df)],[prev_height, now_height+440], color = 'black')
+    ax[ax_nr].scatter( np.abs(df), now_height+440,  marker = 'v', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(f, max(p4[0]) * gauss3 ** 2, facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot( eodf, np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale,color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( abs(df),np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,
+               color=color_df, marker='o', linestyle='None')
+    #ax[ax_nr].plot(df + eodf,
+    #           np.max(df + eodf,p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,
+    #           color=color_stim, marker='o', linestyle='None')
+    #ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_filter(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot( p4[array_nr],f, color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    #ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    #ax[ax_nr].scatter( now_height+440,  np.abs(df), marker = '>', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f,  facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df),
+               color=color_df, marker='o', linestyle='None')
+    ax[ax_nr].plot(
+               np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf,
+               color=color_stim, marker='o', linestyle='None')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+
+def plot_amp(ax, mean1, dev,name = 'amp',nr = 1):
+    np.unique(mean1['type'])
+    all_means = mean1[mean1['type'] == name +' mean']
+    original = all_means[all_means['dev'] == 'original']
+    #m005 = all_means[all_means['dev'] == '005']
+    m05 = all_means[all_means['dev'] == '05']
+    m2 = all_means[all_means['dev'] == '2']
+    # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True)
+    versions = [original, m05, m2] #m005,
+    for i in range(len(versions)):
+        keys = [k for k in versions[i]][2::]
+        try:
+            data = np.array(versions[i][keys])[0]
+        except:
+            break
+        axis = np.arange(0, len(data), 1)
+        axis_new = axis * 1
+        similarity = [keys, data]
+        sim = np.argsort(similarity[0])
+        # similarity[sim]
+        all_means = mean1[mean1['type'] == name+' std']
+        std = all_means[all_means['dev'] == dev[i]]
+        std = np.array(std[keys])[0]
+        #ax[1, 1].set_ylabel('Modulation depth')
+        #ax[nr,i].set_title(dev[i] + ' ms')
+        all_means = mean1[mean1['type'] == name+' 95']
+        std95 = all_means[all_means['dev'] == dev[i]]
+        std95 = np.array(std95[keys])[0]
+        all_means = mean1[mean1['type'] == name+' 05']
+        std05 = all_means[all_means['dev'] == dev[i]]
+        std05 = np.array(std05[keys])[0]
+        ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]),
+                              color='gainsboro')
+        ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]),
+                              color='darkgrey')
+
+        # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf')
+        ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black')
+        # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %')
+    #embed()
+    return ax
+
+def create_beat_corr(hz_range, eod_fr):
+    beat_corr = hz_range%eod_fr
+    beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
+    return beat_corr
+
+def plot_mean_cells_modul( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,1,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+        ax[0] = plt.subplot(plots[0])
+        eod = d['eodf'].iloc[0]
+
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        ax[0].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        if e == 'original':
+            whole_height = height
+            whole_height = np.max(y2)
+        if (e != 'whole') and (e!= 'original'):
+            ax[0].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            ax[0].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            color = 'black')
+        ax[0].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        if i == len(sigma)-1:
+            ax[0].set_xlabel('stimulus frequency [EODf]')
+        ax[0].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        #ax[1].spines['top'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        #ax[0].spines['left'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        #ax[0].set_xlim([0, 5])
+        #embed()
+        #ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+
+
+def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    half_page_width = 7.9 / inch_factor
+    intermediate_width = 12 / inch_factor
+    whole_page_width = 16 * 2 / inch_factor
+
+    small_length = 6 / inch_factor
+    intermediate_length = 12 * 1.5 / inch_factor
+    max_length = 25 / inch_factor
+    whole_page_width = 6.7
+    intermediate_length = 3.7
+
+    #plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    # load data for plot
+
+    # data1 = pd.read_csv('ma_allcells_unsmoothed.csv')
+    # data2 = pd.read_csv('ma_allcells_05.csv')
+    # data3 = pd.read_csv('ma_allcells_2.csv')
+    # data_tob = pd.read_csv('ma_toblerone.csv')
+
+    # smothed = df_beat[df_beat['dev'] == 'original']
+    # data1 = smothed[smothed['type'] == 'amp']
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,1,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+    #fig.suptitle(set)
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        if e != 'whole':
+            ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
+            ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
+        ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
+        #embed()
+        eod = d['eodf'].iloc[0]
+
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        #embed()
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+        #ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
+
+        #ax[0].set_ylabel('MPF [EODf]')
+        #ax[0].set_ylabel('Modulation ')
+        #ax[0, dd].set_title(e + ' ms')
+        ax[0].set_xlim([0, 4])
+
+        #ax[1] = plt.subplot(plots[1])
+        #if e != sigma[-1]:
+        #    ax[1] = remove_tick_marks(ax[1])
+        #ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        #height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        #if e == 'original':
+        #    whole_height = height
+        #    whole_height = np.max(y2)
+        #if (e != 'whole') and (e!= 'original'):
+            #ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            #ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            #color = 'black')
+        #ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        #if i == len(sigma)-1:
+        #    ax[1].set_xlabel('stimulus frequency [EODf]')
+        #ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        #ax[1].spines['top'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        #ax[1].spines['top'].set_visible(False)
+        ax[0].set_xlim([0, 5])
+        #embed()
+        ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+    #fig.tight_layout()
+    # fig.label_axes()
+
+
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    #fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
+    trans = False
+    sigma = [0.0005, 0.002]  # 0.00005,0.00025,
+    if trans == True:
+        col = 1
+        row = 2
+        col_small = len(sigma)+2
+        row_small = 1
+        l = 6
+        t = 'horizontal'
+        wd = [1]
+        hd = [1,2.5]
+        left = 1
+        right = 0
+    else:
+        col = 2
+        row = 1
+        row_small = len(sigma)+2
+        col_small = 1
+        t = 'vertical'
+        l = 9
+        wd = [3, 4]
+        hd = [1]
+        left = 0
+        right = 1
+    default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
+    grid1 = gridspec.GridSpec(1, 1, left=0.15, right=0.95, wspace=0.25, height_ratios=hd,
+                             hspace=0.2)  # ,
+
+    grid3 = gridspec.GridSpecFromSubplotSpec(row, col,subplot_spec=grid1[0],  wspace=0.35, height_ratios=hd,
+                                             width_ratios=wd, hspace=0.2)  # ,
+
+
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
+                                            subplot_spec=grid3[0, 0], wspace=0.15, hspace=hs)
+    wish_df = 240  # 324
+    color_eod = 'darkgreen'  # 'orange'
+    color_stim = 'navy'
+    color_df = 'orange'  # 'green'
+    colors = ['brown']
+    colors_mpf = ['brown']
+    colors_mod = ['brown']
+    line_col = 'black'
+    plot_example_ps_trans(axis, input=['2019-10-21-aa-invivo-1'], line_col=line_col, colors=colors, sigma=sigma,
+                    wish_df=wish_df, color_eod=color_eod, color_stim=color_stim, color_df=color_df)
+    # plt.show()
+    # embed()
+    # fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
+                                            subplot_spec=grid3[left, right], wspace=0.25, hspace=hs)
+    # embed()
+    plot_mean_cells_modul(axis, data=['2019-10-21-aa-invivo-1'], line_col=line_col, lw=1.23, size=22,
+                    sigma=['whole', 'original', '05', '2'], colors=colors_mpf, wish_df=wish_df, color_eod=color_eod,
+                    color_df=color_df, color_modul=colors_mod)  # '005','025',
+
+    plt.savefig('rotatedps_singleamodul.pdf')
+    plt.savefig('../highbeats_pdf/rotatedps_singleamodul.pdf')
+    #plt.savefig('.pdf')
+    plt.show()
+    # plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
+    # plt.show()
+    # plt.close()
+    # embed()
+    # plot_single_tublerones() # original beat_activity
\ No newline at end of file
diff --git a/rotatedps_singlesingle.py b/rotatedps_singlesingle.py
new file mode 100644
index 0000000..96a4842
--- /dev/null
+++ b/rotatedps_singlesingle.py
@@ -0,0 +1,521 @@
+import nixio as nix
+import os
+from IPython import embed
+#from utility import *
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.mlab as ml
+import scipy.integrate as si
+from scipy.ndimage import gaussian_filter
+from IPython import embed
+from myfunctions import *
+from myfunctions import auto_rows
+from functionssimulation import  default_settings
+import matplotlib.gridspec as gridspec
+from myfunctions import remove_tick_marks
+
+def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000):
+
+    #nfft = 4096
+    #trial_cut = 0.1
+    #freq_step = sampling_rate / nfft
+    data_cell = data[data['dataset'] == d]#
+    dfs = np.unique(data_cell['df'])
+    df_here = dfs[np.argmin(np.abs(dfs - wish_df))]
+    dfs310 = data_cell[data_cell['df'] == df_here]
+    #pp = [[]]*len(dfs310)
+    pp = []
+    ppp = []
+    trial_cut = 0.1
+    for i in range(len(dfs310)):
+        duration = dfs310.iloc[i]['durations']
+        #cut_vec = np.arange(0, duration, trial_cut)
+        cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        #for j, cut in enumerate(cut_vec):
+        #    # print(j)
+        #    spike_times = dfs310.iloc[i]['spike_times']
+        #    spikes = spike_times - spike_times[0]
+        #    spikes_cut = spikes[(spikes > cut) & (spikes < cut_vec[j + 1])]
+        #    if cut == cut_vec[-2]:
+        #        #counter_cut += 1
+        #        break
+        #    if len(spikes_cut) < 10:
+        #        #counter_spikes += 1
+        #        break
+        #    spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #    spikes_idx = np.round((spikes_cut - trial_cut * j) * sampling_rate)
+        #    for spike in spikes_idx:
+        #        spikes_mat[int(spike)] = 1#
+        #
+        #    #spikes_mat = np.zeros(int(spikes[-1]* sampling_rate + 5))
+        #    #spikes_idx = np.round((spikes) * sampling_rate)
+        #    #for spike in spikes_idx:
+        #    #    spikes_mat[int(spike)] = 1
+        #    spikes_mat = spikes_mat * sampling_rate
+        #    if type(window) != str:
+        #        spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+        #        # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+        #        # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        #    else:
+        #        smoothened = spikes_mat * 1
+        #    nfft = 4096
+        #    p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #    pp.append(p)
+        spike_times = dfs310.iloc[i]['spike_times']
+        if len(spike_times) < 3:
+            counter_spikes += 1
+            break
+
+        spikes = spike_times - spike_times[0]
+        spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+        if len(spikes_cut) < 3:
+            counter_cut += 1
+            break
+        spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5))
+        spikes_idx = np.round((spikes) * sampling_rate)
+        for spike in spikes_idx:
+            spikes_mat[int(spike)] = 1
+        spikes_mat = spikes_mat * sampling_rate
+        if type(window) != str:
+            spikes_mat = gaussian_filter(spikes_mat, sigma=window)
+            # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate
+            # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate
+        else:
+            spikes_mat = spikes_mat*1
+        nfft = 4096
+        p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        ppp.append(p)
+        #spike_times = data_cell.iloc[i]['spike_times']#
+
+        #if len(spike_times) < 3:
+        #    counter_spikes += 1
+        #    break
+
+        #spikes = spike_times - spike_times[0]
+
+        # cut trial into snippets of 100 ms
+        #cut_vec = np.arange(0, duration, trial_cut)
+        #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)]
+
+        #if len(spikes_cut) < 3:
+        #    counter_cut += 1
+        #    break
+        #spikes_new = spikes_cut - spikes_cut[0]
+
+        #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
+        # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1)
+        #spikes_idx = np.round((spikes_new) * sampling_rate)
+        #for spike in spikes_idx:
+        #    spikes_mat[int(spike)] = 1
+        #spikes_mat = spikes_mat * sampling_rate
+
+        #nfft = 4096
+        #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
+        #ppp.append(p)
+    #p_mean = np.mean(pp,axis = 0)
+    p_mean2 = np.mean(ppp, axis=0)
+    #ref = (np.max(p_mean2))
+    #
+    db = 10 * np.log10(p_mean2 / np.max(p_mean2))
+    #ref = (np.max(p_mean2))
+    #db2 = 10 * np.log10(p_mean2 / ref)
+    #embed()
+    return df_here,p_mean2,f,db
+
+
+
+def plot_example_ps(grid,colors = ['brown'],line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    sampling_rate = 40000
+    #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
+
+
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 6
+    #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl')
+    #iter = np.unique(data['dataset'])
+    iter = ['2019-05-07-by-invivo-1']
+    iter = ['2019-09-23-ad-invivo-1']
+    iter = input
+    for cell in iter:
+        data = pd.read_pickle('data_beat.pkl')
+        beat_results = pd.read_pickle('beat_results_smoothed.pkl')
+        #embed()
+        eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0])
+        df = [[]] * (len(sigma) + 1)
+        p = [[]] * (len(sigma) + 1)
+        f = [[]] * (len(sigma) + 1)
+        db = [[]] * (len(sigma) + 1)
+        sigmaf = [[]] * (len(sigma) + 1)
+        gauss = [[]] * (len(sigma) + 1)
+
+        df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate)
+        for i in range(len(sigma)):
+            df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window =  sigma[i]*sampling_rate,sampling_rate = sampling_rate)
+            sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i])
+            gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2)))
+        db = 'no'
+        stepsize = f[0][1] - f[0][0]
+        if db == 'db':
+            p = db
+
+
+        # fig.suptitle(d, labelpad = 25)
+        #print(d)
+        ax = {}
+        fc = 'lightgrey'
+        ec = 'grey'
+        #fc = 'moccasin'
+        #ec = 'wheat'
+        scale = 1
+        ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[0].legend( loc=(0,1),
+                         ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1),
+
+        ax[0] = remove_tick_marks(ax[0])
+        ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale,  ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+        ax[1] = remove_tick_marks(ax[1])
+        #ax[0].set_ylim([0, 2000])
+
+        wide = 2
+        #embed()
+        for i in range(len(sigma)):
+            ax[i+2] = plt.subplot(grid[i+2])
+            plot_filter(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
+            ax[i+2].set_ylim([0, eodf*1.5])
+        ax[2] = remove_tick_marks(ax[2])
+        #embed()
+        #if db == 'db':
+        #    ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000]
+        #else:
+        #    ax[0].set_ylim([ 0,np.max([p])])
+        ax[int(len(df))-1].set_ylabel('frequency [Hz]')
+        # ax[1].set_ylabel(r'power [Hz$^2$/Hz]')
+        #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0])
+
+
+        #print(df[3])
+        for i in range(len(df)+1):
+            ax[i].spines['right'].set_visible(False)
+            ax[i].spines['top'].set_visible(False)
+        cols = grid.ncols
+        rows = grid.nrows
+        ax[int(len(df))].set_xlabel(' power spectral density [Hz²/Hz]')
+        #ax[2].set_ylabel('Hz²/Hz')
+        #ax[3].set_ylabel('Hz²/Hz')
+        #ax[0].set_ylabel('Hz²/Hz')
+        for i in range(1,len(df)+1):
+            ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
+        plt.tight_layout()
+        #embed()
+        #fig.label_axes()
+
+def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, ax_nr = 0,nr = 0, filter = 'original',  scale = 1,  color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',):
+    ax[ax_nr] = plt.subplot(grid[ax_nr])
+
+    if filter == 'whole':
+        #ax[nr].set_facecolor('lightgrey')
+        ax[ax_nr].plot(p[nr], f[nr], color=colors[0])
+        ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o', linestyle='None', label='Df')
+        ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus')
+        ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf')  # = '+str(int(eodf))+' Hz')
+
+    elif filter == 'original':
+        #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey')
+        #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey')
+        ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
+        ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
+        #embed()
+        ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec)
+        ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0],
+                       color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
+        ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o',
+                       linestyle='None',
+                       label='stimulus',zorder = 2)#,edgecolors = 'black'
+        ax[ax_nr].plot(0, eodf - 1, color=color_eod,
+                   marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
+
+    #plt.plot([np.min(p),np.max(p)],[eodf,eodf], color = 'red')
+    #embed()
+    #ax[nr].plot([0]*5)
+    #ax[nr].plot([1000]*5)
+    # ax[0].fill_between( [max(p[0])]*len(f[1]),f[0], facecolor='lightgrey', edgecolor='grey')
+
+
+    ax[ax_nr].set_ylim([0, eodf * 1.5])
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+def plot_filter(ax, ax_nr,  f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df,  fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'):
+    ax[ax_nr].plot( p4[array_nr],f, color=colors[0])
+    prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
+    now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
+
+    #ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
+    #ax[ax_nr].scatter( now_height+440,  np.abs(df), marker = '>', color='black', zorder = 2)
+    #embed()
+
+
+    ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f,  facecolor=fc, edgecolor=ec)
+    ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o',
+               linestyle='None')
+    ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df),
+               color=color_df, marker='o', linestyle='None')
+    ax[ax_nr].plot(
+               np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf,
+               color=color_stim, marker='o', linestyle='None')
+    ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
+    return ax
+
+
+
+def plot_amp(ax, mean1, dev,name = 'amp',nr = 1):
+    np.unique(mean1['type'])
+    all_means = mean1[mean1['type'] == name +' mean']
+    original = all_means[all_means['dev'] == 'original']
+    #m005 = all_means[all_means['dev'] == '005']
+    m05 = all_means[all_means['dev'] == '05']
+    m2 = all_means[all_means['dev'] == '2']
+    # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True)
+    versions = [original, m05, m2] #m005,
+    for i in range(len(versions)):
+        keys = [k for k in versions[i]][2::]
+        try:
+            data = np.array(versions[i][keys])[0]
+        except:
+            break
+        axis = np.arange(0, len(data), 1)
+        axis_new = axis * 1
+        similarity = [keys, data]
+        sim = np.argsort(similarity[0])
+        # similarity[sim]
+        all_means = mean1[mean1['type'] == name+' std']
+        std = all_means[all_means['dev'] == dev[i]]
+        std = np.array(std[keys])[0]
+        #ax[1, 1].set_ylabel('Modulation depth')
+        #ax[nr,i].set_title(dev[i] + ' ms')
+        all_means = mean1[mean1['type'] == name+' 95']
+        std95 = all_means[all_means['dev'] == dev[i]]
+        std95 = np.array(std95[keys])[0]
+        all_means = mean1[mean1['type'] == name+' 05']
+        std05 = all_means[all_means['dev'] == dev[i]]
+        std05 = np.array(std05[keys])[0]
+        ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]),
+                              color='gainsboro')
+        ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]),
+                              color='darkgrey')
+
+        # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf')
+        ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black')
+        # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %')
+    #embed()
+    return ax
+
+def create_beat_corr(hz_range, eod_fr):
+    beat_corr = hz_range%eod_fr
+    beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
+    return beat_corr
+
+
+def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-1'],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']):
+    #mean1 = pd.read_pickle('mean.pkl')
+    data_all = pd.read_pickle('beat_results_smoothed.pkl')
+    d = data_all[data_all['dataset'] == data[0]]
+
+    #embed()
+
+    inch_factor = 2.54
+
+    half_page_width = 7.9 / inch_factor
+    intermediate_width = 12 / inch_factor
+    whole_page_width = 16 * 2 / inch_factor
+
+    small_length = 6 / inch_factor
+    intermediate_length = 12 * 1.5 / inch_factor
+    max_length = 25 / inch_factor
+    whole_page_width = 6.7
+    intermediate_length = 3.7
+
+    #plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
+    plt.rcParams['font.size'] = 11
+    plt.rcParams['axes.titlesize'] = 12
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['lines.linewidth'] = 1.5
+    plt.rcParams['lines.markersize'] = 8
+    plt.rcParams['legend.loc'] = 'upper right'
+    plt.rcParams["legend.frameon"] = False
+
+    # load data for plot
+
+    # data1 = pd.read_csv('ma_allcells_unsmoothed.csv')
+    # data2 = pd.read_csv('ma_allcells_05.csv')
+    # data3 = pd.read_csv('ma_allcells_2.csv')
+    # data_tob = pd.read_csv('ma_toblerone.csv')
+
+    # smothed = df_beat[df_beat['dev'] == 'original']
+    # data1 = smothed[smothed['type'] == 'amp']
+
+    x = np.arange(0, 2550, 50)
+    corr = create_beat_corr(x, np.array([500] * len(x)))
+
+    #np.unique(mean1['type'])
+    #all_means = mean1[mean1['type'] == 'max mean']
+
+    #versions = [[]]*len(dev)
+    #for i in range(len(dev)):
+    version =[[]]*len(sigma)
+    version2 = [[]] * len(sigma)
+    dev = [[]] * len(sigma)
+    limits = [[]]*len(sigma)
+    minimum = [[]] * len(sigma)
+    y_max = [[]] * len(sigma)
+    y_min = [[]] * len(sigma)
+    ax ={}
+
+    for i, e in enumerate(sigma):
+        y2 = d['result_amplitude_max_' + e]
+        y_max[i] = np.max(y2)
+        y_min[i] = np.min(y2)
+    for i,e in enumerate(sigma):
+        dev[i] = sigma[i]
+        plots = gridspec.GridSpecFromSubplotSpec( 1,1,
+                                                subplot_spec=grid[i], wspace=0.4, hspace=0.5)
+        d = data_all[data_all['dataset'] == data[0]]
+        x = d['delta_f'] / d['eodf'] + 1
+        #embed()
+        data = ['2019-10-21-aa-invivo-1']
+        #end = ['original', '005', '05', '2']
+        y = d['result_frequency_' + e]
+        #embed()
+        y2 = d['result_amplitude_max_' + e]
+        #y_sum[i] = np.nanmax(y)
+        ff = d['delta_f'] / d['eodf'] + 1
+        fe = d['beat_corr']
+    #fig.suptitle(set)
+        ax[0] = plt.subplot(plots[0])
+        if e != sigma[-1]:
+            ax[0] = remove_tick_marks(ax[0])
+        if e != 'whole':
+            ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
+            ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
+        ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
+        #embed()
+        eod = d['eodf'].iloc[0]
+
+        if np.max(y)<d['eodf'].iloc[0]*0.6:
+            color_chosen = color_df
+        else:
+            color_chosen = color_eod
+        #embed()
+        x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
+        #ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
+
+        #ax[0].set_ylabel('MPF [EODf]')
+        #ax[0].set_ylabel('Modulation ')
+        #ax[0, dd].set_title(e + ' ms')
+        ax[0].set_xlim([0, 4])
+
+        #ax[1] = plt.subplot(plots[1])
+        #if e != sigma[-1]:
+        #    ax[1] = remove_tick_marks(ax[1])
+        #ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
+        #height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
+        #if e == 'original':
+        #    whole_height = height
+        #    whole_height = np.max(y2)
+        #if (e != 'whole') and (e!= 'original'):
+            #ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
+            #embed()
+            #ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
+            #color = 'black')
+        #ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
+        #y_all[i] = np.max(y2)
+        #if i == len(sigma)-1:
+        #    ax[1].set_xlabel('stimulus frequency [EODf]')
+        #ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
+        #ax[1].spines['top'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        ax[0].spines['right'].set_visible(False)
+        ax[0].spines['top'].set_visible(False)
+        ax[0].spines['left'].set_visible(False)
+        #ax[1].spines['right'].set_visible(False)
+        #ax[1].spines['top'].set_visible(False)
+        ax[0].set_xlim([0, 5])
+        #embed()
+        ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
+        ax[0].set_yticks([])
+        # fig.tight_layout()
+        # fig.label_axes()
+        print(sigma[i])
+    #embed()
+    plt.subplots_adjust(bottom = 0.13)
+    #fig.tight_layout()
+    # fig.label_axes()
+
+
+
+
+
+if __name__ == "__main__":
+    data = ['2019-10-21-aa-invivo-1']
+    #fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
+    trans = False
+    sigma = [0.0005, 0.002]  # 0.00005,0.00025,
+    if trans == True:
+        col = 1
+        row = 2
+        col_small = len(sigma)+2
+        row_small = 1
+        l = 6
+        t = 'horizontal'
+        wd = [1]
+        hd = [1,2.5]
+        left = 1
+        right = 0
+    else:
+        col = 2
+        row = 1
+        row_small = len(sigma)+2
+        col_small = 1
+        t = 'vertical'
+        l = 9
+        wd = [3, 4]
+        hd = [1]
+        left = 0
+        right = 1
+    default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
+    grid = gridspec.GridSpec(row, col, left = 0.15,right = 0.95,wspace=0.02,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
+    hs = 0.1
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
+    wish_df = 240#324
+    color_eod = 'darkgreen'#'orange'
+    color_stim = 'navy'
+    color_df = 'orange'#'green'
+    colors = ['brown']
+    colors_mpf = ['brown']
+    colors_mod = ['steelblue']
+    line_col = 'black'
+    plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
+    #plt.show()
+    #embed()
+    #fig.savefig()
+    axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
+                                             subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
+    #embed()
+    plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
+    plt.savefig('rotatedps_singlesingle.pdf')
+    plt.savefig('../highbeats_pdf/rotatedps_singlesingle.pdf')
+    #plt.savefig('.pdf')
+    plt.show()
+    # plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
+    # plt.show()
+    # plt.close()
+    # embed()
+    # plot_single_tublerones() # original beat_activity
\ No newline at end of file
diff --git a/singlecellexample5.pdf b/singlecellexample5.pdf
index 06d7fc3..7e3f416 100644
Binary files a/singlecellexample5.pdf and b/singlecellexample5.pdf differ