From feb48211fc7b8ee77d104c380e0b5b2c6be3d0d2 Mon Sep 17 00:00:00 2001 From: weygoldt <88969563+weygoldt@users.noreply.github.com> Date: Fri, 13 Jan 2023 15:37:52 +0100 Subject: [PATCH] started new overlap --- code/chirpdetection.py | 429 ++++++++++++++++++------------------ code/chirpdetector_conf.yml | 3 +- 2 files changed, 218 insertions(+), 214 deletions(-) diff --git a/code/chirpdetection.py b/code/chirpdetection.py index 900e251..c15e1e0 100644 --- a/code/chirpdetection.py +++ b/code/chirpdetection.py @@ -187,22 +187,9 @@ def main(datapath: str) -> None: # start_index = t0 * data.samplerate # stop_index = (t0 + dt) * data.samplerate - fig, axs = plt.subplots( - 7, - 2, - figsize=(20 / 2.54, 12 / 2.54), - constrained_layout=True, - sharex=True, - sharey='row', - ) - # iterate through all fish for i, track_id in enumerate(np.unique(ident[~np.isnan(ident)])[:2]): - # load region of interest of raw data file - data_oi = data[start_index:stop_index, :] - time_oi = raw_time[start_index:stop_index] - # get indices for time array in time window window_index = np.arange(len(idx))[ (ident == track_id) & (time[idx] >= t0) & ( @@ -220,213 +207,229 @@ def main(datapath: str) -> None: if len(freq_temp) < expected_duration * 0.9: continue + fig, axs = plt.subplots( + 7, + config.electrodes, + figsize=(20 / 2.54, 12 / 2.54), + constrained_layout=True, + sharex=True, + sharey='row', + ) # get best electrode - electrode = np.argsort(np.nanmean(powers_temp, axis=0))[-1] + best_electrodes = np.argsort(np.nanmean( + powers_temp, axis=0))[-config.electrodes:] # <------------------------------------------ Iterate through electrodes - # plot wavetracker tracks to spectrogram - # for track_id in np.unique(ident): # <---------- Find freq gaps later - # here - - # # get indices 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]] - - # axs[0].plot(time_temp-t0, freq_temp, lw=2) - # axs[0].set_ylim(500, 1000) - - # track_id = ids - - # frequency where second filter filters - search_freq = 50 - - # filter baseline and above - baseline, search = double_bandpass( - data_oi[:, electrode], data.samplerate, freq_temp, search_freq - ) - - # compute instantaneous frequency on broad signal - broad_baseline = bandpass_filter( - data_oi[:, electrode], - data.samplerate, - lowf=np.mean(freq_temp)-5, - highf=np.mean(freq_temp)+100 - ) - - # compute instantaneous frequency on narrow signal - baseline_freq_time, baseline_freq = instantaneos_frequency( - baseline, data.samplerate - ) - - # compute envelopes - baseline_envelope = envelope( - baseline, data.samplerate, config.envelope_cutoff) - search_envelope = envelope( - search, data.samplerate, config.envelope_cutoff) - - # highpass filter envelopes - baseline_envelope = highpass_filter( - baseline_envelope, - data.samplerate, - config.envelope_highpass_cutoff - ) - - baseline_envelope = np.abs(baseline_envelope) - # search_envelope = highpass_filter( - # search_envelope, - # data.samplerate, - # config.envelope_highpass_cutoff - # ) - - # envelopes of filtered envelope of filtered baseline - baseline_envelope = envelope( - np.abs(baseline_envelope), - data.samplerate, - config.envelope_envelope_cutoff - ) + for i, electrode in enumerate(best_electrodes): + + # load region of interest of raw data file + data_oi = data[start_index:stop_index, :] + time_oi = raw_time[start_index:stop_index] + + # plot wavetracker tracks to spectrogram + # for track_id in np.unique(ident): # <---------- Find freq gaps later + # here + + # # get indices 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]] + + # axs[0].plot(time_temp-t0, freq_temp, lw=2) + # axs[0].set_ylim(500, 1000) + + # track_id = ids + + # frequency where second filter filters + search_freq = 50 + + # filter baseline and above + baseline, search = double_bandpass( + data_oi[:, electrode], data.samplerate, freq_temp, search_freq + ) + + # compute instantaneous frequency on broad signal + broad_baseline = bandpass_filter( + data_oi[:, electrode], + data.samplerate, + lowf=np.mean(freq_temp)-5, + highf=np.mean(freq_temp)+100 + ) + + # compute instantaneous frequency on narrow signal + baseline_freq_time, baseline_freq = instantaneos_frequency( + baseline, data.samplerate + ) + + # compute envelopes + baseline_envelope = envelope( + baseline, data.samplerate, config.envelope_cutoff) + search_envelope = envelope( + search, data.samplerate, config.envelope_cutoff) + + # highpass filter envelopes + baseline_envelope = highpass_filter( + baseline_envelope, + data.samplerate, + config.envelope_highpass_cutoff + ) + + baseline_envelope = np.abs(baseline_envelope) + # search_envelope = highpass_filter( + # search_envelope, + # data.samplerate, + # config.envelope_highpass_cutoff + # ) + + # envelopes of filtered envelope of filtered baseline + baseline_envelope = envelope( + np.abs(baseline_envelope), + data.samplerate, + config.envelope_envelope_cutoff + ) # search_envelope = bandpass_filter( # search_envelope, data.samplerate, lowf=lowf, highf=highf) - # bandpass filter the instantaneous - inst_freq_filtered = bandpass_filter( - baseline_freq, - data.samplerate, - lowf=config.instantaneous_lowf, - highf=config.instantaneous_highf - ) - - # test taking the log of the envelopes - # baseline_envelope = np.log(baseline_envelope) - # search_envelope = np.log(search_envelope) - - # CUT OFF OVERLAP ------------------------------------------------- - - # cut off first and last 0.5 * overlap at start and end - valid = np.arange( - int(window_overlap / 2), len(baseline_envelope) - - int(window_overlap / 2) - ) - baseline_envelope = baseline_envelope[valid] - search_envelope = search_envelope[valid] - - # get inst freq valid snippet - valid_t0 = int(window_overlap / 2) / data.samplerate - valid_t1 = baseline_freq_time[-1] - \ - (int(window_overlap / 2) / data.samplerate) - - inst_freq_filtered = inst_freq_filtered[(baseline_freq_time >= valid_t0) & ( - baseline_freq_time <= valid_t1)] - - baseline_freq = baseline_freq[(baseline_freq_time >= valid_t0) & ( - baseline_freq_time <= valid_t1)] - - baseline_freq_time = baseline_freq_time[(baseline_freq_time >= valid_t0) & ( - baseline_freq_time <= valid_t1)] + t0 - - # overwrite raw time to valid region - time_oi = time_oi[valid] - baseline = baseline[valid] - broad_baseline = broad_baseline[valid] - search = search[valid] - - # PEAK DETECTION -------------------------------------------------- - - # detect peaks baseline_enelope - prominence = np.percentile( - baseline_envelope, config.baseline_prominence_percentile) - baseline_peaks, _ = find_peaks( - np.abs(baseline_envelope), prominence=prominence) - - # detect peaks search_envelope - prominence = np.percentile( - search_envelope, config.search_prominence_percentile) - search_peaks, _ = find_peaks( - search_envelope, prominence=prominence) - - # detect peaks inst_freq_filtered - prominence = np.percentile( - inst_freq_filtered, config.instantaneous_prominence_percentile) - inst_freq_peaks, _ = find_peaks( - np.abs(inst_freq_filtered), prominence=prominence) - - # PLOT ------------------------------------------------------------ - - # plot spectrogram - plot_spectrogram( - axs[0, i], data_oi[:, electrode], data.samplerate, t0) - - # plot baseline instantaneos frequency - axs[1, i].plot(baseline_freq_time, baseline_freq - - np.median(baseline_freq), marker=".") - - # plot waveform of filtered signal - axs[2, i].plot(time_oi, baseline, c="k") - - # plot narrow filtered baseline - axs[2, i].plot( - time_oi, - baseline_envelope, - c="orange", - ) - - # plot broad filtered baseline - axs[2, i].plot( - time_oi, - broad_baseline, - c="green", - ) - - # plot waveform of filtered search signal - axs[3, i].plot(time_oi, search) - - # plot envelope of search signal - axs[3, i].plot( - time_oi, - search_envelope, - c="orange", - ) - - # plot filtered and rectified envelope - axs[4, i].plot(time_oi, baseline_envelope) - axs[4, i].scatter( - (time_oi)[baseline_peaks], - baseline_envelope[baseline_peaks], - c="red", - ) - - # plot envelope of search signal - axs[5, i].plot(time_oi, search_envelope) - axs[5, i].scatter( - (time_oi)[search_peaks], - search_envelope[search_peaks], - c="red", - ) - - # plot filtered instantaneous frequency - axs[6, i].plot(baseline_freq_time, np.abs(inst_freq_filtered)) - axs[6, i].scatter( - baseline_freq_time[inst_freq_peaks], - np.abs(inst_freq_filtered)[inst_freq_peaks], - c="red", - ) - - axs[6, i].set_xlabel("Time [s]") - axs[0, i].set_title("Spectrogram") - axs[1, i].set_title("Fitered baseline instanenous frequency") - axs[2, i].set_title("Fitered baseline") - axs[3, i].set_title("Fitered above") - axs[4, i].set_title("Filtered envelope of baseline envelope") - axs[5, i].set_title("Search envelope") - axs[6, i].set_title("Filtered absolute instantaneous frequency") - - plt.show() + # bandpass filter the instantaneous + inst_freq_filtered = bandpass_filter( + baseline_freq, + data.samplerate, + lowf=config.instantaneous_lowf, + highf=config.instantaneous_highf + ) + + # test taking the log of the envelopes + # baseline_envelope = np.log(baseline_envelope) + # search_envelope = np.log(search_envelope) + + # CUT OFF OVERLAP ------------------------------------------------- + + # cut off first and last 0.5 * overlap at start and end + valid = np.arange( + int(window_overlap / 2), len(baseline_envelope) - + int(window_overlap / 2) + ) + baseline_envelope = baseline_envelope[valid] + search_envelope = search_envelope[valid] + + # get inst freq valid snippet + valid_t0 = int(window_overlap / 2) / data.samplerate + valid_t1 = baseline_freq_time[-1] - \ + (int(window_overlap / 2) / data.samplerate) + + inst_freq_filtered = inst_freq_filtered[(baseline_freq_time >= valid_t0) & ( + baseline_freq_time <= valid_t1)] + + baseline_freq = baseline_freq[(baseline_freq_time >= valid_t0) & ( + baseline_freq_time <= valid_t1)] + + baseline_freq_time = baseline_freq_time[(baseline_freq_time >= valid_t0) & ( + baseline_freq_time <= valid_t1)] + t0 + + # overwrite raw time to valid region + time_oi = time_oi[valid] + baseline = baseline[valid] + broad_baseline = broad_baseline[valid] + search = search[valid] + + # PEAK DETECTION -------------------------------------------------- + + # detect peaks baseline_enelope + prominence = np.percentile( + baseline_envelope, config.baseline_prominence_percentile) + baseline_peaks, _ = find_peaks( + np.abs(baseline_envelope), prominence=prominence) + + # detect peaks search_envelope + prominence = np.percentile( + search_envelope, config.search_prominence_percentile) + search_peaks, _ = find_peaks( + search_envelope, prominence=prominence) + + # detect peaks inst_freq_filtered + prominence = np.percentile( + inst_freq_filtered, config.instantaneous_prominence_percentile) + inst_freq_peaks, _ = find_peaks( + np.abs(inst_freq_filtered), prominence=prominence) + + # PLOT ------------------------------------------------------------ + + # plot spectrogram + plot_spectrogram( + axs[0, i], data_oi[:, electrode], data.samplerate, t0) + + # plot baseline instantaneos frequency + axs[1, i].plot(baseline_freq_time, baseline_freq - + np.median(baseline_freq), marker=".") + + # plot waveform of filtered signal + axs[2, i].plot(time_oi, baseline, c="k") + + # plot narrow filtered baseline + axs[2, i].plot( + time_oi, + baseline_envelope, + c="orange", + ) + + # plot broad filtered baseline + axs[2, i].plot( + time_oi, + broad_baseline, + c="green", + ) + + # plot waveform of filtered search signal + axs[3, i].plot(time_oi, search) + + # plot envelope of search signal + axs[3, i].plot( + time_oi, + search_envelope, + c="orange", + ) + + # plot filtered and rectified envelope + axs[4, i].plot(time_oi, baseline_envelope) + axs[4, i].scatter( + (time_oi)[baseline_peaks], + baseline_envelope[baseline_peaks], + c="red", + ) + + # plot envelope of search signal + axs[5, i].plot(time_oi, search_envelope) + axs[5, i].scatter( + (time_oi)[search_peaks], + search_envelope[search_peaks], + c="red", + ) + + # plot filtered instantaneous frequency + axs[6, i].plot(baseline_freq_time, np.abs(inst_freq_filtered)) + axs[6, i].scatter( + baseline_freq_time[inst_freq_peaks], + np.abs(inst_freq_filtered)[inst_freq_peaks], + c="red", + ) + + axs[6, i].set_xlabel("Time [s]") + axs[0, i].set_title("Spectrogram") + axs[1, i].set_title("Fitered baseline instanenous frequency") + axs[2, i].set_title("Fitered baseline") + axs[3, i].set_title("Fitered above") + axs[4, i].set_title("Filtered envelope of baseline envelope") + axs[5, i].set_title("Search envelope") + axs[6, i].set_title( + "Filtered absolute instantaneous frequency") + + plt.show() if __name__ == "__main__": diff --git a/code/chirpdetector_conf.yml b/code/chirpdetector_conf.yml index 55a6d08..81a048c 100644 --- a/code/chirpdetector_conf.yml +++ b/code/chirpdetector_conf.yml @@ -1,6 +1,7 @@ # Duration and overlap of the analysis window in seconds window: 5 -overlap: 0.5 +overlap: 1 +edges: 0.25 # Number of electrodes to go over electrodes: 3