diff --git a/jar_functions.py b/jar_functions.py index 82a1433..df625bf 100644 --- a/jar_functions.py +++ b/jar_functions.py @@ -107,8 +107,6 @@ def mean_noise_cut(frequencies, time, n): if k == 0: cutf[:kkk] = f cutf[kkk] = f - #t = time[kk] - #cutt[kkk] = t cutf[kkk:] = f return cutf diff --git a/jar_spectogram.png b/jar_spectogram.png deleted file mode 100644 index ca4a533..0000000 Binary files a/jar_spectogram.png and /dev/null differ diff --git a/sin_response_fit.py b/sin_response_fit.py index fc86cf1..45e0200 100644 --- a/sin_response_fit.py +++ b/sin_response_fit.py @@ -10,118 +10,158 @@ from jar_functions import mean_noise_cut def take_second(elem): # function for taking the names out of files return elem[1] -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 = [] - -identifier = '2020lepto06' - -data = sorted(np.load('%s files.npy' %identifier), 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, time, n = n) - cutt = time - plt.plot(time, jm-cutf, label='cut amfreq') - plt.plot(time, jm, label='spec') - #plt.legend() - #plt.show() - - b, a = signal.butter(4, (float(d[1])) / (0.5/dt), 'high', analog=True) # high pass filtering so our fit gets a bit better - y = signal.filtfilt(b, a, jm) - #plt.plot(time, y) - #plt.plot(time, jm) - - sinv, sinc = curve_fit(sin_response, time, y, [float(d[1]), 2, 0.5]) # fitting - print('frequency, phaseshift, amplitude:', sinv) - p = np.sqrt(sinv[1]**2) - A = np.sqrt(sinv[2] ** 2) - f = float(d[1]) - phaseshift.append(p) - gain.append(A) - amfreq.append(f) - ''' - # root mean square - RMS = np.sqrt(np.mean((cutf_arr - sin_response(cutt_arr, sinv[0], sinv[1], sinv[2]))**2)) - rootmeansquare.append(RMS) - thresh = A / np.sqrt(2) - threshold.append(thresh) - ''' - plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv)) +identifier = ['2018lepto1', + '2018lepto4', + '2018lepto5', + '2018lepto76', + '2018lepto98', + '2019lepto03', + '2019lepto24', + '2019lepto27', + '2019lepto30', + '2020lepto04', + '2020lepto06', + '2020lepto16', + '2020lepto19', + '2020lepto20' + ] +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, time, n = n) + cutt = time + # 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) + p = np.sqrt(sinv[1]**2) + A = np.sqrt(sinv[2] ** 2) + f = float(d[1]) + 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) + + # 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]: + 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) + + # predict of gain + for f in amf: + G = np.max(mgain) / np.sqrt(1 + (2*((np.pi*f*3.14)**2))) + predict.append(G) + + 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 + +# phase in degree \ No newline at end of file diff --git a/sin_response_specto.py b/sin_response_specto.py index 85dc9be..3aa099e 100644 --- a/sin_response_specto.py +++ b/sin_response_specto.py @@ -21,18 +21,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\\2020lepto04' - -datasets = ['2020-07-22-ab', - '2020-07-22-ac', - '2020-07-22-ad', - '2020-07-22-ae', - '2020-07-22-af', - '2020-07-22-ag', - '2020-07-23-ab', - '2020-07-23-ac', - '2020-07-23-ad', - '2020-07-23-ae'] +base_path = 'D:\\jar_project\\JAR\\sin\\2019lepto03' time_all = [] freq_all = [] @@ -41,26 +30,28 @@ amfrequencies = [] gains = [] files = [] -ID = [] - -for idx, dataset in enumerate(datasets): +for idx, dataset in enumerate(os.listdir(base_path)): + if dataset == 'prerecordings': + continue datapath = os.path.join(base_path, dataset, '%s.nix' % dataset) print(datapath) data, pre_data, dt = import_data(datapath) - nfft = 2**17 for d, dat in enumerate(data): - file_name = [] + if len(dat) == 1: + continue - for infodataset in datasets: - infodataset = os.path.join(base_path, infodataset, 'info.dat') + file_name = [] + ID = [] - i = parse_infodataset(infodataset) - identifier = i[0] - if not identifier[1:-2] in ID: - ID.append(identifier[1:-1]) + # identifier for file_name + infodatapath = os.path.join(base_path, dataset, 'info.dat') + i = parse_infodataset(infodatapath) + identifier = i[0] + if not identifier[1:-2] in ID: + ID.append(identifier[1:-1]) # file_name file_name.append(ID[0]) @@ -109,6 +100,8 @@ for idx, dataset in enumerate(datasets): # save filenames for this fish np.save('%s files' %ID[0], files) +print(ID) embed() -# running average over on AM-period? \ No newline at end of file +# running average over on AM-period? +