From ed3791bc0d2ced310dedb1a1caa35b9a48a4b604 Mon Sep 17 00:00:00 2001 From: wendtalexander Date: Wed, 11 Jan 2023 13:58:25 +0100 Subject: [PATCH] adding first script --- code/analaysis.py | 106 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 code/analaysis.py diff --git a/code/analaysis.py b/code/analaysis.py new file mode 100644 index 0000000..cc6d5c8 --- /dev/null +++ b/code/analaysis.py @@ -0,0 +1,106 @@ +from thunderfish.dataloader import DataLoader as open_data +from thunderfish.powerspectrum import spectrogram, decibel +from IPython import embed +from audioio import play +import matplotlib.pyplot as plt +import numpy as np +import os +from scipy.ndimage import gaussian_filter1d +from modules.filters import bandpass_filter + + +def main(folder): + file = os.path.join(folder, 'traces-grid.raw') + data = open_data(folder, 60.0, 0, channel=-1) + time = np.load(folder + 'times.npy', allow_pickle=True) + freq = np.load(folder + 'fund_v.npy', allow_pickle=True) + ident = np.load(folder + 'ident_v.npy', allow_pickle=True) + idx = np.load(folder + 'idx_v.npy', allow_pickle=True) + + t0 = 3*60*60 + 6*60 + 43.5 + dt = 60 + data_oi = data[t0 * data.samplerate: (t0+dt)*data.samplerate, :] + + for i in [10]: + # getting the spectogramm + spec_power, spec_freqs, spec_times = spectrogram( + data_oi[:, i], ratetime=data.samplerate, freq_resolution=50, overlap_frac=0.0) + fig, ax = plt.subplots(figsize=(20/2.54, 12/2.54)) + ax.pcolormesh(spec_times, spec_freqs, decibel( + spec_power), vmin=-100, vmax=-50) + + for track_id in np.unique(ident): + # window_index for time array in time window + window_index = np.arange(len(idx))[(ident == track_id) & + (time[idx] >= t0) & + (time[idx] <= (t0+dt))] + freq_temp = freq[window_index] + time_temp = time[idx[window_index]] + #mean_freq = np.mean(freq_temp) + #fdata = bandpass_filter(data_oi[:, track_id], data.samplerate, mean_freq-5, mean_freq+200) + ax.plot(time_temp - t0, freq_temp) + + ax.set_ylim(500, 1000) + plt.show() + # filter plot + id = 10. + i = 10 + window_index = np.arange(len(idx))[(ident == id) & + (time[idx] >= t0) & + (time[idx] <= (t0+dt))] + freq_temp = freq[window_index] + time_temp = time[idx[window_index]] + mean_freq = np.mean(freq_temp) + fdata = bandpass_filter(data_oi[:, i], rate=data.samplerate, lowf=mean_freq-5, highf=mean_freq+200) + fig, ax = plt.subplots() + ax.plot(np.arange(len(fdata))/data.samplerate, fdata, marker='*') + #plt.show() + #freqency analyis of filtered data + + time_fdata = np.arange(len(fdata))/data.samplerate + roll_fdata = np.roll(fdata, shift=1) + period_index = np.arange(len(fdata))[(roll_fdata < 0) & (fdata>=0)] + + plt.plot(time_fdata, fdata) + plt.scatter(time_fdata[period_index], fdata[period_index], c='r') + plt.scatter(time_fdata[period_index-1], fdata[period_index-1], c='r') + + upper_bound = np.abs(fdata[period_index]) + lower_bound = np.abs(fdata[period_index-1]) + + upper_times = np.abs(time_fdata[period_index]) + lower_times = np.abs(time_fdata[period_index-1]) + + lower_ratio = lower_bound/(lower_bound+upper_bound) + upper_ratio = upper_bound/(lower_bound+upper_bound) + + time_delta = upper_times-lower_times + true_zero = lower_times + time_delta*lower_ratio + + plt.scatter(true_zero, np.zeros(len(true_zero))) + + # calculate the frequency + inst_freq = 1 / np.diff(true_zero) + filtered_inst_freq = gaussian_filter1d(inst_freq, 0.005) + fig, ax = plt.subplots() + ax.plot(filtered_inst_freq, marker='.') + # in 5 sekunden welcher fisch auf einer elektrode am + + + + + embed() + exit() + + # data of intrests + + # first look at the raw data, channel 11 is important + # fig, ax = plt.subplots(figsize=(20/2.54, 12/2.54)) + # ax.plot(np.arange(len(data_oi[:, i])), data_oi[:, i]) + + + pass + + +if __name__ == '__main__': + main('/Users/acfw/Documents/uni_tuebingen/chirpdetection/gp_benda/data/2022-06-02-10_00/')