29.09

apteronotus_code/base_eodf.py

@@ -10,8 +10,8 @@ from jar_functions import mean_noise_cut_eigen
 
 base_path = 'D:\\jar_project\\JAR\\sin'
 
-identifier = ['2018lepto1',
+identifier = [#'2018lepto1',
-              '2018lepto4',
+              #'2018lepto4',
               '2018lepto5',
               '2018lepto76',
               '2018lepto98',
@@ -62,4 +62,4 @@ for ID in identifier:
     print(ID)
     print(base_eod)
 
-embed()
+    embed()

apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py

@@ -0,0 +1,50 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+from jar_functions import gain_curve_fit
+
+identifier = ['2018lepto1', '2018lepto4', '2018lepto5', '2018lepto76']
+
+tau = []
+f_c = []
+for ID in identifier:
+    predict = []
+
+    print(ID)
+    amf = np.load('amf_%s.npy' %ID)
+    gain = np.load('gain_%s.npy' %ID)
+    print(gain)
+
+    sinv, sinc = curve_fit(gain_curve_fit, amf, gain)
+    print('tau:', sinv[0])
+    tau.append(sinv[0])
+    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
+    print('f_cutoff:', f_cutoff)
+    f_c.append(f_cutoff)
+
+    # predict of gain
+    for f in amf:
+        G = np.max(gain) / np.sqrt(1 + (2 * ((np.pi * f * sinv[0]) ** 2)))
+        predict.append(G)
+    print(np.max(gain))
+
+    fig = plt.figure()
+    ax = fig.add_subplot()
+    ax.plot(amf, gain,'o' , label = 'gain')
+    ax.plot(amf, predict, label = 'fit')
+    ax.axvline(x=f_cutoff, ymin=0, ymax=5, ls='-', alpha=0.5, label = 'cutoff frequency')
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.set_ylabel('gain [Hz/(mV/cm)]')
+    ax.set_xlabel('envelope_frequency [Hz]')
+    ax.set_title('gaincurve %s' %ID)
+    plt.legend()
+    plt.show()
+
+
+#np.save('f_c', f_c)
+#np.save('tau', tau)

apteronotus_code/figure_apteronotus_jar_filter_fit.py

@@ -0,0 +1,184 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+
+from jar_functions import import_data
+from jar_functions import import_amfreq
+
+from scipy.optimize import curve_fit
+from jar_functions import sin_response
+from jar_functions import mean_noise_cut
+from jar_functions import gain_curve_fit
+
+def take_second(elem):      # function for taking the names out of files
+    return elem[1]
+
+identifier = ['2018lepto1']
+for ident in identifier:
+
+    predict = []
+
+    rootmeansquare = []
+    threshold = []
+
+    gain = []
+    mgain = []
+    phaseshift = []
+    mphaseshift = []
+    amfreq = []
+    amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+
+    currf = None
+    idxlist = []
+
+    data = sorted(np.load('%s files.npy' %ident), key = take_second)      # list with filenames in it
+
+    for i, d in enumerate(data):
+        dd = list(d)
+        if dd[1] == '0.005':
+            jar = np.load('%s.npy' %dd)     # load data for every file name
+            jm = jar - np.mean(jar)         # low-pass filtering by subtracting mean
+            print(dd)
+
+            time = np.load('%s time.npy' %dd)       # time file
+            dt = time[1] - time[0]
+
+            n = int(1/float(d[1])/dt)
+            cutf = mean_noise_cut(jm, n = n)
+            cutt = time
+
+            sinv, sinc = curve_fit(sin_response, time, jm - cutf, [float(d[1]), 2, 0.5])        # fitting
+            print('frequency, phaseshift, amplitude:', sinv)
+            p = sinv[1]
+            A = np.sqrt(sinv[2] ** 2)
+            f = float(d[1])
+            if sinv[2] < 0:
+                p = p + np.pi
+            phaseshift.append(p)
+            gain.append(A)
+            if f not in amfreq:
+                amfreq.append(f)
+
+            # jar trace
+            plt.plot(time, jar, color = 'C0')
+            #plt.hlines(y=np.min(jar) - 2, xmin=0, xmax=400, lw=2.5, color='r', label='stimulus duration')
+            plt.title('JAR trace 2018lepto1, AM-frequency:%sHz' % float(d[1]))
+            plt.xlabel('time[s]')
+            plt.ylabel('frequency[Hz]')
+            plt.show()
+
+            # low pass filter by mean subtraction
+            # plt.plot(time, jm)
+            # plt.title('JAR trace: filtered by mean subtraction')
+            # plt.xlabel('time[s]')
+            # plt.ylabel('frequency[Hz]')
+            # plt.show()
+
+            # filter by running average
+            plt.plot(time, jm, color = 'C0', label = 'JAR: subtracted by mean')
+            plt.plot(time, jm - cutf, color = 'darkorange', label = 'JAR: subtracted by mean and step response')
+            plt.title('JAR trace spectogram 2018lepto1: subtraction of mean and step response')
+            plt.xlabel('time[s]')
+            plt.ylabel('frequency[Hz]')
+            plt.legend()
+            plt.show()
+
+            # jar trace and fit
+            plt.plot(time, jm - cutf, color = 'darkorange', label = 'JAR: subtracted by mean and step response')
+            phase_gain = [(((sinv[1] % (2 * np.pi)) * 360) / (2 * np.pi)), sinv[2]]
+
+            plt.plot(time, sin_response(time, *sinv), color = 'limegreen', label='fit: phaseshift=%.2f°, gain=%.2f[Hz/(mV/cm)]' % tuple(phase_gain))
+            plt.title('JAR trace spectogram 2018lepto1 with fit')
+            plt.xlabel('time[s]')
+            plt.ylabel('frequency[Hz]')
+            plt.legend()
+            plt.show()
+
+            # root mean square
+            RMS = np.sqrt(np.mean(((jm - cutf) - sin_response(cutt, sinv[0], sinv[1], sinv[2]))**2))
+            thresh = A / np.sqrt(2)
+
+            # mean over same amfreqs for phase and gain
+            if currf is None or currf == d[1]:
+                currf = d[1]
+                idxlist.append(i)
+
+            else:  # currf != f
+                meanf = []  # lists to make mean of
+                meanp = []
+                meanrms = []
+                meanthresh = []
+                for x in idxlist:
+                    meanf.append(gain[x])
+                    meanp.append(phaseshift[x])
+                    meanrms.append(RMS)
+                    meanthresh.append(thresh)
+                meanedf = np.mean(meanf)
+                meanedp = np.mean(meanp)
+                meanedrms = np.mean(meanrms)
+                meanedthresh = np.mean(meanthresh)
+
+                mgain.append(meanedf)
+                mphaseshift.append(meanedp)
+                rootmeansquare.append(meanedrms)
+                threshold.append(meanedthresh)
+                currf = d[1]    # set back for next loop
+                idxlist = [i]
+    meanf = []
+    meanp = []
+    meanrms = []
+    meanthresh = []
+    for y in idxlist:
+        meanf.append(gain[y])
+        meanp.append(phaseshift[y])
+        meanrms.append(RMS)
+        meanthresh.append(thresh)
+    meanedf = np.mean(meanf)
+    meanedp = np.mean(meanp)
+    meanedrms = np.mean(meanrms)
+    meanedthresh = np.mean(meanthresh)
+
+    mgain.append(meanedf)
+    mphaseshift.append(meanedp)
+    rootmeansquare.append(meanedrms)
+    threshold.append(meanedthresh)
+
+    # as arrays
+    mgain_arr = np.array(mgain)
+    mphaseshift_arr = np.array(mphaseshift)
+    amfreq_arr = np.array(amfreq)
+    rootmeansquare_arr = np.array(rootmeansquare)
+    threshold_arr = np.array(threshold)
+
+    # condition needed to be fulfilled: RMS < threshold or RMS < mean(RMS)
+    idx_arr = (rootmeansquare_arr < threshold_arr) | (rootmeansquare_arr < np.mean(rootmeansquare_arr))
+
+    fig = plt.figure()
+    ax0 = fig.add_subplot(2, 1, 1)
+    ax0.plot(amfreq_arr[idx_arr], mgain_arr[idx_arr], 'o')
+    ax0.set_yscale('log')
+    ax0.set_xscale('log')
+    ax0.set_title('%s' % data[0][0])
+    ax0.set_ylabel('gain [Hz/(mV/cm)]')
+    ax0.set_xlabel('envelope_frequency [Hz]')
+    #plt.savefig('%s gain' % data[0][0])
+
+    ax1 = fig.add_subplot(2, 1, 2, sharex = ax0)
+    ax1.plot(amfreq, threshold, 'o-', label = 'threshold', color = 'b')
+    ax1.set_xscale('log')
+    ax1.plot(amfreq, rootmeansquare, 'o-', label = 'RMS', color ='orange')
+    ax1.set_xscale('log')
+    ax1.set_xlabel('envelope_frequency [Hz]')
+    ax1.set_ylabel('RMS [Hz]')
+    plt.legend()
+    pylab.show()
+
+    #np.save('phaseshift_%s' % ident, mphaseshift_arr[idx_arr])
+    #np.save('gain_%s' %ident, mgain_arr[idx_arr])
+    #np.save('amf_%s' %ident, amfreq_arr[idx_arr])
+
+embed()

apteronotus_code/figure_apteronotus_rms_gaincurve.py

@@ -0,0 +1,149 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pylab
+from IPython import embed
+from scipy.optimize import curve_fit
+from matplotlib.mlab import specgram
+import os
+
+from jar_functions import import_data
+from jar_functions import import_amfreq
+
+from scipy.optimize import curve_fit
+from jar_functions import sin_response
+from jar_functions import mean_noise_cut
+from jar_functions import gain_curve_fit
+
+def take_second(elem):      # function for taking the names out of files
+    return elem[1]
+
+identifier = ['2018lepto1']
+for ident in identifier:
+
+    predict = []
+
+    rootmeansquare = []
+    threshold = []
+
+    gain = []
+    mgain = []
+    phaseshift = []
+    mphaseshift = []
+    amfreq = []
+    amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+
+    currf = None
+    idxlist = []
+
+    data = sorted(np.load('%s files.npy' %ident), key = take_second)      # list with filenames in it
+
+    for i, d in enumerate(data):
+        dd = list(d)
+
+        jar = np.load('%s.npy' %dd)     # load data for every file name
+        jm = jar - np.mean(jar)         # low-pass filtering by subtracting mean
+        print(dd)
+
+        time = np.load('%s time.npy' %dd)       # time file
+        dt = time[1] - time[0]
+
+        n = int(1/float(d[1])/dt)
+        cutf = mean_noise_cut(jm, n = n)
+        cutt = time
+
+        sinv, sinc = curve_fit(sin_response, time, jm - cutf, [float(d[1]), 2, 0.5])        # fitting
+        print('frequency, phaseshift, amplitude:', sinv)
+        p = sinv[1]
+        A = np.sqrt(sinv[2] ** 2)
+        f = float(d[1])
+        if sinv[2] < 0:
+            p = p + np.pi
+        phaseshift.append(p)
+        gain.append(A)
+        if f not in amfreq:
+            amfreq.append(f)
+
+        # root mean square
+        RMS = np.sqrt(np.mean(((jm - cutf) - sin_response(cutt, sinv[0], sinv[1], sinv[2]))**2))
+        thresh = A / np.sqrt(2)
+
+        # mean over same amfreqs for phase and gain
+        if currf is None or currf == d[1]:
+            currf = d[1]
+            idxlist.append(i)
+
+        else:  # currf != f
+            meanf = []  # lists to make mean of
+            meanp = []
+            meanrms = []
+            meanthresh = []
+            for x in idxlist:
+                meanf.append(gain[x])
+                meanp.append(phaseshift[x])
+                meanrms.append(RMS)
+                meanthresh.append(thresh)
+            meanedf = np.mean(meanf)
+            meanedp = np.mean(meanp)
+            meanedrms = np.mean(meanrms)
+            meanedthresh = np.mean(meanthresh)
+
+            mgain.append(meanedf)
+            mphaseshift.append(meanedp)
+            rootmeansquare.append(meanedrms)
+            threshold.append(meanedthresh)
+            currf = d[1]    # set back for next loop
+            idxlist = [i]
+    meanf = []
+    meanp = []
+    meanrms = []
+    meanthresh = []
+    for y in idxlist:
+        meanf.append(gain[y])
+        meanp.append(phaseshift[y])
+        meanrms.append(RMS)
+        meanthresh.append(thresh)
+    meanedf = np.mean(meanf)
+    meanedp = np.mean(meanp)
+    meanedrms = np.mean(meanrms)
+    meanedthresh = np.mean(meanthresh)
+
+    mgain.append(meanedf)
+    mphaseshift.append(meanedp)
+    rootmeansquare.append(meanedrms)
+    threshold.append(meanedthresh)
+
+    # as arrays
+    mgain_arr = np.array(mgain)
+    mphaseshift_arr = np.array(mphaseshift)
+    amfreq_arr = np.array(amfreq)
+    rootmeansquare_arr = np.array(rootmeansquare)
+    threshold_arr = np.array(threshold)
+
+    # condition needed to be fulfilled: RMS < threshold or RMS < mean(RMS)
+    idx_arr = (rootmeansquare_arr < threshold_arr) | (rootmeansquare_arr < np.mean(rootmeansquare_arr))
+
+    fig = plt.figure()
+    ax0 = fig.add_subplot(2, 1, 1)
+    ax0.plot(amfreq_arr[idx_arr], mgain_arr[idx_arr], 'o')
+    ax0.set_yscale('log')
+    ax0.set_xscale('log')
+    ax0.set_title('gaincurve 2018lepto1')
+    ax0.set_ylabel('gain [Hz/(mV/cm)]')
+    ax0.set_xlabel('envelope_frequency [Hz]')
+    #plt.savefig('%s gain' % data[0][0])
+
+    ax1 = fig.add_subplot(2, 1, 2, sharex = ax0)
+    ax1.plot(amfreq, threshold, 'o-', label = 'threshold', color = 'b')
+    ax1.set_xscale('log')
+    ax1.plot(amfreq, rootmeansquare, 'o-', label = 'RMS', color ='orange')
+    ax1.set_xscale('log')
+    ax1.set_xlabel('envelope_frequency [Hz]')
+    ax1.set_ylabel('RMS [Hz]')
+    plt.legend()
+    pylab.show()
+
+    #np.save('phaseshift_%s' % ident, mphaseshift_arr[idx_arr])
+    #np.save('gain_%s' %ident, mgain_arr[idx_arr])
+    #np.save('amf_%s' %ident, amfreq_arr[idx_arr])
+
+embed()

apteronotus_code/sin_all_normal.py

@@ -7,43 +7,28 @@ from jar_functions import gain_curve_fit
 from jar_functions import avgNestedLists
 
 
-identifier = [#'2018lepto1',
+identifier = ['2018lepto1',
-              #'2018lepto4',
+              '2018lepto4',
-              #'2018lepto5',
+              '2018lepto5',
-              #'2018lepto76',
+              '2018lepto76',
               '2018lepto98',
               '2019lepto03',
-              #'2019lepto24',
+              '2019lepto24',
-              #'2019lepto27',
+              '2019lepto27',
-              #'2019lepto30',
+              '2019lepto30',
-              #'2020lepto04',
+              '2020lepto04',
-              #'2020lepto06',
+              '2020lepto06',
               '2020lepto16',
               '2020lepto19',
               '2020lepto20'
               ]
 
-tau = []
-f_c = []
-for ID in identifier:
-    print(ID)
-    amf = np.load('5Hz_amf_%s.npy' %ID)
-    gain = np.load('5Hz_gain_%s.npy' %ID)
-
-    sinv, sinc = curve_fit(gain_curve_fit, amf, gain)
-    #print('tau:', sinv[0])
-    tau.append(sinv[0])
-    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
-    print('f_cutoff:', f_cutoff)
-    f_c.append(f_cutoff)
-
-
 amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
 
 all = []
 
 for ident in identifier:
-    data = np.load('5Hz_gain_%s.npy' %ident)
+    data = np.load('gain_%s.npy' %ident)
     all.append(data)
 
 av = avgNestedLists(all)
@@ -51,14 +36,39 @@ av = avgNestedLists(all)
 fig = plt.figure()
 ax = fig.add_subplot(111)
 ax.plot(amf, av, 'o')
+
+#plt.show()
+
+tau = []
+f_c = []
+fit = []
+fit_amf = []
+for ID in identifier:
+    print(ID)
+    amf = np.load('amf_%s.npy' %ID)
+    gain = np.load('gain_%s.npy' %ID)
+
+    sinv, sinc = curve_fit(gain_curve_fit, amf, gain)
+    #print('tau:', sinv[0])
+    tau.append(sinv[0])
+    f_cutoff = abs(1 / (2*np.pi*sinv[0]))
+    print('f_cutoff:', f_cutoff)
+    f_c.append(f_cutoff)
+    fit.append(gain_curve_fit(amf, *sinv))
+    fit_amf.append(amf)
+
+#for ff ,f in enumerate(fit):
+#    ax.plot(fit_amf[ff], fit[ff])
 ax.set_xscale('log')
 ax.set_yscale('log')
-ax.set_title('gaincurve_average_allfish_5Hz')
+ax.set_title('gaincurve_average_allfish')
 ax.set_ylabel('gain [Hz/(mV/cm)]')
 ax.set_xlabel('envelope_frequency [Hz]')
 ax.set_ylim(0.0008, )
-ax.plot(f_c, np.full((len(identifier)), 0.0015), 'o', label = 'cutoff frequencies')
+ax.plot(f_c, np.full(len(identifier), 0.0015), 'o', label = 'cutoff frequencies')
 ax.legend()
 
 plt.show()
+
 embed()
+

apteronotus_code/sin_all_uniform.py

@@ -5,7 +5,8 @@ from IPython import embed
 from scipy.optimize import curve_fit
 from jar_functions import gain_curve_fit
 from jar_functions import avgNestedLists
-
+import matplotlib as mpl
+from matplotlib import cm
 
 identifier_uniform = ['2018lepto1',
               # '2018lepto4',
@@ -38,10 +39,32 @@ identifier = ['2018lepto1',
               '2020lepto20'
               ]
 
+amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+#colors = ['dimgray', 'dimgrey', 'gray', 'grey', 'darkgray', 'darkgrey', 'silver', 'lightgray', 'lightgrey', 'gainsboro', 'whitesmoke']
+colorss = ['g', 'b', 'r', 'y', 'c', 'm', 'k']
+all = []
+new_all = []
+for ident in identifier:
+    data = np.load('gain_%s.npy' %ident)
+    all.append(data)
+for ident in identifier_uniform:
+    data = np.load('gain_%s.npy' % ident)
+    new_all.append(data)
+
+av = avgNestedLists(all)
+new_av = avgNestedLists(new_all)
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+#ax.plot(amf, av, 'o', color = 'orange', label = 'normal')
+ax.plot(amf, new_av, 'o', label = 'uniformed')
+"""
 tau = []
 f_c = []
+fit = []
+fit_amf = []
 for ID in identifier:
-    print(ID)
+    #print(ID)
     amf = np.load('amf_%s.npy' %ID)
     gain = np.load('gain_%s.npy' %ID)
 
@@ -49,11 +72,15 @@ for ID in identifier:
     #print('tau:', sinv[0])
     tau.append(sinv[0])
     f_cutoff = abs(1 / (2*np.pi*sinv[0]))
-    print('f_cutoff:', f_cutoff)
+    #print('f_cutoff:', f_cutoff)
     f_c.append(f_cutoff)
-
+    fit.append(gain_curve_fit(amf, *sinv))
+    fit_amf.append(amf)
+"""
 tau_uniform = []
 f_c_uniform = []
+fit_uniform = []
+fit_amf_uniform = []
 for ID in identifier_uniform:
     #print(ID)
     amf = np.load('amf_%s.npy' %ID)
@@ -65,32 +92,26 @@ for ID in identifier_uniform:
     f_cutoff = abs(1 / (2*np.pi*sinv[0]))
     #print('f_cutoff:', f_cutoff)
     f_c_uniform.append(f_cutoff)
-amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+    fit_uniform.append(gain_curve_fit(amf, *sinv))
+    fit_amf_uniform.append(amf)
 
-all = []
+colors_uniform = plt.cm.flag(np.linspace(0.2,0.8,len(fit_uniform)))
-new_all = []
+#colors = plt.cm.flag(np.linspace(0.2,0.8,len(fit)))
-for ident in identifier:
+
-    data = np.load('gain_%s.npy' %ident)
+# for ff ,f in enumerate(fit):
-    all.append(data)
+#     ax.plot(fit_amf[ff], fit[ff],color =  colors[ff])
-for ident in identifier_uniform:
+#     ax.axvline(x=f_c[ff], ymin=0, ymax=5, ls = '-', alpha = 0.5, color= colors[ff])#colors_uniform[ff])
-    data = np.load('gain_%s.npy' % ident)
+
-    new_all.append(data)
+for ff, f in enumerate(fit_uniform):
+    ax.plot(fit_amf_uniform[ff], fit_uniform[ff], color =  colorss[ff]) #colors_uniform[ff])
+    ax.axvline(x=f_c_uniform[ff], ymin=0, ymax=5, ls = '-', alpha = 0.5, color= colorss[ff])#colors_uniform[ff])
 
-av = avgNestedLists(all)
-new_av = avgNestedLists(new_all)
-lim = 0.001
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.plot(amf, av, 'o', color = 'orange', label = 'normal')
-ax.plot(amf, new_av, 'o', color = 'blue', label = 'uniformed')
 ax.set_xscale('log')
 ax.set_yscale('log')
 ax.set_title('gaincurve_average_allfish')
 ax.set_ylabel('gain [Hz/(mV/cm)]')
 ax.set_xlabel('envelope_frequency [Hz]')
 ax.set_ylim(0.0008, )
-ax.plot(f_c, np.full((len(identifier)), 0.0015), 'o', color = 'orange', label = 'all cutoff frequencies')
-ax.plot(f_c_uniform, np.full((len(identifier_uniform)), 0.001), 'o', color = 'blue', label = 'uniformed cutoff frequencies')
 ax.legend()
 
 plt.show()

eigenmannia_code/eigenmannia_jar_stacked.py

@@ -0,0 +1,86 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import nix_helpers as nh
+from IPython import embed
+from matplotlib.mlab import specgram
+#from tqdm import tqdm
+from jar_functions import parse_stimuli_dat
+from jar_functions import norm_function_eigen
+from jar_functions import mean_noise_cut_eigen
+from jar_functions import get_time_zeros
+from jar_functions import import_data_eigen
+from scipy.signal import savgol_filter
+
+base_path = 'D:\\jar_project\\JAR\\eigenmannia\\deltaf'
+
+identifier = ['2013eigen13','2015eigen16', '2015eigen17', '2015eigen19', '2020eigen22','2020eigen32']
+
+response = []
+deltaf = []
+for ID in identifier:
+    for dataset in os.listdir(os.path.join(base_path, ID)):
+        datapath = os.path.join(base_path, ID, dataset, '%s.nix' % dataset)
+        print(datapath)
+        stimuli_dat = os.path.join(base_path, ID, dataset, 'manualjar-eod.dat')
+        #print(stimuli_dat)
+        delta_f, duration = parse_stimuli_dat(stimuli_dat)
+        dur = int(duration[0][0:2])
+        print(delta_f)
+        if delta_f ==[-2.0]:
+            print('HANDLE WITH CARE -2Hz:', datapath)
+        data, pre_data, dt = import_data_eigen(datapath)
+
+        #hstack concatenate: 'glue' pre_data and data
+        dat = np.hstack((pre_data, data))
+
+        # data
+        nfft = 2**17
+        spec, freqs, times = specgram(dat[0], Fs=1 / dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.95)
+        dbspec = 10.0 * np.log10(spec)  # in dB
+        power = dbspec[:, 25]
+
+        fish_p = power[(freqs > 200) & (freqs < 1000)]
+        fish_f = freqs[(freqs > 200) & (freqs < 1000)]
+
+        index = np.argmax(fish_p)
+        eodf = fish_f[index]
+        eodf4 = eodf * 4
+
+        lim0 = eodf4 - 50
+        lim1 = eodf4 + 50
+
+        df = freqs[1] - freqs[0]
+        ix0 = int(np.floor(lim0/df))    # back to index
+        ix1 = int(np.ceil(lim1/df))    # back to index
+        spec4= dbspec[ix0:ix1, :]
+        freq4 = freqs[ix0:ix1]
+        jar4 = freq4[np.argmax(spec4, axis=0)]      # all freqs at max specs over axis 0
+
+        cut_time_jar = times[:len(jar4)]
+
+        #plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10)
+        #plt.plot(cut_time_jar, jar4)
+        #plt.show()
+
+        b = []
+        for idx, i in enumerate(times):
+            if i > 0 and i < 10:
+                b.append(jar4[idx])
+        j = []
+        for idx, i in enumerate(times):
+            if i > 15 and i < 55:
+                j.append(jar4[idx])
+
+        r = np.median(j) - np.median(b)
+        print(r)
+        deltaf.append(delta_f[0])
+        response.append(r)
+
+    res_df = sorted(zip(deltaf,response))
+
+    np.save('res_df_%s_new' %ID, res_df)
+
+# problem: rohdaten(data, pre_data) lassen sich auf grund ihrer 1D-array struktur nicht savgol filtern
+# diese bekomm ich nur über specgram in form von freq / time auftragen, was nicht mehr savgol gefiltert werden kann
+# jedoch könnte ich trotzdem einfach aus jar4 response herauslesen wobei dies dann weniger gefiltert wäre

eigenmannia_code/step_response_eigen.py

@@ -19,7 +19,7 @@ from jar_functions import average
 base_path = 'D:\\jar_project\\JAR\\eigen\\step'
 
 identifier = ['step_2015eigen8',
-              'step_2015eigen15',
+              'step_2015eigen15\\+15Hz',
               'step_2015eigen16',
               'step_2015eigen17',
               'step_2015eigen19']
@@ -46,17 +46,19 @@ for ID in identifier:
     response = []
     stim_ampl = []
     for idx, dataset in enumerate(os.listdir(base_path)):
-        dataset = os.path.join(base_path, dataset, 'beats-eod.dat')
+        data = os.path.join(base_path, dataset, 'beats-eod.dat')
-        print(dataset)
+
+        if dataset == 'prerecordings':
+            continue
         #input of the function
-        frequency, time, amplitude, eodf, deltaf, stimulusf, stimulusamplitude, duration, pause = parse_dataset(dataset)
+        frequency, time, amplitude, eodf, deltaf, stimulusf, stimulusamplitude, duration, pause = parse_dataset(data)
         dm = np.mean(duration)
         pm = np.mean(pause)
         timespan = dm + pm
         start = np.mean([t[0] for t in time])
         stop = np.mean([t[-1] for t in time])
-        if len(frequency) == 5:
+
-            continue
+        print(dataset)
 
         mf, tnew = mean_traces(start, stop, timespan, frequency, time) # maybe fixed timespan/sampling rate
 
@@ -72,49 +74,26 @@ for ID in identifier:
         for index, i in enumerate(ct):
             if i > -45 and i < -5:
                 b.append(cf[index])
-
         j = []
         for indexx, h in enumerate(ct):
-            if h > 195 and h < 145:
+            if h < 195 and h > 145:
                 j.append(cf[indexx])
-                print(h)
-                print(indexx)
-                print(cf[indexx])
-
-        ''' sounds good, doesnt work somehow: in norm devision by 0 (jar) or index doesnt fit
-                norm, base, jar = norm_function(frequency, time, onset_point=dm - dm,
-                                                offset_point=dm)  # dm-dm funktioniert nur wenn onset = 0 sec
-                b = []
-                for index, i in enumerate(ct):
-                    if i > -45 and i < -5:
-                        b.append(cf[index])
-
-                j = []
-                for indexx, h in enumerate(ct):
-                    if h > 195 and h < 145:
-                        j.append(cf[indexx])
-                        print(h)
-                        print(indexx)
-                        print(cf[indexx])
-        b = np.median(cf[(ct >= onset_end) & (ct < onset_point)])
-
-        j = np.median(cf[(ct >= offset_start) & (ct < offset_point)])
-
-        '''
 
         r = np.median(j) - np.median(b)
         response.append(r)
-        stim_ampl.append(stimulusamplitude)
+        stim_ampl.append(float(stimulusamplitude[0]))
+
     res_ampl = sorted(zip(stim_ampl, response))
-    base_line = plt.axhline(y = 0, color = 'black', ls = 'dotted', linewidth = '1')
 
+    plt.plot(stim_ampl, response, 'o')
     plt.xlabel('Stimulusamplitude')
     plt.ylabel('absolute JAR magnitude')
     plt.title('absolute JAR')
-    plt.savefig('relative JAR')
+    plt.xticks(np.arange(0.0, 0.3, step=0.05))
-    plt.legend(loc = 'lower right')
+    #plt.savefig('relative JAR')
+    #plt.legend(loc = 'lower right')
     plt.show()
-    embed()
+embed()
 
 
 # natalie fragen ob sie bei verschiedenen Amplituden messen kann (siehe tim)

notes

@@ -1,13 +1,42 @@
-+ sin_all_uniform - sin_all_normal (also 5Hz, let away 0.001Hz?, gain_fit): fit als spur reinlegen damit klar wird aus was gerade besteht
++ figures:
-+ eigenmannia_jar:
+    apteronotus: fundament by tims bachelor thesis, important that apteronotus only shifts his frequency up (as eigenmannia doesnt --> natalies measurements)
-    - pre_data und data aneinanderlegen damit nur noch ein specgram und keine lücke, absolute differenz lassen
+    + spectogram
+    + jar trace out of specgram
+    + filtering of jar trace: mean noise cut --> subtracting jar response over whole stimulus
+    + fit and jar trace --> gain and phaseshift
+!!! + this for different am-frequencies and delta f (-15/-5Hz) --> compare gain for them
+    + gain curve for one or more single fish
+    + fit of gain curve for cutoff frequency and tau
+    + gain curve for all fish taken together
+    + single gain curves inside gain curve for all fish --> different cutoff frequencies --> comparison to metzen/chacron
+    --> fig_apt_specgram,
+        fig_apt_jar_filter_fit,
+        fig_apt_rms_gaincurve,
+        fig_apt_gaincurve_cutoff_tau,
+        sin_all_normal (without single gaincurves),
+        sin_all_uniform (with gaincurves for 5Hz)
+
+    eigenmannia:
+    + deltaf / response: -2Hz different, show it
+    + spectogram
+    + direct to fit and jar trace --> gain and phaseshift   DURCH SIN RESPONSE SPEC JAGEN!
+    + gain curve for one or more single fish
+    + gain curve for all fish taken together
+    - (step response eigen)
+
+    fish properties:
+    + parameters
+    + cutoff frequency - dominance score
++ eigenmannia deltaf response over all fish mean
++ phaseshift_all: wenn negativer gain in fit --> +pi rechnen, dann modulo
+- plot_eigenmannia_jar(compare res_df_%s / res_df_%s_new)
+- eigenmannia_jar:
     - specgram auch zeigen, vorallem was auch die ausreißer bei -2 Hz betreffen
-+ plot_eigenmannia_jar(compare res_df_%s / res_df_%s_new)
+- fish_properties:
-+ fish_properties:
-    - step_response eigen: hier für fit relative JAR mit Normierung, bei Normierung einfach wenn j < 1Hz raus oderso
     - hauptsächlich auf f_c und tau konzentrieren, vor allem auch beides auftragen, gewicht/größe noch nehmen
-+ phaseshift_all: wenn negativer gain in fit --> +pi rechnen, dann modulo
+- step_response eigen: absolute response
-+ Q10 Wert aus Formel von Jan auf base_frequenz rechnen (adjust-eodf in jar_functions)
+- Q10 Wert aus Formel von Jan auf base_frequenz rechnen (adjust-eodf in jar_functions)
+- sin_all_uniform - sin_all_normal (also 5Hz, let away 0.001Hz?, gain_fit): fit als spur reinlegen damit klar wird aus was gerade besteht
 
 long term:
 - extra datei mit script drin um fertige daten darzustellen, den fit-code nur zur datenverarbeitung verwenden