Merge branch 'master' into behaviour

code/chirpdetector_conf.yml

@@ -1,3 +1,4 @@
+# directory setup
 dataroot: "../data/"
 outputdir: "../output/"
 
@@ -8,41 +9,38 @@ edge: 0.25
 
 # Number of electrodes to go over
 number_electrodes: 3
+minimum_electrodes: 2
 
-# Boundary for search frequency in Hz
+# Search window bandwidth and minimal baseline bandwidth
-search_boundary: 100
+minimal_bandwidth: 10
 
-# Cutoff frequency for envelope estimation by lowpass filter
+# Instantaneous frequency smoothing usint a gaussian kernel of this width
-envelope_cutoff: 25
+baseline_frequency_smoothing: 5
 
-# Cutoff frequency for envelope highpass filter
+# Baseline processing parameters
-envelope_highpass_cutoff: 3
+baseline_envelope_cutoff: 25
+baseline_envelope_bandpass_lowf: 4
+baseline_envelope_bandpass_highf: 100
+baseline_envelope_envelope_cutoff: 4
 
-# Cutoff frequency for envelope of envelope
+# search envelope processing parameters
-envelope_envelope_cutoff: 5
+search_envelope_cutoff: 5
 
 # Instantaneous frequency bandpass filter cutoff frequencies
-instantaneous_lowf: 15
+baseline_frequency_highpass_cutoff: 0.000005
-instantaneous_highf: 8000
+baseline_frequency_envelope_cutoff: 0.000005
 
-# Baseline envelope peak detection parameters
+# peak detecion parameters
-baseline_prominence_percentile: 90
+prominence: 0.005
-
-# Search envelope peak detection parameters
-search_prominence_percentile: 90
-
-# Instantaneous frequency peak detection parameters
-instantaneous_prominence_percentile: 90
 
 # search freq parameter
-search_df_lower: 25
+search_df_lower: 20
 search_df_upper: 100
 search_res: 1
-search_freq_percentiles:
+search_bandwidth: 10
-  - 5
-  - 95
 default_search_freq: 50
 
+# Classify events as chirps if they are less than this time apart
 chirp_window_threshold: 0.05
 
 

code/modules/datahandling.py

@@ -1,5 +1,59 @@
 import numpy as np
 from typing import List, Any
+from scipy.ndimage import gaussian_filter1d
+
+
+def instantaneous_frequency(
+    signal: np.ndarray,
+    samplerate: int,
+    smoothing_window: int,
+) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Compute the instantaneous frequency of a signal that is approximately
+    sinusoidal and symmetric around 0.
+
+    Parameters
+    ----------
+    signal : np.ndarray
+        Signal to compute the instantaneous frequency from.
+    samplerate : int
+        Samplerate of the signal.
+    smoothing_window : int
+        Window size for the gaussian filter.
+
+    Returns
+    -------
+    tuple[np.ndarray, np.ndarray]
+
+    """
+    # calculate instantaneous 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)][
+        1:-1
+    ]
+
+    upper_bound = np.abs(signal[period_index])
+    lower_bound = np.abs(signal[period_index - 1])
+    upper_time = np.abs(time_signal[period_index])
+    lower_time = np.abs(time_signal[period_index - 1])
+
+    # create ratio
+    lower_ratio = lower_bound / (lower_bound + upper_bound)
+
+    # appy to time delta
+    time_delta = upper_time - lower_time
+    true_zero = lower_time + lower_ratio * time_delta
+
+    # create new time array
+    instantaneous_frequency_time = true_zero[:-1] + 0.5 * np.diff(true_zero)
+
+    # compute frequency
+    instantaneous_frequency = gaussian_filter1d(
+        1 / np.diff(true_zero), smoothing_window
+    )
+
+    return instantaneous_frequency_time, instantaneous_frequency
 
 
 def purge_duplicates(
@@ -64,7 +118,7 @@ def purge_duplicates(
 
 
 def group_timestamps(
-    sublists: List[List[float]], n: int, threshold: float
+    sublists: List[List[float]], at_least_in: int, difference_threshold: float
 ) -> List[float]:
     """
     Groups timestamps that are less than `threshold` milliseconds apart from
