Merge branch 'master' of https://whale.am28.uni-tuebingen.de/git/raab/GP2023_chirp_detection

code/chirpdetection.py

@@ -29,8 +29,8 @@ class LoadData:
     def __init__(self, datapath: str) -> None:
 
         # load raw data
-        file = os.path.join(datapath, "traces-grid1.raw")
+        self.file = os.path.join(datapath, "traces-grid1.raw")
-        self.data = DataLoader(file, 60.0, 0, channel=-1)
+        self.data = DataLoader(self.file, 60.0, 0, channel=-1)
         self.samplerate = self.data.samplerate
 
         # load wavetracker files
@@ -40,6 +40,12 @@ class LoadData:
         self.ident = np.load(datapath + "ident_v.npy", allow_pickle=True)
         self.ids = np.unique(self.ident[~np.isnan(self.ident)])
 
+    def __repr__(self) -> str:
+        return f"LoadData({self.file})"
+
+    def __str__(self) -> str:
+        return f"LoadData({self.file})"
+
 
 def instantaneos_frequency(
     signal: np.ndarray, samplerate: int
@@ -62,7 +68,8 @@ def instantaneos_frequency(
     # calculate instantaneos frequency with zero crossings
     roll_signal = np.roll(signal, shift=1)
     time_signal = np.arange(len(signal)) / samplerate
-    period_index = np.arange(len(signal))[(roll_signal < 0) & (signal >= 0)]
+    period_index = np.arange(len(signal))[(
+        roll_signal < 0) & (signal >= 0)][1:-1]
 
     upper_bound = np.abs(signal[period_index])
     lower_bound = np.abs(signal[period_index - 1])
@@ -77,7 +84,7 @@ def instantaneos_frequency(
     true_zero = lower_time + lower_ratio * time_delta
 
     # create new time array
-    inst_freq_time = true_zero[1:] + 0.5 * np.diff(true_zero)
+    inst_freq_time = true_zero[:-1] + 0.5 * np.diff(true_zero)
 
     # compute frequency
     inst_freq = gaussian_filter1d(1 / np.diff(true_zero), 5)
@@ -167,52 +174,74 @@ def main(datapath: str) -> None:
     # load wavetracker files
     time = np.load(datapath + "times.npy", allow_pickle=True)
     freq = np.load(datapath + "fund_v.npy", allow_pickle=True)
+    powers = np.load(datapath + "sign_v.npy", allow_pickle=True)
     idx = np.load(datapath + "idx_v.npy", allow_pickle=True)
     ident = np.load(datapath + "ident_v.npy", allow_pickle=True)
 
     # set time window # <------------------------ Iterate through windows here
-    window_duration = 60 * 5 * data.samplerate    # 5 minutes window
+    window_duration = 5 * data.samplerate          # 5 seconds window
-    window_overlap = 30 * data.samplerate         # 30 seconds overlap
+    window_overlap = 0.5 * data.samplerate         # 30 seconds overlap
     raw_time = np.arange(data.shape[0])
+    t0 = (3 * 60 * 60 + 6 * 60 + 43.5) * data.samplerate
+    dt = 60 * data.samplerate
+
     window_starts = np.arange(
-        raw_time[0], raw_time[-1], window_duration - window_overlap / 2)
+        t0, t0 + dt, window_duration - window_overlap, dtype=int)
+
+    for start_index in window_starts:
+
+        # make t0 and dt
+        t0 = start_index / data.samplerate
+        dt = window_duration / data.samplerate
 
+        # set index window
+        stop_index = start_index + window_duration
 
-    t0 = 3 * 60 * 60 + 6 * 60 + 43.5
+        # t0 = 3 * 60 * 60 + 6 * 60 + 43.5
-    dt = 60
+#     dt = 60
-    start_index = t0 * data.samplerate
+#     start_index = t0 * data.samplerate
-    stop_index = (t0 + dt) * data.samplerate
+#     stop_index = (t0 + dt) * data.samplerate
 
         # load region of interest of raw data file
         data_oi = data[start_index:stop_index, :]
 
+        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 track_id in np.unique(ident[~np.isnan(ident)])[:2]:
+        for i, track_id in enumerate(np.unique(ident[~np.isnan(ident)])[:2]):
 
             # <------------------------------------------ Find best electrodes here
             # <------------------------------------------ Iterate through electrodes
             # get indices for time array in time window
             window_index = np.arange(len(idx))[
-            (ident == track_id) & (time[idx] >= t0) & (time[idx] <= (t0 + dt))
+                (ident == track_id) & (time[idx] >= t0) & (
+                    time[idx] <= (t0 + dt))
             ]
 
-        # filter baseline and above
+            # get tracked frequencies and their times
             freq_temp = freq[window_index]
-        time_temp = time[idx[window_index]]
+            powers_temp = powers[window_index, :]
+            # time_temp = time[idx[window_index]]
+            track_samplerate = np.mean(1 / np.diff(time))
+            expected_duration = ((t0 + dt) - t0) * track_samplerate
 
-        electrode = 10
+            # check if tracked data available in this window
+            if len(freq_temp) < expected_duration * 0.9:
+                continue
 
-        # initialize plot
+            # get best electrode
-        fig, axs = plt.subplots(
+            electrode = np.argsort(np.nanmean(powers_temp, axis=0))[-1]
-            7,
+            # electrode = best_electrodes[0]
-            1,
-            figsize=(20 / 2.54, 12 / 2.54),
-            constrained_layout=True,
-            sharex=True,
-        )
 
             # plot spectrogram
-        plot_spectrogram(axs[0], data_oi[:, electrode], data.samplerate)
+            plot_spectrogram(axs[0, i], data_oi[:, electrode], data.samplerate)
 
             # plot wavetracker tracks to spectrogram
             # for track_id in np.unique(ident):  # <---------- Find freq gaps later
@@ -242,30 +271,37 @@ def main(datapath: str) -> None:
             )
 
             # plot waveform of filtered signal
-        axs[2].plot(np.arange(len(baseline)) / data.samplerate, baseline)
+            axs[2, i].plot(np.arange(len(baseline)) /
+                           data.samplerate, baseline, c="k")
 
             # plot instatneous frequency
-        # broad_baseline = bandpass_filter(data_oi[:, electrode], data.samplerate, lowf=np.mean(
+            broad_baseline = bandpass_filter(data_oi[:, electrode], data.samplerate, lowf=np.mean(
-        #     freq_temp)-5, highf=np.mean(freq_temp)+200)
+                freq_temp)-5, highf=np.mean(freq_temp)+100)
 
             baseline_freq_time, baseline_freq = instantaneos_frequency(
                 baseline, data.samplerate
             )
-        axs[1].plot(baseline_freq_time, baseline_freq)
+            axs[1, i].plot(baseline_freq_time, baseline_freq -
+                           np.median(baseline_freq), marker=".")
 
             # plot waveform of filtered search signal
-        axs[3].plot(np.arange(len(baseline)) / data.samplerate, search)
+            axs[3, i].plot(np.arange(len(baseline)) / data.samplerate, search)
 
             # compute envelopes
             cutoff = 25
             baseline_envelope = envelope(baseline, data.samplerate, cutoff)
-        axs[2].plot(
+            axs[2, i].plot(
                 np.arange(len(baseline)) / data.samplerate,
                 baseline_envelope,
                 c="orange",
             )
+            axs[2, i].plot(
+                np.arange(len(baseline)) / data.samplerate,
+                broad_baseline,
+                c="green",
+            )
             search_envelope = envelope(search, data.samplerate, cutoff)
-        axs[3].plot(
+            axs[3, i].plot(
                 np.arange(len(baseline)) / data.samplerate,
                 search_envelope,
                 c="orange",
@@ -293,21 +329,20 @@ def main(datapath: str) -> None:
             )
 
             # plot filtered and rectified envelope
-        axs[4].plot(
+            axs[4, i].plot(
                 np.arange(len(baseline)) / data.samplerate, baseline_envelope
             )
 
-        axs[5].plot(np.arange(len(baseline)) /
+            axs[5, i].plot(np.arange(len(baseline)) /
                            data.samplerate, search_envelope)
 
-        axs[6].plot(baseline_freq_time, np.abs(inst_freq_filtered))
+            axs[6, i].plot(baseline_freq_time, np.abs(inst_freq_filtered))
 
             # detect peaks baseline_enelope
-        embed()
             prominence = iqr(baseline_envelope)
             baseline_peaks, _ = find_peaks(
                 baseline_envelope, prominence=prominence)
-        axs[4].scatter(
+            axs[4, i].scatter(
                 (np.arange(len(baseline)) / data.samplerate)[baseline_peaks],
                 baseline_envelope[baseline_peaks],
                 c="red",
@@ -315,27 +350,28 @@ def main(datapath: str) -> None:
 
             # detect peaks search_envelope
             search_peaks, _ = find_peaks(search_envelope, height=0.0001)
-        axs[5].scatter(
+            axs[5, i].scatter(
                 (np.arange(len(baseline)) / data.samplerate)[search_peaks],
                 search_envelope[search_peaks],
                 c="red",
             )
 
             # detect peaks inst_freq_filtered
-        inst_freq_peaks, _ = find_peaks(np.abs(inst_freq_filtered), height=2)
+            inst_freq_peaks, _ = find_peaks(
-        axs[6].scatter(
+                np.abs(inst_freq_filtered), height=2)
+            axs[6, i].scatter(
                 baseline_freq_time[inst_freq_peaks],
                 np.abs(inst_freq_filtered)[inst_freq_peaks],
                 c="red",
             )
 
-        axs[0].set_title("Spectrogram")
+            axs[0, i].set_title("Spectrogram")
-        axs[1].set_title("Fitered baseline instanenous frequency")
+            axs[1, i].set_title("Fitered baseline instanenous frequency")
-        axs[2].set_title("Fitered baseline")
+            axs[2, i].set_title("Fitered baseline")
-        axs[3].set_title("Fitered above")
+            axs[3, i].set_title("Fitered above")
-        axs[4].set_title("Filtered envelope of baseline envelope")
+            axs[4, i].set_title("Filtered envelope of baseline envelope")
-        axs[5].set_title("Search envelope")
+            axs[5, i].set_title("Search envelope")
-        axs[6].set_title("Filtered absolute instantaneous frequency")
+            axs[6, i].set_title("Filtered absolute instantaneous frequency")
 
         plt.show()