From 4a58f1382f2fd710a8b0efbcfd6eb94d49f1f214 Mon Sep 17 00:00:00 2001 From: xaver Date: Thu, 20 Aug 2020 18:47:45 +0200 Subject: [PATCH] 20.08 --- jar_functions.py | 20 +++++++++------- sin_response_fit.py | 54 ++++++++++++++++++++---------------------- sin_response_specto.py | 49 +++++++++++++++++++------------------- 3 files changed, 62 insertions(+), 61 deletions(-) diff --git a/jar_functions.py b/jar_functions.py index 09f4589..82a1433 100644 --- a/jar_functions.py +++ b/jar_functions.py @@ -99,14 +99,18 @@ def mean_traces(start, stop, timespan, frequencies, time): return mf, tnew def mean_noise_cut(frequencies, time, n): - cutf = [] - cutt = [] - for k in np.arange(0, len(frequencies), n): - f = np.mean(frequencies[k:k+n]) - t = time[k] - cutf.append(f) - cutt.append(t) - return cutf, cutt + cutf = np.zeros(len(frequencies)) + for k in range(0, len(frequencies) - n): + kk = int(k) + f = np.mean(frequencies[kk:kk+n]) + kkk = int(kk+n/2) + if k == 0: + cutf[:kkk] = f + cutf[kkk] = f + #t = time[kk] + #cutt[kkk] = t + cutf[kkk:] = f + return cutf def norm_function(f, t, onset_point, offset_point): onset_end = onset_point - 10 diff --git a/sin_response_fit.py b/sin_response_fit.py index 54fa7c6..fc86cf1 100644 --- a/sin_response_fit.py +++ b/sin_response_fit.py @@ -25,7 +25,9 @@ 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('files.npy'), key = take_second) # list with filenames in it +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) @@ -34,41 +36,40 @@ for i, d in enumerate(data): print(dd) time = np.load('%s time.npy' %dd) # time file - for float(d[1]) = 0.001: #ähnlich wie das, einfach passenstens n verwenden? - n = int((len(jar) / time[-1]) * (1 / 15)) - print(n) - cutf, cutt = mean_noise_cut(jar, time, n = 2) - plt.plot(cutt, cutf, label='cut amfreq') - plt.plot(time, jar, label='spec') - plt.legend() - plt.show() - #embed() + 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]) / 2) / 10000, 'high', analog=True) # high pass filtering so our fit gets a bit better + 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, jar) + #plt.plot(time, jm) sinv, sinc = curve_fit(sin_response, time, y, [float(d[1]), 2, 0.5]) # fitting print('frequency, phaseshift, amplitude:', sinv) - plt.show() 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((jm - sin_response(time, sinv[0], sinv[1], sinv[2]))**2)) + 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)) - - #mean over same amfreqs for phase and gain + ''' + 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) @@ -85,7 +86,6 @@ for i, d in enumerate(data): mphaseshift.append(meanedp) currf = d[1] # set back for next loop idxlist = [i] - meanf = [] meanp = [] for y in idxlist: @@ -101,8 +101,6 @@ for f in amf: G = np.max(mgain) / np.sqrt(1 + (2*((np.pi*f*3.14)**2))) predict.append(G) - - fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.plot(amf, mgain, 'o') @@ -111,19 +109,19 @@ ax.set_yscale('log') ax.set_xscale('log') ax.set_title('%s' % data[0][0]) ax.set_ylabel('gain [Hz/(mV/cm)]') -ax.set_xlabel('AM-frequency [Hz]') +ax.set_xlabel('envelope_frequency [Hz]') #plt.savefig('%s gain' % data[0][0]) pylab.show() + plt.plot(threshold, label = 'threshold') plt.plot(rootmeansquare, label = 'RMS') plt.legend() plt.show() + embed() -#phase in degree +# phase in degree # Q10 / conductivity # AM-frequency / envelope-frequency scale title? -# bevor fit noch filtern mit 15Hz damit AM-Modulation rausgefiltert wird und nur noch envelope übrig bleibt. -# Dazu running average mit n wobei n über samplingrate von spectogram und delta f bestimmt wird -# samplingrate über overlap muss dabei aber größer sein als samplingrate die noch übrig bleibt wenn ich mit delta f frequency gefiltert hab \ No newline at end of file +# einfach passendes n verwenden um AM-beats rauszufiltern? ...Menge Datenpunkte zu gering diff --git a/sin_response_specto.py b/sin_response_specto.py index da36942..85dc9be 100644 --- a/sin_response_specto.py +++ b/sin_response_specto.py @@ -21,18 +21,18 @@ from jar_functions import average from jar_functions import import_data from jar_functions import import_amfreq -base_path = 'D:\\jar_project\\JAR\\sin\\2020lepto16' - -datasets = ['2020-08-04-ab', - '2020-08-04-ac', - '2020-08-04-ad', - '2020-08-04-ae', - '2020-08-04-af', - '2020-08-05-ab', - '2020-08-05-ac', - '2020-08-05-ad', - '2020-08-05-ae', - '2020-08-05-af'] +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'] time_all = [] freq_all = [] @@ -41,11 +41,7 @@ amfrequencies = [] gains = [] files = [] - ID = [] -col = ['dimgrey', 'grey', 'darkgrey', 'silver', 'lightgrey', 'gainsboro', 'whitesmoke'] -labels = zip(ID, datasets) - for idx, dataset in enumerate(datasets): datapath = os.path.join(base_path, dataset, '%s.nix' % dataset) @@ -53,13 +49,11 @@ for idx, dataset in enumerate(datasets): data, pre_data, dt = import_data(datapath) - nfft = 2**17 for d, dat in enumerate(data): file_name = [] - for infodataset in datasets: infodataset = os.path.join(base_path, infodataset, 'info.dat') @@ -67,6 +61,8 @@ for idx, dataset in enumerate(datasets): identifier = i[0] if not identifier[1:-2] in ID: ID.append(identifier[1:-1]) + + # file_name file_name.append(ID[0]) amfreq = import_amfreq(datapath) @@ -75,10 +71,12 @@ for idx, dataset in enumerate(datasets): file_name.append(str(d)) files.append(file_name) + + # spectogram if float(amfreq) < 0.01: spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.8) else: - spec, freqs, times = specgram(dat, Fs=1 / dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.96) + spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.95) dbspec = 10.0*np.log10(spec) # in dB power = dbspec[:, 50] @@ -93,7 +91,6 @@ for idx, dataset in enumerate(datasets): lim0 = eodf4-10 lim1 = eodf4+15 - # control of plt.imshow df = freqs[1] - freqs[0] ix0 = int(np.floor(lim0/df)) # back to index ix1 = int(np.ceil(lim1/df)) # back to index @@ -102,14 +99,16 @@ for idx, dataset in enumerate(datasets): jar4 = freq4[np.argmax(spec4, axis=0)] # all freqs at max specs over axis 0 jar = jar4 / 4 jm = jar4 - np.mean(jar4) # data we take - cut_times = times[:len(jar4)] + + # plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10) + + # save data np.save('%s time' % file_name, cut_times) np.save('%s' % file_name, jar4) - #plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10) - - np.save('files.npy', files) -#embed() +# save filenames for this fish +np.save('%s files' %ID[0], files) +embed() # running average over on AM-period? \ No newline at end of file