@@ -100,7 +154,7 @@ def group_timestamps(
 
     # Group timestamps that are less than threshold milliseconds apart
     for i in range(1, len(timestamps)):
-        if timestamps[i] - timestamps[i - 1] < threshold:
+        if timestamps[i] - timestamps[i - 1] < difference_threshold:
             current_group.append(timestamps[i])
         else:
             groups.append(current_group)
@@ -111,7 +165,7 @@ def group_timestamps(
     # Retain only groups that contain at least n timestamps
     final_groups = []
     for group in groups:
-        if len(group) >= n:
+        if len(group) >= at_least_in:
             final_groups.append(group)
 
     # Calculate the mean of each group

code/modules/filters.py

@@ -3,8 +3,8 @@ import numpy as np
 
 
 def bandpass_filter(
-        data: np.ndarray,
+        signal: np.ndarray,
-        rate: float,
+        samplerate: float,
         lowf: float,
         highf: float,
 ) -> np.ndarray:
@@ -12,7 +12,7 @@ def bandpass_filter(
 
     Parameters
     ----------
-    data : np.ndarray
+    signal : np.ndarray
         The data to be filtered
     rate : float
         The sampling rate
@@ -26,21 +26,22 @@ def bandpass_filter(
     np.ndarray
         The filtered data
     """
-    sos = butter(2, (lowf, highf), "bandpass", fs=rate, output="sos")
+    sos = butter(2, (lowf, highf), "bandpass", fs=samplerate, output="sos")
-    fdata = sosfiltfilt(sos, data)
+    filtered_signal = sosfiltfilt(sos, signal)
-    return fdata
+
+    return filtered_signal
 
 
 def highpass_filter(
-    data: np.ndarray,
+    signal: np.ndarray,
-    rate: float,
+    samplerate: float,
     cutoff: float,
 ) -> np.ndarray:
     """Highpass filter a signal.
 
     Parameters
     ----------
-    data : np.ndarray
+    signal : np.ndarray
         The data to be filtered
     rate : float
         The sampling rate
@@ -52,14 +53,15 @@ def highpass_filter(
     np.ndarray
         The filtered data
     """
-    sos = butter(2, cutoff, "highpass", fs=rate, output="sos")
+    sos = butter(2, cutoff, "highpass", fs=samplerate, output="sos")
-    fdata = sosfiltfilt(sos, data)
+    filtered_signal = sosfiltfilt(sos, signal)
-    return fdata
+
+    return filtered_signal
 
 
 def lowpass_filter(
-    data: np.ndarray,
+    signal: np.ndarray,
-    rate: float,
+    samplerate: float,
     cutoff: float
 ) -> np.ndarray:
     """Lowpass filter a signal.
@@ -78,21 +80,25 @@ def lowpass_filter(
     np.ndarray
         The filtered data
     """
-    sos = butter(2, cutoff, "lowpass", fs=rate, output="sos")
+    sos = butter(2, cutoff, "lowpass", fs=samplerate, output="sos")
-    fdata = sosfiltfilt(sos, data)
+    filtered_signal = sosfiltfilt(sos, signal)
-    return fdata
+
+    return filtered_signal
 
 
-def envelope(data: np.ndarray, rate: float, freq: float) -> np.ndarray:
+def envelope(signal: np.ndarray,
+             samplerate: float,
+             cutoff_frequency: float
+             ) -> np.ndarray:
     """Calculate the envelope of a signal using a lowpass filter.
 
     Parameters
     ----------
-    data : np.ndarray
+    signal : np.ndarray
         The signal to calculate the envelope of
-    rate : float
+    samplingrate : float
         The sampling rate of the signal
-    freq : float
+    cutoff_frequency : float
         The cutoff frequency of the lowpass filter
 
     Returns
@@ -100,6 +106,7 @@ def envelope(data: np.ndarray, rate: float, freq: float) -> np.ndarray:
     np.ndarray
         The envelope of the signal
     """
-    sos = butter(2, freq, "lowpass", fs=rate, output="sos")
+    sos = butter(2, cutoff_frequency, "lowpass", fs=samplerate, output="sos")
-    envelope = np.sqrt(2) * sosfiltfilt(sos, np.abs(data))
+    envelope = np.sqrt(2) * sosfiltfilt(sos, np.abs(signal))
+
     return envelope