removing the inst freq as a feature, and fintuning

.python-version

@@ -1 +0,0 @@
-chirpdetection

chirp_instantaneous_freq/chirp_exploration.py

@@ -0,0 +1,241 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# # Why can the instantaneous frequency of a band-pass filtered chirp recording go down ...
+# ... if a chirp is an up-modulation of the frequency? 
+# 
+# This is the question we try to answer in this notebook
+
+# In[14]:
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+import thunderfish.fakefish as ff 
+from filters import instantaneous_frequency, bandpass_filter
+#get_ipython().run_line_magic('matplotlib', 'qt')
+
+# parameters that stay the same
+samplerate = 20000
+duration = 0.2
+chirp_freq = 5
+smooth = 3
+
+
+# In[7]:
+
+
+def make_chirp(eodf, size, width, kurtosis, contrast, phase0):
+
+    chirp_trace, amp = ff.chirps(
+        eodf = eodf,
+        samplerate = samplerate,
+        duration = duration,
+        chirp_freq = chirp_freq,
+        chirp_size = size,
+        chirp_width = width,
+        chirp_kurtosis = kurtosis,
+        chirp_contrast = contrast,
+    )
+
+    chirp = ff.wavefish_eods(
+        fish = 'Alepto',
+        frequency = chirp_trace,
+        samplerate = samplerate,
+        duration = duration,
+        phase0 = phase0,
+        noise_std = 0,
+    )
+
+    chirp *= amp
+
+    return chirp_trace, chirp
+
+def filtered_chirp(eodf, size, width, kurtosis, contrast, phase0):
+
+    time = np.arange(0, duration, 1/samplerate)
+    chirp_trace, chirp = make_chirp(
+        eodf = eodf, 
+        size = size, 
+        width = width, 
+        kurtosis = kurtosis, 
+        contrast = contrast, 
+        phase0 = phase0,
+    )
+
+    # apply filters
+    narrow_filtered = bandpass_filter(chirp, samplerate, eodf-10, eodf+10)
+    narrow_freqtime, narrow_freq = instantaneous_frequency(narrow_filtered, samplerate, smooth)
+    broad_filtered = bandpass_filter(chirp, samplerate, eodf-300, eodf+300)
+    broad_freqtime, broad_freq = instantaneous_frequency(broad_filtered, samplerate, smooth)
+
+    original = (time, chirp_trace, chirp)
+    broad = (broad_freqtime, broad_freq, broad_filtered)
+    narrow = (narrow_freqtime, narrow_freq, narrow_filtered)
+
+    return original, broad, narrow
+
+def plot(original, broad, narrow, axs):
+
+    axs[0].plot(original[0], original[1], label='chirp trace')
+    axs[0].plot(broad[0], broad[1], label='broad filtered')
+    axs[0].plot(narrow[0], narrow[1], label='narrow filtered')
+    axs[1].plot(original[0], original[2], label='unfiltered')
+    axs[1].plot(original[0], broad[2], label='broad filtered')
+    axs[1].plot(original[0], narrow[2], label='narrow filtered')
+
+original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, 0.1, 0)
+fig, axs = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
+plot(original, broad, narrow, axs)
+fig.align_labels()
+plt.show()
+
+
+# ## Chirp size
+# now that we have established an easy way to simulate and plot the chirps, lets change the chirp size and see how the narrow-filtered instantaneous frequency changes.
+
+# In[17]:
+
+
+sizes = np.arange(10, 600, 20)[:10]
+fig, axs = plt.subplots(2, len(sizes), figsize=(20, 5), sharex=True, sharey='row')
+integrals = []
+
+for i, size in enumerate(sizes):
+    original, broad, narrow = filtered_chirp(600, size, 0.02, 1, 0.1, 0)
+
+    integral = np.sum(original[1]-600)/(20000)
+    integrals.append(integral)
+
+    plot(original, broad, narrow, axs[:, i])
+    axs[:, i][0].set_xlim(0.06, 0.14)
+    axs[0, i].set_title(np.round(integral, 3))
+    print(f'size {size} Hz; Integral {np.round(integral,3)}')
+ 
+fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)
+axs[0,0].set_ylabel('frequency [Hz]')
+axs[1,0].set_ylabel('amplitude [a.u.]')
+fig.supxlabel('time [s]')
+fig.align_labels()
+plt.show()
+
+
+# ## Chirp width
+
+# In[9]:
+
+
+widths = np.arange(0.02, 0.08, 0.005)
+fig, axs = plt.subplots(2, len(widths), figsize=(10, 5), sharex=True, sharey='row')
+integrals = []
+
+for i, width in enumerate(widths):
+    if i > 9:
+        break
+
+    original, broad, narrow = filtered_chirp(600, 100, width, 1, 0.1, 0)
+
+    integral = np.sum(original[1]-600)/(20000)
+
+    plot(original, broad, narrow, axs[:, i])
+    axs[:, i][0].set_xlim(0.06, 0.14)
+    axs[0, i].set_title(f'width {np.round(width, 2)} s')
+    print(f'width {width} s; Integral {np.round(integral, 3)}')
+ 
+fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)
+axs[0,0].set_ylabel('frequency [Hz]')
+axs[1,0].set_ylabel('amplitude [a.u.]')
+fig.supxlabel('time [s]')
+fig.align_labels()
+plt.show()
+
+
+# ## Chirp kurtosis
+
+# In[10]:
+
+
+kurtosiss = np.arange(0, 20, 1.6)
+fig, axs = plt.subplots(2, len(kurtosiss), figsize=(10, 5), sharex=True, sharey='row')
+integrals = []
+
+for i, kurtosis in enumerate(kurtosiss):
+
+    original, broad, narrow = filtered_chirp(600, 100, 0.02, kurtosis, 0.1, 0)
+
+    integral = np.sum(original[1]-600)/(20000)
+
+    plot(original, broad, narrow, axs[:, i])
+    axs[:, i][0].set_xlim(0.06, 0.14)
+    axs[0, i].set_title(f'kurt {np.round(kurtosis, 2)}')
+    print(f'kurt {kurtosis}; Integral {np.round(integral, 3)}')
+ 
+fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)
+axs[0,0].set_ylabel('frequency [Hz]')
+axs[1,0].set_ylabel('amplitude [a.u.]')
+fig.supxlabel('time [s]')
+fig.align_labels()
+plt.show()
+
+
+# ## Chirp contrast
+
+# In[11]:
+
+
+contrasts = np.arange(0.0, 1.1, 0.1)
+fig, axs = plt.subplots(2, len(sizes), figsize=(10, 5), sharex=True, sharey='row')
+integrals = []
+
+for i, contrast in enumerate(contrasts):
+    if i > 9:
+        break
+    original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, contrast, 0)
+
+    integral = np.trapz(original[2], original[0])
+    integrals.append(integral)
+
+    plot(original, broad, narrow, axs[:, i])
+    axs[:, i][0].set_xlim(0.06, 0.14)
+    axs[0, i].set_title(f'contr {np.round(contrast, 2)}')
+ 
+fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)
+axs[0,0].set_ylabel('frequency [Hz]')
+axs[1,0].set_ylabel('amplitude [a.u.]')
+fig.supxlabel('time [s]')
+fig.align_labels()
+plt.show()
+
+
+# ## Chirp phase 
+
+# In[12]:
+
+
+phases = np.arange(0.0, 2 * np.pi, 0.2)
+fig, axs = plt.subplots(2, len(sizes), figsize=(10, 5), sharex=True, sharey='row')
+integrals = []
+for i, phase in enumerate(phases):
+    if i > 9:
+        break
+
+    original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, 0.1, phase)
+
+    integral = np.trapz(original[2], original[0])
+    integrals.append(integral)
+
+    plot(original, broad, narrow, axs[:, i])
+    axs[:, i][0].set_xlim(0.06, 0.14)
+    axs[0, i].set_title(f'phase {np.round(phase, 2)}')
+
+ 
+fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)
+axs[0,0].set_ylabel('frequency [Hz]')
+axs[1,0].set_ylabel('amplitude [a.u.]')
+fig.supxlabel('time [s]')
+fig.align_labels()
+plt.show()
+
+
+# These experiments show, that the narrow filtered instantaneous freuqency only switches its sign, when the integral of the instantaneous frequency (that was used to make the signal)
+# changes. Specifically, when the instantaneous frequency is 0.57, 1.57, 2.57 etc., the sign swithes. 

