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.0 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/" )