From 90d9b19d9acaf41050935db88b17f1be66de39ea Mon Sep 17 00:00:00 2001 From: xaver Date: Tue, 29 Sep 2020 20:27:26 +0200 Subject: [PATCH] 29.09 --- apteronotus_code/base_eodf.py | 6 +- ...igure_apteronotus_gaincurve_cutofff_tau.py | 50 +++++ .../figure_apteronotus_jar_filter_fit.py | 184 ++++++++++++++++++ .../figure_apteronotus_rms_gaincurve.py | 149 ++++++++++++++ apteronotus_code/sin_all_normal.py | 64 +++--- apteronotus_code/sin_all_uniform.py | 65 ++++--- eigenmannia_code/eigenmannia_jar_stacked.py | 86 ++++++++ eigenmannia_code/step_response_eigen.py | 53 ++--- notes | 45 ++++- 9 files changed, 605 insertions(+), 97 deletions(-) create mode 100644 apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py create mode 100644 apteronotus_code/figure_apteronotus_jar_filter_fit.py create mode 100644 apteronotus_code/figure_apteronotus_rms_gaincurve.py create mode 100644 eigenmannia_code/eigenmannia_jar_stacked.py diff --git a/apteronotus_code/base_eodf.py b/apteronotus_code/base_eodf.py index e856b30..46b914a 100644 --- a/apteronotus_code/base_eodf.py +++ b/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', - '2018lepto4', +identifier = [#'2018lepto1', + #'2018lepto4', '2018lepto5', '2018lepto76', '2018lepto98', @@ -62,4 +62,4 @@ for ID in identifier: print(ID) print(base_eod) -embed() + embed() diff --git a/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py b/apteronotus_code/figure_apteronotus_gaincurve_cutofff_tau.py new file mode 100644 index 0000000..f734834 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/apteronotus_code/figure_apteronotus_jar_filter_fit.py b/apteronotus_code/figure_apteronotus_jar_filter_fit.py new file mode 100644 index 0000000..288169e --- /dev/null +++ b/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() \ No newline at end of file diff --git a/apteronotus_code/figure_apteronotus_rms_gaincurve.py b/apteronotus_code/figure_apteronotus_rms_gaincurve.py new file mode 100644 index 0000000..927fb57 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/apteronotus_code/sin_all_normal.py b/apteronotus_code/sin_all_normal.py index b4a03a0..f364311 100644 --- a/apteronotus_code/sin_all_normal.py +++ b/apteronotus_code/sin_all_normal.py @@ -7,43 +7,28 @@ from jar_functions import gain_curve_fit from jar_functions import avgNestedLists -identifier = [#'2018lepto1', - #'2018lepto4', - #'2018lepto5', - #'2018lepto76', +identifier = ['2018lepto1', + '2018lepto4', + '2018lepto5', + '2018lepto76', '2018lepto98', '2019lepto03', - #'2019lepto24', - #'2019lepto27', - #'2019lepto30', - #'2020lepto04', - #'2020lepto06', + '2019lepto24', + '2019lepto27', + '2019lepto30', + '2020lepto04', + '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() + diff --git a/apteronotus_code/sin_all_uniform.py b/apteronotus_code/sin_all_uniform.py index ef8f4c6..53bf32f 100644 --- a/apteronotus_code/sin_all_uniform.py +++ b/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 = [] -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) +colors_uniform = plt.cm.flag(np.linspace(0.2,0.8,len(fit_uniform))) +#colors = plt.cm.flag(np.linspace(0.2,0.8,len(fit))) + +# for ff ,f in enumerate(fit): +# ax.plot(fit_amf[ff], fit[ff],color = colors[ff]) +# ax.axvline(x=f_c[ff], ymin=0, ymax=5, ls = '-', alpha = 0.5, color= colors[ff])#colors_uniform[ff]) + +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() diff --git a/eigenmannia_code/eigenmannia_jar_stacked.py b/eigenmannia_code/eigenmannia_jar_stacked.py new file mode 100644 index 0000000..ce76220 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/eigenmannia_code/step_response_eigen.py b/eigenmannia_code/step_response_eigen.py index a3da159..bb76733 100644 --- a/eigenmannia_code/step_response_eigen.py +++ b/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') - print(dataset) + data = os.path.join(base_path, dataset, 'beats-eod.dat') + + 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.legend(loc = 'lower right') + plt.xticks(np.arange(0.0, 0.3, step=0.05)) + #plt.savefig('relative JAR') + #plt.legend(loc = 'lower right') plt.show() - embed() +embed() # natalie fragen ob sie bei verschiedenen Amplituden messen kann (siehe tim) diff --git a/notes b/notes index 342ec55..ed5ba7e 100644 --- a/notes +++ b/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 -+ eigenmannia_jar: - - pre_data und data aneinanderlegen damit nur noch ein specgram und keine lücke, absolute differenz lassen ++ figures: + apteronotus: fundament by tims bachelor thesis, important that apteronotus only shifts his frequency up (as eigenmannia doesnt --> natalies measurements) + + 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: - - step_response eigen: hier für fit relative JAR mit Normierung, bei Normierung einfach wenn j < 1Hz raus oderso +- fish_properties: - 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 -+ Q10 Wert aus Formel von Jan auf base_frequenz rechnen (adjust-eodf in jar_functions) +- step_response eigen: absolute response +- 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