diff --git a/jar_functions.py b/jar_functions.py index 8e61ef9..36477a9 100644 --- a/jar_functions.py +++ b/jar_functions.py @@ -171,117 +171,29 @@ def average(freq_all, time_all, start, stop, timespan, dm): return mf_all, tnew_all, values_all -def iload_traces(basedir, repro='', before=0.0, after=0.0): - """ - returns: - info : metadata from stimuli.dat - key : key from stimuli.dat - time : an array for the time axis - data : the data of all traces of a single trial - """ - p = re.compile('([-+]?\d*\.\d+|\d+)\s*(\w+)') - reproid = 'RePro' - deltat = None - - # open traces files: - sf = [] - for trace in range(1, 1000000): - if path.isfile('%s/trace-%i.raw' % (basedir, trace)): - sf.append(open('%s/trace-%i.raw' % (basedir, trace), 'rb')) - else: - break - - basecols = None - baserp = True - - for info, key, dat in iload('%s\\stimuli.dat' % (basedir,)): - if deltat is None: - deltat, tunit = p.match(info[0]['sample interval%i' % 1]).groups() - deltat = float(deltat) - if tunit == 'ms': - deltat *= 0.001 - - if 'repro' in info[-1] or 'RePro' in info[-1]: - if not reproid in info[-1]: - reproid = 'repro' - if len(repro) > 0 and repro != info[-1][reproid] and \ - basecols is None: - continue - baserp = (info[-1][reproid] == 'BaselineActivity') - duration_indices = [i for i, x in enumerate(key[2]) if x == "duration"] - else: - baserp = True - duration_indices = [] - - if dat.shape == (1, 1) and dat[0, 0] == 0: - warnings.warn("iload_traces: Encountered incomplete '-0' trial.") - yield info, key, array([]) - continue - - if baserp: - basecols = [] - basekey = key - baseinfo = info - if len(dat) == 0: - for trace in range(len(sf)): - basecols.append(0) - for d in dat: - if not baserp and not basecols is None: - x = [] - xl = [] - for trace in range(len(sf)): - col = int(d[trace]) - sf[trace].seek(basecols[trace] * 4) - buffer = sf[trace].read((col - basecols[trace]) * 4) - tmp = fromstring(buffer, float32) - x.append(tmp) - xl.append(len(tmp)) - ml = min(xl) - for k in range(len(x)): - if len(x[k]) > ml: - warnings.warn("trunkated trace %d to %d" % (k, ml)) - x[k] = x[k][:ml] - xtime = arange(0.0, len(x[0])) * deltat - before - yield baseinfo, basekey, xtime, asarray(x) - basecols = None - if len(repro) > 0 and repro != info[-1][reproid]: - break - - durations = [d[i] for i in duration_indices if not isnan(d[i])] - if not baserp: - if len(durations) > 0: - duration = max(durations) - if duration < 0.001: # if the duration is less than 1ms - warnings.warn( - "iload_traces: Skipping one trial because its duration is <1ms and therefore it is probably rubbish") - continue - l = int(before / deltat) - r = int((duration + after) / deltat) - else: - continue - x = [] - xl = [] - for trace in range(len(sf)): - col = int(d[trace]) - if baserp: - if col < 0: - col = 0 - basecols.append(col) - continue - sf[trace].seek((col - l) * 4) - buffer = sf[trace].read((l + r) * 4) - tmp = fromstring(buffer, float32) - x.append(tmp) - xl.append(len(tmp)) - if baserp: - break - ml = min(xl) - for k in range(len(x)): - if len(x[k]) > ml: - warnings.warn("trunkated trace %d to %d" % (k, ml)) - x[k] = x[k][:ml] - time = arange(0.0, len(x[0])) * deltat - before - yield info, key, time, asarray(x) - - for trace in range(len(sf)): - sf[trace].close() +def import_data(dataset): + print(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'] + dat = [] + pre_dat = [] + for mt in b.multi_tags: + data = mt.retrieve_data(0, 'EOD-1')[:] # data[0] + dat.append(data) + i0 = int(mt.positions[0][0]/dt) + pre_data = eod[i0-di:i0] + pre_dat.append(pre_data) + + return dat, pre_dat, dt + #nf.close() + + +if __name__ == '__main__': + import_data(os.path.join('JAR', '2020-07-21-ak', '2020-07-21-ak.nix')) + diff --git a/jar_spectogram.png b/jar_spectogram.png new file mode 100644 index 0000000..ca4a533 Binary files /dev/null and b/jar_spectogram.png differ diff --git a/scratch.py b/scratch.py index 49d22d3..60aee9a 100644 --- a/scratch.py +++ b/scratch.py @@ -2,47 +2,42 @@ import os import numpy as np from IPython import embed from jar_functions import mean_noise_cut +from matplotlib.mlab import specgram +import DataLoader as dl -datasets = [(os.path.join('D:\\jar_project\\JAR\\2020-06-22-ac\\beats-eod.dat'))] - -""" -#second_try scratch -minimum = min(len(frequency[0]), len(frequency[1])) -f1 = frequency[0][:minimum] -f2 = frequency[1][:minimum] - -frequency = np.array(frequency) - -mean = np.mean(frequency, axis=0) - -for i in range(len(minimum)): - mean(frequency[0][i], frequency[1][i]) - -for f in frequency: - print(np.mean(f)) +#print(np.logspace(-3, 1, 10)) +'''for idx, dataset in enumerate(datasets): + datapath = os.path.join(base_path, dataset) + for info, key, time, data in dl.iload_traces(datapath, repro='Beats', before=0.0, after=0.0): +''' -# mean_f = np.mean(x) for x in zip(freqeuncies1, frequencies2) -""" +base_path = 'D:\\jar_project\\JAR' -''' -g = [1, 2] +#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 + ] +for idx, dataset in enumerate(datasets): + datapath = os.path.join(base_path, dataset) + for info, key, time, data in dl.iload_traces(datapath, repro='Beats', before=0.0, after=0.0): + print(data[0]) -h = [3, 4] +dat = np.arange(100) -z = np.array([[g], [h]]) +for d in range(int(len(data)/10)): + nfft = 2 -mean0 = np.mean(z, axis=0) -mean1 = np.mean(z, axis=1) + spec, freqs, times = specgram(data[0][d*10:(d+1)*10], NFFT=nfft, noverlap=nfft*0.5) -print(mean0) -print(mean1) -'''''' -t = [600, 650] + #print(freqs) + #print(times) + embed() -x = 1 - np.exp(t / 11) -print(x) -a, b ,c,d = normalized_JAR(fre) -''' -print(np.logspace(-3, 1, 10)) -embed() \ No newline at end of file diff --git a/step_response_specto.py b/step_response_specto.py index 3203d2c..2ed32df 100644 --- a/step_response_specto.py +++ b/step_response_specto.py @@ -17,6 +17,7 @@ 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' @@ -50,44 +51,49 @@ for infodataset in datasets: for idx, dataset in enumerate(datasets): - datapath = os.path.join(base_path, dataset) - for info, key, time, data in dl.iload_traces(datapath, repro='Beats', before=0.0, after=0.0): - - #plt.plot(time, data[0]) # V(t) - #plt.show() - nfft = 2**18 - spec, freqs, times = specgram(data[0], Fs=1.0/(time[1]-time[0]), detrend='mean', NFFT=nfft, noverlap=nfft*0.8) - - dbspec = 10.0*np.log10(spec) # in dB - - power = dbspec[:, 10] - - 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', extent=(times[0], times[-1], freqs[0], freqs[-1]), aspect='auto', vmin=-80, vmax=-30) - - # 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 - - plt.plot(times, jar4, '-k') - - #plt.plot(freqs, dbspec[:,100], label = '0') - #plt.legend() - plt.ylim(lim0, lim1) - plt.show() - embed() + 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()