code/chirpdetection.py

@@ -749,8 +749,9 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
     multiwindow_chirps = []
     multiwindow_ids = []
 
-    for st, window_start_index in enumerate(window_start_indices):
+    for st, window_start_index in enumerate(window_start_indices[1853:]):
-        logger.info(f"Processing window {st+1} of {len(window_start_indices)}")
+
+        logger.info(f"Processing window {st} of {len(window_start_indices)}")
 
         window_start_seconds = window_start_index / data.raw_rate
         window_duration_seconds = window_duration / data.raw_rate
@@ -914,17 +915,19 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
                     baseline_frequency - np.median(baseline_frequency)
                 )
 
-                # check if there is at least one superthreshold peak on the
+                # # check if there is at least one superthreshold peak on the
-                # instantaneous and exit the loop if not. This is used to
+                # # instantaneous and exit the loop if not. This is used to
-                # prevent windows that do definetely not include a chirp
+                # # prevent windows that do definetely not include a chirp
-                # to enter normalization, where small changes due to noise
+                # # to enter normalization, where small changes due to noise
-                # would be amplified
+                # # would be amplified
 
-                if not has_chirp(
+                # if not has_chirp(
-                    baseline_frequency_filtered[amplitude_mask],
+                #     baseline_frequency_filtered[amplitude_mask],
-                    config.baseline_frequency_peakheight,
+                #     config.baseline_frequency_peakheight,
-                ):
+                # ):                
-                    continue
+                #     logger.warning(
+                #     f"Amplitude to small for the chirp detection of track {track_id} window {st},")
+                #     continue
 
                 # CUT OFF OVERLAP ---------------------------------------------
 
