diff --git a/jar_functions.py b/jar_functions.py index 36477a9..20ee141 100644 --- a/jar_functions.py +++ b/jar_functions.py @@ -10,7 +10,7 @@ def step_response(t, a1, a2, tau1, tau2): return r_step def sin_response(t, f, p, A): - r_sin = A*sin(2*np.pi*t*f + p) + r_sin = A*np.sin(2*np.pi*t*f + p) return r_sin def parse_dataset(dataset_name): @@ -172,7 +172,6 @@ def average(freq_all, time_all, start, stop, timespan, dm): def import_data(dataset): - print(dataset) import nixio as nix nf = nix.File.open(dataset, nix.FileMode.ReadOnly) b = nf.blocks[0] @@ -193,6 +192,17 @@ def import_data(dataset): return dat, pre_dat, dt #nf.close() +def import_amfreq(dataset): + import nixio as nix + nf = nix.File.open(dataset, nix.FileMode.ReadOnly) + b = nf.blocks[0] + eod = b.data_arrays['EOD-1'] + dt = eod.dimensions[0].sampling_interval + di = int(50.0/dt) + t = b.tags['Beats_1'] + amfreq = t.metadata['RePro-Info']['settings']['amfreq'] + + return amfreq if __name__ == '__main__': import_data(os.path.join('JAR', '2020-07-21-ak', '2020-07-21-ak.nix')) diff --git a/scratch.py b/scratch.py index 60aee9a..1e715c2 100644 --- a/scratch.py +++ b/scratch.py @@ -1,5 +1,6 @@ import os import numpy as np +import matplotlib.pyplot as plt from IPython import embed from jar_functions import mean_noise_cut from matplotlib.mlab import specgram @@ -11,7 +12,7 @@ import DataLoader as dl for info, key, time, data in dl.iload_traces(datapath, repro='Beats', before=0.0, after=0.0): ''' -base_path = 'D:\\jar_project\\JAR' +'''base_path = 'D:\\jar_project\\JAR' #nicht: -5Hz delta f, 19-aa, 22-ae, 22-ad (?) datasets = [#'2020-06-19-aa', #-5Hz delta f, horrible fit @@ -38,6 +39,10 @@ for d in range(int(len(data)/10)): #print(freqs) #print(times) - embed() + embed()''' +g = [1.2917623576698833, -5.479055166593157, -2.689492238578325, -0.11604244418416806, -0.05353823781665627] +a = [0.2, 0.002, 0.02, 0.5, 1.0] +plt.plot(a, g) +plt.show() diff --git a/sin_response_fit.py b/sin_response_fit.py new file mode 100644 index 0000000..e69de29 diff --git a/sin_response_specto.py b/sin_response_specto.py new file mode 100644 index 0000000..8e40c87 --- /dev/null +++ b/sin_response_specto.py @@ -0,0 +1,135 @@ +import matplotlib.pyplot as plt +import matplotlib as cm +from matplotlib.colors import ListedColormap, LinearSegmentedColormap +from matplotlib.mlab import specgram +import os +import glob +import IPython +import numpy as np +import DataLoader as dl +from IPython import embed +from scipy.optimize import curve_fit +from jar_functions import step_response +from jar_functions import sin_response +from jar_functions import parse_dataset +from jar_functions import parse_infodataset +from jar_functions import mean_traces +from jar_functions import mean_noise_cut +from jar_functions import norm_function +from jar_functions import sort_values +from jar_functions import average +from jar_functions import import_data +from jar_functions import import_amfreq + +base_path = 'D:\\jar_project\\JAR\\sin\\2018lepto98' + +#nicht: -5Hz delta f, 19-aa, 22-ae, 22-ad (?) + +#dat = glob.glob('D:\\jar_project\\JAR\\2020*\\beats-eod.dat') +#infodat = glob.glob('D:\\jar_project\\JAR\\2020*\\info.dat') +datasets = [#'2020-06-19-aa', #-5Hz delta f, horrible fit + #'2020-06-19-ab', #-5Hz delta f, bad fit + #'2020-06-22-aa', #-5Hz delta f, bad fit + #'2020-06-22-ab', #-5Hz delta f, bad fit + #'2020-06-22-ac', #-15Hz delta f, good fit + #'2020-06-22-ad', #-15Hz delta f, horrible fit + #'2020-06-22-ae', #-15Hz delta f, horrible fit + #'2020-06-22-af', #-15Hz delta f, good fit + '2020-07-21-al', #sin + '2020-07-21-am', + '2020-07-21-ak', + #'2020-07-21-an', + #'2020-07-21-ao' + ] + +time_all = [] +freq_all = [] + +amfrequencies = [] +gains = [] + +ID = [] +col = ['dimgrey', 'grey', 'darkgrey', 'silver', 'lightgrey', 'gainsboro', 'whitesmoke'] +labels = zip(ID, datasets) + +for infodataset in datasets: + infodataset = os.path.join(base_path, infodataset, 'info.dat') + + i = parse_infodataset(infodataset) + identifier = i[0] + if not identifier[1:-2] in ID: + ID.append(identifier[1:-2]) + + +for idx, dataset in enumerate(datasets): + 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): + amfreq = float(import_amfreq(datapath)) + print(amfreq) + amfrequencies.append(amfreq) + + file_name = [] + file_name.append(ID) + file_name.append(amfreq) + file_name.append(d) + + spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft*0.6) + + dbspec = 10.0*np.log10(spec) # in dB + power = dbspec[:, 50] + + fish_p = power[(freqs > 400) & (freqs < 1000)] + fish_f = freqs[(freqs > 400) & (freqs < 1000)] + + index = np.argmax(fish_p) + eodf = fish_f[index] + eodf4 = eodf * 4 + + 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 + spec4 = dbspec[ix0:ix1, :] + freq4 = freqs[ix0:ix1] + jar4 = freq4[np.argmax(spec4, axis=0)] # all freqs at max specs over axis 0 + jar = jar4 / 4 + jm = jar4 - np.mean(jar4) + + cut_times = times[:len(jar4)] + #np.save('%s.npy' % file_name, jar4) + + plt.plot(cut_times, jm, '-k') + + cf, ct = mean_noise_cut(jar4, cut_times, n = int(round(len(jar4)/((times[-1] - times [0]) * amfreq)))) + + #plt.plot(ct, cf, '-k') + + #plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10) + #np.save( , spec4) + + embed() + b, a = signal.butter(4, 0.01 / 10000, 'high', analog=True) + y = signal.filtfilt(b, a, jm) + sinv, sinc = curve_fit(sin_response, cut_times, jm, [amfreq, 2, 0.5]) + print('frequency, phaseshift, amplitude:', sinv) + gains.append(sinv[2]) + + plt.plot(cut_times, sin_response(cut_times, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv)) + + #plt.legend() + #plt.ylim(lim0, lim1) + plt.legend() + plt.show() + + #print(np.load('%s.npy' % file_name)) +embed() \ No newline at end of file diff --git a/step_response_specto.py b/step_response_specto.py deleted file mode 100644 index 2ed32df..0000000 --- a/step_response_specto.py +++ /dev/null @@ -1,99 +0,0 @@ -import matplotlib.pyplot as plt -import matplotlib as cm -from matplotlib.colors import ListedColormap, LinearSegmentedColormap -from matplotlib.mlab import specgram -import os -import glob -import IPython -import numpy as np -import DataLoader as dl -from IPython import embed -from scipy.optimize import curve_fit -from jar_functions import parse_dataset -from jar_functions import parse_infodataset -from jar_functions import mean_traces -from jar_functions import mean_noise_cut -from jar_functions import norm_function -from jar_functions import step_response -from jar_functions import sort_values -from jar_functions import average -from jar_functions import import_data - -base_path = 'D:\\jar_project\\JAR' - -#nicht: -5Hz delta f, 19-aa, 22-ae, 22-ad (?) -datasets = [#'2020-06-19-aa', #-5Hz delta f, horrible fit - #'2020-06-19-ab', #-5Hz delta f, bad fit - #'2020-06-22-aa', #-5Hz delta f, bad fit - #'2020-06-22-ab', #-5Hz delta f, bad fit - #'2020-06-22-ac', #-15Hz delta f, good fit - #'2020-06-22-ad', #-15Hz delta f, horrible fit - #'2020-06-22-ae', #-15Hz delta f, horrible fit - #'2020-06-22-af', #-15Hz delta f, good fit - '2020-07-21-ak' #sin - ] - -#dat = glob.glob('D:\\jar_project\\JAR\\2020*\\beats-eod.dat') -#infodat = glob.glob('D:\\jar_project\\JAR\\2020*\\info.dat') - -time_all = [] -freq_all = [] - -ID = [] -col = ['dimgrey', 'grey', 'darkgrey', 'silver', 'lightgrey', 'gainsboro', 'whitesmoke'] -labels = zip(ID, datasets) - -for infodataset in datasets: - infodataset = os.path.join(base_path, infodataset, 'info.dat') - i = parse_infodataset(infodataset) - identifier = i[0] - ID.append(identifier) - - -for idx, dataset in enumerate(datasets): - datapath = os.path.join(base_path, dataset, '%s.nix' % dataset) - - data, pre_data, dt = import_data(datapath) - - - nfft = 2**15 - - spec, freqs, times = specgram(data[0], Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft*0.95) - - dbspec = 10.0*np.log10(spec) # in dB - - power = dbspec[:, 50] - - fish_p = power[(freqs > 400) & (freqs < 1000)] - fish_f = freqs[(freqs > 400) & (freqs < 1000)] - - index = np.argmax(fish_p) - eodf = fish_f[index] - eodf4 = eodf * 4 - print(eodf4) - lim0 = eodf4-20 - lim1 = eodf4+20 - - - #plt.imshow(dbspec, cmap='jet', origin='lower', aspect='auto', vmin=-80, vmax=-10, extent=(times[0], times[-1], freqs[0], freqs[-1])) - - # control of plt.imshow - df = freqs[1] - freqs[0] - ix = int(np.round(eodf4/df)) - ix0 = int(np.round(lim0/df)) - ix1 = int(np.round(lim1/df)) - spec4 = dbspec[ix0:ix1, :] - freq4 = freqs[ix0:ix1] - jar4 = freq4[np.argmax(spec4, axis=0)] - jar = jar4 / 4 - #np.save - plt.plot(times[:4862], jar4, '-k') - #plt.plot(times, jar) - - plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], freq4[0], freq4[-1]), aspect='auto', - vmin=-80, vmax=-10) - -#plt.plot(freqs, dbspec[:,10], label = '0') -#plt.legend() -plt.ylim(lim0, lim1) -plt.show()