From fb14104325ff6ae675887f54181999bfb3c7a3e3 Mon Sep 17 00:00:00 2001 From: xaver Date: Wed, 2 Sep 2020 10:47:31 +0200 Subject: [PATCH] 02.09 --- eigenmannia_jar.py | 70 +++++++++++++++++++++++++++++++++++++ jar_functions.py | 73 ++++++++++++++++++++++++++++++++++++++- notes | 13 +++++++ plot_eigenmannia_jar.py | 50 +++++++++++++++++++++++++++ sin_all.py | 74 +++++++++++++++++++++++++++++++++++++++ sin_response_fit.py | 76 ++++++++++++++++++++--------------------- sin_response_specto.py | 22 +++++++----- step_response.py | 2 +- 8 files changed, 331 insertions(+), 49 deletions(-) create mode 100644 eigenmannia_jar.py create mode 100644 notes create mode 100644 plot_eigenmannia_jar.py create mode 100644 sin_all.py diff --git a/eigenmannia_jar.py b/eigenmannia_jar.py new file mode 100644 index 0000000..e7d3a43 --- /dev/null +++ b/eigenmannia_jar.py @@ -0,0 +1,70 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +import nix_helpers as nh +from IPython import embed +#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 + +base_path = 'D:\\jar_project\\JAR\\eigenmannia\\2015eigen16' + +datasets = ['2020-07-08-aa', '2020-07-08-as'] + +response = [] +deltaf = [] + +for dataset in os.listdir(base_path): + datapath = os.path.join(base_path, dataset, '%s.nix' % dataset) + print(datapath) + stimuli_dat = os.path.join(base_path, dataset, 'manualjar-eod.dat') + + df, duration = parse_stimuli_dat(stimuli_dat) + dur = int(duration[0][0:2]) + print(df) + + time, eod = nh.read_eod(datapath, duration = 2000) + + zeropoints = get_time_zeros(time, eod, threshold = np.max(eod)*0.1) + + frequencies = 1 / np.diff(zeropoints) + + # norm, base, jar = norm_function_eigen(frequencies, zeropoints[:-1], onset_point=(dur - dur)+10, offset_point=dur+10) # dm-dm funktioniert nur wenn onset = 0 sec + + # cf, ct = mean_noise_cut_eigen(frequencies, zeropoints[:-1], n=200) + + window = np.ones(101) / 101 + freq = np.convolve(frequencies, window, mode='same') + + ''' # plt.plot(ct, cf) + plt.plot(zeropoints[:-1], freq) + plt.ylabel('EOD_frequency [Hz]') + plt.xlabel('time [s]') + plt.xlim(1, 140) + plt.ylim(np.median(freq) - 10, np.median(freq) + 10) + plt.title('JAR_deltaf_%s' % deltaf) + plt.show() + ''' + j = [] + for idx, i in enumerate(zeropoints): + if i > 20 and i < 80: + j.append(freq[idx]) + + b = [] + for idx, i in enumerate(zeropoints): + if i < 20: + b.append(freq[idx]) + + r = np.median(j) - np.median(b) + response.append(r) + deltaf.append(df[0]) + +res_df1 = sorted(zip(deltaf[:32],response[:32])) +res_df2 = sorted(zip(deltaf[33:],response[33:])) +res_df = sorted(zip(deltaf,response)) + +np.save('res_df1', res_df1) +np.save('res_df2', res_df2) +np.save('res_df', res_df) \ No newline at end of file diff --git a/jar_functions.py b/jar_functions.py index df625bf..2fb3cf8 100644 --- a/jar_functions.py +++ b/jar_functions.py @@ -84,6 +84,25 @@ def parse_infodataset(dataset_name): identifier.append((l.split(':')[-1].strip())) return identifier +def parse_stimuli_dat(dataset_name): + assert (os.path.exists(dataset_name)) # see if data exists + f = open(dataset_name, 'r') # open data we gave in + lines = f.readlines() # read data + f.close() # ? + + deltaf = [] + duration = [] + + for i in range(len(lines)): + l = lines[i].strip() # all lines of textdata, exclude all empty lines (empty () default for spacebar) + if "#" in l and "Delta f" in l: + ll = (l.split(':')[-1].strip()) + deltaf.append(float(ll.split('.')[0])) + if '#' in l and 'duration' in l: + duration.append((l.split(':')[-1].strip())) + + return deltaf, duration + def mean_traces(start, stop, timespan, frequencies, time): minimumt = min([len(time[k]) for k in range(len(time))]) @@ -98,7 +117,17 @@ def mean_traces(start, stop, timespan, frequencies, time): mf = np.mean(frequency, axis=0) return mf, tnew -def mean_noise_cut(frequencies, time, n): +def mean_noise_cut_eigen(frequencies, time, n): + cutf = [] + cutt = [] + for k in np.arange(0, len(frequencies), n): + t = time[k] + f = np.mean(frequencies[k:k+n]) + cutf.append(f) + cutt.append(t) + return cutf, cutt + +def mean_noise_cut(frequencies, n): cutf = np.zeros(len(frequencies)) for k in range(0, len(frequencies) - n): kk = int(k) @@ -125,6 +154,20 @@ def norm_function(f, t, onset_point, offset_point): normed = ground / jar norm.append(normed) + return norm, base, jar + +def norm_function_eigen(f, t, onset_point, offset_point): + onset_end = onset_point - 10 + offset_start = offset_point - 10 + + + base = np.median(f[(t >= onset_end) & (t < onset_point)]) + + ground = f - base + + jar = np.median(ground[(t >= offset_start) & (t < offset_point)]) + + norm = ground / jar return norm def base_eod(frequencies, time, onset_point): @@ -147,6 +190,34 @@ def JAR_eod(frequencies, time, offset_point): return jar_eod + +def get_time_zeros (time, ampl, threshold = 0.0): + + """ + Ermittelt die Zeitpunkte der Nullpunkte der EOD-Kurve + param time: Zeitachse der Datei + param eod: EOD-Kurve aus Datei + return zeropoints: Liste mit Nullpunkten der EOD-Kurve + """ +#Xavers gedöns + new_time = time[:-1] + if len(new_time) != (len(ampl[:-1]) | len(ampl[1:])): + new_time = time [:-2] + zeropoints = new_time[(ampl[:-1] >= threshold) & (ampl[1:] < threshold)] + dx = np.mean(np.diff(new_time)) + + for index in range(len(zeropoints)): # Daten glätten + zeit_index = int(zeropoints[index] / dx) + if ampl[zeit_index] < threshold: + dy = ampl[zeit_index + 1] - ampl[zeit_index] + else: + dy = ampl[zeit_index] - ampl[zeit_index - 1] + m = (dy / dx) + x = (threshold - ampl[zeit_index]) / m + zeropoints[index] += x + return zeropoints + + def sort_values(values): a = values[:2] tau = np.array(sorted(values[2:], reverse=False)) diff --git a/notes b/notes new file mode 100644 index 0000000..7dccc12 --- /dev/null +++ b/notes @@ -0,0 +1,13 @@ +- mit zu hohem RMS rauskicken: evtl nur ein trace rauskicken wenn nur da RMS zu hoch +- 2019lepto27/30 nochmal anschauen, gainpunkt fehlt +- fit für gain kurven: über tau = 1/cutoff_f * 2 * pi --> wie bestimm ich cutoff_f? +- daten von natalie zu eigenmannia mit + / - delta f anschauen ob unterschiede +- unterschiedliche nffts auf anderem rechner laufen lassen evtl um unterschiede zu sehen + +long term: +- extra datei mit script drin um fertige daten darzustellen, den fit-code als datenverarbeitung allein verwenden +- darstellung: specgram --> rausgezogene jarspur darüber --> filterung --> fit und daten zusammen dargestellt, das ganze für verschiedene frequenzen +- größe/evtl. gewicht nachtragen, eod basefrequenz rausziehen und zuweisen +- liste mit eigenschaften der fische (dominanz/größe), messvariablen (temp/conductivity), eodf und evtl ampl machen um diese plotten zu können + +- phase in degree \ No newline at end of file diff --git a/plot_eigenmannia_jar.py b/plot_eigenmannia_jar.py new file mode 100644 index 0000000..f2c2f55 --- /dev/null +++ b/plot_eigenmannia_jar.py @@ -0,0 +1,50 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +import nix_helpers as nh +from IPython import embed + +res_df1 = np.load('res_df1.npy') +res_df2 = np.load('res_df2.npy') +res_df = np.load('res_df.npy') + +mres = [] +mdf = [] + +currf = None +idxlist = [] + +for i, d in enumerate(res_df): + if currf is None or currf == d[0]: + currf = d[0] + idxlist.append(i) + + else: # currf != f + meanres = [] # lists to make mean of + meandf = [] + for x in idxlist: + meanres.append(res_df[x][1]) + meandf.append(res_df[x][0]) + meanedres = np.mean(meanres) + meaneddf = np.mean(meandf) + mres.append(meanedres) + mdf.append(meaneddf) + currf = d[0] # set back for next loop + idxlist = [i] +meanres = [] # lists to make mean of +meandf = [] +for y in idxlist: + meanres.append(res_df[y][1]) + meandf.append(res_df[y][0]) +meanedres = np.mean(meanres) +meaneddf = np.mean(meandf) +mres.append(meanedres) +mdf.append(meaneddf) + +plt.plot(mdf, mres, 'o') +plt.xlabel('deltaf [Hz]') +plt.ylabel('JAR_respones [Hz]') +plt.axhline(0, color='grey', lw =1) +plt.axvline(0, color='grey', lw = 1) +plt.title('JAR_response_to_deltaf_eigenmannia') +plt.show() \ No newline at end of file diff --git a/sin_all.py b/sin_all.py new file mode 100644 index 0000000..44e5340 --- /dev/null +++ b/sin_all.py @@ -0,0 +1,74 @@ +import matplotlib.pyplot as plt +import numpy as np +import pylab +from IPython import embed + +def avgNestedLists(nested_vals): + """ + Averages a 2-D array and returns a 1-D array of all of the columns + averaged together, regardless of their dimensions. + """ + output = [] + maximum = 0 + for lst in nested_vals: + if len(lst) > maximum: + maximum = len(lst) + for index in range(maximum): # Go through each index of longest list + temp = [] + for lst in nested_vals: # Go through each list + if index < len(lst): # If not an index error + temp.append(lst[index]) + output.append(np.nanmean(temp)) + return output + +identifier = ['2018lepto4', + '2018lepto1', + '2018lepto5', + '2018lepto76', + '2018lepto98', + '2019lepto03', + '2019lepto24', + '2019lepto27', + '2019lepto30', + '2020lepto04', + '2020lepto06', + '2020lepto16', + '2020lepto19', + '2020lepto20' + ] + +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('gain_%s.npy' %ident) + all.append(data) + +av = avgNestedLists(all) + +fig = plt.figure() +ax = fig.add_subplot(111) +ax.plot(amf, av, 'o') +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]') +plt.show() +embed() + +'''len_arr = [] +for a in all: + len_arr.append(len(a)) +max_a = np.max(len_arr) + +arr = np.ma.empty((1,len(all),max_a)) +arr.mask = True + +for x, a in enumerate(all): + arr[:a.shape[0],x] = arr[0][x] +embed() + +print(arr.mean(axis = 2)) +embed()''' \ No newline at end of file diff --git a/sin_response_fit.py b/sin_response_fit.py index 45e0200..960e70b 100644 --- a/sin_response_fit.py +++ b/sin_response_fit.py @@ -10,20 +10,20 @@ from jar_functions import mean_noise_cut def take_second(elem): # function for taking the names out of files return elem[1] -identifier = ['2018lepto1', - '2018lepto4', - '2018lepto5', - '2018lepto76', - '2018lepto98', - '2019lepto03', - '2019lepto24', - '2019lepto27', +identifier = [#'2018lepto1', + #'2018lepto4', + #'2018lepto5', + #'2018lepto76', + #'2018lepto98', + #'2019lepto03', + #'2019lepto24', + #'2019lepto27', '2019lepto30', - '2020lepto04', - '2020lepto06', - '2020lepto16', - '2020lepto19', - '2020lepto20' + #'2020lepto04', + #'2020lepto06', + #'2020lepto16', + #'2020lepto19', + #'2020lepto20' ] for ident in identifier: @@ -56,10 +56,10 @@ for ident in identifier: n = int(1/float(d[1])/dt) cutf = mean_noise_cut(jm, time, n = n) cutt = time - # plt.plot(time, jm-cutf, label='cut amfreq') - # plt.plot(time, jm, label='spec') - # plt.legend() - # plt.show() + #plt.plot(time, jm-cutf, label='cut amfreq') + #plt.plot(time, jm, label='spec') + #plt.legend() + #plt.show() sinv, sinc = curve_fit(sin_response, time, jm - cutf, [float(d[1]), 2, 0.5]) # fitting print('frequency, phaseshift, amplitude:', sinv) @@ -73,12 +73,11 @@ for ident in identifier: # 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) - # plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv)) - # plt.legend() - # plt.show() + #plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv)) + #plt.legend() + #plt.show() # mean over same amfreqs for phase and gain if currf is None or currf == d[1]: @@ -99,6 +98,7 @@ for ident in identifier: meanedp = np.mean(meanp) meanedrms = np.mean(meanrms) meanedthresh = np.mean(meanthresh) + mgain.append(meanedf) mphaseshift.append(meanedp) rootmeansquare.append(meanedrms) @@ -118,6 +118,7 @@ for ident in identifier: meanedp = np.mean(meanp) meanedrms = np.mean(meanrms) meanedthresh = np.mean(meanthresh) + mgain.append(meanedf) mphaseshift.append(meanedp) rootmeansquare.append(meanedrms) @@ -128,9 +129,18 @@ for ident in identifier: G = np.max(mgain) / np.sqrt(1 + (2*((np.pi*f*3.14)**2))) predict.append(G) + # as arrays + mgain_arr = np.array(mgain) + 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, mgain(RMS rausgezogene jarspur darüber --> filterung --> fit und daten zusammen dargestellt, das ganze für verschiedene frequenzen -# liste mit eigenschaften der fische (dominanz/größe) und messvariablen (temp/conductivity) machen um diese plotten zu können + np.save('gain_%s' %ident, mgain_arr[idx_arr]) + np.save('amf%s' %ident, amfreq_arr[idx_arr]) -# phase in degree \ No newline at end of file +embed() \ No newline at end of file diff --git a/sin_response_specto.py b/sin_response_specto.py index 3aa099e..da9417e 100644 --- a/sin_response_specto.py +++ b/sin_response_specto.py @@ -6,8 +6,9 @@ import os import glob import IPython import numpy as np -import DataLoader as dl +#import DataLoader as dl from IPython import embed +from tqdm import tqdm from scipy.optimize import curve_fit from jar_functions import step_response from jar_functions import sin_response @@ -21,7 +22,7 @@ from jar_functions import average from jar_functions import import_data from jar_functions import import_amfreq -base_path = 'D:\\jar_project\\JAR\\sin\\2019lepto03' +base_path = 'D:\\jar_project\\JAR\\sin\\2019lepto30' time_all = [] freq_all = [] @@ -30,7 +31,7 @@ amfrequencies = [] gains = [] files = [] -for idx, dataset in enumerate(os.listdir(base_path)): +for idx, dataset in tqdm(enumerate(os.listdir(base_path))): if dataset == 'prerecordings': continue datapath = os.path.join(base_path, dataset, '%s.nix' % dataset) @@ -42,7 +43,15 @@ for idx, dataset in enumerate(os.listdir(base_path)): for d, dat in enumerate(data): if len(dat) == 1: continue + '''if len(data) > 0: + print(datapath) + amfreq = import_amfreq(datapath) + print(amfreq) + continue + else: + embed() + ''' file_name = [] ID = [] @@ -57,7 +66,7 @@ for idx, dataset in enumerate(os.listdir(base_path)): file_name.append(ID[0]) amfreq = import_amfreq(datapath) - print(amfreq) + #print(amfreq) file_name.append(str(amfreq)) file_name.append(str(d)) @@ -100,8 +109,5 @@ for idx, dataset in enumerate(os.listdir(base_path)): # save filenames for this fish np.save('%s files' %ID[0], files) -print(ID) -embed() - -# running average over on AM-period? +embed() \ No newline at end of file diff --git a/step_response.py b/step_response.py index 5b4eebc..3f940b1 100644 --- a/step_response.py +++ b/step_response.py @@ -72,7 +72,7 @@ for idx, dataset in enumerate(datasets): mf, tnew = mean_traces(start, stop, timespan, norm, time) # maybe fixed timespan/sampling rate - cf, ct = mean_noise_cut(mf, tnew, n=1250) + cf, ct = mean_noise_cut(mf, n=1250) cf_arr = np.array(cf) ct_arr = np.array(ct)