@@ -974,11 +977,11 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
                 # normalize all three feature arrays to the same range to make
                 # peak detection simpler
 
-                baseline_envelope = minmaxnorm([baseline_envelope])[0]
+                # baseline_envelope = minmaxnorm([baseline_envelope])[0]
-                search_envelope = minmaxnorm([search_envelope])[0]
+                # search_envelope = minmaxnorm([search_envelope])[0]
-                baseline_frequency_filtered = minmaxnorm(
+                # baseline_frequency_filtered = minmaxnorm(
-                    [baseline_frequency_filtered]
+                #     [baseline_frequency_filtered]
-                )[0]
+                # )[0]
 
                 # PEAK DETECTION ----------------------------------------------
 
@@ -996,6 +999,7 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
                     prominence=config.frequency_prominence,
                 )
             
+
                 # DETECT CHIRPS IN SEARCH WINDOW ------------------------------
 
                 # get the peak timestamps from the peak indices
@@ -1010,11 +1014,10 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
 
                 # check if one list is empty and if so, skip to the next
                 # electrode because a chirp cannot be detected if one is empty
-
                 one_feature_empty = (
                     len(baseline_peak_timestamps) == 0
                     or len(search_peak_timestamps) == 0
-                    or len(frequency_peak_timestamps) == 0
+                    #or len(frequency_peak_timestamps) == 0
                 )
 
                 if one_feature_empty and (debug == "false"):
@@ -1026,15 +1029,16 @@ def chirpdetection(datapath: str, plot: str, debug: str = "false") -> None:
                 sublists = [
                     list(baseline_peak_timestamps),
                     list(search_peak_timestamps),
-                    list(frequency_peak_timestamps),
+                    #list(frequency_peak_timestamps),
                 ]
 
                 singleelectrode_chirps = group_timestamps(
                     sublists=sublists,
-                    at_least_in=3,
+                    at_least_in=2,
                     difference_threshold=config.chirp_window_threshold,
                 )
 
+
                 # check it there are chirps detected after grouping, continue
                 # with the loop if not
 
@@ -1185,4 +1189,4 @@ if __name__ == "__main__":
     # datapath = "/home/weygoldt/Data/uni/efishdata/2016-colombia/fishgrid/2016-04-09-22_25/"
     # datapath = "/home/weygoldt/Data/uni/chirpdetection/GP2023_chirp_detection/data/mount_data/2020-03-13-10_00/"
     datapath = "../data/2022-06-02-10_00/"
-    chirpdetection(datapath, plot="save", debug="false")
+    chirpdetection(datapath, plot="show", debug="false")

code/chirpdetector_conf.yml

@@ -38,9 +38,9 @@ search_envelope_cutoff: 10  # search envelope estimation cufoff
 # search_prominence: 0.000004   # peak prominence threshold for search envelope
 # frequency_prominence: 2       # peak prominence threshold for baseline freq
 
-baseline_prominence: 0.3  # peak prominence threshold for baseline envelope
+baseline_prominence: 0.00005  # peak prominence threshold for baseline envelope
-search_prominence: 0.3  # peak prominence threshold for search envelope
+search_prominence: 0.000005  # peak prominence threshold for search envelope
-frequency_prominence: 0.3  # peak prominence threshold for baseline freq
+frequency_prominence: 1  # peak prominence threshold for baseline freq
 
 # Classify events as chirps if they are less than this time apart
 chirp_window_threshold: 0.02